## Abstract

We present a rigorous derivation of the weak-form formulation of the Helmholtz equation for electromagnetic structures incorporating general nonreciprocal, dispersive materials such as magnetized ferrites or magnetized free-carrier plasmas. This formulation allows an efficient self-consistent treatment using finite elements of a variety of problems involving magnetic or magneto-optical materials biased by an external DC field where the eigenmodes become nonreciprocal or even unidirectional. The possibilities of this method are illustrated with several examples of TE-polarized modes at microwave frequencies and TM-polarized modes at optical and infrared wavelengths.

© 2017 Optical Society of America

## 1. Introduction

Nonreciprocal electromagnetic components such as isolators, circulators and phase shifters often rely on biased magnetic or magneto-optic materials to operate—either through direct coupling between the oscillating field and the spins of the magnetic medium, as is typically the case at microwave and radio frequencies, or through spin-orbit coupling in the optical regime. A rigorous description of these interactions is an extraordinarily difficult task because it implies solving for two distinct sets of equations in a self-consistent way—those describing the evolution of the spins subject to the bias and time-varying fields and those that characterize the evolution of the time-varying field itself. Solving such coupled equations is a central task of spintronics [1]. Because the magnetic layers used in spintronic architectures are typically only a few nanometers thick, it is possible in many cases to slightly reduce the complexity of the problem by treating the oscillating field as a static quantity [2]. In electromagnetism, an even stronger approximation can be made because one is less concerned about the actual motion of the spins and more concerned by the behavior of the oscillating field. Moreover, the wavevectors involved in electromagnetism are often smaller than those for which the interactions between the spins themselves (e.g. through exchange and anisotropy…) become non-negligible. In most practical cases, then, it is reasonable to decouple the two sets of equations in the following manner [3]. First, the precession of biased and non-interacting spins subject to circularly polarized waves is calculated for an infinite volume of material. Then, the resulting response is homogenized and turned into a macroscopic quantity (a permeability tensor when the spins directly interact with the fields or a permittivity tensor in the case of spin-orbit coupling). Finally, this macroscopic parameter is assumed to be independent of the actual size of the material (but it can be corrected by depolarization factors if the bias fields inside and outside the sample are different). The initial problem is thus reduced to solving Maxwell’s equations in a system where the response of the magnetic or magneto-optic material is entirely defined by a tensorial non-symmetric bulk permeability *µ* or permittivity *ε*, and respectively termed gyromagnetic and gyroelectric materials. Note that this nonsymmetric character of *µ ^{T}*≠

*µ*and/or

*ε*≠

^{T}*ε*is a necessary (but not sufficient) condition for nonreciprocal behavior in linear, time-invariant problems [4].

Despite this crude simplification, nonreciprocal electromagnetic calculations are still challenging. There exists a number of systems where analytical expressions have been derived, as is the case for instance for the bulk and surface waves guided by 2D ferrite slabs at microwave frequencies [5–9] or the magnetoplasma wave phenomena in magnetized free-carrier systems [10]. Most configurations, however, can only be investigated numerically and a variety of home-made, open-source and proprietary electromagnetic codes have been written and/or used for this purpose. Among them, commercial packages based on the finite element method (FEM) have become particularly widespread in the community but typically lack a crucial functionality—the possibility of solving for the eigenmodes of a nonreciprocal system at an arbitrary frequency. Indeed, in a standard frequency domain FEM formulation of eigenmode problems, the eigenvalue solved for is the frequency *ω* for a chosen (parametric) spatial wavector $\overrightarrow{k}$ — a Bloch phase in a periodic structure or a propagation constant *k _{0}n_{eff}* in an uniform waveguide. When the materials present become dispersive, the parametric vectorial eigenvalue problem, $\widehat{L}(\overrightarrow{k};\omega )\overrightarrow{E}={(\omega /c)}^{2}\epsilon (\omega )\cdot \overrightarrow{E}$(with $\widehat{L}=\nabla \times {\mu}^{-1}(\omega ).\nabla \times $) becomes strongly nonlinear and increasingly hard to solve in a self-consistent manner. In particular, for nonreciprocal gyrotropic materials, dispersion cannot be neglected, as the origin of the gyrotropy is almost always connected to a resonant transition in the light-matter interaction – resonant RF absorption near the ferromagnetic resonance frequency in Larmor precession [3] or the Zeeman splitting of the circularly polarized light absorption at infrared and visible frequencies [11]. This limitation of FEM methods for dispersive materials has already been tackled for isotropic dispersive structures such as plasmonic crystals or metamaterials by transforming the nonlinear ω-eigenvalue problem into a quadratic

*k*-eigenvalue problem with ω as a parameter [12]. As a result, material dispersion is readily taken into account.

It is the goal of this paper to extend such an approach for two-dimensional periodic systems containing dispersive nonreciprocal media using the weak form of the finite element method. The resulting *k*-quadratic eigenvalue equation will be shown to contain linear terms in *k*, thereby breaking the inversion symmetry around the *Γ*-center of the Brillouin zone and thus also breaking time reversal symmetry. In section 2, the geometry and eigenvalue equation to be solved are presented. In section 3, a practical implementation in Comsol Multiphysics is described for modeling the nonreciprocal waves supported by a grounded ferrite slab at microwave frequencies. This is a well-known problem that admits an analytical solution and we show that the numerical results produced by our model are in quantitative agreement with the theory. We then study a more complex problem in the same frequency range to demonstrate the possibilities of the model for configurations with no tractable theoretical expressions. In section 4, the method is benchmarked on an extension of a recently proposed concept for unidirectional plasmonic waveguide [13,14]. Magnetized small-gap semiconductors, such as InSb, allow strong cyclotron resonances at reasonably low magnetic fields due to the high mobility and small effective mass of their free carriers. Using the numerical method developed here, we show how this feature can be exploited to create a periodic surface guide sustaining only one-way bands at mid- and far-infrared wavelengths. Concluding remarks appear in section 5.

## 2. Method

Two classes of electromagnetic phenomena involving magnetic or magneto-optic materials can be defined—those that involve a bias field H_{0} parallel to the direction of propagation,such as the Faraday effect, and those for which H_{0} is normal to the propagation. The present article specifically focuses on the second category of problems which comprises a wealth of effects ranging from nonreciprocal surface waves at microwave frequencies [4–8] to the transverse magneto-optical Kerr effect (TMOKE) at optical and infrared wavelengths in guided and free-space configurations [15–19].

Figure 1 shows the generic problem under investigation: the two-dimensional domain is flanked by periodic boundary conditions along the vertical walls while the top and bottom boundaries may be either periodic or perfectly conducting depending on the problem considered (if needed, perfectly matched layers may be inserted before these outer boundaries to emulate an open system). It is also assumed that the system contains at least one magnetic or one magneto-optic element subject to a bias field along the *z* direction.

We first consider the modes that propagate along the horizontal *x* axis and that directly interact with the spins of the magnetic material (a generalization for arbitrary propagation directions and magneto-optical materials will be given later). This situation typically occurs for ferrites and garnets near their gyromagnetic resonance, i.e. in the microwave regime. Following the prescriptions discussed in the introduction, the relative permeability *µ* of the magnetic material can be approximated by the following tensor [3]:

*x*-

*y*plane are sensitive to the nonreciprocal anisotropy of

*µ*. Said otherwise, we only need to consider TE-polarized modes. Note in passing that for

*µ*and

_{r}*κ*purely real, the permeability tensor is Hermitian and the system is lossless. The anti-symmetric off-diagonal elements

*κ*are proportional to the external magnetic field. Due to Onsager's principle, this proportionality has odd character and

*κ*therefore changes sign upon reversal of the magnetization [20,21]. The wave equation satisfied by the electric field $\overrightarrow{E}=\widehat{z}.{E}_{z}$at all points of the domain is:

Note that the values of *ε* and$\tilde{\mu}$depend on the coordinates *x* and *y*. In the regions occupied by the magnetic material, the elements of$\tilde{\mu}$are obtained by inverting Eq. (1). In all other regions, the non-diagonal components *ξ* and *κ* vanish while *η* and *µ _{r}* become unity. Developing Eq. (2) leads to the following scalar equation:

Because we apply periodic boundary conditions on the vertical walls of the computational domain (Fig. 1), the solutions of Eq. (4) can be expressed as:

In this equation, $\overrightarrow{k}$is the in-plane momentum of the mode. We shall now transform Eq. (4) into a weak-form expression in which the in-plane momentum amplitude *k* is the eigenvalue by generalizing the method outlined by Davanco et al. for isotropic and non-magnetic media [12]. Note that the normal, *y*-component of the momentum is implicitly contained in *u(x,y).* This latter component is not a conserved quantity of the system and therefore cannot be used as an eigenvalue. Inserting Eq. (5) into Eq. (4), and multiplying all the terms by a test function *w(x,y)* according to the Galerkin procedure [12,22], the scalar wave equation becomes:

Here it can be noted that, as expected, the presence of the terms containing the gyrotropy *ξ* breaks the reversal symmetry of the eigenvalue problem. Only upon simultaneous sign reversal of *ξ* → –*ξ*, would the reversed system, $\tilde{k}=-k$and$\tilde{x}=-x$, have an eigenvalue that is equal to the original one. As the sign of *κ,* and thus *ξ*, is governed by the sense of the external field, one can conclude that the backward eigenmode in these systems is given by the forward eigenmode under reversed magnetization.

The values of *η* and *ξ* are known everywhere. However, we generally do not know the values of their spatial derivatives. A straightforward way to remove these unknown quantities is to integrate by parts over the domain *Ω*. If we consider the second to last term of Eq. (6), for example, we can formally write:

*y*and

_{min}*y*are the

_{max}*y*coordinates of the bottom and top outer boundaries, respectively. In our model, these boundaries are either periodic or perfectly conducting (Fig. 1). In both cases, the first term of the integration by parts vanishes, either because the integrals over

*x*are equal at

*y*and

_{min}*y*(periodicity), or because

_{max}*w*(

*x*,

*y*) is set to zero at

*y*and

_{min}*y*(perfect conductor/Dirichlet boundary conditions, which may be placed behind perfectly matched layers to emulate open systems). A similar reasoning applies for the spatial derivatives of

_{max}*η*and

*ξ*with respect to

*x*; in this case also, the first term of the integration by parts vanishes because the system is periodic along

*x*. After expressing

*η*and

*ξ*as a function of the material parameters

*μ*and

_{r}*κ*, we finally recast Eq. (6) as:

*k*is the eigenvalue. Equation (8) is our final weak-form quadratic eigenvalue expression, valid for TE-polarized waves interacting with magnetic media characterized by the

*µ*tensor given by Eq. (1). Following the same steps, this equation can be generalized to the case of a domain with periodic conditions applied to all boundaries (e.g. a 2D photonic crystal) and for an arbitrary propagation direction in the

*x*-

*y*plane:

Furthermore, Eqs. (8) and (9) can also be formally used in optics to describe the magnetic field $\overrightarrow{H}=\widehat{z}.{H}_{z}$ of TM-polarized modes interacting with transversely magnetized 2D magneto-optical media having a relative gyro-electric permittivity of the form:

This fundamental duality is obviously a direct consequence of the symmetry in Maxwell's equations, $\overrightarrow{E}\to \overrightarrow{H},\overrightarrow{H}\to -\overrightarrow{E}$and *µ$\to $* *ε*. Apart from a global phase factor *π* on either the electric or the magnetic field, any solution to Maxwell's equations is therefore also a solution for the dual problem with the permeability and the permittivity interchanged. In this precise case, it amounts at replacing *µ _{r}* and

*κ*by ${\epsilon}_{\perp}$ and

*g*and also permuting

*ε*and

*µ*in Eqs. (8) and (9). To summarize, the method outlined here can be used to solve for TE-polarized modes at frequencies close to the gyromagnetic resonance of biased magnetic media and TM-polarized modes in magneto-optics. Equations (8) and (9) can be readily discretized and solved using the finite-element method, as will be shown in the next sections where practical examples are discussed and implemented in the commercial package Comsol Multiphysics.

## 3. Nonreciprocal waves at microwave frequencies

Here we examine for benchmarking purposes the surface and volume waves of biased ferrimagnetic slabs at microwave frequencies. Such modes have been extensively studied since the late 1960s due to their rich properties (existence of several families of solutions, possibility of subwavelength confinement, nonreciprocal propagation) which are relevant for a variety of nonreciprocal microwave components.

The first example considered in this section is the grounded and lossless ferrimagnetic slab. When biased by an external field normal to the *x*-*y* plane, the TE-polarized modes propagating along its top and bottom interfaces obey the following dispersion relation [8]:

*β*and

_{mag}*β*are the normal wavevector components in the magnetic material and in air, respectively,

_{air}*d*is the thickness of the slab, and all the others parameters have been defined previously.

To test the validity of our approach, we have modeled this structure and verified that the eigenmodes follow the dispersion predicted by Eq. (11). A schematic of the computational domain is shown in Fig. 2(a). To utilize the weak-form expressions derived above, a periodic unit cell must be defined along the propagation (*x*-) direction even though the problem is invariant along this axis. This requirement does not constitute a problem as long as the period is sufficiently small to prevent any artificial band folding in the frequency range of interest. The YIG slab is directly in contact with a Dirichlet boundary condition to emulate a ground plane while the opposite boundary can be of any type provided that it is sufficiently far away from the YIG interfaces along which the (non-leaky) modes propagate. The material parameters are specified by defining a wave equation for each subdomain. Here we use Eq. (8) with the specific values of *ε, µ _{r}*, and

*κ*satisfied by the YIG slab and the air region. Because several terms of Eq. (8) vanish for nonmagnetic media, the wave equation in the latter subdomain reduces to the expression already derived by Davanco et al. for isotropic and reciprocal media in [12]. The discretization and numerical calculations can be implemented with a variety of finite-element solvers; in this study, we utilize the weak-form PDE environment of Comsol Multiphysics to simulate the structure. In Comsol vernacular, the default field name is u, its partial derivatives are ux and uy while the test function

*w*(

*x*,

*y*) is test(u). We also make use of the operators ut and utt respectively defined in terms of the eigenvalue

*k*as ut = -

*k*u and utt =

*k*

^{2}u. Moreover, the integral notation is implicit in the weak-form PDE environment and therefore only the integrands of Eq. (8) need to be entered in the model. The mesh generation procedure is the same as for a non-magnetic model because all the gyromagnetic effects are captured by material parameters characterized by local constitutive relations.

Figure 2(b) compares the dispersion relations given by Eq. (11) with our numerical simulations for a YIG slab with a thickness of 1 mm (the other parameters are given in the caption). An excellent agreement is found between the two sets of data. The modes supported by the structure are of two types. The low-frequency branches correspond to the so-called magnetostatic (M) modes and have originally been discovered in the magnetostatic approximation [5]. These surface waves result from the coupling between the electromagnetic field and a collective oscillation of the spins, implying that they can only form in a frequency range close to the gyromagnetic resonance where the YIG permeability components acting on TE-polarized waves have negative values. The high-frequency branches are known as dynamic, or ferrite-dielectric (FD) modes [6–8]. They can be seen as guided modes within the ferrite slab, except that their properties are heavily modified by the anisotropic permeability—in fact, the FD modes can be either surface waves or volume waves depending on the geometry, the operation frequency and the YIG material dispersion [7,8]. In this specific example, they are volume waves since their field pattern within the YIG region does not reach their maximum at the interface and does not decay exponentially, see Fig. 2(c). In contrast, the field pattern of the M modes is exponential on both sides of the interface, as shown on Fig. 2(d).

The different plots of Fig. 2 are clearly nonreciprocal: the dispersion diagram is asymmetric with respect to the origin *k* = 0 and the forward and reverse solutions have distinct field profiles. This behavior can be understood by noting that, had the YIG slab been bounded by air on both sides, the forward and reverse modes would have traveled along opposite interfaces because of the anisotropy of the permeability tensor *µ* [8]. Here, in contrast, the ground plane breaks the spatial symmetry of the problem and all solutions are pushed on the bottom air side, considerably perturbing those that would have otherwise been located around the top interface and turning the system nonreciprocal as a result.

Before leaving this first example, we note that the field patterns of Figs. 2(c) and 2(d) have been calculated both analytically and numerically. One can see that the curves generated by the two methods are essentially indistinguishable, further validating the accuracy of our FEM model. We have verified that this good agreement was not coincidental by performing systematic benchmarking tests on this system as well as on other configurations that admitted theoretical expressions (result not shown here).

We now turn to a system that cannot be treated analytically. It is a variant of the previous configuration where the perfect ground plane has been replaced by a periodic chain of Cu cylinders with a finite conductivity [Fig. 3(a)]. Moreover, magnetic losses typical of polycrystalline YIG have been incorporated in the ferromagnetic layer so that the overall structure is far more realistic than the previous example. The presence of the periodic cylinders induces diffraction effects that couple the surface modes of the YIG slab to propagating waves in free space. It is thus necessary to apply perfectly matched layers on the top and bottom walls of the domain to properly take into account any leakage in the air regions, as schematically represented on Fig. 3(a).

The dispersion diagram of the different eigenmodes of the structure is represented in Fig. 3(b). Here again, the modes propagating in the forward and reverse directions are nonreciprocal because they are not equally perturbed by the presence of the cylinders.

Apart from this similarity with the previous example, the band diagram exhibits much more complex features. In particular, the branches of the M modes are now folded due to two distinct effects. A first folding occurs because of the periodicity at *k* = ± π/*P*. Note, however, that the forward and reverse modes do not reach these points at the same frequency: as a result, the diagram does not exhibit the familiar standing wave pattern typical of reciprocal structures and the branch is translated to the other edge of the Brillouin zone rather than appearing merely flipped. The second type of branch folding is directly related to the losses of the ferrimagnetic slab (these features disappear if losses are removed from the model) and occur for M branches that approach the frequency given by (ω_{0} + 0.5ω_{m})/2π ≈3.85 GHz, which corresponds to the asymptote of the lossless magnetostatic mode in the absence of the metallic cylinders [7]. For example, the red curve experiences such an upturn at *k* = 252 m^{−1}. Above this frequency, the mode experiences an anomalous dispersion, in much the same way as lossy surface plasmons in optics [23]. In practice, the modes in the anomalous dispersion regime are extremely lossy, rendering their observation virtually impossible. This fact can be seen on Fig. 3(c) where the propagation lengths of the different modes are represented.

Interestingly, the two FD branches of Fig. 3b occur at much lower frequencies than those of the grounded ferrite slab of Fig. 2. More generally, it is possible to tune their lower cutoff frequency by playing with the geometrical parameters such as the wire radius and distance of the Cu wires to the YIG slab. In certain cases, they can even coexist with the M modes, opening up interesting opportunities to design highly nonreciprocal subwavelength structures as shown in [24] where we have demonstrated a one-way glass metamaterial constructed from the same periodic unit cell (but with the periodicity perpendicular to the propagation direction).

## 4. Infrared gyro-electric one-way band structures

A second benchmarking example is concerned with an extreme example of nonreciprocity, namely one-way guiding. Such a behavior is characterized by the non-existence of real energy propagating solutions (neither guided nor radiating) to Helmholtz' equation in one direction, as if the system possessed a one-way photonic band gap. In solid-state physics, such a behavior for electrons is known as topological insulation and is characterized by the appearance of a surface electrical conductivity without dissipation or backscattering, even in the presence of large impurities, on a system that is otherwise insulating in its bulk [25,26]. About a decade ago Haldane and Raghu predicted a photonic analogue in certain classes of transversely magnetized photonic crystals (PhCs) [27,28] and its subsequent successful observation has spurred ever since an intense study into topological photonic systems [29,30]. In such systems, the phase diagram (or thus the Brillouin zone) contains gapless zones where topological phase transition boundaries occur between bands of different integer Chern number. The most common example of such phase boundaries are pairs of degenerate Dirac cones in the reciprocal space of 2D PhCs [27]. Proper symmetry breaking (*in casu* time reversal breaking *T* by magneto-optics) opens a gap between these bands. The integer quantization of the topological invariant of the bands of the PhC forces the system to create unidirectional edge modes at the physical boundary of such a crystal in order to transition between the integer jumps of the Chern number. These edge modes are topologically protected and are therefore backscattering immune against edge imperfections or perturbations.

The necessary key feature in such topological unidirectional systems is the availability of strong *T*-breaking to create a bulk gap that itself is sufficiently immune against disorders. Indeed, while the edge states are topologically protected, the bulk gap itself is not governed by a topology invariant. Up till now such a strong *T*-breaking has not been realized at optical or infrared frequencies.

Other *T*-breaking one-way waveguides and mirrors have therefore been proposed that however are not protected by topology. One such an appealing system is formed by a surface plasmon-polariton mode (SPP) on a transversely magnetized free-electron plasma that has a dielectric PhC in its close vicinity. The transverse magnetization breaks the time-reversal symmetry and thus lifts the degeneracy of the forward and backward SPP bands. As a result, if the SPP modes are guided at frequencies within the PhC gap or thus if the surface plasma frequency lies within this gap, the splitting of the forward-backward SPP's leads to a band of unidirectional surface modes. The robustness of this one-way guiding (even though not topologically protected) is governed by the strength of the gap of the dielectric PC, which is independent of *T*-breaking. It can be made very strong by optimal geometrical tuning.

This system has first been proposed by Yu et al. [13] exploiting the nonreciprocity of the cyclotron effect on the free carriers in a transversely magnetized plasma. A classical Drude-Lorentz description of the carrier motion in such a magnetoplasma [10] leads to a nonreciprocal gyro-electric permittivity where the off-diagonal elements express the Hall conductivity of the plasma induced by the Lorentz force on the free carriers:

*ω*, Γ and

_{P}*ω*=

_{B}*eB*/

*m**(with

*m**the effective mass of the carriers) are respectively the plasma frequency, the Drude collision frequency and the cyclotron frequency. The diagonal term ${\epsilon}_{\infty}\overleftrightarrow{I}$collects possible high-frequency interband transitions in the background of the plasma. The prefactor of the tensor is the standard Drude term and the tensor contains the gyro-electric corrections of the Drude-Lorentz description of the cyclotron effects. It evaluates correctly to the unity tensor if the magnetic field is switched off (

*ω*= 0). The transverse magnetic field is as always applied in the

_{B}*z*-direction as in the generic layout of Fig. 2(a) [with the garnet replaced by the magnetoplasma described by Eq. (12)]. The permittivity in the

*z*direction is therefore unaffected.

The dispersion relation of the TM-polarized surface plasmon polariton (*E _{z}* = 0) guided at the surface between a dielectric with permittivity

*ε*and this nonreciprocal magnetoplasma is given by [31]:

_{d}*ik*). Note that this configuration is again that of Fig. 2(a), except that the YIG medium is replaced by the nonreciprocal magnetoplasma. The time reversal symmetry in this dispersion equation is effectively broken by the presence of a term that is odd both in

_{0}n_{eff}x*n*and the gyro-electric gyrotropy

_{eff}*g*. This implies correctly that the backward propagating SPP is degenerate with the forward one if and only if the external magnetic field

*B*is simultaneously switched, i.e.

*ω*→ –

_{B}*ω*and thus

_{B}*g*→ –

*g*. The asymptotic surface plasmon frequencies of these nonreciprocal SPPs are found by the limits

*n*→ ± ∞ of Eq. (13), or thus by solving for

_{eff}*ω*. Using the definitions of Eq. (12) for

*g*and ${\epsilon}_{\perp},$this leads to (assuming infinitesimal losses Γ = 0) ${\epsilon}_{d}+{\epsilon}_{\perp}\mp g=0$and:

For small cyclotron frequency with respect to the plasma frequency, this implies that a frequency band of Δ*ω*_{oneway} ≈*ω _{B}* is opened around the reduced plasma frequency ${\omega}_{P}/\sqrt{{\epsilon}_{d}+{\epsilon}_{\infty}}$where the magneto-surface plasmons can only propagate in the backward sense (if

*g*> 0 and thus

*B*> 0). Due to the magnetic field induced reciprocity breaking, it is important to note that in the above form, Eq. (13) is only applicable for the configuration of Fig. 2(a). A spatial inversion of the geometry in Fig. 2(a) would place the magnetoplasma on the bottom but with unchanged permittivity tensor because of the axial character of the magnetic field, $\widehat{P}\overrightarrow{B}\to \overrightarrow{B}$ (and therefore unchanged

*ω*). As $\widehat{P}$changes the sign of the square roots in Eq. (13), a forward guided SPP in Fig. 2(a) will have the same dispersion as the backward guided SPP in the inverted configuration. In other words, if Fig. 2(a) does not allow a forward SPP between ${\omega}_{P}/\sqrt{{\epsilon}_{\infty}+\epsilon}-{\omega}_{B}/2$ and ${\omega}_{P}/\sqrt{{\epsilon}_{\infty}+\epsilon}+{\omega}_{B}/2$, then the inverted configuration will not allow a backward SPP.

_{B}The simple analytical treatment of this nonreciprocal magnetoplasma problem provides an ideal benchmark test for the nonreciprocal FEM recipe of Section 2. Figure 4 compares the solution of the analytical dispersion relation given by Eq. (13) with the one obtained by our numerical finite element method for the case of a lossless magnetoplasma described by a Drude-Lorentz dispersion (12) with a high background permittivity *ε _{∞}* = 15.4 and exhibiting cyclotron resonances at

*ω*0.01

_{B}=*ω*. The numerical solution (continuous line) of the quadratic finite element eigenvalue equation in the Bloch vector

_{P}*k*[Eq. (8), with

*µ*,

_{r}*κ*and

*ε*replaced by ${\epsilon}_{\perp},$

*g*and

*µ*of the corresponding materials] is in perfect agreement with the analytical solution (markers)

*k*of Eq. (12). The FEM model assumes a Bloch mode. In order to avoid artificial folding of the modes of this uniform guide, we have used a period$P\approx {\lambda}_{P}\sqrt{{\epsilon}_{\infty}+1}/100,$leaving even strongly confined SPPs with

_{0}n_{eff}*n*up to 10 or more unperturbed by the artificial periodicity of the FEM model. The predicted one-way gap with a width equal to the cyclotron frequency and the theoretically predicted broken forward and backward surface plasmon frequencies ${\omega}_{fw,bw}^{SP}$are correctly observed. The calculation is done in units normalized with respect to

_{eff}*ω*. The results can therefore be transposed to any surface magnetoplasmonic system that can be described by the above Drude-Lorentz gyroelectric dispersion.

_{P}The chosen values for the Lorentz-Drude model parameters correspond to experimentally fitted values obtained from MO-FTIR measurements (between 50 and 7500 cm^{−1}) on undoped InSb wafers [32]. This small-gap high-mobility III-V semiconductor is known to be described in very good approximation by a plasma of free electrons. The contribution of the holes can in first instance be neglected due to their much higher effective mass. A plasma frequency of *ω _{P}* = 296 cm

^{−1}and a Drude lifetime 1/

*Γ*= 1.9 ps has been obtained. Fitting the cyclotron frequency

*ω*= 23.4 cm

_{c}^{−1}in a known magnetic induction of B = 0.42 T has led to an experimental value of

*m**= 0.0168

*m*for the electron effective mass. It is precisely this very low effective mass (and high collision lifetime) that makes these small-gap semiconductors ideal for THz and mid-infrared nonreciprocal applications: as compared to standard plasmonic materials such as Ag and Au (where the effective electron mass approaches 1), significantly lower magnetic fields are needed to obtain the same gyro-electric strength.

_{0}The surface magnetoplasmon polariton guide is not strictly one-way because of the presence of the bidirectional radiation continuum in the one-way range. The original idea of Yu and associates [13] has been to suppress this by operating the one-way SPP guide in the band gap of a photonic crystal cladding. This will suppress the air light cone, while the InSb bulk plasmon cone only starts at${\omega}_{P}/\sqrt{{\epsilon}_{\infty}}+{\omega}_{c}/2$. Contrary to the previous paragraph, such a system cannot be described analytically and therefore represents an ideal example to check the validity of the proposed numerical model. We consider the configuration of Fig. 5(a). Five periods of a square photonic lattice of dielectric beams (*ε _{d}* = 8.9) are sandwiched between two parallelly magnetized InSb claddings. A thin air gap,

*t*= 0.2

_{g}*λ*, separates the PC from the InSb magnetoplasmonic cladding. The periodicity

_{P}*a*= 1.492

*λ*of the crystal and the thickness

_{P}*t*= 0.8

*a*of the beams is optimized to create a maximal photonic bandgap around the central frequency of the one-way regime, ${\omega}_{P}/\sqrt{{\epsilon}_{\infty}+1}$. For the chosen values, a 20% TE gap is obtained. As explained above, the reduced one-way regime is roughly given by Δ

*ω*=

_{oneway}*ω*/

_{B}*ω*. For reasonable values for the applied magnetic field (B ≈0.125 T), this will not exceed more than 1% of the plasma frequency. The chosen crystal parameters therefore sufficiently isolate the one-way SPP modes from the Bloch radiation continuum of the PhC cladding. The air gap has not been particularly optimized, but chosen to allow sufficient overlap for the SPP modes to feel the presence of the bandgap cladding. Note in passing that without airgap, the boundary conditions of the SPP dispersion would profoundly change, pushing the one-way regime to too low frequencies due to the substitution 1 ↔

_{P}*ε*.

_{d}Using the gyrotropic FEM recipe of Section 2 on the considered hybrid PhC/ magnetoplasmonic guiding structure between 0.241ω_{P} ≤ ω ≤ 0.254ω_{P}, Fig. 5(b) plots the obtained non-evanescent TM Bloch bands in the reduced first Brillouin zone. All other (nonplotted) bands are highly lossy and evanescent modes corresponding to the gap modes of the PhC. The carrier scattering losses in the InSb plasma have been neglected (*Γ* = 0) in Eq. (12). The plotted Bloch bands are therefore purely real and their propagation sense is given by the sign of the bands' slope. The dashed (resp. continuous) bands should therefore correspond to modes that have a negative (resp. positive) Poynting flux in the *x*-direction over any arbitrary *yz*-plane in the periodic guide.

The peculiar behavior of the bands is a direct proof of the numerical validity of the here developed calculation method. With the notations of Fig. 4, the mode calculation of the single magnetoplasmonic InSb/air interface of Fig. 4 predicts the absence of a forward guided SPP mode at the bottom interface above *ω*_{SPP,-} (and likewise, by the given symmetry considerations, the absence of a backward SPP at the top interface above ω_{SPP,-}). The *H _{z}* [ =

*u(x,y)*exp(-

*jkx*)] field profiles calculated for the two bands at

*ω*= 0.25

*ω*[Fig. 5(d)] in this one-way frequency range, confirm that the backward (resp. forward) band is guided at the bottom (resp. top interface). Deriving the frequency domain Poynting vector in the presence of gyrotropic media [34] leads to:

_{P}*u*is the periodic part of the Bloch function (see Eq. (5) and

*k*is the Bloch vector. Evaluating this expression additionally proves that the bottom (resp. top) SPP propagates energy in the backward (resp. forward) direction. The evaluation of this Poynting map along the vertical dashed cut-line [Fig. 5(c)] shows for both cases that its net flux is indeed either negative or positive [35]. Both interfaces are correctly decoupled by the PhC, as the calculated bands cross without interaction. Finally, it is worth noting how the proposed recipe also correctly leads to the appearance in the Brillouin zone of an accumulation of increasingly flat bands when approaching the surface plasma frequencies ω

_{B}_{SPP, ±}. This is directly due to band folding over 2π/

*a*of the increasingly confined SPPs.

Besides being a convincing benchmark of the nonreciprocal gyro-electric FEM recipe, the above example also proposes a very compact design for an integrated isolating and even circulating circuit at an operation frequency that can be tuned by “choosing” a suitable magnetoplasma. As proven in [13] and [14], adding a cavity inside the PhC couples the upper and lower one-way MO SPP guides and by suitably choosing the magnetic field to be parallel or antiparallel, the circulation sense can be set. The precise example given here (with an undoped InSb plasma frequency *f _{P}* ≈9 THz) could serve as such a type of circulator for far infrared and THz devices (needing only limited magnetic field strengths).

## 5. Conclusion

In conclusion, we have presented a numerical approach to solve for the nonreciprocal eigenmodes of a variety of gyromagnetic and gyroelectric problems. Since we have focused on the electromagnetic aspects of these phenomena, we have made the usual approximation of decoupling the equation of motion of the spins from Maxwell’s equations. It is worth mentioning in this regard that the full, self-consistent resolution of the problem has already been successfully tackled with the finite element method for ferromagnetic thin films, but within the magnetostatic limit [36]. Thus, our work contributes to enrich an emerging finite-element toolbox dedicated to the study of nonreciprocal waves and other spin-related effects. In particular, this recipe is straightforwardly applicable to the variety of photonic topological structures that are currently attracting huge attention.

## References and links

**1. **I. Žutić, J. Fabian, and S. Das Sarma, “Spintronics: Fundamentals and applications,” Rev. Mod. Phys. **76**(2), 323–410 (2004). [CrossRef]

**2. **D. D. Stancil, *Theory of Magnetostatic Waves* (Springer-Verlag, 1993).

**3. **D. M. Pozar, *Microwave Engineering* (John Wiley & Sons, 2012).

**4. **D. Jalas, A. Petrov, M. Eich, W. Freude, S. Fan, Z. Yu, R. Baets, M. Popović, A. Melloni, J. D. Joannopoulos, M. Vanwolleghem, C. R. Doerr, and H. Renner, “what is – and what is not – an optical isolator,” Nat. Photonics **7**(8), 579–582 (2013). [CrossRef]

**5. **R. W. Damon and J. R. Eshbach, “Magnetostatic modes of a ferromagnet slab,” J. Phys. Chem. Solids **19**(3-4), 308–320 (1961). [CrossRef]

**6. **L. Courtois, G. Declercq, and M. Peurichard, “On the non-reciprocal aspect of gyromagnetic surface waves,” AIP Conf. Proc. **5**, 1541–1545 (1972). [CrossRef]

**7. **P. de Santis, “Dispersion Characteristics for a Ferrimagnetic Plate,” Appl. Phys. (Berl.) **2**(4), 197–200 (1973). [CrossRef]

**8. **T. J. Gerson and J. S. Nadan, “Surface Electromagnetic Modes of a Ferrite Slab,” IEEE Trans. Microw. Theory Tech. **22**(8), 757–763 (1974). [CrossRef]

**9. **R. E. Camley, “Nonreciprocal surface waves,” Surf. Sci. Rep. **7**(3-4), 103–187 (1987). [CrossRef]

**10. **E. D. Palik and J. K. Furdyna, “Infrared and microwave magnetoplasma effects in semiconductors,” Rep. Prog. Phys. **33**(3), 1193–1322 (1970). [CrossRef]

**11. **S. Visnovsky, *Optics in Magnetic Multilayers and Nanostructures* (CRC, 2010).

**12. **M. Davanco, Y. Urzhumov, and G. Shvets, “The complex Bloch bands of a 2D plasmonic crystal displaying isotropic negative refraction,” Opt. Express **15**(15), 9681–9691 (2007). [CrossRef] [PubMed]

**13. **Z. Yu, G. Veronis, Z. Wang, and S. Fan, “One-way electromagnetic waveguide formed at the interface between a plasmonic metal under a static magnetic field and a photonic crystal,” Phys. Rev. Lett. **100**(2), 023902 (2008). [CrossRef] [PubMed]

**14. **V. Kuzmiak, S. Eyderman, and M. Vanwolleghem, “Controlling surface plasmon polaritons by a static and/or time-dependent external magnetic field,” Phys. Rev. B **86**(4), 045403 (2012). [CrossRef]

**15. **V. I. Belotelov, I. A. Akimov, M. Pohl, V. A. Kotov, S. Kasture, A. S. Vengurlekar, A. V. Gopal, D. R. Yakovlev, A. K. Zvezdin, and M. Bayer, “Enhanced magneto-optical effects in magnetoplasmonic crystals,” Nat. Nanotechnol. **6**(6), 370–376 (2011). [CrossRef] [PubMed]

**16. **L. Halagačka, K. Postava, M. Vanwolleghem, F. Vaurette, J. Ben Youssef, B. Dagens, and J. Pištora, “Mueller matrix optical and magneto-optical characterization of Bi-substituted gadolinium iron garnet for application in magnetoplasmonic structures,” Opt. Mater. Express **4**(9), 1903–1919 (2014). [CrossRef]

**17. **W. Śmigaj, L. Magdenko, J. Romero-Vivas, S. Guenneau, B. Dagens, B. Gralak, and M. Vanwolleghem, “Compact optical circulator based on a uniformly magnetized ring cavity,” Photonics Nano. Fund. Appl. **10**(1), 83–101 (2012). [CrossRef]

**18. **L. Bi, J. Hu, P. Jiang, D. H. Kim, G. F. Dionne, L. C. Kimerling, and C. A. Ross, “On-chip optical isolation in monolithically integrated non-reciprocal optical resonators,” Nat. Photonics **5**(12), 758–762 (2011). [CrossRef]

**19. **S. Ghosh, S. Keyvaninia, W. Van Roy, T. Mizumoto, G. Roelkens, and R. Baets, “Adhesively bonded Ce:YIG/SOI integrated optical circulator,” Opt. Lett. **38**(6), 965–967 (2013). [CrossRef] [PubMed]

**20. **L. Onsager, “Irreversible processes,” Phys. Rev. **37**(4), 237–241 (1931). [CrossRef]

**21. **L. D. Landau and E. M. Liffschitz, “The symmetry of the kinetic coefficients” in S*tatistical Physics, vol. 5 Course in Theoretical Physics* (Pergamon, 1980).

**22. **J. Jin, *The Finite Element Method in Electromagnetics* (Wiley, 2002).

**23. **H. J. Lezec, J. A. Dionne, and H. A. Atwater, “Negative refraction at visible frequencies,” Science **316**(5823), 430–432 (2007). [CrossRef] [PubMed]

**24. **A. Degiron and D. R. Smith, “One-way glass for microwaves using nonreciprocal metamaterials,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **89**(5), 053203 (2014). [CrossRef] [PubMed]

**25. **L. Landau, “Diamagnetismus der Metalle,” Z. Phys. **64**(9-10), 629–637 (1930). [CrossRef]

**26. **K. V. Klitzing, G. Dorda, and M. Pepper, “New method for high-accuracy determination of the fine-structure constant based on quantized hall resistance,” Phys. Rev. Lett. **45**(6), 494–497 (1980). [CrossRef]

**27. **F. D. M. Haldane and S. Raghu, “Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry,” Phys. Rev. Lett. **100**(1), 013904 (2008). [CrossRef] [PubMed]

**28. **S. Raghu and F. D. M. Haldane, “Analogs of quantum-Hall-effect edge states in photonic crystals,” Phys. Rev. A **78**(3), 033834 (2008). [CrossRef]

**29. **Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljacić, “Observation of unidirectional backscattering-immune topological electromagnetic states,” Nature **461**(7265), 772–775 (2009). [CrossRef] [PubMed]

**30. **L. Lu, J. D. Joannopoulos, and M. Soljacic, “Topological photonics,” Nat. Photonics **8**(11), 821–829 (2014). [CrossRef]

**31. **L. Halagačka, M. Vanwolleghem, K. Postava, B. Dagens, and J. Pištora, “Coupled mode enhanced giant magnetoplasmonics transverse Kerr effect,” Opt. Express **21**(19), 21741–21755 (2013). [CrossRef] [PubMed]

**32. **J. Chochol, K. Postava, M. Čada, M. Vanwolleghem, L. Halagačka, J.-F. Lampin, and J. Pištora, “Magneto-optical properties of InSb for terahertz applications,” AIP Adv. **6**(11), 115021 (2016). [CrossRef]

**33. **A. K. Zvezdin and V. A. Kotov, *Modern magnetooptics and magnetooptical materials* (CRC, 1997).

**34. **J. D. Jackson, *Classical Electrodynamics* (Wiley,1999).

**35. **Note also how inside the plasma the energy flux is locally reversed as is known from plasmonic systems.

**36. **M. Mruczkiewicz and M. Krawczyk, “Nonreciprocal dispersion of spin waves in ferromagnetic thin films covered with a finite-conductivity metal,” J. Appl. Phys. **115**(11), 113909 (2014). [CrossRef]