Axion Mie Theory of Electron Energy Loss Spectroscopy in Topological Insulators

Electronic topological states of matter exhibit novel types of responses to electromagnetic fields. The response of strong topological insulators, for instance, is characterized by a so-called axion term in the electromagnetic Lagrangian which is ultimately due to the presence of topological surface states. Here we develop the axion Mie theory for the electromagnetic response of spherical particles including arbitrary sources of fields, i.e., charge and current distributions. We derive an axion induced mixing of transverse magnetic and transverse electric modes which are experimentally detectable through small induced rotations of the field vectors. Our results extend upon previous analyses of the problem. Our main focus is on the experimentally relevant problem of electron energy loss spectroscopy in topological insulators, a technique that has so far not yet been used to detect the axion electromagnetic response in these materials.


Introduction
It has been long known theoretically that in CP-violating theories topological defects become electrically charged, one of the most prominent examples being the induced fractional electric charge in 't Hooft-Polyakov monopoles [1]. This result follows immediately from the presence of a CP-violating interaction in the Lagrangian of the gauge theory, where F µν is the usual field strength tensor of the gauge sector of the theory and g is the gauge coupling. The significance of the parameter Θ varies depending on the context, but we will refer to it here simply as the axion, in reference to the chiral (or axial) anomaly [2]. Within an Abelian gauge theory context, observable consequences of axion electrodynamics have been discussed for some time, with references to materials already being made in early papers. For instance, Nielsen and Ninomiya [3] mentioned HgCdTe and Wilczek [4] discusses an application of axion electrodynamics to PbTe, following a suggestion by Fradkin et al. [5]. However, it was not until more recently with the prediction and discovery of several topological materials, notably topological insulators (TIs) [6,7] and Weyl semimetals [8], that the many possible interesting experimental consequences of axion electrodynamics came closer to spotlight. In this work we mainly focus on TIs, but the theory described here applies with minor modifications to Weyl semimetals as well. Three-dimensional TIs feature a quantum Hall electromagnetic response at the surface, which emerges as a consequence of an insulator behavior in the bulk characterized effectively by the following axion electrodynamics Lagrangian [9], where E and B are the electric and magnetic field, respectively, α = e 2 /( c) is the finestructure constant and we have assumed a paramagnetic bulk with magnetic permeability µ = 1. For TIs where either time-reversal or inversion symmetry holds, Θ = π. Thus, a TI sample in vacuum may be considered as the problem of a topological dielectric satisfying the Maxwell equations, where we have used E → Ee −ickt , B → Be −ickt , current j → je −ickt , charge density ρ → ρe −ickt , and assumed a frequency-dependent dielectric function, (ω, r) = (ck, r). The resulting constitutive relations then read, When j = 0 the above constitutive equations can be rewritten as, where µ = /[ + (αΘ/π) 2 ]. As noted in Ref. [10], in this case the problem reduces to the one of light scattering by an optically active sphere, solved by Bohren long time ago [11,12]. A closely related problem arises also in the study of a so called chiral dielectric-magnetic medium [13]. Using Eq. (5), we can rewrite Eq. (4) as, Note that since Θ is uniform inside the TI and vanishes outside it, ∇Θ is proportional to a delta function on the TI surface, Σ T I . From the above equation we can easily read off the surface quantum Hall current, where dS(u, v) = dudv[∂r S (u, v)/∂u × ∂r S (u, v)/∂v] and r S (u, v) ∈ Σ T I , with (u, v) ∈ R 2 being associated to a parametrization of Σ T I . We have that σ H = e 2 Θ/(2πh) is the Hall conductivity. For Θ = π or, more generally, Θ = 2π(n + 1/2), n ∈ Z, we obtain a halfquantized Hall conductivity [9]. Note that a non-vanishing Hall conductivity is equivalent to a non-null surface admittance employed in other studies to describe the conducting surface states of a TI [10].
Given the possibility of a more general quantization profile for Θ, it is natural to wonder whether values significantly larger than the canonical Θ = π one are really possible. This would enhance the sensitivity of experimental probes, since this would allow to compensate against factors of α, especially in those experimental responses leading to effects ∼ O(α 2 ) However, values of Θ beyond n = 0 are not only experimentally unrealistic, they can actually be excluded at a more fundamental level. Since very large values of Θ (∼ 1000) have been considered in the recent literature [10,14], let us briefly comment on this point by recalling some well known facts related to the parity anomaly [15] in this context [7,9].
First we note that the axion term in the Lagrangian (2) can be expressed in terms of a total derivative of a Chern-Simons term, which is more easily seen within a covariant formulation as given in Eq. (1), which can be rewritten as, where we have replaced g 2 by α. Hence, assuming for simplicity a flat TI surface perpendicular to the z-axis separating the TI bulk from vacuum, we obtain the following surface Chern-Simons (CS) contribution to the action, On the other hand, it is well known that integrating out gapped Dirac fermions in 2+1 dimensions coupled to an electromagnetic field generate the CS action above with Θ = π, a result reflecting the so-called parity anomaly [15,16]. Precisely the same gapped Dirac modes arise on a TI surface when time-reversal (TR) symmetry is broken. Therefore, other values of Θ, especially ones significantly larger than π, seem to be excluded on the basis of this argument. As a final remark concerning the value of Θ in the bulk when TR is broken, recall that Θ = π is still enforced provided inversion symmetry is not broken [9]. In the past several years a number of predictions of electromagnetic phenomena related to TIs have been made based on the above equations of axion electrodynamics [9,[17][18][19][20][21][22][23]. Experimentally, such phenomena are often difficult to probe, since some effects manifest themselves only through corrections ∼ O(α 2 ). There are effects occurring at O(α), which are more accessible to experimental probes. This is precisely the case with, for instance, the Faraday and Kerr rotation in TIs [9,18], which has been measured using Bi 2 Se 3 [24] and strained HgTe [25] samples. The derivation of the effect is elementary and follows directly from the Maxwell equations above with the boundary conditions modified by the discontinuity of Θ across the TI surface [9]. The experimental verification of this result using samples of inversion symmetric, epitaxially grown Bi 2 Se 3 [24] showed the robustness of the value Θ = π, as we have discussed above.
Another O(α) effect that has been predicted in the literature concerns surface plasmon polaritons (SPPs) on the planar surface of a TI [19]. In this work the standard analysis of SPPs on a planar geometry has been extended to TIs in the same geometry. As a matter of fact, the same Maxwell equations can be used, except for a nontrivial change in the boundary conditions due to the presence of the axion term in Eq. (2) [19]. Indeed, the TI surface introduces a discontinuity in Θ mixing electric and magnetic boundary conditions even in the static case [17,21,23]. As a concequence, while there is only a transverse magnetic and no transverse electric component of the electric field in the usual theory for SPPs on planar surfaces, such a component does not vanish in the TI case.
One of the most important effects of light scattering beyond the planar geometry originates from the so called Mie theory [26]. This theory deals with plane wave scattering by a dielectric sphere and amounts to solving the Maxwell equations with appropriate boundary conditions [27,28]. Developing the axion Mie theory is severely complicated by the disparity between the spherical symmetry of the target relative to the incoming plane waves; including the axion modifies the boundary conditions relative to the standard Mie theory calculations. Even if certain aspects of the Mie theory for spherical TIs has been considered recently [14], there are many relevant issues that remain to be resolved. This includes an accurate treatment of the modified boundary conditions beyond the first order perturbation level as well as the incorporation of arbitrary sources of fields (i.e., charges and currents) and the identification of suitable experimental setups allowing to measure the implications of the axion term. Another important distinction between the standard Mie theory for a dielectric sphere and the corresponding extension to the spherical TI case is that the latter has a metallic surface where induced (surface) currents feature electrons with the spin-momentum locking property. The latter is actually responsible for the peculiar electromagnetic response whose content is captured by the Lagrangian (Eq. (2)).
In the following we elaborate on these general considerations restricting ourselves to spherical symmetry. This allows to derive analytical solutions for the electromagnetic response of TI spheres, notably including Hall currents, magnetic fluxes, and scattering cross-sections. Our considerations essentially amount to an extension of classical Mie theory to include the ramification of the axion term. While the homogeneous case was already treated using different approaches [10,14], we do not restrict ourselves to the charge-and current-free case and allow for arbitrary external sources of fields. In particular, external charges within the simulation volume as present in charge particle scattering experiments such as electron energy loss spectroscopy (EELS) are not covered by previous studies. Note that exploiting spherical symmetry requires isotropic dielectric functions which is in generally (especially for TIs) not the case. To overcome this constraint the problem needs to be tackled numerically. We will, however, focus somewhat on particular case of localized surface plasmons (LSPs), i.e., negative dielectric response bands. This serves as an archetypal model system for the LSPs on TI nanoparticles of more complicated shape, noting that in particular topological characteristics (i.e. the value of Θ) do not depend on the particular particle shape. For other regimes of positive dielectric response, which may be treated with the same formalism, we refer to the literature.

Solution of the field equations
In order to describe the dielectric response of a TI the curl of the Maxwell Eq. (4) is computed inserting Eq. (5) to replace the magnetic flux density For spatially constant and Θ the standard form holds, The general dielectric response of the sphere (i.e. Mie theory) is obtained by solving Eq. (15) exploring the spherical symmetry. In the following we tackle the problem by a piecewise solution in regions of spatially constant and Θ (i.e., within and outside of the sphere) and fixing the missing integration constants through appropriate boundary conditions. Apart from assuming spherical symmetry, there are no approximations in the expressions above and in this spirit we will continue below when presenting the exact solution of the problem. It is furthermore informative to perform a Helmholtz decomposition of the dielectric response equation at this state as it allows to distinguish a longitudinal and transverse part of the solution on this very fundamental level. We will make use of these results further below. We begin with expanding the electric field into vector spherical harmonics [29][30][31] (following the definition by Barrera et al. [31]), which form a complete basis for vector fields in three dimensions and hence give rise to a general representation of electric and magnetic fields adapted to spherical coordinates. That notably includes free vacuum solutions, which are typically considered in the literature of Mie scattering, but also fields occurring in the presence of non-zero charges and currents. While free solutions with ∇ · D = 0 everywhere may be represented by a restricted set of two-dimensional vector spherical harmonics, which typically correspond to the poloidal and toroidal fields [10,12,27], the general inhomogeneous case requires a complete threedimensional representation as the one employed in this work. Note furthermore that the l = 0 case is peculiar in that both Ψ and Φ vanish identically (ultimately a consequence of the hairy ball theorem). As a consequence the l = 0 modes play a special role throughout. Taking into account the completeness relations of the vector spherical harmonics Eqs., (15) and (16) reduce to a system of three ordinary differential equations for each vector spherical harmonics. After setting E ⊥ lm = E (⊥) lm /r we explicitely obtain for the first component E (⊥) . The second component E (1) is directly linked to the first through the first Maxwell equation Note that the apparent divergence in case of l = 0 dissolves under closer inspection as the Ψ component vanishes. Indeed, l = 0 solutions do not occur in the homogeneous case as the previous two differential equations are not compatibel in this case. The third differential equation for E (2) has again the same mathematical structure as the first Both, the first and third equation, are inhomogeneous ordinary second order differential equation in E (⊥) lm of (modified) spherical Bessel type (after absorbing k √ = kn into r). The fate of the equation being of modified type or not is decided by the sign of of (if the latter is real). While in vacuum is positive corresponding to a spherical Bessel equation, within the sphere both positive and negative sign (depending on the frequency) can occur. The general solution, notably including whispering gallery modes when sphere > vac (positive sign of ) as well as evanescent excitations i.e., surface plasmons (negative sign of ) leads to Bessel functions of complex arguments. Note furthermore that the solution space can be further restricted to those, which do not diverge at the origin. The differential equations above do not contain any modifications due to the topological term. Indeed, being a divergence term, the axion coupling only affects the boundary conditions applying at the interface between vacuum and TI sphere, as anticipated. These internal boundary conditions are derived from partial derivatives normal to the surface inserting Maxwell's equations as usual.
From the homogeneous Maxwell equations (second and third) we obtain and The BCs obtained from the inhomogeneous Maxwell equations (first and fourth) read and Similarly to the field equations only 2 of the 4 BCs are affected by the topological term (i.e., those pertaining to the inhomogeneous Maxwell equations).

Plasmonic Regime
The solutions of the differential Eqs. (18) - (20) allows to treat several problems e.g., scattering of electromagnetic waves or the electromagnetic response of a sphere to external charges and currents. However, we focus on plasmonic excitations which can be excited by plane waves or other external sources. The fundamental condition for plasmonic behavior is an opposite sign of the dielectric function of the nanoparticle with respect to the surrounding medium. As a consequence, we analyse the case of negative dielectric function of the nanoparticle, since we assume the vacuum ( = 1) as a surrounding medium. This condition typically holds for optical and near-infrared frequencies. Representatives of plasmonic TIs at optical frequencies are Bi 2 Se 3 [32], [33] and Bi 2 Te 3 [34].

Homogeneous Case
The general solutions of the above differential equations split into a homogeneous and inhomogeneous part. The homogeneous solution, which includes the important problem of light scattering on a sphere (i.e. Mie scattering), is considered first (for details see Appendix A.1). Even though the homogeneous axion Mie problem has already been treated comprehensively [10,14], it is worth to revisit it here, since we use an alternative set of vector spherical harmonics [31] that is more suitable to solve the full inhomogeneous problem (including arbitrary external charges and currents), and introduce the scattering matrix formalism employed later on. General solutions for the above spherical Bessel differential equations may be generally expanded into spherical Bessel functions j l and y l (first and second kind, respectively). For the radial component (Eq. (18)) this expansion explicitely reads Here the indices i and o correspond to the inner and outer areas of the sphere. Note that the modified Bessel function of second kind does not appear in the interior of the sphere as it diverges at the origin. From Eq. (19) the first tangential component can be derived The solution to Eq. (20) for the second tangential component reads Again, we observe that the radial and first tangential component form a coupled solution, whereas the second tangential component is independent. In the first case the magnetic field (i.e., the curl of the electric field) is strictly tangential to the sphere, whereas in the second case the electric field is tangential. Following the notation of [35] we will refer to the latter as magnetic (m) and the former as electric (e) modes here. Keep in mind that we set lm /r in order to have less headaches solving Eq. (18). For further calculations we need to transform E (⊥) lm and in consequence a Indeed, the coupled radial and first tangential component form the poloidal component of the electrical field (more precisely: of the the electric displacement for which ∇ · D = 0 everywhere), whereas the second tangential component is the toroidal field. Both modes are completely decoupled (and orthogonal), which is exploited in the classical Mie solution by expanding the fields directly in poloidal and toroidal fields from the outset. Note furthermore that in both cases the homogeneous solutions are degenerate with respect to m.
In order to determine the six expansion coefficients (remember that E (1) can be computed once E (⊥) is known) of the homogeneous solutions a corresponding number of boundary conditions has to be provided. Four of them stem from the internal boundaries at the surface of the sphere, noted previously. The remaining two are given either by normalization conditions or external boundary conditions. Indeed, it is only at this stage, where deviations from the classical Mie theory are introduced by the axion term, as anticipated. We first consider the impact of the topological term on the so-called normal modes, i.e., modes decaying to zero at infinity. They are obtained by solving the system of 4 internal BCs (Eq. system (79) for the 4 coefficients b l , c l and fixing the a l by normalization conditions (see Appendix A. From the equations above it is clear that indeed the modifications due to the axion term amount to a mixing of the originally decoupled electric and magnetic modes, which has been previously noted for the half-plane boundary case by Karch [19]. Noting that fields of electric and magnetic modes are perpendicular, their weak mixing is equivalent to a rotation of the electric or magnetic field vectors linear in α, which has been successfully detected in the half plane geometry and may be also amenable to experimental detection in case of the sphere. Fig. 1 shows a comparison between classical electric and magnetic normal modes for l = 1 and the corresponding TI solutions (normalized by setting a ⊥ l or a (2) l to zero, respectively). In the topologically trivial setting (left columns of Fig. 1 respectively) we readily observe that electric and magnetic normal modes are completely decoupled. In case of the TI sphere on the right columns of Fig. 1 respectively, on the other hand, originally purely electric and magnetic modes mix due to the axion term, i.e., they acquire a small magnetic mode and electric mode character respectively and become hybrid modes. Similarly, the 3D representation of the l = 1, m = 1 modes (Fig. 2) reveals a tilting out of the tangential plane of the B-field (Efield) in case of the electric (magnetic) mode, which corresponds to the above noted mixing of electric and magnetic modes. |/α in a. u. R Figure 1: Radial dependency of electric (Eq. (25 -27)) and magnetic (derived using Eq. (5)) field components (absolute value) of l = 1 normal modes for topologically trivial (R = 50 nm, = −1, Θ = 0) and TI sphere (R = 50 nm, = −1, Θ = π) embedded in vacuum respectively at ω = 2 eV. Note that the y-axis was always scaled in the same arbitrary units.
More insight into the origin of that behavior may be obtained from the analysis of the Hall charges (Eq. 23) and Hall currents (Eq. (12) and Eq.(24)) associated to the axion term. Fig. 3 shows the Hall charge and current of the l = 1, m = 1 electric mode. Accordingly, the Hall current has a source and sink due to oscillating charges and produces magnetic fields with components normal to the surface, i.e. a magnetic mode. Moreover, the number of sources and sinks increases with the mode order, e.g., a quadrupolar structure is visible for the l = 2, m = 2 electric mode (see Fig. 3). Slightly reformulating the dependencies between the expansion coefficients and fixing the missing 2 BCs by an external excitation allows the computation of any electromagnetic response of the TI in the absence to charges and currents, i.e., the description of (resonant) photon scattering. As an example we note the solution to the scattering of a linearly po-larized plane wave. Note that due to conservation of angular momentum the whole problem separates in l. In order to obtain an analytic expression of the scattering and extinction cross section (σ s , Eq. (38) and σ e , Eq. (39)) we first express the expansion in the vacuum regions in terms of spherical Hankel functions of first and second kind with excitation coefficients d l = b l −ic l 2 and e l = b l +ic l 2 respectively, which allow a separation of outgoing and incoming spherical waves, respectively. In a second step the outgoing waves are expressed as a linear function of the incoming one (see Eq. system (79)): Noting that we are ultimately interested in the computation of the induced outgoing wave (i.e., total outgoing wave subtracted by outgoing Hankel functions contained in the external plane wave), we finally define a 2x2 scattering matrix t l by d l − e l = t l e l . Terms of the kind xp l (x) can be identified as spherical Riccati-Bessel functions and will be denoted byp l (x) = xp l (x) from now on (here p l represents j i l , j o l , y l , h 1 l and h 2 l ). For clarity we furthermore define the arguments of the Bessel and Hankel functions as ρ i = nkr and ρ o = kr for the inner and outer region of the sphere respectively. The components of the scattering matrix t l than read: respectively. Accordingly, the classical Mie theory result for the scattering matrix, which only contains diagonal entries (see, e.g., [27]), is recovered when dropping the topological term. Note furthermore that we get an additional factor of − 1 2 in the definition of t l in comparison to the well-known classical result, which is a consequence of our choice of the vector spherical harmonics defined by Barrera et al. [31], deviating from those used conventionally. Note, that we already correct for the factor − 1 2 in the Eqs. The scattering cross section σ s and extinction cross section σ e then reads: and (see Eq. (81) and Eq. (83) for the full analytic solution). Accordingly, we observe that the well-known classical Mie theory result (Θ = 0, see, e.g., [27]) is modified by a topological term of the order α 2 in both the nominator and denominator. While the former slightly shifts the energy of the resonances at zeros of the nominator the latter leads to a modification of the scattering intensity. Due to their smallness both are beyond current experimental detection.

Inhomogeneous Case
We finally turn our attention to the general inhomogeneous case. We first note that the solution of the pertaining differential Eqs. (18) and (20) for arbitrary external currents and charges generally does not admit analytical solutions and ultimately requires numerical methods. These consist of (A) expanding the external charge gradients and current sources into vector spherical harmonics, (B) imposing suitable boundary conditions (e.g., Dirichlet boundary conditions at r = 0 and r → ∞), (C) converting the ordinary differential equations into an algebraic system of equations by imposing a suitable grid, and (D) inverting this system to obtain a numerical solution for each E lm . Highly symmetric charge and current distributions (such as the homogeneously charged sphere) may be treated analytically.
Having in mind an experimental setups facilitating a direct measurement of the axion response, we focus on the dielectric response to an external charge moving along a straight line in the following. This setup corresponds to probing the dielectric response of nanostructures by measuring the energy loss and deflection of a highly-focused and highly-accelerated electron beam in the Transmission Electron Microscope (TEM). This Electron Energy Loss Spectroscopy (EELS) technique is widely employed for studying dielectric phenomena such as surface plasmons at the nanoscale and we will come up with a novel detection setup suited to analyze the axion term as a result of the following considerations.
The charge density and current (in frequency space) of an electron moving with velocity and respectively. Here, r 0⊥ denotes the impact parameter in the x − y plane andẑ is the unit vector in z-direction. It is obvious that such an external source defies any spherical symmetry, hence would require the tedious numerical treatment as outlined above. In the following we will therefore adopt several approximations as described in Ref. [35], ultimately leading to analytical expressions for the energy loss and the deflection of the electron beam. These approximations consist of neglecting (A) possible longitudinal solutions (only the projection on the transverse components is considered) and (B) certain boundary terms in the poloidaltoroidal expansion of the transverse solutions [35], which, however, are typically justified for energy losses of fast electrons in the optical regime. The benefit of these restrictions is that one can directly adopt the results of the homogeneous calculations treating the electron as a source of external poloidal and toroidal fields instead as a moving external charge with ρ = 0. The energy loss ∆E suffered by a fast electron is determined by the path integral of the force along the beam path z where we introduced the loss probability Γ loss (ω) corresponding to the EEL spectrum. Repeating the derivation detailed in Ref. [35] and introducing the axion term the latter reads  43)) we observe that the peaks in the loss probabilities, which correspond to zeros in the denominators of the scattering matrix elements, are shifted by terms of O α 2 due to the axion terms in complete analogy to the photon scattering case. Consequently, directly probing the axion term by measuring peak positions in the loss probability is out of range with current TEM-EELS instrumentation. There is, however, a contribution linear in α to the magnitude of the loss probability, which may be detectable. Note, however, that a measurement of such small modifications on top of the total loss probability requires precise reference measurements on identical non-topological samples, which is typically not possible. Similar to the photon scattering case a direct measurement of the off-diagonal terms of the scattering matrix could solve that problem. To that end we consider the lateral deflection of the electron beam given by p kin We now restrict ourselves to a discussion of the lateral deflection component being tangential to the surface of the sphere (see Fig. 5). The logic behind is that such a deflection is not present in the non-topological case. It is generated by a toroidal electric and a poloidal magnetic field, which are both induced by the off-diagonal t 21 coupling term from the dominant poloidal electric field. The static analogue to this field, referred to as a Pearl vortex, has been discussed at length in Ref. [22]. Seeking the computation of the axion-induced toroidal (poloidal) electric (magnetic) field, we only consider d where We finally calculate the transverse momentum per electron at a certain energy loss ω (experimentally energy loss interval around ω) as recorded by EELS. To that end we normalize the spectral representation of the transverse momentum by the number of electrons being inelastically scattered at a particular energy loss.
At this point we finally need computational support in order to calculate the loss probability and the transverse momentum transfer. Fortunately the classical case (Θ = 0) is already implemented within the MNPBEM toolbox of U. Hohenester [36]. However, we need to add the computation of the off-diagonal elements of the scattering matrix (t 12 l and t 21 l ) vanishing in the topological trivial case as well as the axion contribution of the diagonal elements (t 11 l and t 22 l ). After this generalization of the existing implementation we can calculate the loss probability and the lateral deflection for a real scenario. For an electron with 40 keV primary energy a comparison between a topological trivial and a TI sphere is shown in Fig. 4 (loss probability on the left hand side and lateral deflection angle on the right hand side). Here the smallness of the axion contribution to the loss probability manifests itself in almost identical curves for the classical (Θ = 0) and the topological (Θ = π) case. The angle of lateral deflection β vanishing in the topological trivial case was calculated to a maximum value of ≈ 4.5 × 10 −2 µrad for the TI sphere. This angle is at the edge of current detection limits using advanced TEM setups like inelastic momentum transfer (IMT) (see Ref. [37]). This method detects the deflection angle of fast electrons due to induced electromagnetic fields excited by the same electron (see Fig. 5). The calculated value of β is roughly one order of magnitude smaller than the angular resolution achievable with the IMT setup applied in Ref. [37]. Note, however, that the lateral deflection depends (in strictly difference to the value of Θ) strongly on the geometry of the sample and is potentially larger for non-spherical particles, e.g., cuboids exhibiting long straight faces. To calculate the response of arbitrary shaped particles numerical solvers are indispensable. Furthermore the experimental IMT setup can be further improved in terms of increased angular resolution.

Conclusion
In summary, we have extended the description of the electromagnetic response of spherical particles to include the ramifications of an axion term such as present in topological insulators. We have studied the consequences for the excitation of localized surface plasmons on spherical TI nanoparticles. We found an intrinsic mixing of electric and magnetic modes as leading order effect (linear in α). This mixing manifests itself as a small rotation of the induced electric and magnetic fields, which we anticipate to be detectable by suitable setups. Promising candidates utilizing localized probes (i.e. electron beams) are polarization-resolved cathodoluminescence and inelastic momentum transfer measurements [37]. Corresponding shifts in LSP resonance energies and scattering or extinction coefficients, on the other hand, are of the order of α 2 , hence difficult to detect. We also note that the presented formalism allows to include arbitrary external currents and charges, typically not included in the standard formulation of Mie theory. As an example we discussed the dielectric response of a TI sphere transmitted by a focused electron beam as utilized in electron energy loss spectroscopy. The impact parameter of the exciting electron was set to 1 nm at 40 keV primary energy. Note that the loss probability follows almost the same course for the topological trivial and the TI case due to the tiny contribution of the axion term to the loss probability quadratic in α.

A.1 Field equations in spherical coordinates
The general dielectric response of the sphere (i.e. Mie theory) is obtained by solving the response Eq. (15) exploiting spherical symmetry. In the following we tackle the problem by a piecewise solution in regions of constant and Θ (i.e., within and outside of the sphere) and fixing the missing integration constants through appropriate boundary conditions. Accordingly, the response equation in homogeneous regions reads In the next step we expand the solution into vector spherical harmonics (following the definition of Barrera et al. [31]) which gives us later less headaches, when defining the boundary conditions, compared to the more common Hertz potential expansion. Inserting the expansion into Eq. (51) yields after evaluating one curl and finally ∞ l=0 l m=−l l(l + 1) or, upon expanding Eq. (52), ∞ l=0 l m=−l 2 + l(l + 1) For notional simplicity we did not write the vector harmonics expansion on the RHS (inhomogeneous term) of both equations, which is straightforward, however. Taking into account the closure relations of the vector spherical harmonics both equations reduce to a system of three ordinary differential equations for each vector spherical harmonics. Note furthermore that only the first two are coupled, respectively, whereas the last one can be solved independently. We proceed by first solving the coupled system of the first two. Subtracting the Ψ lm expansion coefficients of both equations gives 1 r Inserting into the Y lm expansion coefficient of Eq. (55) then leads to l(l + 1) We finally compute the derivative of the E lm term by inserting Eq. (58) obtaining (62) This is an inhomogeneous second order differential equation in E (⊥) lm of (modified) spherical Bessel equation type (after absorbing k | | = kn into r). The third equation has the same structure lm .
Both sets of solutions are independent. Note furthermore that the solution space can be further restricted by only admitting solutions, which do not diverge at the origin (i.e. modified spherical Bessel functions of the first kind in the interior). The homogeneous solution to the first expansion then explicitely reads l y l (kr), and from Eq. (58) The second solution reads l y l (kr), In the first case the magnetic field is strictly tangential to the sphere, whereas in the second case the electric field is tangential. In accordance with the convention employed for the plane boundary (i.e. sphere radius → ∞) they are called magnetic (m) and electric (e) modes. Note furthermore that in both cases the homogeneous solutions are degenerate with respect to m.

A.2 Boundary Conditions
In order to determine the 6 expansion coefficients of the homogeneous solutions above we require 4 BCs (two for E (⊥) and two for E (2) ) for each homogeneous domain, E (1) can be computed once E (⊥) is known. The two remaining degrees of freedom are fixed by normalization of the electric and magnetic modes respectively. We begin with the internal boundaries at the surface of the sphere. Indeed, it is only at this stage, where deviations from the classical Mie theory are introduced by the axion term as we anticipated. These boundary conditions are derived from partial derivatives normal to the surface as usual: 1. Boundary condition from the second Maxwell equation and hence 2. Boundary condition from the third Maxwell equation 3. Boundary condition from the first Maxwell equation (r = outward unit normal) where the last line is obtained after inserting the continuity of the E lm components derived from the second boundary condition Eq. (70) and using Θ 2 = 0 (topological trivial material outside of the sphere) and Θ 1 = Θ. 4. Boundary condition from the fourth Maxwell equation (t 1,2 = tangent unit vector into Ψ lm and Φ lm direction) and hence The solution of the latter equation system (Eq. (32) and Eq. (33)) directly leads to the components of the scattering matrix t l defined by d l − e l = t l e l (Eqs. (35)- (37). The resulting complete expressions for the scattering cross section σ s and extinction cross section In the second step we inserted the radial dependency (outside of the sphere) of the outgoing toroidal solution and inserted the definition of the vector spherical harmonic field Φ = r × ∇Y lm . Moreover, the excitation coefficient has been denoted by d (2) lm as in the main text. Now, the partial derivatives can be carried out yielding The magnetic deflection into the y-direction is due to the poloidal magnetic field, which is obtained from the toroidal E m field according to and hence B m The second term in the bracket can be directly derived from the electric deflection by partial integration of (87) The definition of the A lm prefactors from Ref. [35] read which may be solved analytically [35].