Layered optomagnonic structures: Time Floquet scattering-matrix approach

TL;DR

This paper develops a dynamic theoretical framework using a time Floquet scattering-matrix method to analyze layered optomagnonic structures, revealing nonlinear interactions and enhanced scattering effects in a dual photonic-magnonic cavity.

Contribution

It introduces a comprehensive time-dependent scattering-matrix approach for layered optomagnonic systems, enabling analysis of nonlinear and strong-coupling phenomena.

Findings

Identification of nonlinear photon-magnon interaction effects.

Observation of enhanced inelastic light scattering in strong-coupling regime.

Limitations of the quasistatic adiabatic approximation are established.

Abstract

A fully dynamic theoretical approach to layered optomagnonic structures, based on a time Floquet scattering-matrix method, is developed. Its applicability is demonstrated on a simple design of a dual photonic-magnonic cavity, formed by sandwiching a magnetic garnet thin film between two dielectric Bragg mirrors, subject to continuous excitation of a perpendicular standing spin wave. Some remarkable phenomena, including nonlinear photon-magnon interaction effects and enhanced inelastic light scattering in the strong-coupling regime, fulfilling a triple-resonance condition, are analyzed and the limitations of the quasistatic adiabatic approximation are established.

Tables

Table 1. Table 1: Intensities of the elastic and the first few inelastic outgoing beams upon excitation of the higher optical resonance at q ∥ ​ a = 1.2 subscript 𝑞 parallel-to 𝑎 1.2 q_{\parallel}a=1.2 in Fig.  1 (c), subject to the action of a second-order perpendicular standing spin wave with a relative amplitude η = 0.1 𝜂 0.1 \eta=0.1 .

$n$	$-4$	$-3$	$-2$	$-1$	$0$	$1$	$2$	$3$	$4$
${ℐ}_n$	$0.015$	$0.025$	$0.042$	$0.070$	$0.648$	$0.071$	$0.043$	$0.025$	$0.015$

Equations

$$\boldsymbol{\epsilon}=\left(\begin{array}[]{ccc}\epsilon&if&0\\ -if&\epsilon&0\\ 0&0&\epsilon\\ \end{array}\right)\;$$

$${\Omega }_{\kappa }=\gamma {\mu }_0{M}_{\mathrm{s}}\left[\frac{H}{M_{\mathrm{s}}}-1+\beta +\alpha {\left(\frac{\kappa \pi }{d}\right)}^2\right],\kappa =1,2,...,$$

$$\mathbf{M}(\mathbf{r},t)/M_{\mathrm{s}}$$

$$+\eta \sin (\kappa \pi z/d)\sin ({\Omega }_{\kappa }t)\hat{\mathbf{y}}+\hat{\mathbf{z}},$$

$$\delta{\boldsymbol{\epsilon}}(z,t)=\frac{1}{2}\big{[}\delta{\boldsymbol{\epsilon}}(z)\exp{(-i\Omega_{\kappa}t)}+\delta{\boldsymbol{\epsilon}}^{\dagger}(z)\exp{(i\Omega_{\kappa}t)}\big{]}\;$$

$$\delta{\boldsymbol{\epsilon}}(z)=f\eta\sin(\kappa\pi z/d)\left(\begin{array}[]{ccc}0&0&1\\ 0&0&i\\ -1&-i&0\end{array}\right)\;$$

$$D_{\alpha }(\mathbf{r},t)\cong {ϵ}_0\sum _{\beta }\int d{t}^{\prime }{\overline{ϵ}}_{\alpha \beta }(t-t^{\prime })m_{\alpha \beta }(\mathbf{r},t)E_{\beta }(\mathbf{r},t^{\prime })$$

$$\cong {ϵ}_0\mathrm{Re}\{\sum _{\beta }{ϵ}_{\alpha \beta }(\omega )m_{\alpha \beta }(\mathbf{r},t){\mathcal{E}}_{\beta }(\mathbf{r},t)\exp (-i\omega t)\},$$

$${\mathbf{E}}(\mathbf{r},t)=\mathrm{Re}\Big{\{}\sum_{n=0,\pm 1,\ldots}\mathbf{E}_{n}(\mathbf{r})\exp[-i(\omega-n\Omega_{\kappa})t]\Big{\}}\;.$$

$$G=if\eta\delta(\mathbf{q}_{\parallel}-\mathbf{q}_{\parallel}^{\prime})\big{[}\delta(\omega-\omega^{\prime}-\Omega_{\kappa})g_{-}+\delta(\omega-\omega^{\prime}+\Omega_{\kappa})g_{+}\big{]}$$

$$g_{\pm }=\int$$

$$\displaystyle\pm i\big{[}E_{x}^{\mathrm{out\star}}(z)E_{z}^{\mathrm{in}}(z)-E_{z}^{\mathrm{out\star}}(z)E_{x}^{\mathrm{in}}(z)\big{]}\Big{\}}~{}.$$

$$\nabla \times \mathcal{E}(\mathbf{r},t)\exp (-i\omega t)$$

$$\nabla \times \mathcal{H}(\mathbf{r},t)\exp (-i\omega t)$$

$$\mathbf{ϵ}(t)$$

$$\mathcal{E}(\mathbf{r},t)$$

$$\mathcal{H}(\mathbf{r},t)$$

$$\mathbf{q}\times \mathbf{h}(n)+\frac{(\omega -n\Omega )}{c}\sum _{n^{\prime }=-N}^N\mathbf{ϵ}(n-n^{\prime })\mathbf{e}(n^{\prime })=0$$

$$\mathbf{q}\times \mathbf{e}(n)-\frac{(\omega -n\Omega )}{c}\mu \mathbf{h}(n)=0,$$

$$\displaystyle\left(\begin{array}[]{cc}\underline{\underline{\mathbf{k}}}\,\underline{\underline{\boldsymbol{\epsilon}}}&\underline{\underline{\mathbf{C}}}_{1}\\ -\underline{\underline{\mathbf{C}}}_{1}&\underline{\underline{\mathbf{k}}}\mu\end{array}\right)^{-1}\left(\begin{array}[]{cc}\underline{\underline{\mathbf{0}}}&\underline{\underline{\mathbf{C}}}_{2}\\ -\underline{\underline{\mathbf{C}}}_{2}&\underline{\underline{\mathbf{0}}}\end{array}\right)\left(\begin{array}[]{c}\underline{\mathbf{e}}\\ \underline{\mathbf{h}}\\ \end{array}\!\right)=\dfrac{1}{q_{z}}\left(\begin{array}[]{c}\underline{\mathbf{e}}\\ \underline{\mathbf{h}}\\ \end{array}\!\right)\;,$$

$$\mathbf{C}_{1}=\left(\begin{array}[]{ccc}0&0&q_{y}\\ 0&0&-q_{x}\\ -q_{y}&q_{x}&0\end{array}\right),\mathbf{C}_{2}=\left(\begin{array}[]{ccc}0&1&0\\ -1&0&0\\ 0&0&0\end{array}\right)\;,$$

$$\underline{\underline{\boldsymbol{\epsilon}}}=\left(\begin{array}[]{cccc}\boldsymbol{\epsilon}(0)&\boldsymbol{\epsilon}(-1)&\cdots&\boldsymbol{\epsilon}(-2N)\\ \boldsymbol{\epsilon}(1)&\boldsymbol{\epsilon}(0)&\cdots&\boldsymbol{\epsilon}(-2N+1)\\ \vdots&\vdots&&\vdots\\ \boldsymbol{\epsilon}(2N)&\boldsymbol{\epsilon}(2N-1)&\cdots&\boldsymbol{\epsilon}(0)\end{array}\right)\;.$$

$$\underline{\mathbf{e}}=\left(\begin{array}{c}\mathbf{e}(-N)\\ \mathbf{e}(-N+1)\\ ⋮\\ \mathbf{e}(N)\end{array}\right),\underline{\mathbf{h}}=\left(\begin{array}{c}\mathbf{h}(-N)\\ \mathbf{h}(-N+1)\\ ⋮\\ \mathbf{h}(N)\end{array}\right).$$

$${\mathcal{E}}_{p\nu }^s(\mathbf{r},t)$$

$${\mathcal{H}}_{p\nu }^s(\mathbf{r},t)$$

$$\displaystyle\left(\begin{array}[]{cc}{k}_{n}\boldsymbol{\epsilon}&\mathbf{C}_{1}\\ -\mathbf{C}_{1}&{k}_{n}\mu\mathbf{I}\end{array}\right)^{-1}\left(\begin{array}[]{cc}\mathbf{0}&\mathbf{C}_{2}\\ -\mathbf{C}_{2}&\mathbf{0}\end{array}\right)\left(\begin{array}[]{c}\mathbf{e}\\ \mathbf{h}\\ \end{array}\right)=\dfrac{1}{q_{z}}\left(\begin{array}[]{c}\mathbf{{e}}\\ \mathbf{{h}}\\ \end{array}\right),$$

$$\displaystyle\left(\begin{array}[]{cccccccc}-{\underline{{e}}}_{1-Nx}^{-(1)}&-{\underline{{e}}}_{2-Nx}^{-(1)}&\cdots&-{\underline{{e}}}_{2Nx}^{-(1)}&{\underline{{e}}}_{1-Nx}^{+(2)}&{\underline{{e}}}_{2-Nx}^{+(2)}&\cdots&{\underline{{e}}}_{2Nx}^{+(2)}\\ -{\underline{{e}}}_{1-Ny}^{-(1)}&-{\underline{{e}}}_{2-Ny}^{-(1)}&\cdots&-{\underline{{e}}}_{2Ny}^{-(1)}&{\underline{{e}}}_{1-Ny}^{+(2)}&{\underline{{e}}}_{2-Ny}^{+(2)}&\cdots&{\underline{{e}}}_{2Ny}^{+(2)}\\ -{{\underline{h}}}_{1-Nx}^{-(1)}&-{{\underline{h}}}_{2-Nx}^{-(1)}&\cdots&-{{\underline{h}}}_{2Nx}^{-(1)}&{{\underline{h}}}_{1-Nx}^{+(2)}&{{\underline{h}}}_{2-Nx}^{+(2)}&\cdots&{{\underline{h}}}_{2Nx}^{+(2)}\\ -{{\underline{h}}}_{1-Ny}^{-(1)}&-{{\underline{h}}}_{2-Ny}^{-(1)}&\cdots&-{{\underline{h}}}_{2Ny}^{-(1)}&{\underline{{h}}}_{1-Ny}^{+(2)}&{\underline{{h}}}_{2-Ny}^{+(2)}&\cdots&{{\underline{h}}}_{2Ny}^{+(2)}\\ \end{array}\right)\begin{pmatrix}S_{1-N;p^{\prime}\nu^{\prime}}^{-s^{\prime}}\\ S_{2-N;p^{\prime}\nu^{\prime}}^{-s^{\prime}}\\ \vdots\\ S_{2N;p^{\prime}\nu^{\prime}}^{-s^{\prime}}\\ S_{1-N;p^{\prime}\nu^{\prime}}^{+s^{\prime}}\\ S_{2-N;p^{\prime}\nu^{\prime}}^{+s^{\prime}}\\ \vdots\\ S_{2N;p^{\prime}\nu^{\prime}}^{+s^{\prime}}\end{pmatrix}=s^{\prime}\begin{pmatrix}{\underline{{e}}}_{p^{\prime}\nu^{\prime}x}^{s^{\prime}(\mathrm{m})}\\ {\underline{{e}}}_{p^{\prime}\nu^{\prime}y}^{s^{\prime}(\mathrm{m})}\\ {\underline{{h}}}_{p^{\prime}\nu^{\prime}x}^{s^{\prime}(\mathrm{m})}\\ {\underline{{h}}}_{p^{\prime}\nu^{\prime}y}^{s^{\prime}(\mathrm{m})}\\ \end{pmatrix}\;,$$

$$\begin{array}{rl}& Q_{p\nu ;p^{\prime }{\nu }^{\prime }}^{\mathrm{I}}=S_{p\nu ;p^{\prime }{\nu }^{\prime }}^{++}\exp \{i[{\mathbf{q}}_{p\nu }^{+(2)}\cdot {\mathbf{d}}_2+{\mathbf{q}}_{p^{\prime }{\nu }^{\prime }}^{+(1)}\cdot {\mathbf{d}}_1]\}\\ & Q_{p\nu ;p^{\prime }{\nu }^{\prime }}^{\mathrm{II}}=S_{p\nu ;p^{\prime }{\nu }^{\prime }}^{+-}\exp \{i[{\mathbf{q}}_{p\nu }^{+(2)}\cdot {\mathbf{d}}_2-{\mathbf{q}}_{p^{\prime }{\nu }^{\prime }}^{-(2)}\cdot {\mathbf{d}}_2]\}\\ & Q_{p\nu ;p^{\prime }{\nu }^{\prime }}^{\mathrm{III}}=S_{p\nu ;p^{\prime }{\nu }^{\prime }}^{-+}\exp \{-i[{\mathbf{q}}_{p\nu }^{-(1)}\cdot {\mathbf{d}}_1-{\mathbf{q}}_{p^{\prime }{\nu }^{\prime }}^{+(1)}\cdot {\mathbf{d}}_1]\}\\ & Q_{p\nu ;p^{\prime }{\nu }^{\prime }}^{\mathrm{IV}}=S_{p\nu ;p^{\prime }{\nu }^{\prime }}^{--}\exp \{-i[{\mathbf{q}}_{p\nu }^{-(1)}\cdot {\mathbf{d}}_1+{\mathbf{q}}_{p^{\prime }{\nu }^{\prime }}^{-(2)}\cdot {\mathbf{d}}_2]\}.\end{array}$$

$$\displaystyle\left(\!\!\begin{array}[]{ccccc}-{{{e}}}_{1nx}^{-(1)}&-{{{e}}}_{2nx}^{-(1)}&{{{e}}}_{1nx}^{+(2)}&{{{e}}}_{2nx}^{+(2)}\\ -{{{e}}}_{1ny}^{-(1)}&-{{{e}}}_{2ny}^{-(1)}&{{{e}}}_{1ny}^{+(2)}&{{{e}}}_{2ny}^{+(2)}\\ -{{{h}}}_{1nx}^{-(1)}&-{{{h}}}_{2nx}^{-(1)}&{{{h}}}_{1nx}^{+(2)}&{{{h}}}_{2nx}^{+(2)}\\ -{{{h}}}_{1ny}^{-(1)}&-{{{h}}}_{2ny}^{-(1)}&{{{h}}}_{1ny}^{+(2)}&{{{h}}}_{2ny}^{+(2)}\\ \end{array}\right)\!\!\!\begin{pmatrix}S_{1n;p^{\prime}n}^{-s^{\prime}}\\ S_{2n;p^{\prime}n}^{-s^{\prime}}\\ S_{1n;p^{\prime}n}^{+s^{\prime}}\\ S_{2n;p^{\prime}n}^{+s^{\prime}}\end{pmatrix}\!\!=\!\!s^{\prime}\!\!\begin{pmatrix}{{{e}}}_{p^{\prime}nx}^{s^{\prime}(\mathrm{m})}\\ {{{e}}}_{p^{\prime}ny}^{s^{\prime}(\mathrm{m})}\\ {{{h}}}_{p^{\prime}nx}^{s^{\prime}(\mathrm{m})}\\ {{{h}}}_{p^{\prime}ny}^{s^{\prime}(\mathrm{m})}\\ \end{pmatrix}\;,$$

$${\mathbf{Q}}^{\mathrm{I}}$$

$${\mathbf{Q}}^{\mathrm{II}}$$

$${\mathbf{Q}}^{\mathrm{III}}$$

$${\mathbf{Q}}^{\mathrm{IV}}$$

$${\mathcal{E}}_{\mathrm{tr}}(\mathbf{r},t)$$

$${\mathcal{E}}_{\mathrm{rf}}(\mathbf{r},t)$$

$$\mathcal{T}=\sum _{p,n}{\mathcal{T}}_{pn}=\sum _{p,n}|Q_{pn;p^{\prime }{n}^{\prime }}^{\mathrm{I}}{|}^2\frac{\mathrm{Re}[q_{nz}^{+(\mathrm{B})}]{\mu }_{\mathrm{A}}}{\mathrm{Re}[q_{n^{\prime }z}^{+(\mathrm{A})}]{\mu }_{\mathrm{B}}}$$

Full text

††preprint: Submitted to PRB

I Introduction

Propagation of electromagnetic (EM) waves in dynamic media with a periodic (spatio)temporal variation of their refractive index, realized, for example, by the action of a progressive sinusoidal disturbance, have long been investigated and a number of intriguing effects, such as formation of frequency and wave vector band gaps, wavelength conversion, pulse shaping, and parametric amplification, were properly analyzed. morgenthaler ; simon ; hessel ; holberg ; harfoush ; notomi ; biancalana ; zurita ; koutserimpas The lack of time invariance of the linear medium can also lead to nonreciprocal optical response, thus offering a promising alternative in the design of compact nonreciprocal photonic devices (see, e.g., Ref. sounas, and references therein).

Spatiotemporal modulation of the electric permittivity tensor of magnetic materials is induced, for instance, by the various types of spin waves that can be excited in different confined geometries, stancil_book due to the periodic variation of the magneto-optic coupling coefficients, hu causing diverse interesting and useful effects. For example, diffraction and mode conversion of guided optical waves, induced by magnetostatic waves in magnetic garnet thin films, offer the potential of large time-bandwidth optical signal processing at microwave frequencies of 1 to 20 GHz and higher. fisher ; stancil ; prabhakar ; fetisov On the other hand, inelastic light scattering by thermally excited spin waves is widely used for probing magnetic properties of bulk and microstructured materials through Raman or Brillouin spectroscopy. fleury ; sandercock ; grunberg ; grimsditch ; hillebrands ; cochran ; jorzick ; demokritov ; meloche ; sebastian It should be noted that the magneto-optic interaction is inherently weak and, therefore, the above effects can be in most cases described by first-order perturbation theory.

Greatly increased light–spin-wave interaction phenomena may arise from the simultaneous localization of both fields in the same ultra-small region of space for a long time period. This can be achieved, for example, in(sub)millimeter magnetic garnet spheres, which support optical whispering gallery modes osada ; zhang ; haigh1 ; hisatomi ; sharma ; haigh2 or Mie resonances, almpanis1 as well as in appropriately designed, so-called photomagnonic (or optomagnonic), stratified structures that comprise magnetic dielectric layers and exhibit dual, photonic and magnonic, functionalities klos ; dadoenkova ; pantazopoulos ; pantazopoulos2 combining the properties of photonic joannopoulos and magnonic kruglyak2010 ; zivieri2012 crystals. Strong parametric photon-magnon coupling effects, allowing for coherent manipulation of elementary magnetic excitations in solids by optical means and vice versa, have been anticipated viola ; liu while large dynamic frequency shift and enhanced modulation of the optical field through multi-magnon absorption and emission processes by a photon have also been reported almpanis1 ; pantazopoulos in both of the above (spherical and planar) optomagnonic cavity designs.

Since the frequency of a pump spin wave is typically a few orders of magnitude smaller than that of a probe light beam in the infrared and visible range of the EM spectrum, interaction between the two fields can be in principle described by a quasi-static adiabatic approximation. In practice, this translates in calculating the relevant optical scattering amplitudes at a sequence of frozen snapshots of the spin wave, during a period, and, at the end of the calculation, obtain the frequency-domain response by Fourier transform. almpanis1 ; pantazopoulos A similar approach was also successfully applied to analyze acousto-optic interaction effects in corresponding phoxonic cavities, gantzounis ; psarobas ; almpanis2 ; almpanis3 encompassing both weak- and strong-coupling regimes and recovering the results of the widely employed linear photoelastic model in the weak-coupling limit. almpanis3 However, the quasi-static adiabatic approximation has a number of drawbacks, which may be more likely manifested in the case of a spin instead of an acoustic pump wave because the former can attain higher frequencies without being strongly attenuated. The adiabatic approximation precludes energy transfer from one wave to the other while, even more importantly, the frequency of the pump wave is decoupled from that of the light beam and thus its actual value becomes immaterial. Consequently, this approximation cannot properly describe, e.g., inelastic light scattering processes in a triple-resonance condition, which occur in optomagnonic cavities, osada ; zhang ; haigh1 where the frequency of the magnon matches a photon transition between two resonant modes.

It is therefore tempting to formulate a rigorous, fully dynamic Floquet scattering-matrix method for (spatio)temporal periodic media, adapted to the case of stratified optomagnonic architectures, following, e.g., an approach analogous to that developed in relation to quantum transport in periodically driven systems. wagner ; li ; moskalets ; agarwal ; bilitewski The purpose of the present paper is to develop such a method and apply it to a specific example of layered optomagnonic structure, establishing the limitations of the quasi-static adiabatic approach. The remaining of the paper is structured as follows. In Sec. II we introduce our optomagnonic model structure, we present some approximate methods for its theoretical description, namely the first-order Born approximation and the adiabatic approach , and briefly recall the results of the latter in the case under consideration. In Sec. III we develop a fully dynamic time-Floquet scattering-matrix method and report a detailed study of the given optomagnonic structure, encompassing the quasi-static limit and providing a consistent interpretation of some remarkable effects which cannot be accounted for by the adiabatic approximation. The last section concludes the paper.

II Structure Description and Approximate Treatments

In magnetic materials, the interaction of visible and infrared light with spin waves is inherently weak. Optomagnonic cavities aim at enhancing this interaction by simultaneously confining both light and spin waves in the same ultra-small region of space for a long time period. In the present paper, following our previous work, pantazopoulos we investigate further a planar optomagnonic cavity formed by a versatile trilayer configuration that comprises a middle magnetic dielectric film made of bismuth-substituted yttrium iron garnet (Bi:YIG) sandwiched between two dielectric Bragg mirrors, as schematically depicted in Fig. 1(a). More specifically, we consider a $\mathrm{TiO_{2}}$/Bi:YIG/$\mathrm{TiO_{2}}$ trilayer and two symmetric Bragg mirrors, each consisting of 14 periods, with lattice constant $a$, of alternating $\mathrm{SiO}_{2}$ and $\mathrm{TiO}_{2}$ layers, of thickness $0.6a$ and $0.4a$, respectively, stacked along the $z$ direction. The whole structure is embedded in air. The trilayer consists of a Bi:YIG film, $0.7a$ thick, magnetically saturated perpendicular to the interfaces, placed between two identical $\mathrm{TiO_{2}}$ layers, each of thickness $0.4a$. Assuming $a$ to be of the order of a few hundred nanometers, the operation wavelengths are in the near-infrared region, where the material relative magnetic permeability is equal to unity and Bi:YIG is characterized by a relative electric permittivity tensor of the form

[TABLE]

with $\epsilon=5.5$ and $f=-0.01$.zvezdin_book $\mathrm{TiO}_{2}$ and $\mathrm{SiO}_{2}$ are optically isotropic materials described by a relative permittivity constant equal to 5.35 and 2.13, respectively. Schiefke ; gao

The periodic structuring of the Bragg mirrors gives rise to frequency gaps, where light propagation is forbidden. joannopoulos Figure 1(b) displays the transmission spectrum of the structure described here above, under illumination by p-polarized light, i.e. light linearly polarized in the plane of incidence, with a wave-vector component parallel to the interfaces $q_{\parallel}=({q_{x}^{2}+q_{y}^{2}})^{1/2}$, in a frequency window about the lowest Bragg gap. We note that, due to translation invariance in the $x-y$ plane, $\mathbf{q}_{\parallel}$ remains constant. The presence of the magnetic film affects the optical response of the structure in two ways. Firstly, it breaks periodicity leading to the appearance of two resonant modes in the gap, spatially localized in the defect region. joannopoulos Secondly, it removes reflection symmetry with respect to the plane of incidence, so that the p and s polarization modes, i.e., modes linearly polarized in and normal to the plane of incidence, respectively, are no longer true eigenmodes of the system. The defect modes do not have a well-defined linear polarization character and thus both couple, of course to a different degree, to any linearly polarized incident light beam. Here, p-polarized light is used to excite these modes, as shown in Fig. 1(c). We note in passing that two distinct modes subsist also in the unmagnetized slab. They originate from the splitting of a doubly degenerate mode at normal incidence, where the two polarization degrees of freedom are, obviously, equivalent. The difference is that, in the case of the unmagnetized slab, reflection symmetry with respect to the plane of incidence implies a specific linear polarization character, p or s, to each of these two modes at off-normal incidence.

Apart from the localization of the EM field in the magnetic film, the considered optomagnonic structure serves, also, as a cavity for spin waves. A general theory of dipole-exchange spin waves in a homogeneous magnetic film of thickness $d$, with out-of-plane arbitrary orientation of the magnetization, has been reported by Kalinikos and Slavin. kalinikos1986 Here, we shall be concerned with perpendicular standing spin waves, assuming the easy axis of magnetization perpendicular to the film, along the $z$ direction. If the film is saturated to $M_{\mathrm{s}}$ by an externally applied static magnetic field $H\widehat{\mathbf{z}}$, under the boundary condition of perfect pinning of the spin waves at the film boundaries we obtain a series of modes at frequencies

[TABLE]

where $\gamma$ is the gyromagnetic ratio, $\mu_{0}$ the magnetic permeability of vacuum, $\alpha$ the exchange constant, and $\beta$ the dimensionless anisotropy coefficient. The corresponding magnetization field profiles are given by

[TABLE]

where $\eta$ is the relative spin wave amplitude. pantazopoulos Assuming $\eta=0.1$, the magnetization precession angle is about $5^{\mathrm{o}}$, which is a tolerable value for linear spin waves.

The interaction between light and spin waves enters through the electric permittivity tensor. Specifically, the magnetization field given by Eq. (3) induces a spatiotemporal perturbation pantazopoulos

[TABLE]

in the permittivity tensor of the statically magnetized material, where

[TABLE]

and the dagger denotes hermitian conjugation.

In a simple anisotropic material, the displacement field is related to the electric field through the electric permittivity tensor. The underlying constitutive relation assuming local response, i.e., that at each point of space the displacement field depends only on the value of the electric field at the same point, reads $\mathbf{D}(\mathbf{r},t)=\epsilon_{0}\int dt^{\prime}\overline{\boldsymbol{\epsilon}}(t-t^{\prime})\mathbf{E}(\mathbf{r},t^{\prime})$, where $\epsilon_{0}$ is the electric permittivity of vacuum. The relative electric permittivity tensor in the time domain, $\overline{\boldsymbol{\epsilon}}$, depends on the difference $t-t^{\prime}$, as implied by time homogeneity, and is related to the usual electric permittivity tensor $\boldsymbol{\epsilon}$ (in the frequency domain) though a Fourier transform $\boldsymbol{\epsilon}(\omega)=\int dt\overline{\boldsymbol{\epsilon}}(t)\exp(i\omega t)$. For monochromatic time-harmonic waves of angular frequency $\omega$, of the form $\mathbf{E}(\mathbf{r},t)=\mathrm{Re}[\mathbf{E}(\mathbf{r})\exp(-i\omega t)]$, we obtain $\mathbf{D}(\mathbf{r},t)=\epsilon_{0}\mathrm{Re}[\boldsymbol{\epsilon}(\omega)\mathbf{E}(\mathbf{r})\exp(-i\omega t)]$. In the presence of a relatively slow (with respect to the period of the EM wave) modulation, $m_{\alpha\beta}(\mathbf{r},t)$, of the elements of the electric permittivity tensor, such as that induced by a spin wave, one can write $\mathbf{E}(\mathbf{r},t)=\mathrm{Re}[\boldsymbol{\mathcal{E}}(\mathbf{r},t)\exp(-i\omega t)]$, where $\boldsymbol{\mathcal{E}}(\mathbf{r},t)$ is a slowly varying envelope function instead of a time-independent phasor. Therefore, ignoring the variation of $m_{\alpha\beta}$ and $\boldsymbol{\mathcal{E}}$ in the time integral of the constitutive relation we obtain

[TABLE]

i.e., we can write, approximately, the constitutive relation in a mixed time-frequency domain representation in the form of an instantaneous-response, instead of a time-convolution, equation.

Due to the periodic time variation of the electric permittivity tensor of the structure under consideration in the presence of a perpendicular standing spin wave, $\boldsymbol{\mathcal{E}}(\mathbf{r},t)$ can be expanded in a Fourier series as follows: $\!\boldsymbol{\mathcal{E}}(\mathbf{r},t)=\sum_{n=0,\pm 1,\ldots}\!\!\!\!\!\mathbf{E}_{n}(\mathbf{r})\exp(in\Omega_{\kappa}t)\;,$ which yields

[TABLE]

This equation implies that, for an incident wave of angular frequency $\omega$, the total outgoing (transmitted plus reflected) field consists, in general, of an infinite number of monochromatic beams with angular frequencies $\omega,\omega\pm\Omega_{\kappa},\omega\pm 2\Omega_{\kappa},\ldots$, which are produced by elastic and inelastic photon scattering that involves absorption and/or emission of zero, one, two, etc. magnons. We note that magnon absorption (emission) processes stem from the first (second) term of Eq. (4).

In the weak-coupling regime, one can restrict to the first-order Born approximation. In this approximation, the coupling strength associated to the photon-magnon scattering is proportional to the overlap integral $G=\bra{\mathrm{out}}\delta{\boldsymbol{\epsilon}}\ket{\mathrm{in}}$, where $\braket{\alpha\mathbf{r}t}{{\mathrm{in}}}=E_{\alpha}^{\mathrm{in}}(z)\exp[i(\mathbf{q}_{\parallel}\cdot\mathbf{r}-\omega t)]$ and $\braket{{\mathrm{out}}}{\alpha^{\prime}\mathbf{r}^{\prime}t^{\prime}}=E_{\alpha^{\prime}}^{\mathrm{out}\star}(z)\exp[-i(\mathbf{q}_{\parallel}^{\prime}\cdot\mathbf{r}-\omega^{\prime}t)]$ denote appropriate incoming and outgoing monochromatic time-harmonic waves in the static magnetic layered structure. Using Eq. (4) we obtain

[TABLE]

where

[TABLE]

The delta functions in Eq. (8) express conservation of in-plane momentum and energy in inelastic light scattering processes that involve emission and absorption of one magnon by a photon, as expected in the linear regime. In our case here, a simple selection rule can be deduced by considering the symmetry of the EM field and the spin wave. In a structure such as that described in Fig. 1(a), which has the mirror symmetry with respect to the $x-y$ plane, $\widehat{\sigma}_{z}$, the electic field modes, $\mathbf{E}(z)$, are either even or odd functions of $z$. pantazopoulos On the other hand, since $\widehat{\sigma}_{z}\mathbf{E}(z)=E_{x}(-z)\widehat{\mathbf{x}}+E_{y}(-z)\widehat{\mathbf{y}}-E_{z}(-z)\widehat{\mathbf{z}}$, either $E_{x},E_{y}$ are even function of $z$ and $E_{z}$ is an odd function of $z$ or $E_{x},E_{y}$ are odd function of $z$ and $E_{z}$ is an even function of $z$. Therefore, the integral in Eq. (9) vanishes identically, unless the spatial profile of the perpendicular standing spin wave is an odd (even) function of $z$ for transitions between two optical modes of the same (different) symmetry. In the present work, since both optical resonances have the same symmetry, we choose a second-order ($\kappa=2$) perpendicular standing spin wave which is an odd function of $z$.

On the other hand, since the spin-wave frequencies are, typically, a few orders of magnitude smaller than those of visible and infrared light, the slow gradual change of the electric permittivity tensor, given by Eq. (4), can be looked upon as an adiabatic process. Therefore, the optical response of the driven system can be described by first evaluating it with the external parameters held fixed and then, at the end of the calculation, allow them to change, in discrete steps within a time period of the spin wave. This, so-called quasi-static adiabatic, approach, which is not restricted to the weak-coupling regime, has been undertaken in our previous work pantazopoulos and details about its implementation can be found elsewhere. pantazopoulos ; pantazopoulos2 Here, we briefly summarize our main results when the higher optical resonance of Fig. 1(c) for $q_{\parallel}a=1.2$, at $\omega_{0}a/c=1.88375$, is excited by p-polarized light incident with the given $q_{\parallel}$, i.e., at an angle of about $40^{\circ}$, on the structure of Fig. 1(a) in the presence of a second-order ($\kappa=2$) perpendicular standing spin wave with a relative amplitude $\eta=0.1$. The continuous variation of the permittivity tensor of the magnetic film, as implied by Eq. (4), is taken into consideration in the numerical calculation following a discretization approach. Subdividing the Bi:YIG film into 50 homogeneous elementary sublayers, at each of the 60 in total time steps considered, is sufficient to obtain excellent convergence.

When the dynamic magnetization field is switched on, the permittivity tensor and, consequently, the optical transmission spectrum vary periodically in time with the period of the spin wave. Even though the overall spectral features do not change significantly as time evolves, the oscillation of the optical defect mode with an amplitude $\Delta\omega$ that is of the order of the resonance width is remarkable. pantazopoulos The ratio $\Delta\omega/\Gamma$, where $\Gamma$ is the calculated width of the higher resonance in the optical transmission spectrum for $q_{\parallel}a=1.2$ [see Fig. 1(c)] that defines its inverse lifetime, can be adopted as a quantitative measure of the strength of the photon-magnon interaction. Here we obtain $\Delta\omega a/c=0.2\times 10^{-5}$ and $\Gamma a/c=0.1\times 10^{-5}$. This strong modulation of the optical cavity mode by the perpendicular standing spin wave should be ascribed to the simultaneous concentration of both fields for a long time period in the cavity region. In the wave picture, this is manifested as a relatively large-amplitude periodic oscillation of the position of the sharp optical resonance. Correspondingly, in the particle picture, strong inelastic light scattering occurs, with considerable probabilities for absorption and emission of many magnons by a photon. Indeed, the calculated intensities, $\mathcal{I}_{n}$, of the outgoing (transmitted plus reflected) light beams, reported in Table 1, are significant also beyond the relatively large first-order ($n=\pm 1$) components. It is worth-noting that the quasi-static adiabatic approximation precludes energy transfer between the photon and the magnon fields. pantazopoulos2 Therefore, in the absence of dissipative losses, $\sum_{n=-\infty}^{\infty}\mathcal{I}_{n}$ must be equal to unity, because at each frozen snapshot of the spin wave the electric permittivity tensor is Hermitian [see Eqs. (1), (4), and (5)], which is confirmed by our calculations. However, even if there is no energy transfer, coupling between the two fields arises through the time variation of the electric permittivity tensor, induced by the spin wave, which affects the optical response of the structure. The small asymmetry between the Stokes ($n>0$) and anti-Stokes ($n<0$) components in Table 1 is due to the lack of time-reversal symmetry of the magnetic structure. pantazopoulos

Though the quasi-static adiabatic approximation is a useful tool for the study of optomagnonic cavities, it has a number of drawbacks (see Appendix). Namely, as stated above, it precludes energy transfer between the photon and the magnon fields. Even more importantly, the frequency of the pump spin wave is decoupled from that of the light beam and thus its actual value is immaterial. Consequently, this approximation cannot properly describe, e.g., inelastic light scattering processes in a triple-resonance condition, which occur in optomagnonic whispering gallery microresonators, osada ; zhang ; haigh1 when the frequency of the magnon matches a photon transition between two resonant modes.

III Fully Dynamic Description

III.1 Formalism

Let us first consider a homogeneous medium, characterized by a scalar relative magnetic permeability $\mu$ and a relative time-periodic electric permittivity tensor, $\boldsymbol{\epsilon}(t)=\boldsymbol{\epsilon}(t+T)$. If the time variation of the permittivity is very slow compared with the period of an optical wave, assuming slowly varying envelope functions for the electric and magnetic fields, of the form \mathbf{F}(\mathbf{r},t)=\mathrm{Re}\Big{[}\boldsymbol{\mathcal{F}}(\mathbf{r},t)\exp(-i\omega t)\Big{]} where $\mathbf{F}=\mathbf{E},\mathbf{H}$, with the help of Eq. (6) Maxwell equations read

[TABLE]

The time-periodic variation of the permittivity tensor results in time-periodic envelope functions and, for this reason, Maxwell equations can be solved by expanding these quantities in (truncated) Fourier series. Considering solutions in the form of plane waves with wavevector $\mathbf{q}$, we have

[TABLE]

where $\Omega=2\pi/T$, $Z_{0}$ is the impedance of free space, $\mathbf{e}$ and $\mathbf{h}$ are polarization vectors chosen such that $\sum_{n=-N}^{N}\mathbf{e}(n)\cdot\mathbf{e}^{\star}(n)=1$, and $E_{0}$ is the field amplitude. Substituting Eqs. (11) into Eq. (10), we obtain

[TABLE]

for $n=-N,-N+1,\cdots,N$. For given $x$- and $y$-components of the wavevector $\mathbf{q}$, Eqs. (12) can be cast in the form of a $6(2N+1)\times 6(2N+1)$ linear eigenvalue problem

[TABLE]

which can be solved by standard numerical algorithms. numerical_recipes In Eq. (21), by a double underscore we denote $3(2N+1)\times 3(2N+1)$ matrices; $\underline{\underline{\mathbf{0}}}$ is the zero matrix, $\underline{\underline{\mathbf{k}}}=\mathrm{diag}(\mathbf{k}_{-N},\mathbf{k}_{-N+1},\cdots,\mathbf{k}_{N})$ with $\mathbf{k}_{n}=\mathbf{I}(\omega-n\Omega)/c\equiv\mathbf{I}k_{n}$, and $\underline{\underline{\mathbf{C}}}_{\,i}=\mathrm{diag}(\mathbf{C}_{i},\mathbf{C}_{i},\cdots,\mathbf{C}_{i})$, $i=1,2$, are block diagonal matrices with

[TABLE]

while $\underline{\underline{\boldsymbol{\epsilon}}}$ is the Toeplitz matrix of the Fourier coefficients of the electric permittivity tensor

[TABLE]

Correspondingly $\underline{\mathbf{e}}$, $\underline{\mathbf{h}}$ are $3(2N+1)$-dimensional vectors

[TABLE]

We characterize the $4(2N+1)$ physically acceptable solutions (non-zero eigenvalues) of Eq. (21) by a superscript $s=+(-)$, which denotes waves propagating or decaying in the positive (negative) $z$ direction, a subscript $p=1,2$ indicating the two eigen-polarizations, and another subscript $\nu=-N,-N+1,\cdots,N$ which labels the different eigenmodes. The electric and magnetic field components of these modes, with unit amplitude, wave vector $\mathbf{q}_{p\nu}^{s}=\mathbf{q}_{\parallel}+q_{p\nu z}^{s}\widehat{\mathbf{z}}$, and polarization eigenvectors ${\mathbf{e}}_{p\nu}^{s}(n)$ and ${\mathbf{h}}_{p\nu}^{s}(n)$ are given by

[TABLE]

In the case of static homogeneous media, Eq. (21) is reduced to a set of $2N+1$ independent eigenvalue equations christofi

[TABLE]

one for each value of $n=-N,-N+1,\cdots,N$, where $\boldsymbol{\epsilon}$ are the elements of the corresponding static dielectric tensor. From the physically acceptable eigenvectors ${\mathbf{e}}_{pn}^{s}$ and ${\mathbf{h}}_{pn}^{s}$ of Eqs. (34), we construct the eigenvectors which appear in Eqs. (25) as follows: $\mathbf{e}_{p\nu}^{s}(n)=\delta_{\nu n}{\mathbf{e}}_{pn}^{s}$, $\mathbf{h}_{p\nu}^{s}(n)=\delta_{\nu n}\mathbf{h}_{pn}^{s}$, while the corresponding wave vectors are $\mathbf{q}_{p\nu}^{s}=\mathbf{q}_{pn}^{s}$.

We now consider a planar interface between two different homogeneous, in general time-periodic with the same periodicity, media: (1) on the left and (2) on the right side of the interface. The interface is perpendicular to the $z$ axis, which is directed from left to right, at $z=0$. Let us assume an eigenmode of the EM field given by Eqs. (25), with amplitude equal to unity, wave vector $\mathbf{q}^{+(1)}_{p^{\prime}\nu^{\prime}}$, and polarization described by ${\mathbf{\underline{e}}}^{+(1)}_{p^{\prime}\nu^{\prime}}$, ${\mathbf{\underline{h}}}^{+(1)}_{p^{\prime}\nu^{\prime}}$, incident on the interface from the left. Scattering at the interface gives rise to transmitted and reflected waves with wave vectors $\mathbf{q}^{+(2)}_{p\nu}$ and $\mathbf{q}^{-(1)}_{p\nu}$, $p=1,2$, $\nu=-N,-N+1,\cdots,N$, respectively, characterized by the associated polarization eigenvectors. As a result of the common temporal periodicity of both media on either side of the interface (this comprises the case where the one medium is time-invariant), the Floquet quasi-frequency $\omega$ remains constant, similarly to the Floquet quasi-momentum, or else the Bloch wave vector, when there is spatial periodicity. comphy1 ; comphy2 In addition, translation invariance parallel to the interface implies that $\mathbf{q}_{\parallel}$ remains the same for all of these waves. We denote the relative amplitudes of the transmitted and reflected waves (with respect to the amplitude of the incident wave) by $S_{p\nu;p^{\prime}\nu^{\prime}}^{++}$ and $S_{p\nu;p^{\prime}\nu^{\prime}}^{-+}$, respectively. In a similar manner we can proceed for an EM wave incident on the interface from the right and obtain transmitted and reflected waves of relative amplitudes $S_{p\nu;p^{\prime}\nu^{\prime}}^{--}$ and $S_{p\nu;p^{\prime}\nu^{\prime}}^{+-}$, respectively. Continuity of the tangential components of the EM field at the interface, in each case, yields $4(2N+1)$ linear systems, each of dimension $4(2N+1)\times 4(2N+1)$, which can be written in compact form as

[TABLE]

for $p^{\prime}=1,2$, $\nu^{\prime}=-N,-N+1,\cdots,N$, and $s^{\prime}=+(-)$. We note that $\mathrm{m}=1$ or 2 denotes the appropriate medium, which corresponds to $s^{\prime}=+$ or $-$, respectively.

In order to evaluate the scattering properties of multilayer structures, it is convenient to express the waves on the left (right) side of a given interface about different points, at a distance $-\mathbf{d}_{1}$ ($\mathbf{d}_{2}$) from a center at the interface, so that all backward and forward propagating (or evanescent) waves in medium m between two consecutive interfaces refer to the same (arbitrary) origin at $\mathbf{R}_{\mathrm{m}}$. This requires redefining the transmission and reflection coefficients of the interface by multiplying with the appropriate phase factors, as follows comphy1 ; comphy2

[TABLE]

In the particular case of an interface between two static homogeneous media, Eq. (39) is reduced to a set of $2N+1$ independent linear equations christofi

[TABLE]

one for each value of $n=-N,-N+1,\cdots,N$. Therefore, the transmission and reflection matrices become diagonal in $n$, which reflects frequency conservation.

The transmission and reflection matrices of a pair of consecutive interfaces, $\mathrm{i}$ and $\mathrm{i}+1$, are obtained by properly combining those of the two interfaces so as to describe multiple scattering to any order. This leads to the following expressions after summing the infinite geometric series involved comphy1 ; comphy2

[TABLE]

It is obvious that one can repeat the above process to obtain the transmission and reflection matrices $\mathbf{Q}$ of three consecutive interfaces, by combining those of the pair of the first interfaces with those of the third interface, and so on, by properly combining the ${Q}$ matrices of component units, one can obtain the ${Q}$ matrices of a slab which comprises any finite number of interfaces.

In general, we want to describe a multilayer slab, which comprises time-periodic media, sandwiched between two semi-infinite homogeneous and isotropic static media: (A) on its left and (B) on its right, characterized by scalar relative electric permittivities $\epsilon_{\mathrm{A}}$, $\epsilon_{\mathrm{B}}$ and magnetic permeabilities $\mu_{\mathrm{A}}$, $\mu_{\mathrm{B}}$, respectively. The transmission and reflection coefficients of the structure are obtained using the pairing procedure described by Eq. (46). For a wave incident on the slab from the left, of the form $\boldsymbol{\mathcal{E}}_{\mathrm{in}}(\mathbf{{r}},t)=\exp{\{i[\mathbf{q}_{n^{\prime}}^{+\mathrm{(A)}}\cdot\mathbf{r}_{\mathrm{A}}+n^{\prime}\Omega t)]\}}\mathbf{e}_{p^{\prime}n^{\prime}}^{+(\mathrm{A})}$, the transmitted and reflected fields are given by

[TABLE]

where $\mathbf{{r}}_{\mathrm{m}}=\mathbf{r}-\mathbf{R}_{\mathrm{m}}$, $\mathbf{q}_{n}^{\pm(\mathrm{m})}=\mathbf{q}_{\parallel}+q_{nz}^{\pm(\mathrm{m})}\widehat{\mathbf{z}}$, $q_{nz}^{\pm(\mathrm{m})}=\pm\sqrt{k_{n}^{2}\epsilon_{\mathrm{m}}\mu_{\mathrm{m}}-q_{\parallel}^{2}}$, with $\mathrm{m}=\mathrm{A,B}$ denoting the appropriate medium. We note that, in an isotropic medium such as A or B, the polarization eigenmodes are degenerate; this is why the polarization index $p$ (or $p^{\prime}$) has been dropped from the corresponding wave vectors. Having calculated the transmitted and reflected fields, we obtain the transmittance $\mathcal{T}$ and reflectance $\mathcal{R}$ of the slab. These are defined as the ratio of the flux of the transmitted or of the reflected field to the flux of the incident wave, respectively. Integrating the Poynting vector over the $x-y$ plane, on the each time appropriate side of the slab, and taking the average over a long time interval $\tau\rightarrow\infty$, we obtain

[TABLE]

provided that $\omega\neq n\Omega/2$, zurita a condition that will be always satisfied in our case. Because of the time variation of the permittivity tensor, the EM energy is not conserved even in the absence of dissipative (thermal) losses. In this case, $\mathcal{A}=1-\mathcal{T}-\mathcal{R}>0\,(<0)$ means energy tranfer from (to) the EM to (from) the spin-wave field.

III.2 Results

We shall now demonstrate the validity of the fully dynamic time-Floquet scattering-matrix method described above, and compare with the adiabatic approximation, on the structure of Fig. 1(a). We recall that this structure supports two neighboring defect modes in the lowest Bragg gap. As in Sec. II, we consider p-polarized incident light, with frequency $\omega_{0}a/c=1.88375$ and $q_{\parallel}a=1.2$ that corresponds to the higher-frequency defect mode in Fig. 1(c), and assume excitation of a second-order perpendicular standing spin wave of frequency $\Omega_{2}$ with a relative amplitude $\eta=0.1$. The presence of the spin wave induces a spatial and time-periodic variation in the permittivity tensor [see Eq. (4)]. To deal with spatial inhomogeneity we consider 50 elementary homogeneous sublayers in the magnetic film while the dynamic response is described assuming a truncation order $N=20$ in the Fourier series expansions involved. These parameters ensure excellent convergence of the numerical results.

Our method proceeds as schematically described in Fig. 2. We first calculate the eigenmodes of the EM field in each of the 53 (infinite) homogeneous media involved, i.e., the 50 dynamic media that represent the Bi:YIG film under the action of the spin wave and the 3 static media (air, $\mathrm{SiO}_{2}$, $\mathrm{TiO}_{2}$), solving the appropriate eigenvalue-eigenvector equations, Eq. (21) and Eq. (34), respectively. Then, we evaluate the transmission and reflection matrices $\mathbf{Q}$ of the consecutive interfaces by Eq. (39), or Eq. (45) in the simple case of an interface between two static media, and Eqs. (40). At each step of the iterative procedure, the $Q$ matrices of the entire slab built-up up to the current interface are obtained according to Eqs. (46). Therefore, at the end of the procedure, the $Q$ matrices of the whole structure are readily available and the corresponding transmittance and reflectance are calculated from Eqs. (49) and (50), respectively.

Figure 3 displays the total (transmitted plus reflected) intensity of the elastic and inelastic outgoing light beams, associated with absorption and emission of zero, one, and two magnons by a photon, as a function of the spin-wave frequency. It can be seen that, at spin-wave frequencies smaller than the width of the optical resonance ($\Gamma a/c\sim 10^{-6}$), we recover the results of the quasi-static adiabatic approximation. The calculated intensities of the elastic and inelastic outgoing light beams are in perfect agreement with those reported in Table 1, while the actual value of the spin-wave frequency is, indeed, immaterial. In Figs. 3(a) and (b) we see that, as $\Omega_{2}$ increases, the quasi-static adiabatic approximation breaks down and a fully dynamic approach, which properly takes into the account the initial and final photon states is required. If $\Omega_{2}$ is of the order of the dynamical frequency shift of the optical resonance $(\Delta\omega a/c\sim 10^{-6})$, the intensities of the scattered light beams exhibit an oscillatory behavior. By further increasing $\Omega_{2}$, inelastic light scattering processes are suppressed, since the final states are within the gap, where the optical density of states is very low, and we are essentially left with only elastic scattering. At higher values of $\Omega_{2}$, the magnon frequency can match the frequency difference between the two optical resonances shown in Fig. 1(c) (triple-resonance condition). In this case, one-magnon emission processes are favored, leading to enhanced intensities of the corresponding $n=1$ inelastically transmitted and reflected light beams, for both polarization-conserving and polarization-converting processes. At the same time, elastic scattering is considerably reduced, by about 45%, while the other inelastic processes are also resonantly affected, though to a much lower degree. Overall, there is an excess number of magnons produced, which is manifested as a peak in the total optical absorption spectrum as depicted in Fig. 3(c), thus clearly indicating a resonant energy transfer from light to spin waves. This effect, also, cannot be accounted for by the adiabatic approximation. It is worth-noting that the triple-resonance condition can be accomplished by many-magnon emission processes as well, however the effects are considerably weaker as shown in Fig. 3(d).

So far, quantities are expressed in terms of the lattice parameter $a$. Taking $a=300$ nm, the optical wave considered is tuned to a frequency $\omega_{0}/(2\pi)=300$ THz, i.e., a wavelength $\lambda_{0}=1\;\mu\mathrm{m}$. Correspondingly, the spin-wave frequency $\Omega_{2}/(2\pi)$ in Figs. 3(a) and (b) ranges from about $10^{-2}$ GHz to $10^{2}$ GHz. It should be pointed out that the different features observed can be tuned in frequency by appropriately choosing the geometric parameters of the structure. For example, the triple-resonance condition can be shifted to lower frequencies, say to the few GHz region which is typical for perpendicular standing spin waves, navabi2017 by increasing the lattice parameter $a$ and, accordingly, the optical wavelength. This will also require to use a lower Faraday coefficient, zvezdin_book leading to a smaller separation between the two optical resonances involved. We recall that the amplitude of the multiple-magnon processes increases with the coupling strength $\Delta\omega/\Gamma$. In this respect, it should be noted that the lifetime of a resonant defect mode ($\propto 1/\Gamma$) can be made arbitrarily long by increasing the number of unit cells in the Bragg mirrors, provided that absorptive losses are negligible, while $\Delta\omega$, can be widened by raising the spin-wave amplitude.

IV Conclusion

In summary, we formulated a fully dynamic time-Floquet scattering-matrix method, with arbitrary-order accuracy, in order to describe the interaction of visible and infrared light with spin waves in layered optomagnonic structures. The method is robust, computationally efficient and, contrary to brute-force time-domain methods, is ideally suited for fast and accurate calculations of the optical response of periodically driven systems where two very different time scales are involved. The applicability of the method is demonstrated on a specific model of a dual photonic-magnonic cavity, which operates in the strong-coupling regime. The cavity is formed by sandwiching a magnetic Bi:YIG film between two dielectric Bragg mirrors and is subject to continuous excitation of a perpendicular standing spin wave. For spin-wave frequencies smaller than the width of the optical resonance, we recover the results of the adiabatic approach, namely enhanced modulation of light through multi-magnon absorption and emission mechanisms, beyond the linear-response approximation. For higher spin-wave frequencies where the adiabatic approach breaks down, our fully dynamic calculations provide evidence for the occurrence of remarkable effects, such as enhanced inelastic light scattering in the nonlinear strong-coupling regime when a triple-resonance condition is fulfilled and resonant energy transfer between the photon and the magnon fields. We note that similar effects should occur, also, for perpendicular standing spin waves in in-plane magnetized films as well as for surface Dammon-Eshbach and backward volume waves with an in-plane propagation wave vector. In the latter cases, however, the spin wave induces a periodic variation of the electric permittivity tensor both in time and in the propagation direction. The implementation of a space-time-Floquet scattering-matrix method for a fully dynamic description of such cases is currently in progress.

Acknowledgements.

Petros-Andreas Pantazopoulos is supported from the General Secretariat for Research and Technology and the Hellenic Foundation for Research and Innovation through a PhD scholarship (No. 906). *

Appendix A

In order to gain further insight into the accuracy and limits of validity of the quasi-static adiabatic approximation, let us consider a generic one-dimensional quantum mechanical model: A particle of mass $m$ in a periodically oscillating delta-function potential at the origin, in which case Schrödinger’s equation reads

[TABLE]

We note that this model was studied, though in a different perspective, also by Martinez and Reichl. martinez We want to describe scattering of the particle, when it impinges with (kinetic) energy $E$ on the potential from $x<0$. It is straightforward to show from Eq. (51) that the wave function of the incident particle, with amplitude equal to unity, is

[TABLE]

for $x<0$, where

[TABLE]

We now assume that the variation of the potential is relatively slow, i.e., $\Omega\ll E/\hbar\equiv\omega$, and undertake an adiabatic approach by solving the problem with the potential strength kept fixed and then, at the end of the calculation, allow it to change in successive snapshots within a period $T=2\pi/\Omega$. This implies transmitted and reflected waves of the form

[TABLE]

for $x>0$ and

[TABLE]

for $x<0$, which can be readily obtained by applying the well-known results for the static delta-function potential. griffiths

Expanding the periodically varying transmission and reflection coefficients into Fourier series we have

[TABLE]

for $x>0$ and

[TABLE]

for $x<0$, where

[TABLE]

and $A^{0}_{\mathrm{rf}}(n)=A^{0}_{\mathrm{tr}}(n)-\delta_{n0}$. The Fourier coefficients $A^{0}_{\mathrm{tr}}(n)$ can be evaluated analytically by contour integration, which yields

[TABLE]

where $\alpha_{0}\equiv k_{0}\hbar^{2}/(mg)$, and the actual value of $\Omega$ is of course immaterial.

Equation (56) and (57) tell us that the transmitted and reflected wave functions consist of an infinite number, $n=0,\pm 1,\pm 2,\ldots$, of plane wave beams of the same wave number $k_{0}$, frequencies $\omega-n\Omega$, and amplitudes $A^{0}_{\mathrm{tr}}(n)$ and $A^{0}_{\mathrm{rf}}(n)$, respectively. The corresponding transmittance, $\mathcal{T}$, and reflectance, $\mathcal{R}$, are defined as the time-averaged [$\lim_{\tau\longrightarrow\infty}\frac{1}{\tau}\int_{0}^{\tau}dt(\cdots)$] probability current [$\mathbf{J}=\frac{\hbar}{m}\mathrm{Im}(\Psi^{*}\nabla\Psi)$] associated with the transmitted and reflected waves in the direction of flow, $\widehat{\mathbf{x}}$ and $-\widehat{\mathbf{x}}$, respectively, relative to that of the incident wave. Using Eqs.(57) we obtain

[TABLE]

[TABLE]

We note that, considering the Fourier series expansions $\alpha_{0}/(\alpha_{0}+i\cos\phi)=\sum_{n=0,\pm 1,\pm 2,\ldots}A^{0}_{\mathrm{tr}}(n)\exp(in\phi)$ and $-i\cos\phi/(\alpha_{0}+i\cos\phi)=\sum_{n=0,\pm 1,\pm 2,\ldots}A^{0}_{\mathrm{rf}}(n)\exp(in\phi)$, Parseval’s theorem implies $\mathcal{T}+\mathcal{R}=1$, as expected for a nondissipative adiabatically changing system.

On the other hand, a fully dynamic solution to the scattering problem described above can be obtained by a Fourier expansion of the incident, transmitted, and reflected wave functions, prior to applying the appropriate boundary conditions. We have

[TABLE]

where

[TABLE]

Continuity of the total wave function $\Psi(x,t)$ ($\Psi=\Psi_{\mathrm{in}}+\Psi_{\mathrm{rf}}$ for $x<0$ and $\Psi=\Psi_{\mathrm{tr}}$ for $x>0$) requires

[TABLE]

for $n=0,\pm 1,\pm 2,\ldots$. Moreover, integrating Eq. (51), we obtain a discontinuity in the first derivative of the total wave function

[TABLE]

which, with the help of Eq. (64), leads to

[TABLE]

for $n=0,\pm 1,\pm 2,\ldots$, where $\alpha_{n}\equiv k_{n}\hbar^{2}/(mg)$.

Equation (66) provides the fully dynamic solution for the transmitted wave function while the reflected wave function can be readily obtained from Eq. (64). The corresponding transmittance and reflectance are given by

[TABLE]

We note that here, contrary to the quasi-static adiabatic description, for $n>\omega/\Omega$ we obtain evanescent plane wave beams that do not contribute to the transmittance and the reflectance because $\mathrm{Re}(k_{n})=0$.

It would be now interesting to compare the results of the adiabatic approximation with those of the fully dynamic description in the limit $\Omega\ll\omega$. Taking $\alpha_{n}=\alpha_{0}$, Eq. (66), for $A_{\mathrm{in}}(n)=\delta_{n0}$ as appropriate to the case under consideration, yields directly $A_{\mathrm{tr}}(n)$ equal to the middle ($n=0$) column of the inverse of a tridiagonal Toeplitz matrix with the elements of its main diagonal equal to unity and the elements of its first diagonals below and above the main diagonal equal to $i/(2\alpha_{0})$. Using the explicit expressions for the inverse of a general tridiagonal Toeplitz matrix of dimensions $(2N+1)\times(2N+1)$, fonseca we obtain in our case here

[TABLE]

for $n=0,\pm 1,\pm 2,\ldots,\pm N$, where $U_{p}(z)$ are the Chebyshev polynomials of the second kind abramowitz_book given by

[TABLE]

with $r_{\pm}=z\pm\sqrt{z^{2}-1}$. Interestingly, restricting to $|n|\ll N$, in the limit of matrix dimensions tending to infinity, Eq. (68) leads strictly to Eq. (59): $A_{\mathrm{tr}}(n)=A^{0}_{\mathrm{tr}}(n)$.

For $\alpha_{n}\neq\alpha_{0}$, the matrix to be inverted is a tridiagonal Jacobi (non-Toeplitz) matrix and an elegant, concise formula for its inverse is also available. fonseca Though in this general case no simple analytic solution can be obtained, considering $\Omega/\omega\leq 10^{-3}$, numerical calculations show that the absolute relative difference between the results of the adiabatic approximation and the fully dynamic approach does not exceed $10^{-3}$. The calculated partial transmittance associated with the different plane wave beams, $\mathcal{T}_{n}$, considering sufficiently large tridiagonal matrices to ensure convergence, is depicted in Fig. 4 for some characteristic cases. It can be seen that strong deviations from the results of the adiabatic approximation appear only for simultaneously large-amplitude ($\alpha_{0}=0.05$) and fast ($\Omega/\omega=0.1$) variations. It is worth-noting that, in quantum mechanics, the transmittance and reflectance usually refer to probability currents, as in the example considered here, and not to energy currents as commonly described in electromagnetism by the Poynting vector. jackson Probability current conservation ($\mathcal{T}+\mathcal{R}=1$) is here always satisfied (to machine accuracy $\sim 10^{-16}$) in both adiabatic and fully dynamic calculations, as expected for the Hermitian Hamiltonian system under study. In addition to that, the adiabatic approximation should also preserve energy current conservation, i.e., since particles in the $n$-th beam have an energy $\hbar(\omega-n\Omega)$, the quantity $\sum_{n=0,\pm 1,\pm 2,\ldots}n\Omega(\mathcal{T}_{n}+\mathcal{R}_{n})$ must vanish identically, which is indeed ensured by the symmetry properties $\mathcal{T}_{-n}=\mathcal{T}_{n}$ and $\mathcal{R}_{-n}=\mathcal{R}_{n}$ implied by Eq. (59). Of course, these symmetry properties are not satisfied in the fully dynamic description of the problem, as can be clearly seen in the bottom right diagram of Fig. 4, while small deviations from the perfectly symmetric partial transmittance also exist in the other diagrams but they are not discernible in the scale of the figure. This is translated to the absence of energy conservation, as expected in the actual time-varying system under consideration. It becomes, therefore, clear from the example elaborated above that the quasi-static adiabatic approximation provides a reasonably good description of wave scattering in systems with a relatively slow (but not necessarily weak) periodic variation with time, at least to the extent where the scattering spectral features do not exhibit strong changes over a frequency range of the order of $\Omega$ such as in inelastic scattering between resonant modes.