STRATIFY: a comprehensive and versatile MATLAB code for a multilayered sphere

We present a computer code for calculating near- and far-field electromagnetic properties of multilayered spheres. STRATIFY is one-of-a-kind open-source package which allows for the efficient calculation of electromagnetic near-field, energy density, total electromagnetic energy, radiative and non-radiative decay rates of a dipole emitter located in any (non-absorbing) shell (including a host medium), and fundamental cross-sections of a multilayered sphere, all within a single program. The developed software is typically more than 50 times faster than freely available packages based on boundary-element-method. Because of its speed and broad applicability, our package is a valuable tool for analysis of numerous light scattering problems, including, but not limited to fluorescence enhancement, upconversion, downconversion, second harmonic generation, surface enhanced Raman spectroscopy. The software is available for download from GitLab https://gitlab.com/iliarasskazov/stratify

To date, such matryoshkas with various combinations of dielectric and metal shells have wellestablished protocols for the efficient and controllable synthesis [6,[53][54][55][56][57], greatly improved on initial synthesis attempts [3][4][5].Theoretical and numerical studies in these fields are inevitable for a thoughtful design of experiment and reliable interpretation of measured data.The theory of electromagnetic light scattering from multilayered spheres has a long history since the pioneering work of Aden and Kerker for two concentric spheres [58].An impressive number of closed-form solutions and discussions on coated [59,60] and general multilayered [61] spheres have been reported.There is a thorough understanding of scattering [62] and absorption [63] of light from multilayered spheres, including aspects for various types of illuminations [64][65][66][67], spontaneous decay rates of a dipole emitter in a presence of a multilayered sphere [68,69], electromagnetic energy [70] and near-field [71] distribution within and near the multilayered spheres, and the respective strategies for numerical implementation of the developed theories [72][73][74][75].
For a single shell, a number of features can be qualitatively understood from the quasi-static analysis [2], a full, quantitative understanding of the mechanisms underlying the applications of multilayered spheres is usually quite involved, requiring sophisticated and complex theoretical and numerical studies.The latter are often handled with commercial brute-force finite-difference timedomain (FDTD), finite-element-method (FEM) or boundary-element-method (BEM) solvers, which provide accurate results at relatively high computational price.However, spherical multilayered particles are computationally feasible per se due to the existence of the closed-form analytical solutions inherently convenient for computational implementation.Nonetheless, there are few freely available, user-friendly and comprehensive codes for light scattering from a general arXiv:2006.06512v1[physics.optics]11 Jun 2020 multilayered spheres based on these solutions [76].In most of the cases, a very limited number of properties are calculated within a single freely available package.Fundamental cross-sections (i.e.scattering, absorption and extinction) and/or near fields are the usual choice implemented in freely available codes [77,78], and in a number of so-called "online Mie calculators", which commonly handle only homogeneous spheres with a rare exceptions for core-shells and even more rarely for general multilayered spheres.Computation of the orientation-averaged electric or magnetic field intensities (i.e.averaged over a spherical surface of a given fixed radius r), spontaneous decay rates and other quantities are almost unavailable, yet generally very useful.Here we present free MATLAB code based on compact and easy-to-implement transfer-matrix formalism, which can be conveniently used to handle most of the problems of light scattering from multilayered spheres by using appropriate combinations and manipulations of transfer-matrices.The use of closed-form solutions significantly enhances the performance of the code compared with brute-force solvers without sacrificing the accuracy.
Our recursive transfer-matrix method (RTMM) has been inspired by the success of such a RTMM for planar stratified media developed by Abelès [79][80][81], summarized in Born & Wolf classical textbook [82,Sec.1.6],and can be seen as its analogue for spherical interfaces.From a historical perspective, the RTMM presented here had been developed as early as in 1998 and implemented in the F77 sphere.fcode [76].It was first applied to a number of different theoretical settings with external plane wave source [10,12,14], and successfully tested against experiment [6,13].Later it was extended to incorporate a dipole source [68,69], energy calculations [70], and has enabled exhaustive optimization of plasmon-enhanced fluorescence [30].
The paper is organized as follows.In Sec. 2, we provide an overview of the RTMM for multilayered spheres and formulate the fundamental properties of spheres.In Sec. 3, we describe the developed code and benchmark it with robust BEM [83] implemented in freely available "MNPBEM" MATLAB package [84][85][86][87].In Sec. 4, we discuss possible application of the code in SERS, plasmon-enhanced fluorescence, upconversion, and end with conclusive remarks and propose further developments of the package.Gaussian units are used throughout the paper.

Recursive Transfer-Matrix method
Consider a multilayered sphere with N concentric shells as shown in Fig. 1.The sphere core counts as a shell with number n = 1 and the host medium is the n = N + 1 shell.Occasionally, the host medium will be denoted as n = h.Each shell has the outer radius r n and is assumed to be homogeneous and isotropic with scalar permittivity ε n and permeability µ n .Respective refractive indices are η n = √ ε n µ n .We assume that the multilayered sphere is illuminated with a harmonic electromagnetic wave (either a plane wave or a dipole source) having vacuum wavelength λ.Corresponding wave vector in the n-th shell is k n = η n ω/c = 2πη n /λ, where c is the speed of light in vacuum, and ω is frequency.Electromagnetic fields in any shell are described by the stationary macroscopic Maxwell's equations (with time dependence e −iωt assumed and suppressed throughout the paper): where the permittivity ε and permeability µ are scalars.Following the notation of [68], the basis of normalized (the normalization here refers to angular integration) transverse vector multipole fields ∇ • F pL ≡ 0 that satisfy the vector Helmholtz equation for n-th shell, 1 ≤ n ≤ N + 1, can be formed as: where with fEL = f M L and fM L = f E L .Here L = , m is a composite angular momentum index and f pL is a suitable linear combination of spherical Bessel functions.Provided that the multipole fields in (3) represent E, the respective subscripts p = M and p = E denote the magnetic, or transverse electric (TE), and electric, or transverse magnetic (TM), polarizations [88].The so-called magnetic, Y (m) L , longitudinal, Y (o) L , and electric, Y (e) L , vector spherical harmonics of degree and order m, are defined in spherical (ϕ, ϑ, r) coordinates as [77,[88][89][90]] where êϕ , êϑ , and êr are corresponding unit vectors, and P m (x) are the associated Legendre functions of the first kind [91] of degree and order m: . General solution for the electric field in the n-th shell, 1 ≤ n ≤ N + 1, is [68, Eq. ( 8)]: Fig. 1.A sketch of the multilayered sphere in a homogeneous isotropic host medium with permittivity ε h = ε N +1 and permeability µ h = µ N +1 .Center of the sphere is located at the origin of the reference system.
In order to emphasize and for the sake of tracking the Bessel function dependence, the multipoles F pL in Eq. ( 6) are denoted as J pL and H pL for the respective cases that f p = j and f p = h (1) , where j and h (1) being the spherical Bessel functions of the first and third kind, correspondingly.Given the general solution ( 6), the corresponding expansion of magnetic field H follows from that of the electric field E by the stationary macroscopic Maxwell's equations (1) on using relations (4) [68, Eqs. ( 9)-( 10)]: The H pL 's in Eqs. ( 6)-( 7) with angular-momentum index L refer to a basis multipole which distinguishes them from the magnetic field H p .Very much as in the case of planar stratified media [79][80][81][82], our RTMM for spherical interfaces exploits to the maximum the property that once the coefficients A pL (n + 1) and B pL (n + 1) on one side of a n-th shell interface are known, the coefficients A pL (n) and B pL (n) on the other side of the shell interface can be unambiguously determined, and vice versa.Schematically, where T ± p (n) are 2 × 2 transfer matrices [68].The matrices, which are the fundamental quantities of our approach, can be viewed as a kind of raising and lowering ladder operators of quantum mechanics.In the present case, they lower and raise the argument n of expansion coefficients A pL (n) and B pL (n).The transfer matrices are interrelated by and are unambiguously determined by matching fields across the n-th shell interface, i.e. requiring that the tangential components of E and H are continuous.One finds [68]: where ψ (x) = x j (x) and ζ (x) = xh (1) (x) are the Riccati-Bessel functions, prime denotes the derivative with respect to the argument in parentheses, x n = k n r n and xn = x n / ηn = k n+1 r n are the internal and external dimensionless size parameters, ηn = η n /η n+1 and μn = µ n /µ n+1 are relative refractive indices and permeabilities.For the sake of clarity, the n-subscript has been suppressed on the rhs of Eqs. ( 10) - (13).The above expressions for T ± p (n) are general, independent on the incident field, and valid for any homogeneous and isotropic medium, including gain ( (ε n ) < 0) or magnetic materials (µ n 1).
The formalism becomes compact upon the use of the composite transfer matrices T p (n) and M p (n) defined as ordered (from the left to the right) products of the constituent raising and lowering 2 × 2 matrices from Eqs. ( 10)-( 13): Composite matrices T p (n) and M p (n) transfer expansion coefficients to the n-th shell from the sphere core or from the surrounding medium, respectively.Note that T p (n) are defined for 2 ≤ n ≤ N + 1, while M p (n) are defined for 1 ≤ n ≤ N.They are used in our formalism to chiefly relate the expansion coefficients in the core to those in a source region.For example, T p (N + 1) enables one to obviate all the intermediary shell interfaces and to relate the expansion coefficients in the core directly to those in the surrounding host medium, .
The latter formally reduces the problem of a multilayered sphere with a source in the host medium to that of a homogeneous sphere.Analogously to (9): The reader familiar with the RTMM for planar stratified media (e.g.multilayer coatings and interference filters) will recognize in the relations (8), ( 9), ( 14), ( 15), ( 16) familiar properties of the constituent and composite transfer matrices for planar interfaces [79][80][81][82].
In order to unambiguously determine the expansion coefficients A pL (n) and B pL (n) in each shell, one has to impose two boundary conditions, which is the subject of the following section.

Boundary conditions
The total number of different sets of expansion coefficients comprising all j's from the interval 1 ≤ j ≤ N + 1 is larger by two than the number of corresponding equations.Therefore, boundary conditions have to be imposed to unambiguously determine the expansion coefficients at any shell.They are 1.The regularity condition of the solution at the sphere origin, which eliminates h (1) (0) → ∞ for f p in Eq. ( 3): 2. For a source located outside a sphere, the A pL (N + 1) coefficients at any given frequency ω are equal to the expansion coefficients of an incident electromagnetic field in spherical coordinates.
The regularity condition (17) alone suffices to unambiguously determine the m-independent (due to spherical symmetry of a problem) ratio [88]: which defines the familiar T-matrix of scattering theory, T p , where i, j in T i j;p (n) label the (i, j)-th element of the 2 × 2 matrix T p (n). Corresponding closed-form analytic expressions for applying the second boundary condition are presented in (i) Refs.[90], [70, Eq. ( 10)] for the plane electromagnetic wave, and in (ii) Ref. [68,Eqs. (44), (50)] for the electric dipole source.In the case of an elementary dipole radiating inside a sphere, there is no source outside a sphere, and the second boundary condition reduces to A pL (N + 1) = 0 [68].

Far-field properties
Fundamental cross sections σ (scattering, absorption and extinction) are determined as an infinite sum over polarizations (p = E, M) and all partial -waves: where takes the real part, and T p is found from Eq. ( 18).On using polarized scattering waves parallel and perpendicular to the scattering plane, one can easily build up scattering matrix, find Stokes parameters [77, Eq. (4.77)], and obtain angle-resolved scattering pattern.

Near-field properties
For a given illumination, and, thus, known expansion coefficients A pL (n) and B pL (n), the electromagnetic field can be unambiguously defined by Eqs. ( 6) and (7).Below we consider more sophisticated near-field properties.

Electromagnetic energy
Total electromagnetic energy W stored within a multilayered sphere is a sum of energies, W n , stored within any given n-th shell, which in turn is an integral of an electromagnetic energy density w n (r) over the n-th shell: where r 0 = 0, i.e. the core "inner radius".Here we limit the discussion to a plane wave excitation.For non-dispersive (usually dielectric) shells [92], electric and magnetic coefficients in (21) are while for the dispersive and absorbing shells [93]: where γ D is the free electron damping constant in the Drude formula, and takes the imaginary part.In our code, we use Drude fit parameters from either Blaber et al. [94] or Ordal et al. [95].
Electric and magnetic components from Eq. ( 21) are nothing but the orientation-averaged electric and magnetic field intensities: Here fp , ±1 = Āp (n) j , ±1 (k n r) + Bp (n)h (1)  , ±1 (k n r).Note that in all cases the spherical Bessel functions of either and ± 1 orders are multiplied by the expansions coefficients Āp and Bp with the index .These expansion coefficients are determined via raising T p (n) and lowering M p (n) composite transfer matrices as Because the surface integrals of electric and magnetic field intensities are performed analytically [70], the calculation of average intensity costs the same computational time as determining intensity at a single given point.The radial integrations in ( 21) are performed by using Lommel's integration formulas [70]: Here x = k n r, and purely imaginary functions Fp = 2i x fp +1 (x) f * p (x) are cancelled by purely imaginary x 2 − x * 2 = 4i (x) (x) in the denominator, which results in purely real integrals in (27).For a special case of lossless shells, the denominator x 2 − x * 2 vanishes.It can be eliminated by using l'Hôpital's rule [70] to get the following amendments in Eq. ( 27): . Substitution of ( 24) and ( 27) into ( 21) yields in explicit expressions for the total electromagnetic energy W n stored within each shell of the multilayered sphere and for the electromagnetic energy density w n (r).These relations are general and valid for any shell, including the (N + 1)-th shell being a surrounding medium.For the comprehensive derivation of the equations above, we refer the Reader to the Ref. [70], which generalizes the pioneering work by Bott and Zdunkowski on homogeneous spheres [92] and subsequent follow-ups for magnetic spheres [96] and core-shells [97].

Spontaneous decay rates
Radiative and nonradiative decay rates (normalized with respect to Γ rad;0 , the intrinsic radiative decay rate in the absence of a multilayered sphere) for a dipole emitter located in n d -th shell at r d distance from a center of a sphere (see Fig. 1) are given by [68]: where Coefficients N rad and N nrad depend on whether the decay rates were normalized with respect to the radiative decay rates in infinite homogeneous medium having the refractive index of (i) the host or (ii) having the refractive index of the shell where the dipole emitter is located: Functions F p (x d ) and D p ;a (x d ) depend on the relative position of the emitter with respect to the sphere, and to a n a -th absorbing shell: The summation over nonradiative decay channels in Eqs. ( 28) includes each absorbing shell, i.e. shell with (ε a ) > 0. Respective volume integrals I E ;a and I M ;a are 2 dr , Here k a = 2πη a /λ is the wavevector in the absorbing shell, and coefficients A p ;a and B p ;a are , , n a 1 and Note that A p ;a and B p ;a depend only on the relative position of the dipole with respect to the absorbing shell.We make use of this property to optimize calculations, namely, we get I E ;a and I M ;a once for a set of r d , considering only the relative position of the dipole emitter.
For the detailed discussion of the spontaneous decay rates of a dipole emitter in a presence of a general multilayered sphere, we refer the reader to Ref. [68].

Convergence criteria
Numerical implementation of equations above, which involve infinite summation over , requires truncation at some finite number max .For far-field properties, Eqs. ( 19) and (20), the classic choice is the Wiscombe criterion [98] where x N = k h r N is the usual size parameter.Keeping in mind plasmonics and nanophotonics applications as intended for our package, and recalling finite precision of the MATLAB, we limit the discussion for multilayered spheres with x N < 50 requiring max 70, which is more than sufficient.
For the near-field, Eqs. ( 6), (7), and for the electromagnetic energy, Eqs. ( 21), ( 24), ( 27), slightly different cut-off is recommended [99]: For the dipole source implied in Eq. ( 28), especially when considering non-radiative decay rates [68,100], the accuracy has to be set instead of max , since there are no general guidelines for choosing the latter in this case.

Electron free path correction
For thin metallic shell with thickness (r n −r n−1 ) less than a free electron path, the bulk permittivity ε n,bulk has to be corrected to take into account electron scattering from the shell surface [101]: where υ F is the Fermi velocity and ω p is plasma frequency.This modification of thin metallic shells permittivity may have crucial consequences for their electromagnetic properties [20,21,[102][103][104].

Overview
Fundamental properties are calculated with the following functions: • t_mat.mcalculates transfer matrices with Eqs. ( 10) -( 13) and returns their ordered products (14).Any other function which returns electromagnetic properties (decay rates, electromagnetic energy and fields, far-field properties and etc) makes use of pre-calculated ordered products of transfer matrices for a better performance; • decay.mreturns decay rates calculated with Eqs.(28).Decay rates are normalized (see Eqs. ( 29)) with respect to radiative decay rates of a dipole in a homogeneous medium (without a multilayered sphere) with the refractive index of a shell where the dipole is embedded, η d , or with the refractive index of a host, η h ; • nrg_dens.mreturns normalized by the factor π|E 0 | 2 ε h [70, Sec.3E] electric and magnetic components of electromagnetic energy density, w n (r), given by Eq. ( 21), and orientation-averaged intensities of electric and magnetic fields given by Eq. ( 24); • nrg_tot.m returns normalized by the factor (2/3)π|E 0 | 2 (r n − r 3 n−1 )ε h [70, Sec.3E] electric and magnetic components of total electromagnetic energy W n stored within n-th shell given by Eq. ( 21); • near_fld.mreturns electric, E(r) = E E (r) + E M (r), and magnetic, H(r) = H E (r) + H M (r) near-field distributions given by Eqs. ( 6) and (7).We take the advantage of the spherical symmetry of the problem and pre-calculate computationally expensive r-dependent Bessel functions and cos θ-dependent associated Legendre polynomials only for unique values of r and cos θ.Usually, the amount of these unique r and cos θ is significantly smaller than the respective total number of points in the rectangular mesh, which results in faster execution; • far_fld.m returns polarized scattering waves from Eq. ( 20); • crs_sec.m returns fundamental cross sections from Eq. ( 19); • el_fr_pth.mreturns corrected permittivity according to Eq. (37).
Suitable combinations or minor post-processing of these fundamental characteristics may be used for almost every known application of multilayered spheres briefly discussed in Sec. 4 below.

Verification and performance
We have compared STRATIFY with freely available exact BEM-based [83] solver MNPBEM [84][85][86][87].The latter package has been chosen for a benchmark since it is a freely available comprehensive MATLAB code capable of calculating most of the quantities considered in our code.For testing purposes, we have considered fundamental cross-sections, electric near-field distribution, and spontaneous decay rates of dipole emitter in the presence of matryoshkas composed of different combinations of Au and SiO 2 layers.It can be easily seen from Fig. 2 that our code produces the same results as the exact BEM method, but for sufficiently lower computational price.

Discussion and Conclusions
The developed package is ready to be applied in a number of well-established applications of multilayered spheres in optics and photonics.Below we discuss the most common examples.Due to enormous electric field enhancement in multilayered metal-dielectric nanospheres, they are considered as good candidates for a number of applications.Squared electric field intensities from Eq. (24) or Eq. ( 27) are measure for performance of nanostructures in SERS [25,41,42,[106][107][108] and second harmonic generation [8,9].For fluorescence or upconversion enhancement, where the spontaneous decay rates of the dipole emitters are modified, the generic enhancement factor is nothing but a product of excitation rate enhancement (at excitation wavelength, λ exc ) and quantum yield (at emission wavelength, λ ems ): Here γ exc /γ exc;0 ∝ |E| 2P /|E 0 | 2P , P is the number of photons involved in the excitation process (P = 1 for fluorescence enhancement, P = 2 for the upconversion), q ems;0 is the intrinsic quantum yield of the emitter, and Γ rad,nrad = (2Γ rad,nrad + Γ ⊥ rad,nrad )/3 is an orientationally averaged decay rate determined at a fixed dipole radial position by averaging over all possible orientations of a dipole emitter.Orientationnaly averaged fields in γ exc , and orientationally averaged decay rates, Γ rad,nrad , are reasonable choice for reliable estimates, unless, of course, the emitter is positioned with a controllable orientation of its electric dipole moment at a particular (e.g. a hot spot) location.
Cloaking, Kerker effect [109], super-or optimally tuned scattering [51,[110][111][112][113][114][115] and absorption [20,51,113,[115][116][117][118], embedded photonic eigenvalues [119], spasing [44,47,120,121] and other intriguing phenomena [122] are easily understood from fundamental cross sections (19) and scattering patterns (20).We have summarized a self-consistent and comprehensive RTMM theory reported earlier in our [68,70] for electromagnetic light scattering from general multilayered spheres composed of isotropic shells.Within the framework of RTMM, we have developed an efficient multi-purpose MATLAB package for calculating fundamental properties of multilayered spheres.Our package is one-of-a-kind freely available software which allows for a simultaneous calculation of a wide range of electromagnetic properties and is ready to be used for a broad number of applications in chemistry, optics and photonics, including optimization problems and machine-learning studies.We hope that the generalization presented here and corresponding MATLAB code will serve as a useful tool for photonics, physics, chemistry and other scientific communities and will boost the researches involving various kinds of multilayered spheres.
Note that T p are nothing but familiar expansion coefficients a and b in Bohren and Huffman's representation [77, Eqs.(4.56),(4.57)],but with the opposite sign: T E = −a and T M = −b .