1. Introduction

Magneto-optical (MO) or gyromagnetic materials are a class of rare-earth iron garnets such as yttrium iron garnet (YIG). MO effect, which will happen when electromagnetic (EM) waves propagate through magnetized medium, has been widely studied and explored since Faraday effect was discovered in 1845 [1]. This phenomenon will cause polarization rotation when light passes through a gyromagnetic material, which originates from the non-diagonal tensor form of magnetic permeability in MO materials. Due to the nonreciprocal nature of Faraday effect, it has been utilized for application in optical isolators and circulators [2–4]. Further, magneto-optical Kerr effect reveals the change of reflection light from a magnetic surface [5], which has taken essential role in optical data storage. In addition, the studies of Zeeman effect and magnetic birefringence like Cotton-Mouton effect or Voigt effect [6–8] have enriched the MO effect and all these phenomena have shown the close relationship between magnetic field and light.

With the development of nanofabrication technology, the combination of MO materials with subwavelength photonic devices has becoming a promising field. The nanostructures built by MO materials have shown peculiar properties [9]. Recently, there were reports of applying MO effect in metamaterials, such as magneto-plasmonic [10–13], MO subwavelength gratings [14–16], and MO resonators [17–19], which are promising in fabricating novel nonreciprocal devices. There were also reports of utilizing graphene-based MO plasmonic for designing nonreciprocal device in mid-infrared region [20,21]. These works have shown that the interaction between light and magnetization at nanoscale reveals not only plentiful physics properties but also the potential for different applications. The MO effect has also been used to break time-reversal symmetry for realizing topological states in photonics [22]. Topological waveguide, which has potential in unidirectional propagation and be immune against defects, have received more and more attention recently and could be applied in a wide range of wave propagation and manipulation devices [23–27]. Nonetheless, the physical scheme behind these unique phenomena is not clear and further theoretical investigations are required. In our recent work, we developed a simple theory for calculating the electromagnetic waves transport between the interface of isotropic and MO media [28]. We also demonstrated optical Lorenz force to explain the topological photonic states and robust unidirectional transport behaviors [29].

The MO devices are majorly waveguides or transmission lines containing MO materials. For unidirectional propagation, the most common model relies on the Faraday rotation in MO medium and 45° crossed polarizers. Recently, there were reports of utilizing Faraday effect for designing nonreciprocal waveguide devices, based on nonreciprocal mode conversion [30,31]. The investigation of modes in waveguides containing MO materials is the first step of theoretical understanding of their electromagnetic propagation behavior [32]. However, even under simplified conditions, the accurate solution of anisotropic dielectric waveguide could face formidable boundary value problems. Over the years, people increasingly rely on the use of various numerical simulation methods [33,34], for example, the field distribution in the cross section of waveguide can be accurately calculated by finite element software. However, it is time-consuming to find suitable structure parameters. Besides, perturbation theory can be utilized for simplifying the mode analysis, when the anisotropic or inhomogeneous properties can be treated as the perturbation added on the isotropic medium [35]. It has been shown that under perturbed situation, the approximate analysis can also get results with enough precision. Nonetheless, it is no longer suitable when the external magnetic field is strong or MO effect is obvious enough. Therefore, it is still necessary to develop analytical or semi-analytical theoretical solution method, which has greater advantages for people to deeply understand the interaction between light and matter, conduct mode analysis, and predict the evolution of electromagnetic field, etc., and usually consumes less calculation time.

In this work, we considered a simple planar MO waveguide model. The MO layer is surrounded by isotropic media and the external magnetic field is added in the propagation direction. By setting the suitable form of the electromagnetic components, the analytical solution of MO mode can be solved through a transcendental equation. Then, the analytical solutions of MO modes can give us a complete picture of dispersion curves of modes, transport behavior and energy flows in the waveguide. In Section 2, we present the derivation process of EM waves transport in the MO waveguide from Maxwell’s equations. Then in Section 3, the perturbation theory is introduced for understanding the propagation behaviors of MO modes. It is shown that the MO modes can be treated as a coupling of TE and TM modes. Next, specific MO mode and transport analysis is discussed in Section 4. We have revealed that mode degeneracy and anti-crossing phenomena could happen in MO waveguide, which is unable to be observed in isotropic waveguide. Besides, we have observed nonreciprocal Faraday rotation during the propagation. The analytical solutions can give detailed information of dispersion, and energy transport of MO modes, which can be very useful for setting structure and material parameters in MO waveguide design.

2. Model and formalism for magneto-optical planar waveguide

The planar waveguide model we consider here is schematically displayed in Fig. 1 with three homogeneous medium. The core layer is set as magneto-optical (MO) medium, and it is surrounded by isotropic cladding layer and substrate layer. In this configuration, we set the relative permittivities and relative permeabilities of substrate layer, core layer and cladding layer as ${\varepsilon _1}$, ${\mu _1}$, ${\varepsilon _2}$, ${\hat{\mu }_2}$, ${\varepsilon _3}$, and ${\mu _3}$. For gyromagnetic MO medium, the permeability has tensor form with non-diagonal elements when external magnetic field is applied. Assuming that a magnetic field is added in the z direction, the relative permeability ${\hat{\mu }_2}$ can be written as

 

Fig. 1. The schematic diagram of planar waveguide model. (a) The cladding, core, and substrate layers of planar waveguide, which is infinitely extended in the y and z direction. (b) The x-z cross section of the planar waveguide model. The magnetic field is added in the z direction, as shown by the red arrows.

Download Full Size | PDF

The EM wave transmitting in the core layer can be described by the source-free Maxwell’s curl equations as below,

Here $\mathbf{E}$ and $\mathbf{H}$ represent the electric and magnetic field. We suppose that the planar waveguide extends infinitely both in the y and z direction. Thus, the EM field vectors $\mathbf{E}$ and $\mathbf{H}$ of guided mode are only confined in the x direction. Considering that the mode propagates along the z direction, as displayed in Fig. 1, then $\mathbf{E}$ and $\mathbf{H}$ can be written as

Since that the slab waveguide structure, permittivity, and permeability maintain unchanged along the y-axis, we have $\partial E/\partial y = 0$ and $\partial H/\partial y = 0$. The relation between the six EM components of $\mathbf{E}(x)$ and $\mathbf{H}(x)$ (i.e. ${E_x}$, ${E_y}$, ${E_z}$, ${H_x}$, ${H_y}$, and ${H_z}$) can be obtained by combing Eqs. (2) and (4). When in non-magnetic medium, the EM components could be separated into two isolated groups of equations as transverse electric (TE) mode and transverse magnetic (TM) mode. For example, in the substrate layer, the TE mode satisfies the following wave equation [36],

Similarly, the wave equation of TM mode could be written as below,

However, these six components could not be independent in MO medium. By combining Eq. (1) and Eq. (2), the coupling relationship of the EM components in the core layer can be obtained by the following set of formulas,

Here, TE and TM modes are no longer existing since these six EM components cannot be separated into two groups, i.e., ${H_y}$ and ${H_x}$ can be connected through ${\mu _k}$. By expressing all the rest EM components in Eq. (9) with ${H_y}$ component, the wave equation can be expressed as follow,

Usually, we have ${\mu _z} = 1$ so that Eq. (10) can be simplified as,

The propagation constants and EM fields of the guided modes can be obtained by solving Eq. (10) or Eq. (11). Considering that the EM fields would be mainly confined in the core layer, and decay exponentially into the substrate and cladding layer, the complete form of H_y in these three layers can be written as follows,

To solve the dispersion equation, boundary condition for electromagnetic fields should be considered as constraints. Herein, it means that ${E_y}$, ${E_z}$, ${H_y}$, and ${H_z}$ components are continuous at the boundaries. The rest EM components in the core layer can be obtained from ${H_y}$ through Eq. (9). For example, the relation between ${E_y}$ and ${H_y}$ in the core layer is,

Notice that ${E_y}$ component is isolated from ${H_y}$ in the non-magnetic material, which belong to TE and TM components separately. Expecting that ${E_y}$ would be evanescent in the cladding and substrate layers, the solutions of ${E_y}$ have the form as

Here, the undetermined coefficients ${A_2}$ and ${D_2}$ are different from ${A_1}$ and ${D_1}$. ${E_z}$ and ${H_z}$ components in the core layer can be obtained through ${E_z} ={-} \frac{1}{{i\omega {\varepsilon _0}{\varepsilon _2}}}\frac{{\partial {H_y}}}{{\partial x}}$ and ${H_z} = \frac{1}{{i\omega {\mu _0}{\mu _z}}}\frac{{\partial {E_y}}}{{\partial x}}$. Therefore, we have

According to the boundary conditions, ${E_y}$, ${E_z}$, ${H_y}$, and ${H_z}$ should be continuous at x = 0 and x = d. Then there will be eight equations connecting all eight undetermined coefficients. The undetermined coefficients ${B_1}$, ${C_1}$, ${B_2}$, ${C_2}$, ${D_1}$, and ${D_2}$ could be replaced by ${A_1}$ and ${A_2}$, and the boundary conditions will be simplified as

Finally, the dispersion equation can be obtained by eliminating the undetermined coefficients ${A_1}$ and ${A_2}$. And the result is the following transcendental equation as

The propagation constants or eigenvalues $\beta $ can be determined from the transcendental equation Eq. (22). If the external magnetic field is removed, the relative permeability ${\mu _2}$ in the core layer will return to isotropic state, which means ${\mu _r} = {\mu _z} = 1$ and ${\mu _k} = 0$. In this situation, the EM waves propagating in the core layer will be separated into TE and TM waves. The wavenumbers ${k_{21}}$ and ${k_{22}}$ in core layer will degenerate into ${k_2} = \sqrt {k_0^2{\varepsilon _2}{\mu _2} - {\beta ^2}} $. Equation (20) will be separated into two isolated equations as $f_1^{\textrm{TM}} = 0$ and $f_3^{\textrm{TE}} = 0$, where $f_1^{\textrm{TM}} = (\frac{{{k_1}{\varepsilon _2}}}{{{\varepsilon _1}}} + \frac{{{k_3}{\varepsilon _2}}}{{{\varepsilon _3}}})\cos ({k_2}d) + ( - {k_{2}} + \frac{{{k_3}\varepsilon _2^2{k_1}}}{{{\varepsilon _3}{\varepsilon _1}{k_2}}})\sin ({k_2}d)$, and $f_3^{\textrm{TE}} = (\frac{{{k_3}{\mu _2}}}{{{\mu _3}}} + \frac{{{k_1}{\varepsilon _2}}}{{{\varepsilon _1}}})\cos ({k_2}d) + (\frac{{{k_3}{\varepsilon _2}{\mu _2}{k_1}}}{{{\varepsilon _1}{\mu _3}{k_2}}} - {k_{2}})\sin ({k_2}d)$. Thus, the dispersion equations of TM waves and TE waves can be obtained separately as

The propagation modes in the external magnetic field configuration could not be separated into TE and TM modes as in the case of the homogeneous medium, so we mark these special modes in the MO waveguide as ${m_0}$, ${m_1}$, ${m_2}$, and so on for convenience. Through the dispersion equations, Fig. 2 shows the graphical method to obtain the propagation constant $\beta $ of the symmetrical slab waveguide. Here, we simply set the relevant parameters as $f = 5$GHz, $d = 2$cm, the core layer is consisted of YIG, and the cladding and substrate layers are set as air, where ${\varepsilon _1} = {\mu _1} = {\varepsilon _3} = {\mu _3} = 1$. When no external magnetic field is exerted, the EM waves propagating in the core layer are same with the ordinary TE and TM modes. Figure 2(a) shows the plots for calculation of TE modes, where the blue lines represent $\tan ({k_2}d)$ and red lines represent the right part of Eq. (23b) as ${F_{TE}}$, and their intersections correspond to solutions to the dispersion equation. And in the case of an applied magnetic field ${H_0} = 2000$G, the curves for calculating the MO modes are displayed in Fig. 2(b). The blue lines in Fig. 2(b) are the curves of function ${F_1}({k_{21}}d) = {f_1}/{f_2}$ and the red lines represent ${F_2}({k_{21}}d) = {f_3}/{f_4}$. The black points marked in the figure correspond to different planar modes in Fig. 2. In order to obtain accurate solution, numerical methods such as Newton-Raphson method or the bisection method are necessary. According to the expression of Eq. (12), the wave numbers ${k_1}$, ${k_{21}}$, ${k_{22}}$, and ${k_3}$ should be real numbers, otherwise the EM fields cannot be confined in the core region. And from these, the cut-off condition of the MO guided modes can be determined.

 

Fig. 2. Display of numerical method for the dispersion equation when the frequency is 5 GHz. Black dots where red and blue lines intersect correspond to the eigen guided modes of different orders. (a) The calculation of TE modes, without external field. (b) The calculation of magneto modes in YIG waveguide, with the external field ${H_0} = 2000$G.

Download Full Size | PDF

3. Comparison of perturbed and accurate solutions

The EM field distributions of MO modes are very different from non-magnetic TM and TE modes. For better understanding the propagation behaviors of the MO modes, we first consider the case of weak MO response, that is, when the intensity of applied external magnetic field is weak or the operating frequency is far from the resonance region of the MO material. Under this circumstance, the MO modes can be viewed as the perturbation of basic TE and TM modes supported in the isotropic waveguide without external magnetic field, and sufficiently accurate results can also be obtained through perturbation analysis, compared with the exact solution. Such perturbation analysis theory of weak MO effect can be extended to the design of MO devices in terahertz and infrared bands which have attracted much attention in recent years.

The off-diagonal elements in the permeability ${\hat{\mu }_2}$ dominate the interaction of TE and TM modes. Therefore, we suppose that ${\mu _k}$ is relatively small compared to ${\mu _r}$ and utilize perturbation method for mode analysis. The perturbation system can be handled by using of Rayleigh-Ritz variational technique, which is referred to Ref. [35,37]. Regardless of the physical nature of real materials, we would discuss the perturbation situation in a mathematical sense and compare it with the accurate solutions. Combined with perturbation theory, the Maxwell equations can be written into the form of solving eigenvalue problem [38]:

When written in terms of electromagnetic components, the expression of Eq. (28) will be

The anisotropy of the medium will manifest in the operator $\mathrm{{\cal L}}$. Here, $\mathrm{{\cal L}}$ can be rewritten as $\mathrm{{\cal L}} = {\mathrm{{\cal L}}_0} + L$. The homogeneous part is set as ${\mathrm{{\cal L}}_0}$, and the perturbed operator L contains only the off-diagonal elements of the permeability tensor. The parameter $\sigma $ is introduced to assist the perturbation terms, and will be set as 1 in the final result. Considering only the first order of the perturbed system, we have $\mathrm{{\cal L}} = {\mathrm{{\cal L}}_0} + \sigma L$, $\mathbf{\Phi } = {\mathbf{\Phi }_0} + \sigma {\mathbf{\Phi }_1}$, and $\beta = {\beta _0} + \sigma {\beta _1}$. The variational method in the first order approximation will lead to the following form as

Therefore, the EM field in the perturbation system can be regarded as the coupling of different TE and TM modes, according to Eqs. (29) and (30). And the coupling terms of TE and TM modes contain off-diagonal component ${\mu _k}$. For simplicity, we discuss a case where the perturbation system mainly includes the interaction of two basic TE and TM modes, and their coupling with other modes can be ignored. The coupling between TM (or TE) modes are forbidden since that no off-diagonal components terms in Eq. (29) support it. The perturbation modes can be written as,

Suppose that the coupling is between TE and TM modes, the parameters in Eq. (32) are

Then the propagation constants of perturbation modes can be obtained as

Herein, we would compare the calculated results among isotropic TE and TM modes, the perturbation modes, and accurate solutions. We set ${\varepsilon _2} = 15.26$, ${\mu _r} = 1$, the value of ${\mu _k}$ changes from 0.01 to 0.03, and the cladding and substrate layers are set as air. The element ${\mu _k}$ is relatively small so that the perturbation theory can be applied. The solid lines in Fig. 3 are the dispersion curves of TE_m and TM_m modes where m = 0 and 1. The effective index of mode in the waveguide is defined as ${n_{eff}} = \beta /{k_0}$. The perturbation modes are calculated by the coupling of TE₀ with TM₀ mode or TE₁ with TM₁ mode, which are shown as the dashed lines in Fig. 3. Note that, the coupling of TE_m and TM_m modes with different m-order are negligible. As seen in Fig. 3, the dispersion curves of perturbation modes are close to that of the isotropic TE_m or TM_m modes, so we mark the perturbation solutions as quasi-TE_m (qTE_m) and quasi-TM_m (qTM_m). In Eq. (34), ${\widetilde \beta _ + }$ and ${\widetilde \beta _ - }$ belongs to qTE_m modes and qTM_m modes, respectively. For the qTE_m modes, the TE_m components remain dominant and the TM_m components are the perturbation. And vice versa, the TM_m components still dominate in the qTM_m modes.

 

Fig. 3. Comparison of dispersion diagram for perturbed MO modes, accurate MO modes, nonmagnetic TE and TM modes. The TE and TM modes are solid lines, while the perturbed modes are marked as dash lines. The accurate MO modes are displayed as chain lines. (a) ${\mu _k} = 0.01$. (b) ${\mu _k} = 0.03$.

Download Full Size | PDF

In addition, the introduction of ${\mu _k}$ enlarges the disparities between qTE_m and qTM_m with same m-order, and they will not converge like TE and TM modes as the normalized width $d/\lambda $ grows. With the growing of ${\mu _k}$, the disparities between qTE_m and qTM_m will be larger. When ${\mu _k} = 0.03$, qTM₀ and qTE₁ modes will have same effective index at particular $d/\lambda $, as displayed in Fig. 3(b) at the intersection of red and purple dash lines, which shows the degeneracy of these two modes. The first four accurate MO modes are also plotted as chain lines in Fig. 3 for comparison. When $d/\lambda $ is small, the difference of effective index among isotropic TE and TM modes, perturbation modes, and accurate modes is negligible. In this case, the perturbation theory could simplify the calculation of MO modes and still achieve sufficient precision. When $d/\lambda $ is increased, the disparities of effective indexes between MO modes and perturbation modes will increase. It should be noticed that the accurate MO modes could be cutoff at larger $d/\lambda $, which corresponds to the situation of higher frequency or larger width of core layer. However, the perturbation modes will not have cutoff since isotropic TE and TM modes can still exist in this case. Therefore, the perturbation theory will not be applicable when $d/\lambda $ become larger. Besides, when ${\mu _k}$ is increased from 0.01 to 0.03, the effective index of perturbation modes become more deviated from the accurate ones, showing that perturbation theory has enough precision only in weak MO effect situation.

From the above analysis, when no magnetic field exists, the propagation modes in the core layer remain TE_m and TM_m form. When the external magnetic field is imposed, there will be coupling between TE and TM modes and each MO mode will have six EM components. The dispersion curves of MO modes are close to the isotropic TE and TM modes if the MO effect is weak, where the perturbation theory can be adopted for fairly good precision analysis. In weak MO effect configuration, the TE_m (TM_m) components remain dominant in qTE_m (qTM_m) modes. When the intensity of external magnetic field become larger, the MO modes will differ greatly from the TE_m (TM_m) modes, showing that the coupling between TE and TM modes is stronger and more complex. For general situation, the MO modes could be solved through the dispersion equation Eq. (22).

4. Crossing, avoided crossing, and dispersion curves of MO modes in general situation

The coupling of TE and TM modes in the MO waveguide will show unique phenomena that could not be found in isotropic waveguide. The analytical solutions of MO modes allow us to perform a comprehensive and detailed analysis of the EM transmission in the core layer. Regardless of real materials, we set ${\mu _r}$ and ${\mu _k}$ to have specific value and ignore their dispersion with frequency change, and the relative permittivity of core layer is fixed as ${\varepsilon _2} = 15.26$. The substrate and cladding layers are set as air, and their relative permittivity and relative permeability are 1. All these parameters remain unchanged in the following discussion. The dispersion curves of MO modes are displayed as solid lines in Fig. 4 when ${\mu _r} = 1$ and ${\mu _k} = 0.5$, and TE_m and TM_m modes are presented as dashed lines as comparison. Comparing the dispersion curves of m₀ and m₁ modes with TE₀ and TM₀ modes, we find that the curve of m₀ mode is close to TE₀ mode and m₁ mode is close to TM₀ mode when the effective index remains at a low value. Hence, m₀ and m₁ modes originate from TE₀ and TM₀ modes; similarly, m₂ and m₃ modes come from TE₁ and TM₁; while m₄ and m₅ modes are related with TE₂ and TM₂. However, the strong coupling of TE and TM modes makes the MO mode dispersion much different from the isotropic modes. As seen in Fig. 4, the curves of m₁ and m₂ modes, m₃ and m₄ modes undergo crossing, where the intersection shows the degeneracy of modes. Besides, when the dispersion curves of m₁ and m₄, m₄ and m₅ modes are getting closer, the slope of these curves will change dramatically and keep the curves of these modes away from each other. In other words, the curves of m₁, m₄ and m₅ modes will avoid crossing with each other. Unlike in isotropic medium, in which all the TE and TM modes are orthogonal and the overlap between two different modes is always zero, the crossings and avoided crossings of dispersion curves of MO waveguide are the results of introducing the off-diagonal term ${\mu _k}$.

 

Fig. 4. The dispersion curves of MO modes (solid lines) and nonmagnetic TE and TM modes (dash lines).

Download Full Size | PDF

As one step further, we investigate the formation mechanism and influence of crossings and avoided crossings on the EM wave transport in the waveguide. After the propagation constant β is obtained through Eq. (22), the EM field distribution can be given once the coefficients A₁, A₂, B₁, C₁, B₂, C₂, D₁, and D₂ are determined. These coefficients can all be expressed by constant A₁, and A₁ can be determined when we specify the optical power P carried by the waveguide. The power P is expressed as

This is the power of light in the z-axis propagation direction, for per unit length along the y-axis. Here, we set the total power P = 1 W for determining the constant A₁ of each MO mode and the field distribution in the waveguide can be obtained. In order to investigate the influence of degeneracy when two dispersion curves cross, we mark out the crossing region and show the H_y fields of m₁ and m₂ modes in Fig. 5. The width d of core layer is fixed as 2 cm, and the frequency we calculate is mainly at gigahertz region since that the MO effect is intense in microwave band. As seen in Fig. 5(a), m₁ and m₂ modes will be degenerate around 2.7 GHz, meaning that these two modes will have same propagation constants. In this special case, these two modes will stay synchronized as they propagate in the waveguide along z-axis. In Figs. 5(b) – 5(e), we show the H_y field of m₁ and m₂ when f is 2.6 and 2.8 GHz, respectively, where the core layer region is ranging from 0 to 2 cm. As the frequency increases and passes the degeneracy point, the effective index of m₂ mode will exceed that of m₁ mode, while their field distributions remain similar in this process. This means that the nature of each mode remains unchanged, and no interaction happens between the crossing modes. On the other hand, we explore the change in field distribution when avoided crossing phenomenon occurs. Figure 6(a) shows the avoided crossing dispersion region of m₁ and m₄ modes, and the center point of avoided crossing is marked as a star where the frequency is 5.503 GHz. The avoided crossing will change the slopes of the dispersion curves and expand the difference in propagation constants between these two MO modes. In Figs. 6(b) and 6(e) when f = 4.8 GHz, the m₁ and m₄ modes are away from the avoided crossing region and their H_y field distributions have not been affected by each other. Then when f = 5.3 GHz, the effective index of m₁ and m₄ become closer, and their field distributions start to change, as seen in Figs. 6(c) and 6(f). The crest of H_y field in m₁ mode become trough in the center of core layer, and the trough between two peaks in m₄ mode is lifted up. In the right side of avoided crossing region when f = 5.6 GHz, in addition to the trough in the center of core layer, two small peaks appear at both ends of the m₁ mode, while the trough of m₄ mode is much elevated, as displayed in Figs. 6(d) and 6(g). Obviously, mode coupling, transformation and energy exchange occur between these two MO modes in the avoided crossing region.

 

Fig. 5. (a) The crossing region of MO modes. The red line is the dispersion curve of m₁ mode and the green line belong to m₂ mode. (b)-(c) The H_y field distributions of m₁ mode when (b) f = 2.6 GHz, (c) f = 2.8 GHz. (d)-(e) The H_y field distributions of m₂ mode when (d) f = 2.6 GHz, (e) f = 2.8 GHz.

Download Full Size | PDF

 

Fig. 6. (a) The avoided crossing region of MO modes. The red line is the dispersion curve of m₁ mode and the cyan line belongs to m₄ mode. (b)-(d) The H_y field distributions of m₁ mode when (b) f = 4.8 GHz, (c) f = 5.3 GHz, (d) f = 5.6 GHz. (e)-(g) The H_y field distributions of m₄ mode when (e) f = 4.8 GHz, (f) f = 5.3 GHz, (g) f = 5.6 GHz.

Download Full Size | PDF

Note that, TE_m and TM_m modes with numbers m = 0, 2, 4… can be classified as even modes while the modes with m = 1, 3, 5… can be categorized as odd modes. Since m₁, m₂, and m₄ modes originate from TM₀, TE₁, and TE₂ respectively, the crossing and avoided crossing of these modes can be regarded as the results of the coupling relationship between odd and even modes. In the crossing cases between m₁ and m₂ modes, the coupling between odd TM and even TE is negligible small, thus the modes just degenerate at the frequency of crossing point. On the contrary, the m₁ and m₄ modes both emerges from even modes, TM₀ and TE₂ modes respectively. Hence, there can be strong coupling between them, due to the significant overlap of their mode fields. Therefore, as the frequency increases and the effective index of these two modes approach each other, the mode coupling and energy exchange between them increases, and then the avoided crossing phenomenon occurs. In general situation, the avoided crossings of MO modes in dispersion curves will happen when odd (even) TE modes are strongly coupled with odd (even) TM modes, which leads to the obvious variation of their slope in dispersion curves. It means that for coupled waves in the avoided crossing region, the group velocity will undergo a sharp change.

In the next step, we delve into the effects of parameters ${\mu _r}$ and ${\mu _k}$ on dispersion, and the dispersion curves of modes in waveguide with different ${\mu _r}$ and ${\mu _k}$ are calculated and the results are shown in Fig. 7. In Fig. 7(a), when ${\mu _r} = 1$ and ${\mu _k} = 0$, there is no external magnetic field applied to the waveguide and the EM waves can be separated into TE (blue lines) and TM (red lines) modes. When ${\mu _r} = 2$ and ${\mu _k} = 0$, the core layer could be treated as uniaxial medium, and the dispersion curves are shown in Fig. 7(b). In this case, the EM waves propagating in it can still be separated into TE and TM modes, while the propagation constants of TM_m modes will exceed that of TE_m modes with same m order, which result in degeneracy of TE and TM modes. The dispersion curves of MO modes when ${\mu _r} = 1$ and ${\mu _k} = $0.1, 0.3, and 0.5 are displayed as blue, green, orange and purple lines in Figs. 7(c) – 7(e), respectively. In these cases, the MO effect is intense. The grey lines indicate the upper limit and lower limit of the effective index, namely the cutoff conditions of MO modes. When ${\mu _k} = 0.1$, the external magnetic field is modest and only crossings of MO modes arise. With the increasing of ${\mu _k}$, the avoided crossings in the dispersion curves can be observed, showing that stronger external magnetic field leads to more complex and strong coupling between TE and TM modes. Besides, as ${\mu _k}$ gets closer to ${\mu _r}$, the upper limit of effective index goes down. Note that, when ${{(\mu _r^2 - \mu _k^2)} / {{\mu _r}}} \le 0$, there will be no propagation modes allowed in the MO medium because the propagation constants become imaginary in this case. In Fig. 7(f), when ${\mu _r} = 2$ and ${\mu _k} = 0.5$, the dispersion curves of MO modes are more complicated, which can be viewed as the combinations of dispersion curves in Fig. 7(b) and Fig. 7(e). The analytical calculation of MO modes allows us to give the dispersion and field distribution variations of the modes in the planar waveguide from the case of no magnetism to strong external magnetism. In particular, the crossing and avoided crossing regions can be identified and understood, which is useful in MO waveguide design.

 

Fig. 7. The dispersion curves with different μ_r and μ_k. (a) μ_r = 1, μ_k = 0. (b) μ_r = 2, μ_k = 0. (c) μ_r = 1, μ_k = 0.1. (d) μ_r = 1, μ_k = 0.3. (e) μ_r = 1, μ_k = 0.5. (f) μ_r = 2, μ_k = 0.5.

Download Full Size | PDF

5. Transport behaviors of electromagnetic waves in MO waveguide

Once the propagation constants and undetermined coefficients are settled, the propagation fields and transport behaviors of MO modes can also be obtained. Figures 8(a)-8(f) show the EM components distribution of MO mode m₁ under parameters f = 4.8 GHz, d = 2 cm, μ_r = 1, and μ_k = 0.5, the core layer ranges from 0 to 2 cm in x-axis. The calculations are based on our analytical solutions. The H_y and E_y components can be obtained through Eq. (12) and Eq. (16), and the rest components E_x, H_x, E_z and H_z can also be obtained as

 

Fig. 8. The field distributions of EM components obtained by analytical theory and finite element simulation at f = 4.8 GHz. (a)-(f) The theoretically calculated fields of E_x, E_y, E_z, H_x, H_y, and H_z, respectively (see Visualization 1). (g)-(l) The simulated fields of E_x, E_y, E_z, H_x, H_y, and H_z, respectively.

Download Full Size | PDF

As seen from Figs. 8(a)–8(f), the E_y, H_x, and H_z components are π/2 phase ahead or delay of the H_y, E_x, and E_z components, because these components have additional imaginary unit i in the expression. The π/2 phase difference between these EM components leads to Faraday rotation in the waveguide. In our configuration with external magnetic field pointed in the z-axis, the phenomenon of Faraday rotation is predictable, and in the following analysis, it can be seen the analytical solutions could give comprehensive detail of transport behavior of MO modes.

In order to examine the validity of our analytical solutions, we also simulate the MO modes through finite element software. Figures 8(g)–8(l) show the simulation results of the same m₁ MO mode at 4.8 GHz by COMSOL Multiphysics. The MO mode is excited by applying the Port boundary condition to the planar waveguide in COMSOL, with a total exiting power of 1 W. It can be seen that the simulation results are almost completely consistent with the analytical results in Figs. 8(a)–8(f), which proves the correctness of the analytical theory. Note that, the slight difference in the field patterns between the analytical solution and the simulation results is only due to the simulation accuracy.

Figure 9 shows the changes of E and H vectors through propagation along z-axis in the center of core layer, where x = 1 cm. These vectors will rotate during the propagation in core layer, which confirms the Faraday rotation in the waveguide. When the external magnetic field is imposed in the reversed direction, the rotation will also undergo a switch between clockwise and counterclockwise states, which exhibits obvious nonreciprocal behavior. As displayed in Figs. 9(a) and 9(b), the reversion of magnetic field can be realized by taking the negative of μ_k in our calculation. It should be mentioned that the rotation of EM vectors will happen in every position in the core layer, and the EM waves will rotate around the center.

 

Fig. 9. The E and H vectors along z-axis in the center of the core layer. The red vectors represent E and blue vectors represent H (see Visualization 2). (a) μ_k = 0.5. (b) μ_k = -0.5.

Download Full Size | PDF

The energy flow in the waveguide will give more information about the transport behaviors in the core layer. The Poynting vector S, based on the definition of S = E × H, can be calculated in different positions of waveguide. In isotropic case, the Poynting vectors of TE and TM modes are all pointing towards the propagation direction and are always lying in the planar waveguide plane. Figure 10(a) shows the S vectors of TE₀ mode in the core layer when f = 4.8 GHz and d = 2 cm. Here, the S vectors we calculate are time-averaged Poynting vectors, which are defined as ${S_{av}} = \frac{1}{2}\textrm{Re} [\mathbf{E} \times {\mathbf{H}^\ast }]$. The cyan plane is the core layer of planar waveguide and the S vectors are all towards z-axis in the xoz plane. And the S vectors of the m₁ and m₄ MO modes when μ_r = 1, and μ_k = 0.5 are displayed in Figs. 10(b) and 10(c), respectively. Different from TE and TM modes, the Poynting vectors of MO modes have out-of-plane components in the y axis. As seen in Fig. 10(b), the overall Poynting vectors of m₁ mode have the shape of a helical blade, where the S vectors in the center of core layer keep forward in the z-axis and the out-of-plane components on both sides of the center point towards opposite direction. In other words, the EM vectors not only spin in every position inside the MO waveguide, but also rotate around the center. Since that the S vector in the center of core layer located at x = 1 cm has only z components during the propagation, we can treat it as a node of the energy flow. And the EM waves rotate in the core layer with precession along the energy nodes in the waveguide. By extension, it is reasonable to assume that for a rectangular or circular waveguide, which is the extension of planar waveguide model, the EM waves will also rotate in the waveguide with precession along the center of waveguide at the same time. In Fig. 10(c), there are more energy flow nodes of S vectors that only have z component, and on each side of these nodes, the S vectors have opposite out-of-plane components and rotate around the nodes, showing that the precession of EM waves in high-order modes has a more complex form. If the parameter μ_k is set to its negative value, only the out-of-plane components of the S vectors will be flipped while the in-plane components will remain unchanged, which also manifest the nonreciprocal property. For better visual effect of the EM transport behaviors in Figs. 8, 9, 10(b), and 10(c), the animations of propagation fields, E and H vectors, and S vectors evolving over time can be seen in Visualization 1, Visualization 2, and Visualization 3.

 

Fig. 10. The distributions of Poynting vectors S of (a) the TE₀ mode, (b) the m₁ MO mode and (c) the m₄ MO mode in the core layer (cyan plane) of the planar waveguide when f = 4.8 GHz (see Visualization 3).

Download Full Size | PDF

In the end, we investigate the energy flow change in the avoided crossing region. Figure 11 shows the S vectors change of m₁ and m₄ MO modes, which correspond to the analysis of Fig. 6. When the frequency f is 5 GHz, the S vector distribution of the two modes is obviously different, in which m₄ mode has more nodes than m₁ modes, as seen in Figs. 11(a) and 11(d), and now these two modes start to undergo the avoided crossing. As the frequency increases to f = 5.3 GHz, the energy flows of these two MO modes tend to be similar due to the strong coupling between modes. Finally, when transcending the avoided crossing region and f = 5.6 GHz, the S vector distribution of m₁ mode in Fig. 11(c) becomes similar to that of m₄ mode in Fig. 11(d), while the S vector distribution of m₄ mode in Fig. 11(f) becomes similar to that of m₁ mode in Fig. 11(a). That means, when passing through the avoided crossing region, the two MO modes can be converted into each other through the exchange of energy, which also confirms the strong coupling between modes at the avoided crossing region.

 

Fig. 11. The distributions of Poynting vectors S of m₁ MO mode (red arrows) and m₄ MO mode (blue arrows) in the core layer (cyan plane) of the planar waveguide, with frequencies in the avoided crossing region. (a) and (d) f = 5 GHz. (b) and (e) f = 5.3 GHz. (c) and (f) f = 5.6 GHz.

Download Full Size | PDF

6. Conclusion

In conclusion, we have demonstrated the theoretical methods that allows one to have analytical solutions of modes in planar waveguide. When the external magnetic field is set along the propagation direction, the TE and TM modes will be strongly coupled and form new hybrid MO modes. Based on the analytical theory developed in this work, we find that crossings and avoided crossings phenomena can appear in the dispersion curves of MO planar waveguides, which have not been discovered before. Moreover, we deeply study the complex field distributions, mode coupling and transport behaviors of the MO modes. The Faraday rotation and the unique energy flow precession can be observed during the MO mode transmission in the waveguide. By flipping the magnetic field, the nonreciprocal characteristics of the energy flow of MO modes are also clarified. In addition, the perturbation method can be effectively applied to scenarios with weak MO effect, such as the design of magneto-optical devices in terahertz and infrared bands, which are of great interest in recent years. The analytical solutions can give comprehensive information of MO modes, dispersion, and propagation, which have advantage in MO waveguide design. Compared with finite element simulation, the analytical solution can be time-saving in setting geometric and material parameters of MO waveguide devices. As a basic example, the planar MO model can be extended to ridge or circular MO waveguide, which has potential in designing MO nonreciprocal devices. Moreover, the methods developed in this paper will be of great significance for further understanding of the detailed energy flow processes and light-matter interaction in various nonreciprocal photonic devices such as topological waveguides based on MO materials, and for the future development of subwavelength MO integrated photonic devices.