{'fraction_non_alphanumeric': 0.04880044276535842, 'fraction_numerical': 0.01768655100223789, 'mean_word_length': 4.670350661754673, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Singular perturbations approach to localized surface-plasmon resonance: Nearly touching metal nanospheres

Metallic nano-structures characterised by multiple geometric length scales support low-frequency surface-plasmon modes, which enable strong light localization and field enhancement. We suggest studying such configurations using singular perturbation methods, and demonstrate the efficacy of this approach by considering, in the quasi-static limit, a pair of nearly touching metallic nano-spheres subjected to an incident electromagnetic wave polarized with the electric field along the line of sphere centers. Rather than attempting an exact analytical solution, we construct the pertinent (longitudinal) eigen-modes by matching relatively simple asymptotic expansions valid in overlapping spatial domains. We thereby arrive at an effective boundary eigenvalue problem in a half-space representing the metal region in the vicinity of the gap. Coupling with the gap field gives rise to a mixed-type boundary condition with varying coefficients, whereas coupling with the particle-scale field enters through an integral eigenvalue selection rule involving the electrostatic capacitance of the configuration. By solving the reduced problem we obtain accurate closed-form expressions for the resonance values of the metal dielectric function. Furthermore, together with an energy-like integral relation, the latter eigen-solutions yield also closed-form approximations for the induced-dipole moment and gap-field enhancement under resonance. We demonstrate agreement between the asymptotic formulas and a semi-numerical computation. The analysis, underpinned by asymptotic scaling arguments, elucidates how metal polarization together with geometrical confinement enables a strong plasmon-frequency redshift and amplified near-field at resonance.


I. INTRODUCTION
Advances in micro-fabrication technology along with new innovative ideas in nano-photonics and metamaterials have rekindled interest in the optics of small metal particles and structures [1][2][3]. Nano-metallic configurations exhibit unique optical characteristics, which are largely associated with the excitation of localized surface-plasmon (LSP) modes: conjoint resonant oscillations of electron-density and electromagnetic fields. LSP resonances are manifested in spectral peaks of the optical cross-sections and in significant near-field enhancements. Such phenomena are attractive for sensing applications [4], plasmonic manipulation of light below the diffraction limit [5], enhanced Raman scattering [6,7], plasmonic rulers [8], lasers [9] and nanoantennas [10][11][12], second [13] and higher-order [14] harmonic generation, and photovoltaics [15]. LSP resonance is exemplified by light scattering from a nano-metric metal sphere. In the quasi-static limit, the near-field, which appears to be driven by a uniform time-harmonic field, is derived from an electric potential governed by Laplace's equation. The induced electricdipole moment is readily found to be proportional to the Clausius-Mossotti factor (ǫ − 1)/(ǫ + 2), where ǫ is the frequency-dependent dielectric function of the metal particle relative to its dielectric surrounding. The evident divergence at the Fröhlich value ǫ = −2 is associated with the existence of a solution to the homogeneous quasistatic problem (in which the incident field is absent). At visible and ultraviolet frequencies, a good approximation for the Drude dielectric function is [16] Re where ω is the frequency of light, and ω p and γ are the plasma and collision frequencies, respectively (representative values for noble metals: ω p ≈ 9 eV, γ ≈ 0.13 eV). Owing to the smallness of γ/ω p , a damped resonance occurs when Re[ǫ] ≈ −2. LSP frequencies and the magnitude of resonant enhancement can be "designed" by shape and material variations [17][18][19][20]. Appreciable modification is achievable with "singular" configurations characterised by nonsmoothness or multiple-scale geometries: sharp corners [21,22]; elongated particles [23][24][25]; tips [26]; nearly touching [27][28][29], touching [30,31], and embedded [32] particle pairs; "plasmonic molecules" [20]; and particlewall configurations [33]. A recurring feature encountered in multiple-scale geometries is the existence of modes, that with increasing scale disparity, significantly redshift to low frequencies and large −Re [ǫ] values.
Singular configurations of the latter sort, which enable a particularly strong near-field enhancement [34], have been extensively studied through numerical simulations [35][36][37][38][39][40], analytic solutions in separable coordinate systems [21,30,41], and the powerful method of transformation optics (TO) [42,43]. In two dimensions TO furnishes explicit solutions [22,27,28] that can be generalised to address effects of retardation [44] and non-locality [45]. TO is less straightforward to apply in three dimensions, and typically furnishes an efficient computational scheme [46] rather than simple analytic expressions; nevertheless, the method has been successfully applied in various scenarios, e.g. in the design of broadband plasmonic structures [31] and in the study of van der Waals forces [29]. Yet another approach, motivated by general scaling arguments [47], is based on approximate reductions of implicit infinite-series solutions for nano-wire and nanosphere dimers [48][49][50]; we delay a critical discussion of the accuracy of the expressions in Refs. 48 and 50 to §VI.
The above methodologies are "top-bottom" -necessarily relying on "exact" analytical or numerical solutions, which then need to be reduced or evaluated in singular geometric limits (e.g. small gaps for dimers, increasing aspect ratio for elongated particles). In this paper, we suggest an alternative "bottom-up" approach, which is specifically suitable for analyzing the extreme physical characteristics of LSP resonance in multiple-scale metallic nano-structures. The approach is built upon the paradigm of singular perturbation methods, where an asymptotic rather than an exact solution is sought by carrying out an ab initio singular analysis. In particular, applying asymptotic scaling arguments in conjunction with the method of matched asymptotic expansions [51] allows exploiting the spatial non-uniformity of the singular LSP modes towards systematically decomposing the problem governing these into several considerably simpler ones. The main outcomes are asymptotic expressions for the resonance frequencies, along with a minimalistic description of the corresponding modes. In the context of excitation problems, when losses are low an asymptotic approximation for the field enhancement may be obtained with little additional effort.
To demonstrate this approach, we shall focus in this paper on the prototypical configuration of a pair of nearly touching metal nano-spheres [35,37,39,47,48,50,[52][53][54][55][56][57]. We limit ourselves to a quasi-static treatment of the near field, and focus on the axisymmetric and longitudinal (bonding) modes that singularly redshift with vanishing gap width, and the excitation of these by a long-wave electromagnetic plane wave polarized with the electric field along the line of centers. The paper proceeds as follows. In §II we formulate the problem and note useful generalities. In §III we consider the plasmonic eigenvalue problem and obtain closed-form expressions for the dielectric-function resonance values. In §IV we employ the eigenvalues and eigen-potentials found in §III to obtain formulas for the induced dipole moment and field enhancement in the gap under resonance conditions. In §V we present a validation of our results against a seminumerical solution of the quasi-static problem. Lastly, in §VI we recapitulate our results, compare with related works, and discuss recent developments along with future directions. by their radius a, and their frequency-dependent dielectric function relative to the background material, ǫ 0 ǫ, ǫ 0 being the vacuum permittivity. The minimum separation between the spheres is denoted by 2ha. The pair is subjected to an incident time-harmonic electromagnetic plane wave of frequency ω, polarized such that the electric field is in a directionî parallel to the line of sphere centers, see Fig. 1. We shall focus on the quasistatic limit of sufficiently small particles, a ≪ c/(ω |ǫ|), c being the speed of light in the surrounding dielectric [16,58]. In this description, the pair experiences the incident wave as a uniform time harmonic fieldîERe[e −iωt ], where E is a real constant magnitude, and the electric near-field is approximately curl-free, derived from a timeharmonic potential aERe[ϕe −iωt ]. Henceforth, we normalise all lengths by a. In particular, we denote by x the dimensionless position vector relative to the centre of the gap. It is convenient to employ a cylindrical coordinate system (r, z, φ), with origin at x = 0 and a z-axis directed alongî (see Fig. 2). The problem at hand is to find a continuous potential ϕ that satisfies where E is unity in the dielectric medium and ǫ in the metal spheres, together with the far-field condition The additive freedom of ϕ is removed in (3) by choosing the disturbance from the incident field to attenuate at large distances.

B. General properties
The following generalities can be inferred from (2)-(3), see also Refs. [3,36,[59][60][61]. First, the potential is axisymmetric about theî axis, i.e. ∂ϕ/∂φ = 0. In addition, the potential is antisymmetric about the plane z = 0. It is therefore possible to consider only the halfspace z > 0, with We next note that each of the spheres is overall charge free. That is, where the integration is over an arbitrary closed surface entirely within the dielectric medium, which encloses one or both spheres. Thus ϕ lacks a monopole term in its farfield expansion (3), which can therefore be extended to The complex-valued constant µ is the induced dipole moment of the pair, normalised by 4πa 3 ǫ 0 ǫ d E (ǫ d being the relative dielectric constant of the surrounding medium). Green identities applied to ϕ and its complex conjugate ϕ * , in conjunction with (2) and (6), can be shown to give where the integration is over the volume of the two spheres. Eq. (7) is a work-dissipation relation: The work done by the incident field on the effective induced dipole of the dimer is lost to Ohmic dissipation in the metal. This relation serves to remind us that we are in the quasistatic limit, where radiation losses are negligible. We note that 4πIm[µ] gives the quasi-static approximation for the absorption cross-section normalised by ωa 3 /c [2]. Excitation problems of the above sort are ill posed for a set of discrete, real and negative, ǫ values, for which the homogeneous variant of the problem obtained by omitting the incident field in (3) possesses non-trivial solutions. Since in practice Im[ǫ] = 0, such solutions are often termed "virtual" LSP modes. Consider now an eigen-potential ϕ r corresponding to one such eigenvalue ǫ r . At large distances, where µ r is arbitrary. We note that the integral relation (7) does not apply to ϕ r , since condition (3) is used to derive it. Say, however, that ϕ is the solution to the actual quasi-static problem, with ǫ = ǫ r + ∆ǫ. Applying the Green identities to the pair of potentials (ϕ, ϕ r ), one can easily show that Whereas (7) bears a clear physical significant, we shall find in §IV that (9) is more useful for deriving quantitative estimates for the resonant field enhancement.

C. Near-contact limit
Henceforth we consider the near-contact limit where h → 0. Specifically, our goal will be to derive closed-form asymptotic formulas for the eigenvalues ǫ (n) , and the values at resonance conditions (i.e. when Re[ǫ] = ǫ (n) ) of the induced-dipole moment µ and the field enhancement in the gap, G = − (∂ϕ/∂z) x=0 .

III. LSP MODES
In this section we consider the boundary eigenvalue problem obtained by replacing the far field (3) with the homogeneous condition ϕ → 0 as |x| → ∞. It is clear that the only relevant modes are those that are both longitudinal, i.e. anti-symmetric about the plane z = 0, and axisymmetric.

A. Scaling arguments
How do the eigenvalues ǫ scale with h? In the vicinity of the gap the spherical boundaries can be approximated by two oppositely facing paraboloids. The separation therefore remains O(h) for O(h 1/2 ) radial distances. We shall refer to the subdomain of the dielectric medium where z = O(h) and r = O(h 1/2 ) as the gap domain, and assume without loss of generality that the eigenpotential ϕ corresponding to ǫ is O(1) there; since ϕ is antisymmetric about the z plane, the transverse gap field must be O(1/h). Consider now the particle-scale dielectric "outer" domain, away from the gap, where ϕ varies over distances |x| = O(1) and is at most O(1). It follows from (4) that ϕ = O(h) when z = O(h) and r = O(1). This implies that the gap potential must vary in the radial direction over O h 1/2 distances. Consider next the metal sphere in z > 0, in the vicinity of the gap boundary. Continuity of potential implies that ϕ is O(1) there, and varies rapidly over O h 1/2 distances along the boundary. On the latter length scale, the metallic volume appears transversely semi-infinite. Thus, the symmetry of Laplace's equation implies comparably rapid transverse variations. We shall refer to the metal subdomain where r, z = O(h 1/2 ) as the pole domain.
The scaling of ǫ now follows by considering the continuity of electric displacement across the gap-pole interface: The product of ǫ and the O(h −1/2 ) transverse field in the metal pole must be on the order of the O(1/h) transverse field in the gap. We accordingly find where, as we shall see, the pre-factor λ itself depends weakly upon h, being a function of ln(1/h). This more subtle feature can also be anticipated from scaling arguments, as we now discuss. The largeness of |ǫ| suggests that, on the particlescale, the spheres are approximately equipotential. Thus, the particle-scale potential distribution in the dielectric medium is approximately that about a pair of nearly touching spheres held at uniform potentials, say ±V. But the latter distribution implies a net charge on each sphere, in contradiction with (5). This means that the deviation in the vicinity of the gap of ϕ from the particlescale distribution must be consistent with a localized accumulation of charge such that the overall charge of each sphere is zero. In fact, the contributions to the integral in (5) from the gap and outer dielectric domains are readily seen to be comparable: The former is given by the integral of an O(1/h) transverse field over an O(h) surface area, while the latter is obviously O(1). This implies a connection between the solution in the gap and pole "inner" domains, and the capacitance of a pair of perfectly conducting spheres held at fixed potentials ±V, which is known from electrostatic theory to depend logarithmically upon h as h → 0 [62].

B. Gap and pole domains
We begin constructing the eigen-potentials by considering the gap domain. To this end, we introduce the stretched coordinates The gap domain is where R, Z ∼ O(1), with R > 0 and Z bounded by the spherical surfaces (see Fig. 3) As in §III A, we assume without loss of generality that ϕ = O(1) in the gap domain. We accordingly pose the asymptotic expansion Substituting (13) into Laplace's equation, written in terms of R and Z, yields at leading order ϕ ∼φ The resonant modes are constructed in §III by matching simple asymptotic expansions describing the gap, pole, and particle-scale "outer" domains.
Integration of (14) in conjunction with (4) gives where the function A(R) remains to be determined. Consider next the pole domain of the metal sphere in z > 0, described by the alternative pair of stretched coordinates The pole domain is where R,Z ∼ O(1), with R > 0 andZ bounded from below by the sphere surfaceZ ∼ O(h 1/2 ); under the rescaling (16), the pole domain is to leading order a semi-infinite half space. We expand the poledomain potential as whereΦ is governed by Laplace's equation, Expanding the eigenvalue ǫ according to (10), requiring leading-order continuity of potential and electric displacement yields Eqs. (19) may be combined to yield the mixed-type boundary condition which involvesΦ alone. We also find the relation which provides the gap potential once the pole potential has been determined. As discussed in §III A, the particle-scale metal domain is equipotential to leading order. Asymptotic matching thus furnishes the "far-field" conditionΦ where V is the above-mentioned uniform potential of the sphere in z > 0. For later reference, we note that (21) and (22) The pole potentialΦ is governed by Laplace's equation (18) in the half spaceZ > 0, the boundary condition (20), and the far-field condition (22). We shall see in §III F that this boundary-value problem possesses nontrivial solutions for essentially arbitrary λ > 0 and V, excluding a discrete set of λ values for which a solution exists only if V is set to zero. It is tempting to guess that the latter set furnishes the requisite eigenvalues; we shall see in §III E however that these only give a leading term in an expansion of λ in inverse powers of ln(1/h). Avoiding such gross logarithmic errors, the eigenvalues cannot be solely determined from local considerations, as already anticipated in §III A. We accordingly divert at this stage to consider the particle-scale dielectric domain.

C. Outer domains and global constraint
In the outer limit, h → 0 with x = O(1), the spheres appear to be in contact to leading order. It is therefore convenient to introduce the tangent-sphere coordinates defined through [63] The plane z = 0 is mapped to ξ = 0, the upper-half-space sphere to ξ ∼ 1 + O(h), and the origin to η → ∞ (see Fig. 3). We expand the potential in the outer dielectric domain as and in the outer metal domain as the ± sign corresponding to the sphere in z ≷ 0, respectively.
Forφ we require a solution to Laplace's equation that attenuates at large distances and satisfiesφ = V at ξ = 1 andφ = 0 at ξ = 0. In addition, the solution must asymptotically match the gap and pole potentials. Standard Hankel-Transform methods yield the requisite potential as [62] where J 0 is the zeroth-order Bessel function of the first kind. It may be readily shown thatφ ∼ 2Vz/r 2 as η → ∞, whereby (26) matches with the gap and pole domains. Asymptotic evaluation of (26) for small ξ and η, together with (23), yields the asymptotic behaviour of (26) at large distances [cf.
The leading-order matching of (26) with the solution in the gap and pole regions has been achieved, apparently without introducing any additional constraint on the problem governingΦ. Instead of continuing the asymptotic analysis to the next order in search of the awaited eigenvalue selection rule, it proves simpler to follow the discussion in §III A and consider the charge constraint (5) applied to the sphere in z > 0. In accordance with the discussion there, an attempt to calculate the contributions to the integral in (5) from the gap and outer region separately results in diverging integralsthis is typical in cases where contributions from overlapping asymptotic domains are comparable [51]. Following the analysis in Ref. [62], we write (5) as where η 0 , R 0 ≫ 1 are connected through h 1/2 R 0 = 2/η 0 + O(1/η 3 0 ). The first integral in (28), which represents the excess contribution from the gap region, can be written as The second integral in (28), which represents the contribution of the outer domain, can be worked out to be where γ e ≈ 0.5772 is the Euler-Gamma constant. Noting that ln H(R 0 ) ∼ ln(2/h) − ln η 2 0 , and invoking (21), the singular parts of the two contributions cancel out, and we find the anticipated "global" constraint:

D. Reduced eigenvalue problem
Together with (31), the problem governing the pole potentialΦ constitutes a reduced boundary-eigenvalue problem. It is convenient to recapitulate this problem in terms of the disturbance potential, ψ(R,Z) =Φ − V. We have Laplace's equation in the upper half-spaceZ > 0, the integro-differential varying-coefficient boundary condition (33) atZ = 0 [we will sometimes write −V for the right-handside, see (31)], and the far-field condition The logarithmic dependence of the eigenvalues λ and their corresponding eigen-potentials ψ upon h is evident from (33). We also see that a rough leading-order solution, entailing an O(ln −1 1 h ) relative error, may be obtained by neglecting the right hand side of (33). In the following subsections we proceed in two paths: (i) in §III E we calculate succesive terms in inverse logarithmic powers; (ii) in §III F we tackle the reduced eigenvalue problem as a whole, thereby finding an "algebraically accurate" approximation for the resonance ǫ values.

E. Two terms in a logarithmic expansion
We here expand the potentials ψ and the eigenvalues λ as At leading order, the boundary condition (33) reads We look for a solution in the form where is the zeroth-order Hankel transform [64]. Applying the Hankel transform to Laplace's equation (32), a solution that satisfies (34) isψ where Y 0 (s) depends on the conditions atZ = 0. Transforming the boundary condition (36), we find the s-space differential equation The first term in (40) is obtained through integration by parts twice, which is only applicable if [64] Y 0 ≪ 1 s 2 and In addition, Y 0 (s) must decay sufficiently fast to allow convergence of (37). Defining substitution shows that T (p) satisfies Laguerre's equation Non-singular solutions of the latter, which are consistent with (41) and the required attenuation rate of Y 0 (s), exist only when the factor multiplying T in (43) is a nonnegative integer. In those cases the solutions are proportional to the Laguerre polynomials. Hence, the leadingorder eigenvalues are The corresponding transformed solutions are where K (n) is an arbitrary constant, and L n is the Legendre polynomial of order n. It is straightforward to invert (37) to get the eigen-potentials ψ (n) 0 . For example, noting that L 0 (p) = 1 and L 1 (p) = 1 − p, the first two modes are and ψ (1) Proceeding to the next logarithmic order, we see that the problem governing ψ is similar to the one governing ψ (n) 0 , only that the homogeneous boundary condition (36) is replaced by an inhomogeneous one. Dropping the (n) superscript for the moment, the latter condition follows from (33) as From the preceding analysis it is clear that, when λ 0 is given by (44), the problem governing ψ 1 possesses a nontrivial homogeneous solution. Hence, the Fredholm alternative suggests a solvability condition on the forcing terms in (48). Indeed, by subtracting (48) multiplied by Rψ 0 /H from (36) multiplied by Rψ 1 /H, followed by integration with respect to R and use of Green's second identity, we find The integrals in (49) can be readily evaluated by use of Parseval's identity and the inversion-symmetry of the Hankel transform, together with the orthonormality of the Laguerre polynomials, see (A.1). One finds the simple result λ 1 = −2 √ 2 (λ 0 ) 2 , whereby reinstating the (n) superscript gives the two-term logarithmic expansion

F. Algebraically accurate result
It is clear that the logarithmic error in (50) limits the usefulness of that result to exceedingly small h. To improve on this, we now return to the reduced problem of §III D and attempt a complete solution for the λ eigenvalues, treating logarithmic terms in par with O(1) terms. This is equivalent to calculating all of the terms in expansion (35). We begin by writing the solution as [cf. (37), (38) and (39)] Applying the Hankel transform to (33), we find 1 s d ds s dY ds where δ(s) is the Dirac delta function. Eq. (52) is to be interpreted in the sense of distributions [65], which is easier to do when identifying the Hankel transform as a two-dimensional Fourier transform of a radially symmetric function, s being the radial distance in Fourier-space. Crucially, conditions (41) are now inapplicable; rather, balancing singularities at the origin yields the condition Following ( dT dp +ñT = 0, p > 0, whereñ is defined through We arrive at an eigenvalue problem consisting of (54), the "boundary" condition and the condition that T does not grow too fast as p → ∞. The solutions found in §III E, which correspond to non-negative integerñ and V = 0, are inconsistent with condition (56). Hence, in the present scheme, where logarithmic errors are not tolerated, V = 0. In the latter case, excluding non-negative integerñ, solutions of (54) which are logarithmically singular as p → 0 and have acceptable behaviour at large p are proportional to U (−ñ, 1, p), where U is the confluent hypergeometric function of the second kind [66]. Noting that (1) (57) as p → 0, where Γ(x) is the Gamma function and Ψ(x) = Γ ′ (x)/Γ(x) is the Digamma function, condition (56) is satisfied by It remains to satisfy the global condition, that is, to relate V to the right hand side of (33) [cf. (31)]. It is convenient to first rewrite that condition as The integral on the right can be evaluated by using Parseval's identity. Noting that the Hankel transform of 1/H is 2K 0 (s √ 2), K 0 being the modified Bessel function of the second kind, that integral is Given (57), and since K 0 (x) ∼ − ln(x/2) − γ e as x → 0, the apparent logarithmic singularity of the bracketed terms cancels, whereby the sifting property of the Dirac delta function can be effected to show that (60) = 2 ln 4 + 2Ψ(−ñ) + 2γ e . Substitution of this result into (59) furnishes a transcendental equation forñ: Since Ψ(x) diverges at non-positive integer x, it is easy to see that (61) has multiple non-integer solutionsñ which (very slowly) approach n = 0, 1, 2, as h → 0 [67]. Recapitulating, definition (55) furnishes the eigenvalues as whereñ is the solution of (61) which approaches n = 0, 1, 2, . . . as h → 0. We reiterate that (62) is algebraically accurate, viz., together with (10) it provides the resonance values of the relative dielectric function up to a relative error which scales with some power of h. In particular, (62) is far more accurate then (50). We note that it is straightforward to recover (50) by expanding Ψ in (61) about its poles. Higher order logarithmic terms can be derived in a similar manner. Invertingψ provides the distribution of the eigenpotentials in the pole to algebraic order [cf. (46) and (47)]. The dominantly transverse field in the gap, which is a function of R alone, is then obtained (up to a minus sign) as ∂Φ/∂Z = A(R) = ψ(R, 0)/H(R). The radial gap-field distributions of the first three modes are plotted in Fig. (4).

A. Resonance conditions
We now return to the excitation problem formulated in §II, where the far-field condition (3) applies, and ǫ is a complex-valued function of frequency. Whereas a negative real ǫ is unphysical, the reduced Drude model (1) implies that −Re[ǫ] ≫ Im[ǫ] as long as ω ≫ γ. Resonance frequencies are safe within this range: even for an extremely narrow gap, say h = 10 −4 , Eqs. (1) and (10) suggest ω ≈ 0.1ω p -still quite larger than γ for silver and gold. Thus, towards studying the damped near-field resonance of the system, it is useful to consider a resonance scenario Here, for a given n = 0, 1, 2, . . ., ǫ (n) is a resonance value is given by (62), and iǫ i is the imaginary deviation from the resonance values. Given (63), it turns out that an asymptotic description of the field distribution -and in particular, closedform approximations for the induced-dipole moment µ and field enhancement in the gap G -can be obtained based on the LSP modes calculated in §III. The crucial observation is that the nth LSP mode is "excited" under (63), viz., the leading-order potential is of the form calculated in §III, but with an asymptotically large multiplicative factor. Note that if the potential was not ≫ 1, (63) would imply a contradictory leading-order problem which, on the one hand, possesses a non-trivial homogeneous solution, and on the other, is forced at large distances by the incident field. Thus, the integral relation (9) can be written, asymptotically, as We can now substitute our solutions from §III (allowing for a large pre-factor) into (64), and solve for the amplification factor (say, in terms of µ). We note that, in contrast to (64), the work-dissipation relation (7) is ambiguous with respect to phase and hence cannot be employed towards making quantitative assessments without further information. We once more precede a detailed analysis with scaling arguments.

B. Scaling arguments
Let ϕ g denote the order-of-magnitude of the potential in the gap region. The asymptotic structure of the potential then follows from the description provided in §III. In particular, the transverse field in the gap is G = O(ϕ g /h), the potential in the pole region is O(ϕ g ), and on the particle scale the sphere is approximately equipotential, with V = O(ϕ g / ln(1/h)), see (31). Note that (27) implies µ = O(V). Consider now the integral on the right hand side of (64) [the work-dissipation relation (7) can also be employed]. The dominant contribution arises from the pole domains, where the potential varies rapidly. Indeed, the volume of the pole domains is O(h 3/2 ) and the integrand there is O(ϕ 2 g /h), leading to an O(ϕ 2 g h 1/2 ) estimate. The contribution from the approximately equipotential O(1) volume is O(ϕ 2 g h) at most. We thus find from (64) that ϕ g = O(h −1/2 /(ǫ i ln(1/h))), from which estimates for µ and G readily follow as

C. Asymptotic formulas
To go beyond scaling arguments, we use Green's first identity together with (5) to rewrite (64) as where ψ denotes the deviation of the potential in the (z > 0) pole domain from V, see §III. Applying Parseval's theorem to the right-hand-side of (64), it becomes [cf. (51)] Upon substitution of Y (see below), and noting that µ ∼ π 2 V/3 [cf. (27)], (66) can be solved for V, and hence µ.
The field enhancement we can then be obtained from where Φ is the gap potential [cf. (13)], and the derivative is evaluated at R, Z = 0. Using (19), G can also be written as where the derivative is evaluated at R = 0,Z = 0, or in terms of Y (s), The calculation can be performed at two levels, depending on whether the logarithmically or algebraically accurate expressions for Y and λ are employed in the above formulas. We start with the former case, substituting Y 0 (s) for Y (see (45)) and λ 0 for λ (see (44)). Noting (31), the arbitrary constant K in Y 0 is to within a logarithmic relative error K ∼ V ln 1 h /λ 0 . Together with the integrals (A.1) and (A.2), we find the explicit closedform expressions and Expressions (71) and (72) have the advantage of being simple, but since they entail an O(1/ ln(1/h)) error, they are hardly accurate, as we shall later demonstrate. We note the agreement with the scaling results (65), and that the excited resonant mode is imaginary, i.e. out of phase with the incident field.
To obtain algebraically accurate formulas, we substitute (58) for Y and (62) for λ. Making use of the integrals (A.3) and (A.4), the solution scheme described above yields and where Ψ ′ (x) is the TriGamma function, the derivative of the DiGamma function, 2 F 1 is the hypergeometric function, andñ is as before the solution of (61) that approaches n = 0, 1, 2, . . . as h → 0. We note that (71) and (72) can be derived as leading approximations of (73) and (74), respectively, by noting thatñ − n ∼

V. COMPARISON WITH SEMI-NUMERICAL SOLUTION
We have validated our asymptotic formulas against a semi-numerical solution of the quasi-static problem. We employ a scheme which is based on the separability of Laplace's equation in bi-spherical coordinates. Briefly, it consists of writing down a formal solution in the form of two infinite series, and determining the coefficients by numerically solving a truncated tri-diagonal system of coupled algebraic equations. The relevant formulas are available elsewhere [3,52,69].
An example calculation is shown in Fig. 5, where we have employed the empirical ǫ-data for silver provided by FIG. 5. Absolute magnitude of the dimensionless induceddipole moment µ as a function of energy ω and h, half the dimensionless gap width. The coloured contours represent the semi-numerical solution of the quasi-static problem. The enumerated white lines depict the resonance frequencies of the first three modes as predicted by the asymptotic formula (62), together with (10). We employ the empirical silver data of Johnson and Christy [68]. Johnson and Cristy [68]. Fig. 5 shows the absolute magnitude of the induced-dipole moment µ as a function of the dimensional energy ω in eV units, for a range of h. We observe magnitude peaks which, as expected, redshift as h → 0. The white lines in Fig. 5 depict the first three resonance frequencies as predicted from the asymptotic formula (62) in conjunction with the data of Johnson and Cristy. For small h the latter lines successfully trace the amplitude peaks. Detailed results for the zeroth-mode ǫeigenvalues are shown in Fig. 6. It depicts as a function of h the numerically obtained eigenvalues (symbols) along with the predictions of (62) and the two-term expansion (50). The asymptoticness of (62) as h → 0 is evident, with good agreement for h 0.02 (i.e. radius-to-gap ratio of ≈ 25). As expected, the logarithmically accurate approximations are quite poor. Note the approach at large h of ǫ to −2, the Fröhlich value of isolated spheres.
Consider next the values of µ and G. Focusing again on the zeroth-mode, we have recorded for a range of h the numerical values attained by µ and G at their zerothmode peaks. The absolute magnitude of these values (which are predominantly imaginary at resonance) are depicted by the symbols in Fig. 7. The latter numerical values are compared with the logarithmically accurate asymptotic approximations (71) and (72) -dashed lines, and the algebraically accurate approximations (73) and (74) -thick lines. Again we see that the latter approximations are rather good, whereas the former are poor for realistically small h.

VI. DISCUSSION
We have employed singular perturbation methods towards analyzing the longitudinal surface modes of a pair of metal spheres in near contact, and the quasi-static excitation of these by an electromagnetic plane wave. Whereas the techniques of matched asymptotic expan-  (71) and (72), and the algebraically accurate approximations (73) and (74). We employ the empirical silver data of Johnson and Christy [68].
sions are routine in many areas of research, it appears that their application in the present context is novel. As demonstrated herein, these allow describing limits of physical extremity in a relatively simple and intuitive manner, notably without the prerequisite of an analytically cumbersome or computationally heavy "exact" solution. In particular, the present analysis has furnished accurate asymptotic expressions for the resonance values of the dielectric function (or frequencies) -see (62) and (10), and the values at resonance of the induced dipole moment and field-enhancement in the gap -see (73) and (74). These expressions have been validated in §V against a semi-numerical solution of the quasi-static problem. The above-mentioned expressions are asymptotic in the limit h → 0, with algebraic accuracy. That is, they involve a relative error on the order of some power of h. These results were preceded with simpler expressions -see (50), (71) and (72), which incur however a much larger relative error on the order of 1/ ln(1/h). As the comparison in §V unsurprisingly suggests, the latter logarithmic expressions are far less accurate for reasonably small h. Nevertheless, the analysis leading to these revealed the essential structure of the longitudinal LSP modes, thereby laying a path towards their algebraic counterparts. We note that similar logarithmic expressions have been derived by approximate manipulations of implicit series solutions in bi-spherical coordinates [48,50]. The limited value of these expressions, however, was not emphasised in those works. To be specific, the leading term in (50) agrees with both references, and (72) agrees with the expression given in [50]. The logarithmic correction for ǫ given in Ref. [50], however, appears to posses a factor-2 error when compared with the second term in (50).
In addition to furnishing closed-form approximations, our analysis offers insight into the physical and geometric circumstances enabling a strong frequency redshift and field enhancement. The former has to do largely with the gap morphology, which enables stronger field localization in the dielectric than in the metal -a prerequisite for large-|ǫ| eigenvalues. In fact, the scaling of ǫ with 1/h 1/2 was shown in §III A to follow from the particular imbalance in transverse localization pertinent to a locally parabolic boundary. This scaling should therefore apply quite generally to the large-|ǫ| resonances of essentially arbitrary plasmonic two-and three-dimensional clusters in near contact, or to particles in near contact with a substrate (see Ref. [47] for an alternative intuitive argument). The origin of field enhancement is more subtle. The amplitude of an excited mode, represented e.g. by the configuration-scale induced-dipole moment, is determined by a balance of work and dissipation, which is in turn affected by the geometrically enabled localization. Disregarding a weak logarithmic dependence upon h, the latter enhancement is on the order of |ǫ/ǫ i |, which in the near-contact limit is O(h −1/2 /ǫ i ). The giant O(h −3/2 /ǫ i ) enhancement in the gap represents a further O(1/h) geometric amplification relative to the induced particle-scale field. We note that in the corresponding planar problem of a nanowire dimer, where scaling arguments similar to those employed herein reveal that the longitudinal modes are determined solely from the gap morphology, the enhancement is smaller, O(h −1 /ǫ i ) [49]. This fundamental difference is apparently overlooked in Ref. 3, where in the case of a spherical pair the resonant field enhancement in the gap is estimated incorrectly on pg. 278 [cf. (65)].
For the sake of demonstration, we have focused on the axisymmetric longitudinal modes of the spherical dimer configuration, which strongly redshift with vanishing separations and are the ones excited by a field polarized along the line of sphere centers. In many applications the latter modes play the key role, but this is not always the case, and it is therefore desirable to apply the present approach to LSP modes having different symmetries [48]. For example, it has been shown that nonaxisymmetric modes of high azimuthal order provide the dominant contribution to van der Waals (vdW) energies at small separations [70]. In Ref. 70, approximate expressions for the vdW energies of sphere dimers and a sphere near a wall were derived by reduction of implicit infiniteseries solutions in the double limit of small separations and high azimuthal mode number. While, as already discussed, similar methods gave in Ref. 48 only logarithmic accuracy in the case of longitudinal modes, the dominant modes in the vdW calculation are more highly localized in the gap and accordingly do not involve a logarithmic coupling with the particle-scale field.
Our methodology can be applied to study various other plasmonic nano-structures that are characterised by multiple length scales. In particular we note the related particle-substrate configuration, which can be addressed along the lines of the present analysis, the substantially simpler two-dimensional problem of a nanowire dimer, and elongated particles. Further desirable generalisations are concerned with physical modelling. Our analysis originates from a quasi-static formulation, and focuses on the near-field response of a small nano-metric particle (still large compared with the gap width). Given our expressions for the induced-dipole moment, estimates for the optical cross section can be obtained via the usual connection formulas [2,16,58]. Ad hoc techniques for improving and extending such estimates to larger nano-structures, wherein the quasi-static induced-dipole is renormalised to account for radiation damping, have proven quite useful [44,71]. More rigorous however would be to connect the near and far fields through systematic application of matched asymptotic expansions [65,72], taking special account of the largeness of |ǫ|, the largeness of the near-field at resonance, and the multiple-scale geometry of the metallic structure; strictly speaking, each of these disparities entails slaving a distinct small perturbation parameter to the Helmholtz parameter aω/c.
Probably the most important direction is to incorporate into the present study non-classical physics, which at sub-nanometric separations ultimately suppress the frequency redshift and field enhancement [29,45,73]. In the framework of the hydrodynamic Drude model [74,75], my co-workers and I recently showed [76] that the nearcontact asymptotics predicted by the latter "nonlocal" model constitute a renormalisation of those predicted by the classical "local" model considered herein. This finding greatly extends the utility of the present study. ACKNOWLEDGMENTS I wish to thank Vincenzo Giannini, Richard Craster and Stefan Maier for fruitful discussions on this work.