Plasmon dispersion in a multilayer solid torus in terms of three-term vector recurrence relations and matrix continued fractions

Toroidal confinement, which has played a crucial role in magnetized plasmas and Tokamak physics, is emerging as an effective means to obtain useful electronic and optical response in solids. In particular, excitation of surface plasmons in metal nanorings by photons or electrons finds important applications due to the engendered field distribution and electromagnetic energy confinement. However, in contrast to the case of a plasma, often the solid nanorings are multilayered and/or embedded in a medium. The non-simply connected geometry of the torus results in surface modes that are not linearly independent. A three-term difference equation was recently shown to arise when seeking the nonretarded plasmon dispersion relations for a stratified solid torus (Garapati et al 2017 Phys. Rev. B 95 165422). The reported generalized plasmon dispersion relations are here investigated in terms of the involved matrix continued fractions and their convergence properties including the determinant forms of the dispersion relations obtained for computing the plasmon eigenmodes. We also present the intricacies of the derivation and properties of the Green’s function employed to solve the three term amplitude equation that determines the response of the toroidal structure to arbitrary external excitations.


Introduction
A variety of particles are emerging in response to the needs of nanoscale functionality for modifying existing material properties or for creating new properties, see for example the editorial by Roco and Pinna [1]. Apart from the complexity of the material at the atomic, molecular, and cluster scales, both the local geometry and size of the particles have been shown to affect the response of single particles as well as many particle systems. In plasmonics [2], these experimental observations are frequently quite satisfactorily accounted for by theoretically invoking nonretarded electrodynamics. The scattering properties of sub-wavelength structures can thus be obtained in the quasi-static limit, where the field is primarily given by the scalar electric potential. For geometries that permit analytical solutions, such calculations can reveal the resonance behavior of the nanoparticle surface modes in response to electromagnetic excitation [3]. Therefore, surface plasmon dispersion relations may be obtained for the nanoparticles so that experiments can be expedited or measurement results can be better interpreted.
Recently, toroidal nanoparticles such as metal and dielectric nanorings, have gained considerable attention due to their potential use in trapping cold polar molecules [4], levitating and trapping dielectric nanoparticles [5], metamaterials [6][7][8], soft Coulomb interactions [9], light trapping in energy-harvesting devices [10], and plasmonic nanoantennas [11]. We here suffice by noting that our motivation for studying such structures parallels that for the extensive investigation of cartesian thin film stratified systems for development of optical filters, photonic band gap materials, and metamaterials [12][13][14][15]. With reference to the works of Love on the calculation of oscillatory modes of a cold plasma [16,17], the complete set of the dispersion relations of a single plasma ring in vacuum and its composite configurations were recently [18] obtained in the quasi-static limit from the separation of variables of the Laplace equation in the toroidal coordinate system, where the surface of a torus is obtained by fixing the value of one of the coordinates. The aim of this article is to rigorously investigate the analytical structure, existence, and convergence properties of the dispersion relations of a system of composite solid tori. Specifically, the derivation and properties of equations (28) and (29) in [18], obtained to analyze a k-layered composite toroidal structure with a single or multiple metal-dielectric interfaces, are presented. In doing so, we consider the confocal toroidal multilayers, shown in figure 1. We here note that for concentric multilayers, as graphically shown in [18], the Laplace equation is not separable [19]. This point was not emphasized in [18], although the analytical treatment pertained to confocal multilayers. The objective there was the experimentally realizable nanostructures of current interest to nanoscience, where the fabrication of multilayer nanorings is practically limited to singly or doubly coated rings such that the ratio of the coating thicknesses to the core ring radius remains small. As a result, the deviation of the confocal rings from concentricity remains negligible.
This article is organized as follows: in section 2, closely following the formulation in [18], we derive a canonical vector three-term recurrence relation that implicitly contains the eigenvalue spectrum corresponding to the quasi-static plasmonic modes of a k-layered toroidal structure. Using a formal argument followed by a rigorous mathematical proof, in section 3, we give the explicit forms of the plasmon dispersion relations for such a geometry in terms of matrix continued fractions (MCFs) and present the determinant forms of the dispersion relations to provide a numerical platform for determination of the resonance values of the involved dielectric functions. In section 4, we introduce the Green's function for the solution of the three term difference equation for an arbitrary continuous toroidal charge distribution. As examples of the application of the analytical results, we visualize the potential for some simple charge distributions. Concluding remarks are provided in section 5.

Model of a k-layered torus and the vector three-term recurrence relation
A multilayered torus can be described as a solid torus with toroidal surface m m = , 1 dielectric function ε 1 and minor radius r , 1 together withk 1 sublayers of confocal toroidal shells, each with dielectric function ε i , where = ¼ i k 2, , (see figure 1). Thus, a single solid torus corresponds to k=1, with no sublayer between the torus and the outside medium. With distances typically of the order of a few nanometers between the layers embedded in a medium, one can divide the space into + k 1 regions based on e e + , ..., , k 1 1 where ε k+1 denotes the dielectric function of the outside medium. The first region corresponds to interior of the solid torus and is given by  m m . 1 The remainingk 1 toroidal shells, described by   m m mi i1 ( = i k 2, ..., ), share the same focal length a with the solid torus via the relations [20] Finally, the + ( ) k 1 th region, which lies outside the k-layered torus, is described by  m m < 0 k (see figure 1 for the case k = 3). Recalling [16,18]   with constants C mn and D mn to be determined, as a suitable ansatz for the scalar electric potential m h j F( ) , , of a toroidal structure.
After dividing the space into + k 1 regions, we denote the associated potential in each region by F F + , ..., , implies that one can treat equations (6) and(7) separately for each integer m (the potential field has rotational symmetry with respect to the z-axis). With m maintained fixed, we may therefore suppress its notation in upcoming equations. Moreover, we adopt the notations P , Introducing theḱ 1 coefficient vectors and applying the boundary condition equation (6) implies n n n n with k×k bidiagonal matrices  n and  n (see appendix A, equations (A.3) and (A.4)). Both of these bidiagonal matrices have non-zero diagonal entries and hence are invertible. Thus, one can solve, for example,  n in terms of  n as The application of the second boundary condition equation (7) together with equation (10) gives the vector three-term recurrence For details regarding the algebraic manipulations resulting in equations (9)- (14) and definitions of the k×k matrices  , n  , n  ¢ , n  ¢ , n m D , E , 1 E , 2 see appendix A. Since the normal modes of the system are independent of the chosen set of coefficients, in our case the vector  , n one is naturally led to ask whether equation (11) would remain invariant with respect to a different choice of a coefficient vector. This is a relevant question and needs to be addressed since the equation system (A.2) can be represented using a different set of coefficient vectors than  n and  n given by equation (8). We show this invariance under the extra feasibility assumption, which stems from the physical basis of the problem. To see this, we generate a different set of coefficient vectors by rewriting equation (9) as: with the possibility of a negative sign in case a vector has been moved from one side of equation (15) to the another side, and the coefficient vectors  n and  n consist of the corresponding coefficients for each chosen column according to equation (15). Similar to equation (9), this is an alternative representation of equation (A.2). In general, neither  n nor  n has to be invertible. However, in order to obtain the vector three-term recurrence, and thus the dispersion relations, one should be able to solve one set of k variables, chosen from , in terms of the remaining k variables using the boundary condition equation (6). As a result, we must require that the rearrangement equation (16) is feasible; i.e., at least one of the matrices  n or  , n say  , n is invertible. In this case, we have     = -. It is easily seen that  n can be expressed as: where  1 and  2 are diagonal matrices with only zeros or ones in the diagonal such that   + = I , k 1 2 the identity matrix. Substituting  n in equation (17) gives  16) gives the vectorial three-term recurrence equation for  , n which is easily seen to be reduced to equation (11). Consequently, under the feasibility condition, equations (11) and (14) are the canonical vector threeterm recurrence relations for a k-layered toroidal structure.
Having solidified the form of equation (11), we close this section by considering the limiting case of equation (14), which will be used throughout the rest of the paper. For the proof of theorem 1, see appendix A. Theorem 1. Let R n be defined as in equation (14), then = m Therefore for large n, the vector three-term recurrence relation all roots of the characteristic equation for the vector three-term recurrence equation (19) are distinct with distinct moduli. This observation in the scalar case, i.e., k=1 in equation (11), implies the utilization of the classical results due to Perron [21] and Pincherle [22] (see also [23]), which prove the convergence of the obtained dispersion relations (see also [17]). Theorem 1 plays the key role in providing the application of Perron-Pincherle type theorem for the multidimensional case  k 2 (see [24][25][26] for details).

Dispersion relations for a k-layered torus and MCFs
While much of the discussion on the relevance and applications of the charge density normal modes of nanorings have been covered recently [18], we here continue to provide some of the same equations for convenience. We start with the facts that the toroidal harmonics -( ) P z n m 1 2 and -( ) Q z , n m 1 2 and their derivatives with respect to z are symmetric in n, see the appendix in [16]. Therefore, it follows from equation (14) that = -R R n n holds for all =   ¼ n 1, 2, . The symmetric property of R n together with substitutions: 0, 1, ... , n n n n n n transforms the bi-infinite vector three-term recurrence equation (11) into semi-infinite vector three-term recurrences: 1, 2, ..., 21 n n n n 1 1 and Recall that for dispersion relation, our goal is not to seek solutions  { } n and  { } n for equations (20)-(23), but to obtain a relation from which the appropriate set of dielectrics and hence the normal modes of the system can be derived. In the well-known cases of simply connected regions such as plane, sphere, spheroid, cylinder, etc, the dispersion relations are derived from the requirement that the problem must possess non-trivial solution. Here, we do not consider the effects of any nonlinearities or losses in the system that are known to lead to further diversification of the eigenmodes similar to the case of a lossy metallic nanowire [27]. This approach in the more complicated case of a torus, signified by the above equations, fails to yield the desired dispersion relations. To illustrate this idea and as the first natural attempt, one can express equations (20)-(23) as infinite system of linear equations: At this point, one may be tempted to claim that equations (24) and(25) possess non-trivial solutions if the determinants of the corresponding infinite matrices vanish identically, and thus obtaining the dispersion relations. Statements similar to the above paragraph have appeared in [16,28]. Unlike systems with finite dimensions where existence of nontrivial solutions are equivalent to the vanishing of the determinant of coefficient matrix of the system, such requirements are not necessarily valid in the infinite dimensional cases. More details regarding the problem with this approach is addressed in appendix B. As we shall see, the correct approach is to employ the theory of MCFs, which is closely related to the vector three-term recurrence relations. First, note that equations (20)-(23) can be written as a single vector three-term recurrence: with two initial conditions: The usage of the phrase 'initial conditions' is due to the fact that equation (26) can be cast into a first order nonlinear matrix recurrence relation, where equations (27) and(28) serve as initial conditions, see equations (30)-(32).
Next, suppose that the sequence of k×k nonsingular matrices { } X n is a solution to the matrix three-term recurrence: - Assuming  n constitutes a single column of X , n equations (27) and(28) imply the two initial conditions: where I denotes the k×k identity matrix and the quotient of two k×k matrices A and B with B non-singular will be denoted by 1 throughout the rest of the paper. Multiplying equation (29) from the right by -X n 1 and defining V n by = + -V X X , n n n 1 1 we get the following first order nonlinear matrix recurrence: or written more compactly from now on as: Assuming the limit of the right-hand side of equation (33) exists as  ¥ n , one gets Utilization of the initial conditions equations (30) and(31) leads to the formal expressions: It turns out that one can show the following facts.
converges and its (unique) limit is independent of any specific choice of the dielectric values e e e + , , ..., .
II. The convergence of the MCF(37) is equivalent to the existence of the minimal solution to the matrix threeterm recurrence equation (29).
The proof of the above statements is long and non-trivial and uses many classical results in connection between the theory of continued fractions and linear recurrence relations. For the sake of readability, we therefore have moved the proof to appendix C. Now, in light of the above results, it follows that equations (35) and(36) are indeed the dispersion relations for the k-layered torus. Clearly since these equations can not hold simultaneously, we have two separate dispersion relations. From the mathematical point of view, however, there is no prior knowledge to assure either of the initial conditions equation (27) or (28)  obtained in this way are the eigenmodes of the system from which the normal modes can be calculated.
Since there is no exact analytical expression for the MCF equation (37), the exact solutions must be calculated using numerical analysis. In order to examine the numerical structure of the dispersion relations equations (35) and(36), we rewrite them as: and make the following trivial observations: where R R R , , , ...
are defined by equation (14). In attempting to find those e e + { } , ..., k 1 1 which solve equations (38) and (39), we may also obtain extraneous roots. Nevertheless, for numerical calculation, it is preferable to work with equations (38) and (39) instead of equations (35) and(36) due to the fast convergence of determinant in comparison to any matrixnorm method. To illustrate this, consider a sequence of 2×2 matrices { } F n whose asymptotic behavior for large n is given by One can easily verify the zero limit by a few iterations of the determinants, while it will be utterly timeconsuming to verify this fact utilizing any matrix-norm calculation. The situation described in the above example is also generic in our case by identifying F n with the nth approximants (see equation (C.2) of the MCF (37)).
The obtained dispersion relation equations (38) and (39) are numerically validated to be in exact agreement with the limiting case results regarding a single-layered torus. More precisely, we obtain the results of a singlelayered torus (see figure 3 in [18] or [28] and [16]) using the dispersion relations (38) and (39) for a k-layered torus model (here only k = 2 is provided) with the following limiting case justifications: Case Limits Parameter setting Aspect ratio Thus, having established the properties of the analytical and numerical approaches to obtaining the generalized dispersion relations of the quasi-static normal modes, we now aim to address the properties of the surface modes that arise when the system is subject to electromagnetic perturbation.

The Green's function approach
The Green's function approach was employed by Love [17] to solve for the quasi-static response of a dielectric ring to an external uniform field. Unlike the homogenous difference equation (11), describing the normal modes of the multilayer solid torus, the presence of a nonzero external field leads to a source term for the difference equation that ultimately describes the potential distribution of the ring responding to the external field. For a simple (uniform and isotropic) ring, Love obtained a three term difference equation with a source term corresponding to an external uniform axial field (see Love's equation (2.13)). Love's result was generalized by Garapati et al [18] to obtain the response of a composite ring to arbitrary external fields created by discrete multipolar charge distributions. This generalization however warrants further consideration with respect to its convergence properties. Here, extending this generalization to also account for continuous charge distributions, we begin by considering a toroidal charge distribution with density ρ that generates a potential F ap and expand    All other used notations, except for  n given by equation (46), are the same as those already introduced in section 2.
Despite the fact that the explicit form of homogenous solutions for equation (47) is unknown, it is still possible to find an explicit expression for  n with the aid of the Green's functions. As we shall see, the success of Green's method is mostly due to the symmetric property of R n and that the homogenous part of equation (47) is of Poincaré-type, discussed in section 3. Our method is the generalization of the one described in Love [17], see also [31].
We begin with the simple observation via a direct substitution that is the solution of equation (47) To analyze equation (52), we consider the (backward) matrix three-term recurrence: which is convergent. Thus, with a similar argument presented for the case of equation (53), the three-term recurrence equation (52) has a minimal solution. Letting  n N , denote this minimal solution for n N 2, it follows from equation (63) that Next, we rewrite equations (59) and (64) as: Now, for each fixed n, one may use equations (65)-(69) to obtain  n from equation (50). Thus one can finally calculate theḱ 1 vectors  n and  n from equations (12) and (46), respectively.
As a final remark regarding the utility of the Green's function, it is illuminating to consider the response of a coated ring to monopoles and dipoles. Other basic cases were treated recently [18]. We consider a k-layered In case of an uniform field polarized along the z-axis, the generalized form of the applied external potential may be given by (see [17])  Applying the principle of superposition as before we obtain equation (41) but with m=0 considered. As an example we present figure 5 which shows the potential distribution of a coated nanoring in an uniform field with polarization parallel to the symmetry axis of the ring where in the vectors  n and  n are computed as discussed in the paragraph following equation (69). For comparison, the potential distributions are provided for both confocal and concentric rings, as shown in figure 5. Other relevant cases have been studied in [18].
In closing, regarding the experimental investigation of nanorings, currently within nanophotonics and nanoplasmonics (typical length scale ∼ few tens of nm), several (emerging) approaches have been reported for fabrication of these structures [33][34][35][36][37][38][39][40]. The main methods include nanofabrication (electron beam lithography) and chemical synthesis. Single rings, ring arrays, and in-fluid suspended rings (metallic and dielectric) have been attempted. Considering the photon response of the rings, singly and possibly doubly coated metallo-dielectric rings are most conceivable, while it is unlikely that the number of coatings will exceed much beyond a few. Further, we note that typically, the coating thicknesses are only a fraction of the minor radius of the ring. Therefore, for practical purposes, as also seen in figure 1, the deviation from full concentricity is negligible. This is clearly seen also from the comparison of the computationally determined potential distributions for the concentric rings versus those of the confocal rings as shown in figures 2-4.

Conclusion
While in general it is recognized that oscillations in the electronic charge density at toroidal interfaces occur at frequencies that can be approximately obtained from the eigenvalues associated with the quasi-static boundary value problem, the scalar three term difference equation that arises when the scalar electric potential is required to satisfy the boundary conditions across a single toroidal interface is seen to cause significant algebraic complexity. The 'three-term coupling' is due to the quasi-separation of variables in toroidal coordinates while attempting to solve the Laplace equation. This phenomenon is demonstrated in the coupling of the toroidal coordinates μ and η via the function f in equation (2). As a result, the generalization of the scalar three-term difference equations to vector three-term recurrences was shown to be necessary to obtain the plasmon dispersion relations in a multilayered toroidal structure. The structure of the nonretarded surface plasmon dispersion relations for a general multilayer toroidal system may thus be studied following the convergence properties of the MCFs and infinite determinants associated with the vector three-term recurrences. However, contrary to the cases of finite dimensional matrices representing the underlying system of equations, in the toroidal case, as described here, it is reasonable not to simply assume that the vanishing of the matrix determinant would guarantee nontrivial solutions. The convergence properties of the obtained MCFs associated with the vector three-term recurrence relation equations (20)-(23) were shown to provide a numerical platform for computation of the eigenmodes using the determinant forms of the dispersion relations given by equations (38) and (39). The presented results help facilitating numerical analysis of the plasmon dispersion for specific metals such as gold, silver and aluminum assembled, in conjunction with suitable dielectric media such as silicon and quartz, into a complex nanoring. Equations (38) and (39) can be used to obtain the dependence of the plasmon excitation frequency upon the aspect ratio of various nanorings, which enters via the frequency dependence of the dielectric functions of the involved materials. Thus, plasmon dispersion can be investigated in arbitrarily coated toroidal nanoparticles (see also [18]). In addition to the normal mode calculation, the response of a composite ring to arbitrary external fields was shown to be attainable via convergent MCFs.

Acknowledgments
This work was supported in part by the laboratory directed research and development (LDRD) fund at Oak Ridge National Laboratory (ORNL). ORNL is managed by UT-Battelle, LLC, for the US DOE under contract DE-AC05-00OR22725.

Appendix A
In this section, we supply the details in obtaining the vector three-term recurrence equation (11) and give a proof of theorem 1.
First, we apply the boundary condition equation ( Similarly, applying the boundary condition equation (7) for each = ¼ i k 1, , , we obtain å å e e ¢ + ¢ + ¢ + ¢   The k×k bidiagonal matrices  ¢ n and  ¢ n are defined by where J n is given by equation (13) as n n n n n n 2 1 1 1 Using equation (A.11), one can easily see n n n n n n n n n n 1 2 1 1 1 It follows now that the substitution  This concludes the algebraic derivation of equation (11).
Proof of theorem 1.  2), may not exist but the system still possesses non-trivial solutions, see [18]. On the other hand, one can show the validity of the above statement under certain rather restrictive conditions. One such classic result (see Koch [45] and also [41][42][43]) states that for a certain class  of infinite matrices if the moduli of the off diagonal terms are square summable and the infinite product of the diagonal terms converges absolutely, and if the vector  is also square summable, then ( ) A det exists and the infinite system equation (B.1) has a non-trivial solution if and only if = ( ) A det 0, see Denk [42] for details. As it stands, neither of the infinite block matrices in equations (24) and (25) belongs to the class . It should also be emphasized that we can not utilize the truncated system method, see [46], where the infinite matrix A has to satisfy similar conditions to those of the class . Finally, we note that a formal calculation of the determinants for infinite matrices in equations (24) and(25) based on the definition equation (B.2) and the result regarding the determinant of a finite tridigonal matrix [47] does not imply the obtained dispersion relations. This justification, see [48], is beyond the scope of this paper and is therefore omitted.

Appendix C
Here, we prove the claims (I) and (II) of section 3 with regard to the MCF(37) and its relation to the matrix three-term recurrence equation (29).
We start by assessing the convergence of equation ( Matrix F n is called the nth approximant of the MCF equation (C.1).
The classical approach to continued fractions with the aid of the theory of Möbius or linear fractional transformations is also fruitful in treating MCFs, see Ahlbrandt [24,25] and references therein. In this paper, however, we only consider and thus explain the portion which is relevant to our special case equation (C.1). So let M be a nonsingularḱ where , ,  and  are k×k matrices. The matrix linear fractional transformation   T : Note that we will suppress the symbol M from our notation whenever it is understood. It is easily seen that Before we can show the connection between equations (C. 16) and (C.14), we need some definitions. Two solutions X n and Y n of equation (29) are said to be linearly independent if their Wronskian (or Casoratian) is nonsingular; that is, theḱ is nonsingular (see [24] for details).

Definition 1.
A solution X n of the matrix three-term recurrence equation (29) is called minimal if the following conditions hold. For large values of n, Y n is nonsingular. Therefore, the last equality gives  = + This completes the proof. , We should mention that theorem 2 is a generalization of the classic Pincherle's theorem [22], and the proof provided here essentially follows the one given in [26], see theorem 2.1. In fact the converse of theorem 2 also holds. More precisely, if the MCF equation (C.16) converges, then there exists a solution of equation (C.14), say Z n as in equation (C. 18), such that both Z 0 and X 0 are nonsingular, Y n is nonsingular for large n, and = ¥ lim 0. n X Y n n Such a Z n is called a fundamental system of solutions of equation (C.14). In our case, it is easily seen that the definition of minimality and fundamental system of solutions are equivalent. Therefore, one can state the following extension of theorem 2 (see also [26], theorem 2.1). At this point, we have only been able to show that the convergence of the continued fraction equation (C.16) is equivalent to the existence of a minimal solution for equation (29). The next theorem is a combination of original results given by Perron [51], Máté and Nevai [52], and completes this picture. For the full version of this theorem see [26]. We point out that the recurrence equation (C.21) with the property q q = ¥ lim n n is said to be of Poincarétype and the corresponding MCF is called limit periodic. In view of theorem 4, we can give our last theorem of this section, which concludes what we ought to show.