A method of effective potentials for calculating the frequency spectrum of eccentrically layered spherical cavity resonators
aa r X i v : . [ c ond - m a t . d i s - nn ] M a y A method of effective potentials for calculating the frequencyspectrum of eccentrically layered spherical cavity resonators
Z. E. Eremenko, Yu. V. Tarasov, and I. N. Volovichev
A. Ya. Usikov Institute for Radiophysics and Electronics, NAS of Ukraine,12 Proskura Str., 61085 Kharkov, Ukraine
ARTICLE HISTORY
Compiled May 12, 2020
ABSTRACT
A novel method for the calculation of eigenfrequencies of non-uniformly filled spher-ical cavity resonators is developed. The impact of the system symmetry on theelectromagnetic field distribution as well as on its degrees of freedom (the set ofresonant modes) is examined. It is shown that in the case of angularly symmetriccavity, regardless of its radial non-uniformity, the set of resonator modes is, as antic-ipated, a superposition of TE and TM oscillations which can be described in termsof a single scalar function independently of each other. The spectrum is basicallydetermined through the introduction of effective “dynamic” potentials which en-code the infill inhomogeneity. The violation of polar symmetry in the infill dielectricproperties, the azimuthal symmetry being simultaneously preserved, suppresses allazimuthally non-uniform modes of electric-type (TM) oscillations. In the absence ofangular symmetry of both electric and magnetic properties of the resonator infill,only azimuthally uniform distribution of both TM and TE fields is expected to oc-cur in the resonator. The comparison is made of the results obtained through theproposed method and of the test problem solution obtained with use of commercialsolvers. The method appears to be efficient for computational complex algorithmsfor solving spectral problems, including those for studying the chaotic properties ofelectrodynamic systems’ spectra.
KEYWORDS
Spherical cavity; layered resonator; Debye potentials; mode decoupling; intermodescattering; spectral chaos
1. Introduction
Isolated wave systems of macroscopic and mesoscopic dimensions, in particular, quasi-optical electromagnetic (EM) resonators arouse large interest of many researchersdue to their numerous applications and fundamental physical properties. This re-lates, specifically, to metallic cavities filled with an anisotropic or gyroelectric medium[1–3], metamaterial cavities [4], oscillators and filters [5], the measurement of the fer-rite materials’ dielectric tensor [6,7], etc. EM cavity resonators are also used quiteextensively to simulate quantum mechanical systems of finite dimension, both openand closed. The main tool to examine their chaotic properties, just the same as forclassical dynamic systems, is the statistical analysis of their spectra [8]. Most of theconclusions of the statistical theory regarding the dynamical system spectra settle on
CONTACT Yu. V. Tarasov. Email: [email protected] heir symmetry properties [9,10], since direct calculation of the spectral parametersis normally difficult to perform. Similar methods apply to wave systems as well, andthe conclusions about their chaotic properties are mostly made on the basis of ray(“trajectory”) optics [11]. Such an approach, in view of significant uncertainty in theposition of indicative points on the ray trajectories, if taking account of finiteness ofthe wave length, cannot provide sufficient reliability of the results, especially whenit comes to closed systems such as cavity resonators. Besides, to obtain the reliablespectra of resonators with some inhomogeneities in the bulk it appears to be, as a rule,computationally expensive task, if one considers the number of (Maxwell’s) equationssubject to the solution.The latter circumstance has previously been overcome for homogeneous systemsthrough their examination in terms of Hertz functions [12] and/or Debye potentials[13] rather than directly through EM field components. The opportunity to express allvector-field components through two scalar functions was noted back in a long timepublication by Whittaker [14] (more detailed reference list can be found in paper byNisbet [15]). The approach has proven to be quite productive, and not for uniformsystems only but also for some systems with gyrotropic properties [16–18]. Yet, for along time it was not possible to extend the approach to arbitrary nonhomogeneoussystems, which has substantially reduced the capability of the method. One of the fewsuccessful attempts to overcome this issue was made in Ref. [19], where the method forsolving the Helmholtz equation in heterogeneous closed waveguides was developed byreducing it to the Hamiltonian form through the introduction of four scalar potentials.In our study we propose a novel method to calculate the spectra of inhomogeneouscavity resonators by introducing a couple of Debye-type scalar effective potentials ,with the subsequent solution of the emergent perturbed Helmholtz equations withinthe whole resonator rather than in its partial regions. The potentials originate fromtwo Hertz vectors through which all the EM field components are expressed in thesystem of any geometry [12]. The scalar nature of these potentials relates to the factthat in the spherical and some other coordinate systems one component only of eachof the two Hertz vectors suffices to represent through them the entire set of vectorfield components [13].To solve the Helmholtz equations subject to the perturbation potentials we applythe method of oscillation mode separation applicable for both homo- and heteroge-neous systems [20], which is a strict procedure, in contrast to the approximate sep-aration normally performed within the framework of the conventional coupled-modetheory [21]. The method we suggest for mode decomposition involves the solution ofmatrix equations either to directly find their eigenvalues or to calculate the Greenfunction containing the full information regarding the system spectrum. In both cases,it is far from easy to find the practical analytical solution, so at some stage it is nec-essary to apply also some numerical methods. For this, in its turn, it is necessaryto deliberate some crucial issues. In the theory set forth below, the expansion of theDebye-type scalar potentials is carried out over an infinite orthonormal sets of ba-sis functions. Obviously, in computer calculations one has to limit oneself to a finitenumber of those functions. To do this, it is necessary, for one thing, to determinethe sensitivity of the spectrum to the dimension of the reduced basis set. Secondly,there appears a need to analyze the accuracy and the stability of the computationalscheme when finding the matrix elements of the effective potentials. And finally, thecalculation of the spectrum in a sufficiently large frequency interval implies the con-sumption of quite a lot of processor time. So, it is important to choose the efficientcomputational algorithm, in particular, for numerical integration in 3D. In our paper,2he speculations and the results on this subjects are presented as well.In the present paper we expound the general method for calculating the spectra ofspherical cavity resonators and the results that relate to the specific kind of inhomo-geneity, namely, the inhomogeneity in the form of concentric sharp-bounded layers.We also discuss the effect of the resonator infill symmetry on the number of degreesof freedom of the EM field, and the possibility to represent the resonator eigenmodesin the form of TE and TM components.
2. Statement of the problem and derivation of basic equations
The primary goal of our study is to determine the degree of chaos in the spectrum ofa spherical EM cavity resonator of radius R out filled with a piecewise homogeneousdielectric substance, as shown in Figure 1. The permittivity of the outer dielectriclayer is ε out whereas the permittivity of the inclusion, which we assume to be also Figure 1.:
The eccentrically layered spherical cavity resonator. the dielectric sphere yet of radius R in and the center shifted to arbitrary distance d from the center of outer sphere, is ε in . The permeability of both parts of the resonatorinfill is assumed to be the same, and we put it equal to unity everywhere in the finalformulae.Normally, to find the EM fields in a system with heterogeneously distributed electricand magnetic parameters one has to solve the set of coupled Maxwell equations, whichin the stationary case is represented by the following eight scalar equalities,rot E = ik µ ( r ) H , (1a)rot H = − ik ε ( r ) E , (1b)div (cid:2) ε ( r ) E (cid:3) = 0 , (1c)div (cid:2) µ ( r ) H (cid:3) = 0 . (1d)Here, k = ω/c , ε ( r ) and µ ( r ) are the permittivity and the permeability of the resonatorinfill, both of which at this stage are assumed to be quite arbitrary functions of thecoordinate vector r . The number of equations to be solved may be reduced if one goesover from the vector-valued electric and magnetic fields to the scalar-valued Hertzfunctions. Until now, these functions have been consistently introduced for entirely3omogeneous systems [12,13] and for gyrotropic systems with distinguished axis [16–18]. In our study we show that the analogous Hertzian functions can be introducedfor arbitrarily heterogeneous systems as well, and without any special restrictions onthe degree and the character of inhomogeneity.Following Refs. [12,13], we represent the fields E and H as the sums of the compo-nents of “electric” ( e ) and “magnetic” ( m ) type, E = E e + E m , (2a) H = H e + H m , (2b)both satisfying the set of equations (1) independently. These components, in view ofthe particular symmetry of the set (1), may, in principle, be defined through differ-ent, yet symmetrical, dependencies on the conventional vector potential A ( r ) and thescalar potential ϕ ( r ). In our study we prefer a different definition, namely, troughthe vanishing of the radial components of magnetic field H e and electric field E m . Inaccordance with this definition, the fields marked by indexes e and m will be termedthe TM and TE polarized fields, respectively. In the uniform and unbounded mediumthese fields are uncoupled by definition. In the cavity resonator, however, the boundaryconditions (BCs) couple them, so it may seem not to make much sense to distinguishthe fields by their polarization. Yet, below it will be shown that such a distinctionis quite reasonable since there is the possibility to satisfy the BCs for TM and TEpolarized fields independently of each other.The practicality of the above definition comes to be appreciable if one writes downequations (1) in spherical coordinates associated with the center of outer sphere, seeFigure 1. In this curvilinear system equations (1a) and (1b) take the form ∂E r ∂ϕ − ∂∂r (cid:0) r sin ϑ E ϕ (cid:1) = ik µ ( r ) r sin ϑH ϑ , (3a) ∂∂r (cid:0) rE ϑ (cid:1) − ∂E r ∂ϑ = ik µ ( r ) rH ϕ , (3b) ∂∂ϑ (cid:0) r sin ϑ E ϕ (cid:1) − ∂∂ϕ (cid:0) rE ϑ (cid:1) = ik µ ( r ) r sin ϑ H r , (3c) ∂H r ∂ϕ − ∂∂r (cid:0) r sin ϑ H ϕ (cid:1) = − ik ε ( r ) r sin ϑE ϑ , (3d) ∂∂r (cid:0) rH ϑ (cid:1) − ∂H r ∂ϑ = − ik ε ( r ) rE ϕ , (3e) ∂∂ϑ (cid:0) r sin ϑ H ϕ (cid:1) − ∂∂ϕ (cid:0) rH ϑ (cid:1) = − ik ε ( r ) r sin ϑ E r . (3f)It looks quite appreciable to consider this set of equations separately by choosingeither H r ≡ E r ≡
0. In such a way we would split a general solution into thesum of solutions of two types which are conventionally designated as the TM and TEpolarized fields for a spherical resonator [13].Instead of being written in terms of double six EM field components marked by in-dices ( e ) and ( m ), the set of equations (3) may be rewritten in terms of two scalar fieldsonly, U ( r ) and V ( r ), which are referred to as the electric and magnetic Hertz func-tions, respectively. For e -type fields in Eqs. (2) we will use the substitution formulas4 eϕ = 1 ε ( r ) r sin ϑ · ∂ U∂ϕ∂r , E eϑ = 1 ε ( r ) r · ∂ U∂ϑ∂r ,E er = − ε ( r ) r sin ϑ (cid:20) ∂∂ϑ (cid:18) sin ϑ ∂U∂ϑ (cid:19) + 1sin ϑ ∂ U∂ϕ (cid:21) ,H eϕ = ikr ∂U∂ϑ , H eϑ = − ikr sin ϑ ∂U∂ϕ , H er ≡ , (4a)which define the transverse-magnetic (TM) part of the entire EM field, while for thefield components of m -type (TE field) the substitutions will be applied H mϕ = 1 µ ( r ) r sin ϑ · ∂ V∂ϕ∂r , H mϑ = 1 µ ( r ) r · ∂ V∂ϑ∂r ,H mr = − µ ( r ) r sin ϑ (cid:20) ∂∂ϑ (cid:18) sin ϑ ∂V∂ϑ (cid:19) + 1sin ϑ ∂ V∂ϕ (cid:21) ,E mϕ = − ikr ∂V∂ϑ , E mϑ = ikr sin ϑ ∂V∂ϕ , E mr ≡ . (4b)The motivation for both of these kinds of field representation is mostly taken fromRef. [13], where in the case of a uniform medium, with ε = const and µ = const,substitution of (4a) and (4b) into Maxwell equations (3) has resulted in the identicalHelmholtz equations for Debye potentials (DPs) u ( r ) = U ( r ) /r and v ( r ) = V ( r ) /r ,namely, (cid:16) ∆ + k εµ (cid:17) (cid:18) uv (cid:19) = 0 . (5)The nonuniformity of the medium, if any, may be, in principle, included by putting intothe right-hand-sides of equations (5) the “source” terms in the form of either electricor magnetic polarization [15]. In the present study, however, we suggest a different wayfor considering the medium heterogeneity, namely, by merely replacing the constant ε and µ with variable functions ε ( r ) and µ ( r ), as shown in Eqs. (4a). Certainly, we thusarrive at equations for U and V that differ from simple equations (5) by additionalterms accounting for medium heterogeneity. Yet, regarding these terms as some kindof “dynamic” potentials we can take them into account as a perturbation for thestandard Helmholtz operator. 5ith substitutions (4a), Maxwell equations (3a) and (3b) acquire the form ∂∂r (cid:20) ε ( r ) ∂ U∂ϕ∂r (cid:21) + 1 r ∂∂ϕ (cid:26) ε ( r ) (cid:20) ϑ ∂∂ϑ (cid:18) sin ϑ ∂∂ϑ (cid:19) + 1sin ϑ ∂ ∂ϕ (cid:21) U (cid:27) ++ k µ ( r ) ∂U∂ϕ = 0 , (6a) ∂∂r (cid:20) ε ( r ) ∂ U∂ϑ∂r (cid:21) + 1 r ∂∂ϑ (cid:26) ε ( r ) (cid:20) ϑ ∂∂ϑ (cid:18) sin ϑ ∂∂ϑ (cid:19) + 1sin ϑ ∂ ∂ϕ (cid:21) U (cid:27) ++ k µ ( r ) ∂U∂ϑ = 0 , (6b)while equations (3d) and (3e) are fulfilled identically. Equation (3f) also is fulfilledidentically while Eq. (3c) reduces to ∂ε ( r ) ∂ϑ · ∂ U∂ϕ∂r − ∂ε ( r ) ∂ϕ · ∂ U∂ϑ∂r = 0 . (6c)Analogously, substitution of (4b) result in the following set of equations for the fieldswe recognize as the magnetic-type ones, namely, ∂∂r (cid:20) µ ( r ) ∂ V∂ϕ∂r (cid:21) + 1 r ∂∂ϕ (cid:26) µ ( r ) (cid:20) ϑ ∂∂ϑ (cid:18) sin ϑ ∂∂ϑ (cid:19) + 1sin ϑ ∂ ∂ϕ (cid:21) V (cid:27) ++ k ε ( r ) ∂V∂ϕ = 0 , (7a) ∂∂r (cid:20) µ ( r ) ∂ V∂ϑ∂r (cid:21) + 1 r ∂∂ϑ (cid:26) µ ( r ) (cid:20) ϑ ∂∂ϑ (cid:18) sin ϑ ∂∂ϑ (cid:19) + 1sin ϑ ∂ ∂ϕ (cid:21) V (cid:27) ++ k ε ( r ) ∂V∂ϑ = 0 , (7b) ∂µ ( r ) ∂ϑ · ∂ V∂ϕ∂r − ∂µ ( r ) ∂ϕ · ∂ V∂ϑ∂r = 0 . (7c)Besides curl-based equations (1a) and (1b), which are now written in the form ofEqs. (6) and (7), we must also take into account the divergency-based equations (1c)and (1d). Being written in terms of Hertz functions they look as follows, ∂∂ϑ (cid:20) ε ( r ) ∂V∂ϕ (cid:21) − ∂∂ϕ (cid:20) ε ( r ) ∂V∂ϑ (cid:21) = 0 , (8a) ∂∂ϑ (cid:20) µ ( r ) ∂U∂ϕ (cid:21) − ∂∂ϕ (cid:20) µ ( r ) ∂U∂ϑ (cid:21) = 0 . (8b)For the particular case of axially symmetric inclusion in the spherical resonator, where6 ε/∂ϕ = ∂µ/∂ϕ ≡
0, equations (6c) and (7c) get simplified to ∂ε∂ϑ · ∂ U∂ϕ∂r = 0 , (9a) ∂µ∂ϑ · ∂ V∂ϕ∂r = 0 , (9b)while equations (8) reduce to equalities ∂ε∂ϑ · ∂V∂ϕ = 0 , (10a) ∂µ∂ϑ · ∂U∂ϕ = 0 . (10b)In the case of central-symmetric non-uniformity, when ∂ε/∂ϑ = ∂µ/∂ϑ ≡
0, equa-tions (10) become redundant as they do not impose any restrictions on the EM fieldcomponents. If, however, there appears asymmetry along polar axis, say, if the centralpoints of inner and outer spheres in Figure 1 are shifted by some distance d = 0, thenspecific restrictions appear in the form of equalities H mϕ = E mϑ = 0 , (11a) E eϕ = H eϑ = 0 , (11b)which are valid, respectively, in the regions where ∂ε/∂ϑ = 0 and ∂µ/∂ϑ = 0.Before proceeding with the solution of Eqs. (6) and (7) it is worth noting that func-tions U ( r ) and V ( r ) directly are not so convenient to represent with them EM fields ina spherical resonator. More convenient for this purpose are Debye potentials [13,22],as indicated above in Eq. (5). The equations for these potentials are more advanta-geous since, for one thing, they allow to seek the solution of governing equations in theform of expansion in the Fourier series with respect to well-practical spherical func-tions. One more advantage of using Debye potentials is the possibility to develop forthem a constructive perturbation theory allowing to account for different geometricaland electromagnetic properties of inhomogeneities of the resonator infill. That is whywe shall continue the exposition in terms of precisely the DPs rather than Hertzianfunctions.From here on we will restrict our consideration to the case of non-magnetic resonatorinfill by putting in all the formulas µ ( r ) ≡
1. Using the uniformity of the resonatordielectric properties in azimuth angle ϕ , after transition to Fourier representation withrespect to this variable, with Eq. (10a) taken into account, we obtain the following setof control equations for both electric and magnetic Debye potentials, m h ˆ∆ m + k ε − ˆ V ( ε ) − ˆ V ( r ) i u m ( r, ϑ ) = 0 , (12a) ∂∂ϑ h ˆ∆ m + k ε − ˆ V ( ε ) − ˆ V ( r ) i u m ( r, ϑ )++ (cid:26) ∂∂ϑ h ˆ V ( ε ) + ˆ V ( r ) i − r ˆ W ( ϑ ) m (cid:27) u m ( r, ϑ ) = 0 , (12b)7 h ˆ∆ m + k ε − ˆ V ( ε ) i v m ( r, ϑ ) = 0 , (13a) ∂∂ϑ h ˆ∆ m + k ε − ˆ V ( ε ) i v m ( r, ϑ ) = 0 . (13b)Here, u m ( r, ϑ ) and v m ( r, ϑ ) are the m -th azimuth Fourier components of the DPs,ˆ∆ m = 1 r (cid:20) ∂∂r (cid:18) r ∂∂r (cid:19) + 1sin ϑ ∂∂ϑ (cid:18) sin ϑ ∂∂ϑ (cid:19) − m sin ϑ (cid:21) (14)is the corresponding m -th azimuth component of the full Laplace operator, ε is theaverage (over r and ϑ ) value of the permittivity ε ( r ) ≡ ε ( r, ϑ ) within the resonator,ˆ V ( ε ) = − k ∆ ε ( r, ϑ ) , (15)where ∆ ε ( r, ϑ ) = ε ( r, ϑ ) − ε , is the local amplitude-type potential accounting for thepermittivity heterogeneity, ˆ V ( r ) = ∂ ln ε ( r, ϑ ) ∂r (cid:18) r + ∂∂r (cid:19) (16)is the effective gradient-type potential having the operator nature, whose appearancein equations (12) is due to the radial inhomogeneity of the permittivity, andˆ W ( ϑ ) m = ∂ ln ε ( r, ϑ ) ∂ϑ (cid:20) ϑ ∂∂ϑ (cid:18) sin ϑ ∂∂ϑ (cid:19) − m sin ϑ (cid:21) (17)is the operator accounting for the polar-direction (tangential) non-uniformity of thepermittivity function ε ( r, ϑ ).Equations (12), (13) which describe TM and TE oscillations, respectively, accordingto their origination must be satisfied in pairs. Equations (12a) and (13a) may befulfilled either by m = 0 or by equating to zero the expressions in square brackets. Wewill consider these two possibilities independently. For zeroth azimuth modes the first equations of each of the pairs (12) and (13) arefulfilled automatically. After putting m = 0 into a pair of other equations, (12b) and(13b), and using the easily verified commutation relationship ∂∂ϑ ˆ∆ = ˆ∆ ± ∂∂ϑ , (18)we obtain the following equation pair, h ˆ∆ ± + k ε − ˆ V ( ε ) − ˆ V ( r ) − ˆ V ( ϑ ) i ∂u ( r, ϑ ) ∂ϑ = 0 , (19a) h ˆ∆ ± + k ε − ˆ V ( ε ) i ∂v ( r, ϑ ) ∂ϑ = 0 . (19b)8ere, one more gradient-type effective potential has appeared, viz.,ˆ V ( ϑ ) = ∂ ln ε ( r, ϑ ) ∂ϑ · r (cid:18) cot ϑ + ∂∂ϑ (cid:19) , (20)which reflects the eccentricity property of the resonator infill. Starting, first, from equations (13), one can make sure that for ∀ m = 0 both of themare satisfied if the equation is met h ˆ∆ m + k ε − ˆ V ( ε ) i v m ( r, ϑ ) = 0 . (21)The same does not hold for equations (12). While the fulfillment of the first of themis provided due to the action of the operator in square brackets, the second one,Eq. (12b), reduces to (cid:26) ∂∂ϑ h ˆ V ( ε ) + ˆ V ( r ) i − r ˆ W ( ϑ ) m (cid:27) u m ( r, ϑ ) = 0 . (22)Taking account of definitions (15)–(17) and the differentiation with respect to ϑ inEq. (22) one can observe that this equation is with certainty satisfied provided that thepermittivity is ϑ -uniform, i. e., if the equality holds true ∂ε ( r, ϑ ) /∂ϑ ≡
0. The oppositecase of eccentrically positioned inner sphere requires some extra examination.Consider the actual case where the permittivity is step-wise dependent on spatialradius-vector, namely, when it is represented by formula ε ( r ) ≡ ε ( r, ϑ ) = ε in Θ ( r ∈ Ω in ) + ε out Θ ( r ∈ Ω out \ Ω in ) , (23)where the notations Ω in and Ω out are used for the spatial regions enclosed within theinner and outer spheres, respectively; Θ ( A ) is the logical unit-step function defined as Θ ( A ) = ( , if A = true , , if A = f alse . (24)For the model (23), the earlier introduced average permittivity is given by ε = ε in V in + ε out ( V out − V in ) V out = ε out + ( ε in − ε out ) ( R in /R out ) , (25)where V in / out denote the volumes of the inner and outer spheres in Figure 1.Due to the presence of theta-functions in formula (23), all the terms in Eq. (22)have singular nature. To specify their singularities let us express the derivatives withrespect to ϑ and r as the derivatives in the directions normal and tangential to thesurface of the inner sphere. In Figure 2, the corresponding schematics is shown, where O ′ and O denote the centers of inner and outer sphere, respectively.For the model (23) the derivatives of function ε ( r, ϑ ) with respect to both r and ϑ are everywhere equal to zero except for the boundary of the inclusion. At some point A ′ igure 2.: The diagram for the derivatives calculation. of the boundary these derivatives are related to the normal and tangential derivativesby formulas ∂∂r = cos α ( r, ϑ ) ∂∂l n − sin α ( r, ϑ ) ∂∂l t , (26a)1 r ∂∂ϑ = sin α ( r, ϑ ) ∂∂l n + cos α ( r, ϑ ) ∂∂l t , (26b)where α ( r, ϑ ) is the running angle between the radius-vector of the point A ′ in primedreference frame and the normal to its non-primed counterpart, see Figure 2. Thetangential derivative of ε ( r, ϑ ) in that point equals to zero while the normal derivativeproduces a δ -function. After simple calculus we obtain ∂ε ( r, ϑ ) ∂r (cid:12)(cid:12)(cid:12)(cid:12) A ′ = cos α ( r, ϑ )( ε out − ε in ) δ (cid:2) r ′ ( r, ϑ ) − R in (cid:3) , (27a)1 r ∂ε ( r, ϑ ) ∂ϑ (cid:12)(cid:12)(cid:12)(cid:12) A ′ = sin α ( r, ϑ )( ε out − ε in ) δ (cid:2) r ′ ( r, ϑ ) − R in (cid:3) . (27b)Consider the term in Eq. (22) which contains the “radial-gradient” potential, ˆ V ( r ) .The singularity of that term arisen from the double-derivative of function ε ( r, ϑ ) withrespect to both of its arguments is compensated by the presence of factor (cid:18) r + ∂∂r (cid:19) u m = 1 r ∂ ( ru m ) ∂r = 1 r ∂U m ∂r , (28)which for nonzero azimuth modes disappears according to Maxwell equation (9a). Twoother terms in Eq. (22) result in equation ∂ε∂ϑ (cid:26) k ε − ˆ V ( ε ) + 1 r (cid:20) ϑ ∂∂ϑ (cid:18) sin ϑ ∂∂ϑ (cid:19) − m sin ϑ (cid:21)(cid:27) u m = 0 , (29)which can never be satisfied simultaneously with Eq. (12a) except for the cases of10ither ∂ε/∂ϑ ≡ u m ≡ ∀ m = 0. This makes us conclude that in the case ofeccentrically nonuniform double-spherical resonator only zeroth azimuth mode of TMoscillations is allowed for the existence, while the magnetic-type (TE) oscillations canhave any azimuth-mode index.
3. Separation of TM and TE oscillation modes
Although we have uniquely defined the TM and TE oscillations in the resonator underexamination, as yet we are not sure if these oscillations can be decoupled from oneanother, in spite of the fact that they are not directly intermixed in equation (19a) forTM and equations (19b) and (21) for TE oscillation modes. To clarify this position, wemust account for boundary conditions (BCs) the fields in our resonator should satisfy.The main BC for our problem is the vanishing of tangential component of the electricfield on the boundary of the resonator, i. e., on the outer sphere in Figure 1. For thegeometry chosen this implies the synchronous vanishing of the components E ϑ and E ϕ of the net electric field at r = R out at any angles ϑ and ϕ . From Eqs. (4a) and(4b) it can be seen that the BCs represented in terms of Debye potentials result in thefollowing two equalities, (cid:20) ε out r ∂∂ϑ (cid:18) u + r ∂u∂r (cid:19) + ik sin ϑ ∂v∂ϕ (cid:21) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r = R out = 0 , (30a) (cid:20) ε out r sin ϑ ∂∂ϕ (cid:18) u + r ∂u∂r (cid:19) − ik ∂v∂ϑ (cid:21) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r = R out = 0 , (30b)where the fields of TM and TE polarization are seemingly intermixed.Yet, it is not difficult to observe that equalities (30) would definitely be satisfied ifa couple of more simple relations were met, specifically, (cid:18) u + r ∂u∂r (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) r = R out =0 , (31a)and v ( r = R out ) =0 . (31b)The choice of the boundary conditions for the above-pointed governing equations,(19) and (21), in the form of equalities (31) rather than (30) would allow one to seekthe electric- and magnetic-type fields in the resonator with broken radial and azimuthalsymmetries independently of each other, thus representing the oscillations of a generaltype as a superposition of uncoupled TM and TE oscillation modes.In the next section we show that the BCs in the form of Eqs. (31) are, indeed,admissible for the resonator with inhomogeneous infill, which gives us the opportunityto consider TM and TE oscillations as the uncoupled ones in spite of their being seem-ingly intermixed by BCs (30). The proof can be given, however, not in the coordinaterepresentation but rather in the representation of spherical resonator eigenmodes.11 . Fourier representation of the governing equations
For the subsequent analysis of working equations (19) and (21) we perform theirFourier transformation from coordinate to momentum representation. For this purposewe will use two different complete eigenfunction sets for the Laplace operator. Fora homogeneous sphere of radius R out these functions in spherical coordinates may bechosen in the common form | r ; µ i = D ( l ) n R out r r J l + (cid:16) λ ( l ) n r/R out (cid:17) Y ml ( ϑ, ϕ ) (32) (cid:0) n = 1 , , . . . , ∞ ; l = 0 , , , . . . , ∞ ; m = − l, − l + 1 , . . . l − , l (cid:1) , where µ = { n, l, m } is the vector-valued mode index whose components correspond toradial, polar, and azimuth coordinates, respectively; D ( l ) n is the normalization coeffi-cient which depends on the boundary conditions and is represented by the followingformulae, h D ( l ) n i = h J l + (cid:0) λ ( l ) n (cid:1)i − " − l ( l + 1) λ ( l ) n − for BC (31a), (33a) h D ( l ) n i = h J l + (cid:0) λ ( l ) n (cid:1)i − for BC (31b). (33b) J p (cid:0) u ) is the Bessel function of the first kind, Y ml ( ϑ, ϕ ) is the spherical function; λ ( l ) n isthe set of positive zeros either of the sum J l + ( u ) + 2 uJ ′ l + ( u ), if boundary condition(31a) is applied, or of function J l + ( u ) in the case of BC (31b) (see, e. g., [23], n. 7.10.4).The zeros will be assumed enumerated by natural number n in ascending order. Theeigenvalues of the Laplace operator, which correspond to functions (32), are given by E µ = − k µ = − λ ( l ) n R out ! , (34)where the degeneracy in azimuth index m is apparent, whose order is 2 l + 1.As far as we have already done the Fourier transform of working equations overazimuth variable ϕ , to complete the transformation of Eqs. (19) and (21) we willapply the normalized eigenfunction set other than (32), namely, | r ; µ m i = D ( l µ m ) n µ m R out r r J l µ m + (cid:16) λ ( l µ m ) n µ m r/R out (cid:17) Θ l µ m m ( ϑ ) . (35)Here, µ m is the two-dimensional vector mode index, where in the initially three-dimensional index µ the azimuth component is put equal to the specific value m ,Θ lm ( ϑ ) = ( − m (cid:20) l + 12 · ( l − m )!( l + m )! (cid:21) / P ml (cos ϑ ) (36)is the normalized adjoint Legendre function. The set of functions (35) is orthonormaland complete, with orthonormality and completeness relations having the standard12orm [24], X µ m | r ; µ m ih r ′ ; µ m | = δ ( r − r ′ ) (cid:0) r , r ′ ∈ Σ out (cid:1) , (37a) Z Σ out d r | r ; µ m ih r ; µ ′ m | = δ µ m µ ′ m . (37b)The summation in Eq. (37a) runs over polar and radial mode indices only, Σ out is theouter sphere central cross-section.The set of eigenfunctions (35) is convenient in that it is quite universal in formwith respect to BCs its constituents may fulfill. We may choose for the Fourier seriesexpansion of functions u ( r ) and v ( r ) different basis function sets, being guided byconvenience considerations. In particular, in the series (cid:18) u ( r ) v ( r ) (cid:19) = X µ (cid:18) u µ v µ (cid:19) | r ; µ i (cid:12)(cid:12) u,v (38)as a function | r ; µ i (cid:12)(cid:12) u we can take function (32) that meets BC (31a), whereas basisfunction | r ; µ i (cid:12)(cid:12) v in the expansion of v ( r ) can be chosen to obey BC (31b). In sucha way we, for one thing, fulfill general boundary conditions (30) and, for the other,this permits us to solve the problem for TM and TE oscillations independently ofeach other. This recipe may be considered as a kind of TM and TE mode separationprocedure. The only difference between formulas for electric-type and magnetic-typefields reduces to different sets of zeros λ ( l ) n entering into the basis functions.Equations (19) and (21), which actually define the entire spectrum of the resonatorunder investigation, are identical in the functional structure and differ only in theset of effective potentials entering into them. To describe the technique of subsequentoperations with these equations it is more convenient to work with them not in theinitial form but rather to employ the equations for the corresponding Green functions.Such an approach is even more justified if we consider that the oscillation spectrum ofthe resonator, which in fact is the ultimate goal of our study, is determined by densityof states (DOS) ν ( k ) = π − Im (cid:8) Tr ˆ G (+) (cid:9) , (39)where ˆ G (+) is the Green operator with upper (+) index signifying its retarded nature.By denoting the effective potential, which generally possesses the operator structure,through ˆ V and accomplishing the Fourier transform of Eqs. (19) and (21) with re-spect to coordinate variables r and ϑ we arrive at the following infinite set of coupledalgebraic equations for the Green operator mode matrix elements, viz., (cid:0) k ε − k µ m − V µ m (cid:1) G µ m µ ′ m − X ν m = µ m U µ m ν m G ν m µ ′ m = δ µ m µ ′ m . (40)Matrix element U µ m ν m = Z Σ out d r h r ; µ m | ˆ V | r ; ν m i (41)13f the effective potential, in view of the restriction ν m = µ m in the sum of Eq. (40),will be henceforth referred to as the inter-mode effective potential, while the diagonalelement of the potential matrix, V µ m ≡ U µ m µ m , which simply renormalizes energeticparameter k µ m without change in the mode index, will be termed the intra-mode potential.In the further analysis we will regard the average dielectric constant of the resonatoras either purely real or complex-valued quantity, considering either a dissipation-freetheoretical model or realistic systems with dissipation. In the latter case we will assumethat ε = ε ′ + iε ′′ , where ε ′ and ε ′′ denote the real and imaginary parts of the averagepermittivity, respectively. In general, the specific value of ε ′′ is dependent on bothof the dissipative parts of ε in and ε out , which even may have different signs. Yet,everywhere below we will assume that ε ′′ > ε ′′ , are quitetime consuming for quasioptical resonators, especially in the case where dissipationis taken finite. The latter case, however, can be effectively analyzed by analyticalmeans, if one uses the method of effective disconnection of the resonance modes ofan inhomogeneous system, which is detailed in Ref. [20]. Here we will give a briefsketch of the method and present the result of its implementation.
5. Analytic solution of equation set (40)
Considering equations (40) one can observe that if all the inter-mode potentials wereequal to zero the Green function matrix would be essentially diagonal, G µ m µ ′ m = G ( V ) µ m δ µ m µ ′ m . (42)Here, function G ( V ) µ m , which we term henceforth the trial mode Green function, is givenby G ( V ) µ m ( k ) = (cid:0) k ε − k µ m − V µ m (cid:1) − . (43)This function already contains some information regarding the cavity non-uniformity,yet this information is only partial, concentrated in the intra-mode potential V µ m . Theinter-mode potentials in Eq. (40) result in effective “dynamic” inelasticity of the trialmode scattering, which makes actual states of the field in the cavity ill defined. Thiskind of inelasticity is due to the failure of the wave phase in the course of its propa-gation in the non-uniform medium, and it has little in common with true inelasticitydue to the absorption or gain.By simple algebraic manipulations (see Ref. [20] for details), from the set of equa-tions (40) the linear connection between non-diagonal and diagonal elements of themode Green matrix results, which in the operator form is given by G ν m µ m = ˆ P ν m (ˆ − ˆ R m ) − ˆ R m ˆ P µ m G µ m µ m . (44)Here, ˆ P µ m is the Feshbach-type projection operator [25,26] whose action reduces tothe assignment of the value µ m to the nearest mode index of any adjacent operator,14o matter where it may stand — to the left or to the right of ˆ P µ m ; ˆ R m = ˆ G ( V ) m ˆ U m isthe mode-mixing operator defined on the reduced space of mode indices, M µ m , whichincludes all the indices of the resonator modes save index µ m ; ˆ G ( V ) m is the operatorwith the diagonal matrix consisting of trial Green functions (43), ˆ U m is the operatorwith null-diagonal matrix, which incorporates all the inter-mode potentials enteringinto Eq. (40). Assuming then in equation (40) µ ′ m = µ m and substituting into theachieved equation all the intermode propagators in the form (44) we arrive at thefollowing expression for the exact intramode propagator G µ m µ m , viz., G µ m µ m = (cid:2) k ε − k µ m − V µ m − T µ m ( k ) (cid:3) − . (45)Here, T µ m ( k ) = ˆ P µ m ˆ U m (ˆ − ˆ R m ) − ˆ R m ˆ P µ m (46)is the portion of the mode µ m eigen-energy which is related to the intermode scatteringonly. The validity condition for the expressions (46) and (44) is the non-singularity ofthe inverse operator entering into them. To meet this requirement, it suffices for theoperator ˆ R m not to have the unity among its eigenvalues. The universal method toaccomplish this task is to make the wave operator in Eqs. (19) and (21) non-Hermitian,just like in the papers by Feshbach [25,26]. However, if there it was naturally achieveddue to the openness of the system, in our case the only suitable way is to add to theaverage permittivity of the resonator infill some imaginary part related to dissipationor gain, whatever small. In the final results this imaginary addendum can be directedto zero if there is a need to consider the dissipation-free system.The main feature of expression (45) for the Green matrix diagonal element is thateffective potential T µ m ( k ), which is known as optical potential in the quantum theory ofscattering [27], depends in a quite complicated manner on k , which is the consequenceof the dependence on frequency of the trial Green functions entering into the mode-mixing operator ˆ R m . Typically, optical potentials are complex-valued, if it comes toopen systems. Yet, we study the conservative system, so the complexity necessaryfor the inverse operator in Eq. (46) to be properly defined is introduced through theseed dissipation or gain hidden in the imaginary part of the infill permittivity. Eveninfinitesimal value of ε ′′ can result in difficult-to-control value of Im T µ , of which onecan make sure by rewriting potential (46) as T µ ( k ) = X ν , ν = µ U µν (cid:26)h ˆ − ˆ G ( V ) ˆ U i − (cid:27) ν ν (cid:16) ˆ G ( V ) ˆ U (cid:17) ν µ = X ν , ν = µ U µν (cid:20)(cid:16) △ + k ε − ˆ V − ˆ U (cid:17) − (cid:21) ν ν U ν µ = X ν , ν = µ U µν (cid:20)(cid:16) △ + k ε ∗ − ˆ V − ˆ U (cid:17) (cid:12)(cid:12)(cid:12)(cid:16) △ + k ε − ˆ V − ˆ U (cid:17)(cid:12)(cid:12)(cid:12) − (cid:21) ν ν U ν µ (47)(to be short, we omit the insignificant index m ). It can be seen that equality Im T µ = 0implies the quantity ε to have exactly zero imaginary part. In the opposite case, thevalue of Im T µ is highly dependent on k , which results in different local widths of theresonances. 15f one would go to find spectral values of k through the poles and other exceptionalpoints of Green function (45), in view of extreme complexity of Re T µ m ( k ) the availablevalues of k were not all real. The situation here is reminiscent of that considered byHatano in studies on open quantum-mechanical systems [28]. In our case the effectiveopenness appears if one turns to the trajectory-based interpretation of resonant statesin the cavity under consideration [29]. Thus arising the finite width of the spectrallines can lead to their significant overlap, which may be interpreted as the interactionof resonances declared to be the essential source of spectral chaos [30].
6. Numerical verification of the developed theory
Unfortunately, in the general case the above developed method does not allow obtain-ing the exact analytic expressions for the eigen-frequencies of an eccentrically layeredspherical cavity resonator. To overcome this disadvantage by numerical methods, wehave developed a software that implements the proposed theory. The software usesthe well-known library for scientific computing [31] and is based on the MPI paral-lel computing technology. The calculations were performed on the grid cluster of theInstitute for Radiophysics and Electronics of NAS of Ukraine.To verify the proposed theory and the software developed, the eigen-frequencies ofthe central symmetric layered cavity resonator (whose spectrum can be found by otherknown methods) were calculated. For comparison, we also used the method of solvingthe Maxwell equations in the outer and inner layers with subsequent matching therespective solutions at the inner boundary of the resonator. The results are presentedin Figure 3, where we provide two different parts of the spectrum of spherical res-onator with centralized spherical insert of a given radius. Thin vertical lines depict the N o r m a li ze d d e n s it y o f s t a t e s Wave number k, cm -1 a b N o r m a li ze d d e n s it y o f s t a t e s Wave number k, cm -1 Figure 3.:
The normalized density of states (black curves) and the matrix eigenvalues (redvertical lines) for TE oscillations of the concentrically layered spherical cavity in two particularfrequency domains. The resonator parameters are taken as follows: R in = 2 cm, R out = 4 .
46 cm, ε in = 2 .
08 + 0 . j , ε out = 1 + 0 . j , d = 0; 1000 basis modes were used when solving matrixequation (40). positions of the expected resonance maxima of the Green function, which are foundnumerically from Eq. (40) with ε ′′ from the outset taken exactly zero.By the solid lines in Figure 3, the spectrum obtained from Eq. (39) with Greenmatrix elements taken in the form (45) is shown, where the dissipation is supposed tobe finite but vanishingly small. It is evident that with the presence of dissipation/gainmechanisms in the system the spectrum becomes highly rarefied due to mutual ab-16orption of the smoothed resonances. Moreover, the spectrum, being (quasi-)discretefor the resonator of simple integrable form, tends to develop into the continuous onewhen the ray dynamics comes to be to a certain degree non-integrable. A significantoverlap of the resonances occurs even at infinitesimal value of dissipation, which stemsfrom the fact that the ray trajectories in the resonator with broken symmetry, nomatter how much it would be, after a sufficiently large period of time evolve intochaotic bundles. According to the widespread Gutzwiller’s periodic-orbit theory [29],the bursts on the continuous spectral line correspond to the most stable and long-livedquasi-classic trajectories while the smooth spectral “background” fits with the set ofchaotic trajectory sets. We will discuss in detail the chaotic properties of the spectrumin the subsequent publication, which we intend to devote to the statistical analysis ofthe inter-frequency intervals extracted from the above obtained formulas.For additional verification, we find the limited set of eigen-frequencies of the centralsymmetric layered cavity by means of the dissipation-free version of equation (40) andin parallel got them from the CST Microwave Studio simulations with the use of theEigen mode solver. A comparison of the results of these two independent methods ispresented in Figure 4, wherefrom their coincidence is certain. f, GHz Figure 4.:
The comparison of the frequency spectrum of the layered central symmetric di-electric cavity. Thin lines result from the dissipation-free version of equation set (40) whereasthick lines follow from the CST Microwave Studio simulations with application of the Eigenmode solver. The resonator parameters are as follows: R in = 5 mm, R out = 10 mm, ε in = 2 . ε out = 1, d = 0. Additionally, in Figure 5 we present the distribution of electromagnetic field of thefirst mode of the eigen- mode set obtained using CST Microwave Studio simulations.This resonant mode, as expected, is uniform in the azimuth angle ϕ for both theelectric and magnetic fields.
7. Concluding remarks
In the present paper we present a theory that is capable of the solution to Maxwellequations in spherical cavity resonators with non-homogeneous dielectric infill. Specif-ically, we have examined the resonator whose non-uniformity is due to a sphericaldielectric insert with constant permittivity different from that of the host resonator.17 a) | E | distribution(b) | H | distribution Figure 5.:
Distribution of the absolute values of electric (a) and magnetic (b) fields forthe lowest-frequency mode (11.47 GHz); cut plane through the cavity center. The resonatorparameters are the same as in Figure 4.
The position of the insert inside the host uniform cavity as well as its radius are as-sumed to be optional. To obtain the oscillation spectrum, the vector Maxwell equationsare reduced to the set of two scalar equations for Hertzian functions of electric andmagnetic types, which have the form of Helmholtz equations subject to the perturba-tion. The similarity of these equations to Schr¨odinger equation has made it possible toavoid their solving in different parts of the resonator, with necessary joining of the so-lutions at the boundaries between them. We choose the method to solve the equationswithin the whole resonator by representing its heterogeneity in terms of effective po-tentials, which appear to be of two types — the amplitude-type and the gradient-typepotentials.The EM field in the angularly symmetric cavity can always be represented in termsof TM and TE Fourier components. In this case the general solution is shown tobe a superposition of TM and TE oscillations which are disentangled in spite of theparticular radial non-uniformity. If the resonator is angularly symmetric, whatever thenon-regularity in radial direction, there are present all angular modes in its spectrum,both of TM and TE type. Each type of the oscillations is described by only one scalareffective potential, which substantially reduces the computational resources requiredfor finding the resonator spectrum by means of computer simulations.The violation of angular symmetry in the infill permittivity dramatically affects thespectrum, namely, the structure and the number of degrees of freedom (independent18ariables) of the EM field in the cavity. Misalignment of the insert with regard to thehost resonator center makes all azimuth modes of TM oscillations, except for the zeroth(uniform) mode, become forbidden in the resonator in question, while TE-polarizedoscillations are allowed to have any index of the azimuthal mode. The asymmetry ofthis particular kind is related to the fact that in the present study the permittivityonly is chosen to be non-uniform while the permeability is taken to be everywhereconstant. In case of the both electrically and magnetically non-smooth axially sym-metric resonator with eccentric infill, only azimuthally uniform eigenmodes can berepresented as the independent TE and TM oscillations, i. e., by the pair of uncoupledscalar Hertz potentials. Other azimuthal modes, even if they exist in the spectrum,cannot be found in frameworks of the method presented here and thus should becategorized as the hybrid modes.The spectrum of non-homogeneous resonator, no matter angularly symmetric or not,is shown to acquire some chaotic properties which cannot be meaningfully analyzedwithout invoking the specific methods of the chaos theory. The detailed analysis ofthis issue will be presented in the subsequent publication.
Appendix A. Evaluation and comparison of effective potentials in Eq. (40)
The effective potential ˆ V in Eq. (41) is, in general, a sum of three terms, ˆ V ( ε ) , ˆ V ( r ) and ˆ V ( ϑ ) , of which the first term enters both into Eqs. (19) and (21) while the othertwo into Eq. (19a) only.Consider, first, the potential ˆ V ( ε ) . As opposed to gradient potentials ˆ V ( r ) and ˆ V ( ϑ ) ,the amplitude-type potential ˆ V ( ε ) is not quite perceptive to the coordinate system onechooses for its calculation. For the purposes of numerical calculations it is preferable totake for this potential the coordinate system K located at the center of outer sphere,see the diagram in Figure A1. Taking account of piecewise continuity of the resonator Figure A1.:
To the derivation of formula (A2) in the K frame. Symbols d ≷ in the figurecorrespond to d ≷ R in . infill, the mode matrix element of ˆ V ( ε ) can be represented as U ( ε ) µ m ν m = − k ( ε in − ε out ) h I ( ε ) µ m ν m − (cid:0) R in /R out (cid:1) δ µ m ν m i , (A1)19here I ( ε ) µ m ν m = Z Σ in d r h r ; µ m | r ; ν m i = 2 R D ( l µ m ) n µ m D ( l ν m ) n ν m ×× Z Σ in d r r J l µ m + (cid:16) λ ( l µ m ) l µ m r/R out (cid:17) J l ν m + (cid:16) λ ( l ν m ) l ν m r/R out (cid:17) Θ l µ m m ( ϑ )Θ l ν m m ( ϑ ) , (A2)Σ in is the inner sphere central cross-section.Transform the integral in (A2) for two variants of the inner sphere, encircling andnot encircling the center of the outer sphere. Consider, first, the inner sphere centeredat point d > in Figure A1, when the point K does not fall into it. The integral overregion Σ in reduces to integration with respect to ϑ from zero to ϑ c = arcsin( R in /d > )while the integral with respect to r runs in the interval r min ( ϑ ) < r < r max ( ϑ ),with r min ( ϑ ) and r max ( ϑ ) being the points where the ray going from K intersects thecircle C . As far as this circle is described by equation r − rd > cos ϑ + d > = R , (A3)one can easily find that r maxmin ( ϑ ) (cid:12)(cid:12) C = d > cos ϑ ± q R − d > sin ϑ . (A4)If the coordinate origin K lies within the inner sphere, integration with respect to ϑ in (A2) is done from 0 to π while integration with respect to r runs from 0 to r max ( ϑ ).In the general case the expression for I ( ε ) µ m ν m can be written in the universal form I ( ε ) µ m ν m = B µ m ν m R Z t min dtP ml µ m ( t ) P ml ν m ( t ) ×× r max ( t ) Z r min ( t ) rdrJ l µ m + (cid:16) λ ( l µ m ) n µ m r/R out (cid:17) J l ν m + (cid:16) λ ( l ν m ) n ν m r/R out (cid:17) , (A5)where the notations are used t min = θ ( d − R in ) p − ( R in /d ) − θ ( R in − d ) , (A6a) r min ( t ) = θ ( d − R in ) h dt − q R − d (1 − t ) i , (A6b) r max ( t ) = dt + q R − d (1 − t ) , (A6c) θ ( . . . ) is the Heaviside unit step function, and the coefficient B µ m ν m has the form B µ m ν m = D ( l µ m ) n µ m D ( l ν m ) n ν m (cid:20) (2 l µ m + 1) ( l µ m − m )!( l µ m + m )! (cid:21) / (cid:20) (2 l ν m + 1) ( l ν m − m )!( l ν m + m )! (cid:21) / . (A7)20radient potentials ˆ V ( r ) and ˆ V ( ϑ ) are more suited for calculation of their matrixelements in the spherical coordinate frame located at the center of the inner sphere.This is due to the fact that in the definitions of these potentials, Eqs. (16) and (20),there are present derivatives of function ε ( r, ϑ ) with respect to both radial and polarcoordinates. These derivatives make the potentials ˆ V ( r ) and ˆ V ( ϑ ) singular, namely,proportional to δ -functions arisen when crossing the boundary of the inner sphere.Although singularities of this form are not problematic for matrix element calculations,to concretize them for a specific position of the inner sphere it is advantageous to go tothe reference system related just to this sphere center. For this purpose we can makeuse of the gage preserving coordinate transformation r sin ϑ = r ′ sin ϑ ′ ,r cos ϑ = r ′ cos ϑ ′ + d , (A8)where the primes relate to the inner-sphere coordinate frame. In this system, thelogarithmic derivatives entering into Eqs. (16) and (20) are calculated to (cid:20) ∂ ln ε ( r, ϑ ) ∂r (cid:21) ( ϑ ′ , r ′ ) = 2 ε out − ε in ε out + ε in · R in + d cos ϑ ′ q R + 2 R in d cos ϑ ′ + d δ ( r ′ − R in ) , (A9a) (cid:20) ∂ ln ε ( r, ϑ ) ∂ϑ (cid:21) ( ϑ ′ , r ′ ) = 2 ε out − ε in ε out + ε in d sin ϑ ′ δ ( r ′ − R in ) . (A9b)When deriving the latter two formulas, we have put the permittivity value on theborder between regions with ε in and ε out equal to their half-sum.With the use of (A9), mode matrix elements of both of the gradient potentialsentering into Eq. (19a) are calculated to U ( r ) µ ν = B µ ν R · ε out − ε in ε out + ε in I ( r ) µ ν , (A10a) U ( ϑ ) µ ν =2 M in B µ ν R · ε out − ε in ε out + ε in I ( ϑ ) µ ν , (A10b)where the integrals I ( r ) µ ν and I ( ϑ ) µ ν have the form I ( r ) µ ν = Z − dt M in tξ (cid:0) t ; M in (cid:1) P l µ (cid:20) t + M in ξ (cid:0) t ; M in (cid:1) (cid:21) P l ν (cid:20) t + M in ξ (cid:0) t ; M in (cid:1) (cid:21) J l µ + (cid:18) λ ( l µ ) n µ R in R out ξ (cid:0) t ; M in (cid:1)(cid:19) × ( ξ (cid:0) t ; M in (cid:1) J l ν + (cid:18) λ ( l ν ) n ν R in R out ξ (cid:0) t ; M in (cid:1)(cid:19) ++ λ ( l ν ) n ν R in R out (cid:20) J l ν − (cid:18) λ ( l ν ) n ν R in R out ξ (cid:0) t ; M in (cid:1)(cid:19) −− J l ν + (cid:18) λ ( l ν ) n ν R in R out ξ (cid:0) t ; M in (cid:1)(cid:19) (cid:21)) , (A11a)21 ( ϑ ) µ ν = Z − dt h ξ (cid:0) t ; M in (cid:1)i J l µ + (cid:20) λ ( l µ ) n µ R in R out ξ (cid:0) t ; M in (cid:1)(cid:21) J l ν + (cid:20) λ ( l ν ) n ν R in R out ξ (cid:0) t ; M in (cid:1)(cid:21) ×× P l µ (cid:20) t + M in ξ (cid:0) t ; M in (cid:1) (cid:21)( P l ν (cid:20) t + M in ξ (cid:0) t ; M in (cid:1) (cid:21)(cid:0) t + M in (cid:1) −− (cid:16) P l ν (cid:17) ′ (cid:20) t + M in ξ (cid:0) t ; M in (cid:1) (cid:21) − t ξ (cid:0) t ; M in (cid:1) ) . (A11b)In Eqs. (A11), function ξ (cid:0) t ; M in (cid:1) = q tM in + M is the scaling factor that repre-sents the dimensionless (in units of R in ) length of radius-vector r ∈ Σ out , M in = d/R in is the dimensionless (in the same units) misalignment of the dielectric spheres of theresonator infill. References [1] Zouros GP. Eigenfrequencies and Modal Analysis of Uniaxial, Biaxial, and Gyro-electric Spherical Cavities. IEEE Transactions on Microwave Theory and Techniques.2017;65(1):20–27.[2] Zouros GP. Analysis of Multilayered Gyroelectric Spherical Cavities by WeakForm VIE Formulation. IEEE Transactions on Microwave Theory and Techniques.2017;65(11):4029–4036.[3] Kolezas GD, Zouros GP. Eigenfrequencies in Gyrotropic Metallic Cavities. IEEE Mi-crowave and Wireless Components Letters. 2018;28(3):197–199.[4] Bozza G, Oliveri G, Raffetto M. Cavities Involving Metamaterials With an Uncount-able Set of Resonant Frequencies. IEEE Microwave and Wireless Components Letters.2007;17(8):565–567.[5] Julien A, Guillon P. Electromagnetic Analysis of Spherical Dielectic Shielded Resonators.IEEE Transactions on Microwave Theory and Techniques. 1986;34(6):723–729.[6] Dankov PI. Two-resonator method for measurement of dielectric anisotropy in multilayersamples. IEEE Transactions on Microwave Theory and Techniques. 2006;54(4):1534–1544.[7] Okada F, Tanaka H. Measurement of ferrite tensor permeability using a spherical cavityresonator. IEEE Transactions on Instrumentation and Measurement. 1991;40(2):476–479.[8] St¨ockmann HJ. Quantum Chaos: An Introduction. New York: Cambridge UniversityPress; 2007.[9] Bohigas O, Giannoni MJ. Level density fluctuations and random matrix theory. Annalsof Physics. 1975;89(2):393 – 422.[10] Bohigas O, Giannoni MJ, Schmit C. Characterization of Chaotic Quantum Spectra andUniversality of Level Fluctuation Laws. Phys Rev Lett. 1984;52(1):1–4.[11] Baryakhtar VG, Yanovsky VV, Naydenov SV, Kurilo AV. Chaos in composite billiards.J Exp Theor Phys. 2006;103(2):292 – 302.[12] Stratton JA. Electromagnetic Theory. IEEE Press Series on Electromagnetic WaveTheory. Wiley; 2007.[13] Vainshtein LA. Electromagnetic waves (2nd revised and enlarged edition). Moscow,Izdatel’stvo Radio i Sviaz’; 1988. In Russian.[14] Whittaker ET. On an expression of the electromagnetic field due to electrons by meansof two scalar potential functions. Proc Lond Math Soc. 1904;s2-1(1):367–372.[15] Nisbet A. Hertzian electromagnetic potentials and associated gauge transformations. ProcR Soc Lond A. 1955;231(1185):250–263.
16] Prze´ziecki S, Hurd RA. A note on scalar Hertz potentials for gyrotropic media. Appliedphysics. 1979;20(4):313–317.[17] Prze´ziecki S, Laprus W. On the representation of electromagnetic fiels in gyrotropic mediain terms of scalar Hertz potentials. Journal of Mathematical Physics. 1982;23(9):1708–1712.[18] Weiglhofer WS. Scalar Hertz potentials for nonhomogeneous uniaxial dielectric – mag-netic mediums. International Journal of Applied Electromagnetics and Mechanics.2000;11(3):131–140.[19] Malykh MD, Sevastianov LA, Tiutiunnik AA, Nikolaev NE. On the representation ofelectromagnetic fields in closed waveguides using four scalar potentials. Journal of Elec-tromagnetic Waves and Applications. 2018;32(7):886–898.[20] Ganapolskii EM, Eremenko ZE, Tarasov YV. Influence of random bulk inhomogeneitieson quasioptical cavity resonator spectrum. Phys Rev E. 2007;75(2):026212.[21] Haus HA, Huang W. Coupled-mode theory. Proceedings of the IEEE. 1991;79(10):1505–1518.[22] Oraevsky AN. Whispering-gallery waves. Quantum Electron. 2002;32(5):377–400.[23] Erdelyi A. Higher Transcendental Functions, vol. 2. McGraw-Hill; 1953.[24] Jackson JD. Classical electrodynamics. Wiley; 1999.[25] Feshbach H. Unified theory of nuclear reactions. Annals of Physics. 1958;5(4):357–390.[26] Feshbach H. A unified theory of nuclear reactions. II. Annals of Physics. 1962;19(2):287–313.[27] Wu Ty, Ohmura T. Quantum Theory of Scattering. Dover Publications; 2011.[28] Hatano N. Equivalence of the effective Hamiltonian approach and the Siegert boundarycondition for resonant states. Fortschritte der Physik. 2013;61(2-3):238–249.[29] Gutzwiller MC. Periodic Orbits and Classical Quantization Conditions. Journal of Math-ematical Physics. 1971;12(3):343–358.[30] Zaslavski˘ı GM, Chirikov BV. Stochastic instability of non-linear oscillations. SovietPhysics Uspekhi. 1972;14(5):549–568.[31] Galassi M, Davies J, Theiler J, Gough B, Jungman G, Booth M, et al. GNU ScientificLibrary Reference Manual. 2nd ed. Gough B, editor. Boston (MA): Computer Science;2003.16] Prze´ziecki S, Hurd RA. A note on scalar Hertz potentials for gyrotropic media. Appliedphysics. 1979;20(4):313–317.[17] Prze´ziecki S, Laprus W. On the representation of electromagnetic fiels in gyrotropic mediain terms of scalar Hertz potentials. Journal of Mathematical Physics. 1982;23(9):1708–1712.[18] Weiglhofer WS. Scalar Hertz potentials for nonhomogeneous uniaxial dielectric – mag-netic mediums. International Journal of Applied Electromagnetics and Mechanics.2000;11(3):131–140.[19] Malykh MD, Sevastianov LA, Tiutiunnik AA, Nikolaev NE. On the representation ofelectromagnetic fields in closed waveguides using four scalar potentials. Journal of Elec-tromagnetic Waves and Applications. 2018;32(7):886–898.[20] Ganapolskii EM, Eremenko ZE, Tarasov YV. Influence of random bulk inhomogeneitieson quasioptical cavity resonator spectrum. Phys Rev E. 2007;75(2):026212.[21] Haus HA, Huang W. Coupled-mode theory. Proceedings of the IEEE. 1991;79(10):1505–1518.[22] Oraevsky AN. Whispering-gallery waves. Quantum Electron. 2002;32(5):377–400.[23] Erdelyi A. Higher Transcendental Functions, vol. 2. McGraw-Hill; 1953.[24] Jackson JD. Classical electrodynamics. Wiley; 1999.[25] Feshbach H. Unified theory of nuclear reactions. Annals of Physics. 1958;5(4):357–390.[26] Feshbach H. A unified theory of nuclear reactions. II. Annals of Physics. 1962;19(2):287–313.[27] Wu Ty, Ohmura T. Quantum Theory of Scattering. Dover Publications; 2011.[28] Hatano N. Equivalence of the effective Hamiltonian approach and the Siegert boundarycondition for resonant states. Fortschritte der Physik. 2013;61(2-3):238–249.[29] Gutzwiller MC. Periodic Orbits and Classical Quantization Conditions. Journal of Math-ematical Physics. 1971;12(3):343–358.[30] Zaslavski˘ı GM, Chirikov BV. Stochastic instability of non-linear oscillations. SovietPhysics Uspekhi. 1972;14(5):549–568.[31] Galassi M, Davies J, Theiler J, Gough B, Jungman G, Booth M, et al. GNU ScientificLibrary Reference Manual. 2nd ed. Gough B, editor. Boston (MA): Computer Science;2003.