Design of an optomagnonic crystal: towards optimal magnon-photon mode matching at the microscale

We put forward the concept of an optomagnonic crystal: a periodically patterned structure at the microscale based on a magnetic dielectric, which can co-localize magnon and photon modes. The co-localization in small volumes can result in large values of the photon-magnon coupling at the single quanta level, which opens perspectives for quantum information processing and quantum conversion schemes with these systems. We study theoretically a simple geometry consisting of a one-dimensional array of holes with an abrupt defect, considering the ferrimagnet Yttrium Iron Garnet (YIG) as the basis material. We show that both magnon and photon modes can be localized at the defect, and use symmetry arguments to select an optimal pair of modes in order to maximize the coupling. We show that an optomagnonic coupling in the kHz range is achievable in this geometry, and discuss possible optimization routes in order to improve both coupling strengths and optical losses.


I. INTRODUCTION
Progress in fundamental quantum physics has by now established a basis for developing new technologies in the fields of information processing, secure communication, and quantum enhanced sensing. In order to perform these tasks, physical systems are needed which are capable of processing, storing, and communicating information in a quantum coherent manner and with a high fidelity. Similar to the classical realm, accomplishing this goal requires different degrees of freedom and efficient couplings between them, giving rise to hybrid systems. In this context, systems at the mesoscopic scale (with dimensions ranging from nanometers to microns) are specially interesting since their collective degrees of freedom can be tailored [1]. An important and successful example of these mesoscopic hybrid systems are optomechanical systems [2], where light couples to mechanical motion. Seminal experiments in these systems have demonstrated extra-sensitive optical detection of small forces and displacements [3][4][5][6][7][8][9][10], manipulation and detection of mechanical motion in the quantum regime with light [11][12][13], and the creation of nonclassical light and mechanical motion states [11,14,15].
In the recent years, the family of hybrid quantum systems has been extended by incorporating magnetic materials, where the collective spin degree of freedom can be exploited. For example in spintronics [16], information is carried by spins (as opposed to electrons) in order to remove Ohmic losses and to increase memory and processing capabilities [17,18]. An ultimate form of spintronics is the new field of quantum magnonics [19,20], where superconducting quantum circuits couple, via microwave fields in a cavity, coherently to magnetic collective excitations (magnons) [21,22]. Such systems are promising for generating and characterizing non-classical quantum magnon states [19,23,24] and for developing microwaveto-optical quantum transducers for quantum information processing [25,26]. The coherent coupling of magnons to optical photons has also been demonstrated in recent experiments [27][28][29][30][31], in what have been denominated optomagnonic systems [26,[32][33][34][35][36][37].
In current experiments exploring optomagnonics, the ferrimagnetic dielectric Yttrium Iron Garnet (YIG) is used as the magnetic element, since YIG presents low absorption and a large Faraday constant in the infrared (α = 0.069 cm −1 and θ F = 240 deg/cm @ 1.2 µm). The coupling between spins and optical photons is an sec-ond order processes involving spin-orbit coupling and it is generally small. This can be enhanced using a well polished sphere that acts as an optical cavity, trapping the photons by total internal reflection in order to effectively enhance the spin-photon coupling. The coupling, however, still remains too small for applications. This is due to the large size of the used YIG spheres (the coupling increases as the volume of the cavity decreases [32]), with radius of the order of hundreds of microns, and, concomitantly, the large difference between the optical and the magnetic mode volume (V mag V opt ), by which most of the magnetic mode volume does not participate in the coupling. This can be partially mitigated by making smaller cavities [38], but care has to be taken both to obtain a good mode matching and to retain a good confinement of the optical mode in order to minimize radiation losses. Recent proposals have investigated one dimensional layered structures to this end [35,39].
In order to tackle these issues, we propose an optomagnonic array at the microscale, which acts simultaneously as a photonic crystal [40], determining the optical properties of the structure, and as a magnonic crystal [41][42][43] with tailored magnetostatic modes. Our proposal is inspired in the success of this approach for optomechanical crystals, which can be designed such as to enhance the phonon-photon coupling by many orders of magnitude . In our case, we use similar concepts in order to design the coupling between photonic and magnonic modes. Although similar conceptually, magnetic materials present new challenges for the design, due to the complexity of the magnon modes.
Photonic crystals are the basis for many novel applications in quantum information, and are of high interest due to their ability to guide [68][69][70][71][72] and confine [73][74][75][76][77][78] light, allowing for example to enhance non-linear optical interactions [79][80][81][82]. In turn, magnonic crystals can be designed to create reprogrammable magnetic band structures [83], to act as band-pass or band-stop filters, or to create single-mode and bend waveguides [84][85][86][87]. Additionally these crystals can be used for spin wave computing via logical gates [88][89][90]. An advantage of magnonic crystals is their scalability, their low energy consumption, and possibly faster operation rates [42,90,91]. Together with the state of the art in optomagnonics detailed above, this provides a great incentive to explore the possibility of an optomagnonic crystal, combining both photonic and magnetic degrees of freedom.
Specifically, we consider an optomagnonic crystal consisting of a dielectric magnetic slab (YIG in our simulations) with a periodic array of holes along the slab and with an abrupt defect in the middle. The repeated holes at each side of the defect act as a Bragg mirror for both optical and magnetic modes, localizing them in the region of the defect (see Fig. 1). We show that this structure can co-localize photonic and magnonic modes, and explore how the symmetry of the modes can be used to optimize their coupling. We find that coupling rates in the range of kHz are achievable in these structures, and that optimization of the geometry can lead to higher coupling values, indicating the promise of this approach. Further optimization is nevertheless needed to improve the decay rates, in particular the optical quality factor is low compared to the state of the art in non-magnetic structures (where Silicon is used as the dielectric).
This manuscript is organized as follows. In Sec. II we derive the general expression for the coupling of magnons to optical photons and discuss the normalization of the modes required to find the photon-magnon coupling at the single quanta level, denominated optomagnonic coupling. The remaining sections refer to the numerical method for evaluating this coupling. For our simulations we choose YIG as the magnetic material, in line with the material of choice in experiments. In Sec. III we discuss the properties of the proposed structure as a photonic crystal. In Sec. IV, in turn, we investigate its properties as a magnonic crystal. Sec. V combines the results in order to numerically evaluate the optomagnonic coupling for appropriately chosen confined modes. For concreteness, we focus on the coupling between one single magnon mode and one single optical mode. Sec. VI is devoted to a discussion on how the structure can be optimized and presents results for an optimized geometry. The conclusions and an outlook are presented in Sec. VII. The supplementary material contains further details of the analytic calculations and of the simulations.

II. OPTOMAGNONIC COUPLING
In this section we derive the theoretical expression for the coupling rate between magnons and photons. The instantaneous electromagnetic energy is [92] with D the displacement field, E the electric field, B the magnetic induction, and H the magnetic field. In complex representation D = (D+D * )/2 and E = (E+E * )/2, and similar for the magnetic induction and field. The effect of the magnetization M is to modify the displacement field as D(r, t) =ε[M(r, t)]·E(r, t) where the components of the permittivity tensorε are [93,94] with ε 0 the vacuum permittivity, ε r the relative permittivity, ijk the Levi-Cevita tensor, and {f F , f C } material dependent magneto-optical constants. At optical frequencies the second term in Eq. (1) can be neglected [94,95], being smaller than the first by the fine structure constant squared, and the permeability of the material can be set to the vacuum permeability µ 0 . The magneto-optical constants can be related to the Faraday rotation θ F and the Cotton-Mouton ellipticity θ C per unit length as θ F = ω/(2c √ ε r )f F M s and with c the velocity of light in vacuum and M s the saturation magnetization.
We are interested in how light couples to the fluctuations of the magnetization around the static ground state. We consider norm-preserving small fluctuations, where the ground state satisfies M 0 · M 0 = M 2 s and the fluctuations are perpendicular to local equilibrium magnetization δM · M 0 = 0. In complex notation δM = [M + M * ]/2. The correction to the electromagnetic energy stemming from the interaction between the light field and the magnetization can be rewritten as ignoring E(r, t)·D(r, t) and E * (r, t)·D * (r, t) in the rotating wave approximation. Inserting the relation between the displacement and the electric field along with the permittivity in Eq.
is the Faraday contribution and the Cotton-Mouton one. We have used the dyadic notation and neglected all terms that represent a constant energy shift or that are higher order in δM. Quantizing this expression leads to the optomagnonic coupling Hamiltonian. By assuming that the magnetic material acts as an optical cavity, the electric field of the light can be quantized by using the annihilation (creation) operatorâ with E ( * ) α the mode shape, and α the mode index. We note that we identified E(r, t) with 2 E + (r, t) from the well known quantization expression of the electric field [96] E(r, t) = E + (r, t) + E − (r, t) In order to find the coupling per photon, we normalize the electromagnetic field amplitude to one photon over the electromagnetic vacuum [97] ω α 2ε 0 =ˆdr ε r (r)|E α (r)| 2 .
The spin waves can be quantized as whereb ( †) γ annihilates (creates) one magnon, δm ( * ) γ is the mode shape, and γ the mode index. We note that as in the optical case we identified M(r, t) with 2 M + (r, t) from the magnetic quantization expression In order to normalize the amplitude of the magnetic fluctuations to one magnon, we use the following expression derived in the supplementary material (appendix A) with g the g-factor, µ B the Bohr magneton, and m 0 = M 0 /M s . This expression is valid for arbitrary magnetic textures and it is consistent with the normalization derived previously for a uniform ground state [98]. The quantized optomagnonic energy, neglecting the constant energy shifts, leads to the coupling Hamiltonian with the coupling constant G αβγ = G F αβγ + G C αβγ , where are, respectively, the Faraday and Cotton-Mouton components of the optomagnonic coupling constant, being λ n the light wavelength inside the material.
In this work we focus on the coupling between a given magnon mode and a given optical mode hosted by the 1D optomagnonic crystal shown in Fig. 1. Hence we set α = β and drop the indices in the following. For concreteness, we focus on G F as the analysis for G C is analogous. G F is proportional to the overlap of the magnon's spatial distribution with the electric component of the optical spin density defined as [99] The optical spin density is finite only for fields with certain degree of circular polarization, and points perpendicular to the plane of polarization.

III. PHOTONIC CRYSTAL
Photonic crystals are engineered structures which, by proper shape design, can confine light to a specific region. These are formed by low-loss media exhibiting a periodic dielectric function ε(r), with a discrete translational symmetry ε(r) = ε(r + R) for any R = na with n an integer and a the lattice constant given by the imposed periodicity.
Photonic band gaps arise at the edges of the Brillouin zone (BZ) k = π/a due to the periodicity imposed by the susceptibility of the crystal on the electric field, with wavelength λ = 2a (corresponding to the edge of the BZ). For example, in a 1D photonic crystal (see Fig. 2a) the symmetry of the unit cell around its center implies that the nodes of the standing light wave must be centered either at each low-ε layer or at each high-ε layer. The latter necessarily has lower energy than the former, resulting in a band gap. The position of the photonic band gap is given by the mid-gap frequency at the BZ edge. In the case of two materials with refractive indices n 1 and n 2 and thicknesses d 1 and d 2 = a − d 1 , the normal incidence gap is maximized for n 1 d 1 = n 2 d 2 . In this case the mid-gap frequency is given by [40] with n 1 = √ ε 1 , n 2 = √ ε 2 . The corresponding vacuum wavelength λ mg = (2πc)/ω mg thereby satisfies the relations λ mg /n 1 = 4d 1 and λ mg /n 2 = 4d 2 meaning that the individual layers are a quarter-wavelength thick. An input at frequencies within the photonic band gap is reflected entirely except for an exponentially decaying tail inside the crystal. Thus, two of such crystals can be used to create a Fabry-Perot like cavity. More concretely, as shown in Fig. 2b, a defect in the form of a layer with a different width breaking the symmetry of the crystal may permit localized modes in the band gap by consecutive reflection on both sides. Since the light is localized in a finite region, the modes are quantized into discrete frequencies. We note that the degree of localization is the largest for modes with frequencies near the center of the gap [40].
For our purposes, we consider a geometry in which the variation of the permittivity is attained by holes carved into a dielectric slab (see Fig. 1). The typical material  Figure 2. General structure of a 1D photonic crystal and mode localization at a defect: (a) 1D photonic crystal consisting of periodic layers alternated by the lattice constant a with different dielectric constants ε1 > ε2 and widths d1 and d2. (b) A defect breaks the symmetry and can pull a band-edge mode into the photonic band gap. Since a mode in the band gap cannot propagate into the structure, the light is Bragg-reflected and is thus localized.
used for photonic crystals is silicon due to its high refractive index at optical frequencies, ε r = 12. We use instead YIG for our study, which is a dielectric magnetic material transparent in the infrared range with ε r = 5 [100]. The lower dielectric constant reduces the confinement of the optical modes along the height of the slab, which is reflected in low optical quality factors as discussed below. This structure is a 1D photonic crystal, periodic in one direction (chosen to be thex-direction), with a band gap along this direction and which confines light through index guiding (a generalization of total internal reflection) in the remaining directions. In order to localize an optical mode in this structure we create a defect by increasing the spacing between the two middle holes, which pulls a mode into the band gap. We note that due to the (discrete) periodicity, the crystal only possesses an incomplete band gap and the localized mode can scatter to air modes [40]. We search for a localized mode in the infrared frequency range where YIG is transparent and presents low absorption [101,102]). Thus, the geometrical parameters of the crystal need to be chosen in such a way that the band gap lies in the desired frequency range. We choose a lattice constant of a = 450 nm which gives a mid-gap frequency of ω mg = 2π · 240 THz (corresponding to λ ≈ 1250 nm), using Eq. (17) with refractive indices of YIG n 1 = n YIG = √ 5 and of air n 2 = n air = 1. Note that we choose a lattice constant that allows us to work in the transparency frequency range for the optics, and which at the same time is small enough in order to reduce the computational cost of the micromagnetic simulations of the corresponding magnonic crystal in the next section. Using the relation d air n air = d YIG n YIG for a maximized normal incidence gap we find the optimal radius of the air holes as with d air = 2r air = a − d YIG , from which we obtain r air = 155.25 nm. In order to find the mode with the least losses, the defect width is optimized in order to localize the desired mode most effectively to the defect.
We find the optimal defect size, defined as the centerto-center distance between the two bounding holes [see Fig. 1], to be d = 731 nm, obtained by evaluating the transmission spectra as a function of the defect size (we used the electromagnetic simulation tool MEEP to this end [103]). In order to get a good quality factor of the localized mode we need to insert as many air holes as possible. For creating a compromise between short computational evaluation time (especially important for the magnetic simulations discussed later) and a good quality factor we chose N = 12 holes at each side of the defect. Therefore the investigated crystal is in total l = 11.75 µm long. The overall width of the wave guide is w = 600 nm and its height is h = 60 nm, again to keep the magnetic simulations (which we detail in the next section) feasible. Such a thin slab will not be good at confining the optical modes along its height, since it is much smaller than the light wavelength in the material. In order to increase the optical quality factor without influencing the magnetics we sandwich the crystal with two Si 3 N 4 layers with a height of h Si = 200 nm as proposed in [104]. Si 3 N 4 has an index of refraction similar to YIG (n Si3N4 = √ 4), so that the combined structure acts approximately as a single cavity for the light and its height is roughly half a wavelength, enough to provide a reasonable confinement. The presented simulations include these two extra layers.
We now turn to categorizing the photonic modes of the crystal by using its three mirror symmetry planes (see Fig. 3a). This imposes several restrictions on the mode shape and the mode polarization. We define the three mirror symmetry operationŝ In the following we restrict the discussion to transverse electric modes with an in-plane electric field, which are the modes of interest for the magnetic configuration we consider, as it will be clear from the next section. We note that structures made of a high-ε material with air holes favour a band gap for transverse electric modes [40], which is advantageous for our purposes. Unlike in 2D, in three dimensions we cannot generally distinguish between transverse electric (TE) and transverse magnetic (TM) modes. However, provided that the crystal has a mirror symmetry along its height (underσ z ), and that its thickness is smaller than the mode wavelength, the fields are mostly polarized in TE-like and TM-like modes [40].  Since a TE-like mode has a non-zero electric field in the plane of the crystal (xy-plane), both E x and E y cannot be odd as a function of z. From Eq. (19), this implies that the mode must be even underσ z . Similar symmetry considerations [40] show that a TE-like mode must satisfyσ For evaluating the optical modes we used the two finite element tools MEEP [103] and Comsol Multiphysics [106] (details see appendix C). The simulated band structure for TE-like modes in the considered photonic crystal shows a broad band gap in the infra-red frequency range with a nicely pulled defect band (see Fig. 4a) which is extended in frequency space resulting from the confinement in real space. The defect mode in the gap at the edge of the BZ has a frequency of ω opt ≈ 2π × 246 THz (obtained by Comsol, 205 THz, λ ≈ 1550 nm, according to MEEP) with a damping factor of κ opt ≈ 2π ×0.1 THz which gives a linewidth (FWHM) of γ opt = 2 · κ opt ≈ 2π × 0.2 THz. Using the values obtained by Comsol this gives an optical quality factor of Q = ω opt /(4πκ opt ) = 1250, which is in the expected range for this kind of geometry [107] (note that MEEP gives a roughly three times larger value). The obtained quality factor is however small compared to 1D crystals made of silicon with a smooth defect, where quality factors in the order of 10 4 − 10 6 can be achieved [46,48,50].
The corresponding mode shape in real space is shown in Fig. 4a. We see that the mode is nicely localized at the defect. Furthermore the E x component is even (odd) as a function of x (y), whereas the E y is even as a function of both x and y fulfilling the symmetry requirements for a TE-like mode given in Eqs. (20). Due to this symmetry, we can disregard the E z component here since E z ≈ 0. For the Faraday component of the optomagnonic coupling, the relevant quantity is the electric component of the optical spin density, S opt ∝ E * × E [see Eq. (14)]. S opt points mostly along z-direction, is odd as a function  of x and y, and is even along z (see Fig. 4a and 7a).

IV. MAGNONIC CRYSTAL
As photonic crystals control the flow of light, magnonic crystals can be used to manipulate the spin wave dynamics in magnetic materials. In general a magnonic crystal is made of a magnetic material with a periodic distribution of material parameters. Examples include the modulation of the saturation magnetization or the magnetocrystalline anisotropy, a periodic distribution of different materials, or the modulation by external parameters, such as an applied magnetic field [41][42][43]. Historically, magnonic crystals precede photonic crystals [108,109]. Unlike in photonic or phononic crystals, the band structure in magnonic crystals depends not only on the periodicity of the crystal but also on the spatial arrangement of the ground state magnetization, resulting in an additional degree of freedom. Hence the band structure depends on the applied external magnetic field, the relative direction of the wave vector, the shape of the magnet, and the magnetocrystalline anisotropy of the material [41][42][43]. In this section we study the properties of the crystal presented in Sec. III (see Fig. 3a), as a magnonic crystal.
In the following we consider magnetic excitations which are non homogeneous in space, and we focus only on systems in the presence of an external magnetic field saturating the magnetization in a chosen direction. In this case spin waves can be divided into three classes: if all spins precess uniformly in phase, the mode is homo-geneous and denominated the Kittel mode. If the dispersion is dominated by dipolar interactions (which is usually the case for wavelengths above 100 nm) the excitations are called dipolar spin waves. For wavelengths below 100 nm the exchange interaction dominates instead, giving rise to exchange spin waves. The frequencies of the dipolar spin waves lie typically in the GHz-regime, whereas the exchange spin waves have frequencies in the THz-regime. Since the size of the structure considered in this work is in the micrometre range, we will focus on dipolar spin waves. For this case, the modes can be classified further by their propagation direction with respect to the magnetization. For an in-plane magnetic field, modes with a frequency higher than the frequency of the uniform precession tend to localize at the surface and have a wave vector pointing perpendicular to the static magnetization M 0 and thus the external field, k ⊥ M 0 H ext (see Fig. 5a). These modes are called surface or Damon-Eshbach modes [110,111]. If the wave vector is parallel to the external field such that k M 0 H ext holds, the waves are called backward volume waves and their frequency is smaller than the frequency of the Kittel mode (see Fig. 5a). Finally, if the external field and the magnetization are normal to the crystal's plane and the wave vector lies in plane k ⊥ M 0 H ext these waves are called forward volume waves (see Fig. 5a) [43,112]. In the following we restrict the discussion to external fields which are applied in the plane of the crystal.
Similar to light modes in photonic crystals, magnon modes can also be localized within a certain region in the magnonic crystal. It is well known that the two di-mensional periodic modification of a continuous film, for example by the insertion of holes (denominated antidot arrays) can drastically change the behaviour of the spin waves [113,114]. In this case the modes have either a localized or extended character. The localized mode is a consequence of non-uniform demagnetization fields created by the antidots. These fields change abruptly at the edges of the antidots and act as potential wells for the spin waves [43]. Thus, the above designed crystal, which localizes the optical mode by the insertion of a defect, is also a good candidate for acting as a magnonic crystal localizing magnetic modes via the holes. Although the geometry of the crystal is optimized for the optics, it should be able to host and localize magnetic modes due to its shape and material (YIG). Therefore we do not change the crystal further and use this structure as a proof of principle. This implies that we expect considerable room of improvement with respect to the optomagnonic coupling rates obtained in this structure. YIG is a good choice for magnonics since it has the lowest spin wave damping when compared to other materials commonly used [42]. It is however difficult to pattern at the microscale, but recent advances in fabrication show great promise in this respect [115,116].
For concreteness, in the following we proceed to design the Faraday part of the optomagnonic coupling G F , see Eq. (14). Since G F is proportional to the overlap integral between the optical spin density and the magnon mode, we search for a magnon mode with the same symmetries as the optical spin density, in order to get the highest possible overlap. Like in the optical case, the magnonic crystal has three mirror symmetry planes (z = 0, y = 0, x = 0). However, the external applied magnetic field saturating the magnetization breaks two of these symmetries and thus only the mirror symmetry w.r.t. the plane perpendicular to the external field remains (see Fig. 5b). Note that the magnetization is a pseudo vector and its components perpendicular to the mirror does not change. Thus, the mirror operation is inverted from Eq. (19), Since the optical spin density pointing alongẑ is odd as a function of x, we require δm z to be odd as well and consequently δm to be even underσ M x . Additionally, a π-rotation around thex-axis symmetry remains unbroken Invoking again the symmetries of the optical spin density (odd as a function of y and even with z) we consider modes with even rotational symmetry. We note that due   Fig. 1. Since the external magnetic field breaks two mirror symmetry planes only the mirror symmetry plane normal to the saturation direction remains. Additionally a π-rotation symmetry around the saturation axis is present.
to the different symmetries respected by the photon and magnon modes, we choose the symmetries of the modes in such a way that they preferably match in the xy-plane, which is the most relevant dimension for thin structures. In this case, the symmetries of the optical and the magnetic mode along the height do not necessarily match. For thin films however they do, see Fig. 7b. Since spin waves are excited by an external magnetic pulse which controls the direction of the wave vector k, the pulse also breaks the mirror symmetries of the crystal. Therefore we focus on a setup which conserves the relevant mirror symmetry, and only excite backward volume waves where the external saturation field and the wave vector of the mode are parallel and lie in the plane of the crystal. We note that this configuration is also the most favourable one from an experimental point of view, and additionally the configuration most likely used in magnonic devices [112].
We evaluated the magnetization dynamics numerically by means of the finite difference tool MuMax3 [117] which solves for the Landau-Lifshitz-Gilbert equation of motion for the local magnetization vector (details see appendix B). In order to excite magnon modes with the desired symmetry, we use a 2D antisymmetric sinc-pulse which should moreover avoid spurious effects in the spectrum [118] H pulse =H pulse sin(ω c t) ω c t pointing along theŷ-direction in order to excite backward volume waves [112]. The cut-off frequency was chosen to be ω c = 2π × 16 GHz and the cut-off wave vector to be k c = π/a in order to concentrate all the excitation energy in the first Brillouin zone. Since this pulse is centered in the middle of the crystal, we only excite modes around the crystal's center. The external saturation field was set to H ext = 400 mT (found by hysteresis) and the pulse field to H pulse = 0.4 mT. We note that the pulse strength should be a small perturbation of the saturation field in order to minimize non linear effects. We used the material parameters for YIG, M s = 140 kA/m (saturation magnetization), A ex = 2 pJ/m (exchange constant), K c1 = −610 J/m 3 (anisotropy constant) with the anisotropy axis alongẑ. In order to accelerate the simulations, we used an increased Gilbert damping parameter γ = 0.008 (compare to γ ≈ 10 −5 − 10 −4 for YIG) [119,120].
In the following considerations we focus only on the δm z component of the magnetization dynamics, since the optical spin density of a TE-like mode mostly points into theẑ-direction, rendering δm x and δm y irrelevant for G F (see Eq. (14)). We find that the optical defect also acts as confinement of the magnetic mode, resulting in the defect like dispersion relation presented in Fig. 4c. The obtained band structure shows modes around the edge of the BZ with extended wave vector character, implying that the modes are highly localized in space. The frequency of the highest excited localized mode at the BZ edge is ω mag = 2π × 13.12 GHz with an estimated linewidth (FWHM) of Γ mag = γ ω mag = 2π × 131.2 kHz where we used the real Gilbert damping of YIG γ = 10 −5 . Note that the simulated linewidth shown in Fig. 4c is larger due to the different choice of the Gilbert damping in order to speed up the simulation.
As we see from its mode shape, this mode is nicely localized at the holes attached to the defect and is odd with respect to x = 0 and y = 0, and hence has the same symmetry as the optical spin density as we aimed for (see Fig. 4c).

V. OPTOMAGNONIC CRYSTAL
As shown above, the crystal in Fig. 1 can host both optical and magnetic modes and therefore can be considered an optomagnonic crystal. In this section, we evaluate the optomagnonic coupling G F given in Eq. (14) (G C is briefly discussed at the end of the section) for the modes found in Secs. III and IV shown in Fig. 4.
Numerically evaluating Eq. (14) gives a Faraday contribution to the optomagnonic coupling per magnon per photon of |G F num | = 2π × 0.5 kHz. In order to gauge this value we want to compare it to the analytical estimate derived in [37]. In the optimal case, the magnetic mode volume and the optical mode volume coincide, V mag ≈ V opt . In this case, we estimate the coupling as which evaluates to |G F optimal | = 2π × 0.6 MHz using the material parameters of YIG ((θ F λ n )/(2π) = 4 · 10 −5 , M s = 140 kA/m) and the optical frequency found in Sec. III, ω opt = 2π × 249 THz. The magnetic mode volume is defined as the one where the magnon intensity is above a certain threshold, giving V mag = 2.8 · 10 −2 µm 3 (see appendix D). The coupling is bounded by the magnon mode volume, since in the investigated structure it is smaller than the optical mode volume (see Fig. 4). In order to take the mismatch in the mode volume into account, we introduce the following overlap measure which is also known as filling factor where V overlap represents the volume where the magnon and photon modes overlap. The volumes are estimated similar to the case of magnons to be V overlap = 9.7 · 10 −3 µm 3 and V opt = 0.7 µm 3 (see appendix D). Note that for the optical volume it was taken into account that the mode leaks out of YIG into the Si 3 N 4 layer and air, which is not shown in Fig. 4 a,b. Thus the overlap measure evaluates to O = 0.01, shrinking the optimal coupling to O · |G F optimal | ≈ 2π × 6 kHz. Hence, even though the optomagnonic crystal localizes both modes in the same region, the overlap measure is rather small due to the much larger optical mode volume (see Fig. 4 and Fig 7a), which is detrimental for the coupling strength. Furthermore by looking at the fine structure of the optical spin density and the magnon mode we see that the amplitude peaks of both do not coincide (see Fig. 6): the magnonic peaks are localized nearer to the center than the optical ones. This results in a smaller overlap volume which would be V mag if the peaks of the modes would be at the same position.
Since the coupling also strongly depends on the relative direction between the vectors of the modes, we additionally introduce a 'directionality' measure evaluating to D = 51% using the numerical results presented above. As we see, although the symmetries of the optical spin density and the magnon mode match, the vectors of the modes do not perfectly align in the defect area (see Figs. 4). Taking also this sub-optimal alignment into account the coupling estimate reduces to |G F expected | = O · D · |G F optimal | = 2π × 3 kHz which coincides well with the numerically obtained value. We conclude that the coupling in the investigated structure is mostly affected by the large difference between the optical and the magnetic mode volumes, shrinking the coupling value by two orders of magnitude. We remind the reader that the obtained values are for a proof of principle structure which has been only partially optimized, since we started from a fixed photonic crystal structure. in the next section we discuss a possible optimization from the magnonics side.
We now proceed to briefly discuss the Cotton-Moutton effect for the results found in Secs. III and IV. For YIG, the Cotton-Moutton coefficient (θ C λ n )/(2π) = −2 × 10 −5 [121] is of the same order of magnitude as the corresponding Faraday coefficient, determined by θ F . Since in the Voigt configuration both effects are of leading order in the magnetization fluctuations (see Eqs. (14) and (15)), it is important to take its contribution into account. Moreover, since the coefficients G F and G C are complex, it is difficult to estimate a priori the total coupling |G F + G C |, due to the unknown possible interference effects. Numerically evaluating Eq. (15) gives an interaction value of |G C num | = 2π × 1.6 kHz. This large value can be explained by the symmetry of the integrand which reduces to m x 0 [E x E * y δm y + E * x E y δm y ] due to the backward volume wave setup and the TE-like character of the optical mode. This integrand is fully even since E x E * y has the same symmetry as δm y . The full optomagnonic coupling (27) is found to be G num = 2π × 1.3 kHz.
Compared to the optomechanical coupling in similar 1D crystals, where coupling values (per photon and phonon) up to 2π × 950 kHz can be obtained [46][47][48][49][50][51], the optomagnonic coupling obtained here is still rather small. However, this is large compared to other optomagnonic systems. As we argued above, the coupling is limited by the imperfect spatial matching of magnons and photons with overlap O = 0.01 while it is enhanced due to small volumes, V mag ∼ 0.01 µm 3 and V opt ∼ 1 µm 3 . In the standard setups involving spheres [27][28][29], typically optical volumes are very large ∼ 10 5 µm 3 with low optomagnonic overlap ∼ 10 −3 , resulting in low couplings ∼ 1 Hz. It was theoretically shown that > 75% overlap in such systems is achievable [122] but the couplings would still be ∼ 2π × 500 Hz. The miniaturization of an optical cavity to ∼ 100 µm 3 was demonstrated in Ref. [38], where the coupling is however still small, 2π ×50 Hz, in this case due to the large magnon volume involved.
An important prerequisite for applications in the quantum regime such as magnon cooling, wavelength conversion, and coherent state transfer based on optomagnonics is a high cooperativity. The cooperativity per photon and magnon is an important figure of merit which compares the strength of the coupling to the lifetime of the coupled modes, and is given by where γ opt is the optical linewidth (FWHM), and Γ mag is the magnonic linewidth (FWHM).
To evaluate the theoretical cooperativity of the structure proposed in this manuscript, we use Γ mag = γω mag where γ = 10 −5 is the Gilbert constant and ω mag = 2π × 13.12 GHz. The optical linewidth is found from simulations to be γ opt = 2π × 0.2 THz.
Using the corresponding parameters the cooperativity per photon and magnon of the optomagnonic crystal is C crystal 0 ∼ 2.5 · 10 −10 . The single-particle cooperativity can be enhanced by the photon number in the cavity, C = n ph C 0 . Experimentally, there is a bound on the photon density that can be supported by the cavity without undesired effects due to heating, and it is empirically given by 5 · 10 4 photons per µm 3 [123]. In our structure, considering the effective mode volume V opt this gives an enhanced cooperativity at maximum photon density of C crystal ∼ 1·10 −5 , which is two orders of magnitude larger than the current experimental state of the art [38,123]. Since our model does not account for fabrication imperfections, this number is expected to be lower in a physical implementation, indicating that optimization is needed. Results for similar 1D optomechanical crystals indicate that optimization can lead to larger cooperativity values (at maximum photon density), e.g. ∼ 10 [47]. The small cooperativity obtained in our structure is a combination of a reduced coupling due to modemismatch, plus the very modest quality factor of the optical mode in this simple geometry.
For boosting the coupling strength we investigate briefly in the following the influence of the optomagnonic crystal's height on the coupling, as proposed in Ref. [104]. Therefore we increase the height of the YIG layer from 30 nm to 90 nm without changing the other parameters of the geometry (including the Si 3 N 4 layer in the optical simulations). As we see from the result (see Fig. 8) the coupling exhibits a V mag dependence. We find that the optical mode volume does not change substantially in the modified geometry, and therefore the observed behavior is consistent with the expected V mag /V opt dependence for a constant optical mode volume. The slight decrease for larger heights can be explained by the shrinking of the directionality measure D, stemming from the difference in symmetries obeyed by the magnetization (rotational) and the electric field (mirror).

VI. OPTIMIZATION
So far, we optimized the crystal in order to minimize optical losses for the given geometry. In this section we investigate how to optimize the geometry for magnonics. The optical optimization was achieved by fixing the hole radius and intra-hole distance, which are both along the length of the crystal. In the following we tune instead only the parameters along the width of the crystal (ŷdirection), in order to perturb as little as possible the optical optimization. We found a promising structure by increasing the width of the crystal and considering elliptical holes, see Fig. 9. From a set of trials, we found that a width of w = 900 nm and a radius of the holes along the width of r w = 380 nm give the highest coupling. An increased defect size of d = 1201.5 nm is also beneficial for decreasing the optical losses in this case, it nicely localizes the optical defect mode in the middle of the band gap and thus does not drastically change the localization behaviour of the photonic crystal.
For evaluating the photonic band structure and the optical modes we use the same procedure as described in Sec. III. We obtain a similar band structure for TElike modes and also a similar localized mode with a frequency of ω opt = 2π × 279 THz (obtained by Comsol, 235 THz according to MEEP) and a damping of κ opt = 2π × 3 THz which gives an optical linewidth (FWHM) of γ opt = 2π × 6 THz (see Fig. 10a). Using the values obtained by Comsol this results in a reduced optical quality factor of Q = 93 (note that MEEP gives a twice as large value). This rather low optical quality factor is a trade off for the magnetic optimization achieved by elliptical holes. Moreover, the optical spin density compared to the original crystal is mostly localized within the defect which is advantageous for our purposes (see Fig. 10b). Similarly, for evaluating the magnon modes we used the parameters and procedures presented in Sec. IV. In the following we focus on the Faraday part of the optomagnonic coupling and therefore consider only the δm z -component of the magnon mode due to the structure of the optical spin density. The Cotton-Mouton term is discussed briefly at the end of the section. The simulated band diagram for backward volume waves again shows extended magnon modes but in this case we obtain one broad band, most likely stemming from a fusion of several bands due to the larger width of the crystal (see Fig. 10c). The frequency of the highest excited localized mode is ω mag = 2π × 13.17 GHz with an estimated linewidth of γω mag = 2π × 131.7 kHz where we used the Gilbert damping of YIG γ = 10 −5 . As in the d a 2r x y z Figure 9. Optimization of the geometry: Through increasing the parameters along the width of the crystal we create more space for the modes without touching the optical optimization of the original crystal (dashed line). We note that we also increased the defect size, not shown here.  previous case, the simulated linewidth is larger due to the larger Gilbert damping used in the simulations. As we see from its mode shape, this mode is nicely localized at the holes attached to the defect and has approximately the same shape and symmetry as the optical spin density (see Fig. 10c).
Using the results discussed above, the Faraday component of the optomagnonic coupling of Eq. 14 for the optimized crystal evaluates to |G F num | = 2π × 2.9 kHz. Therefore the optimized coupling is one order of magnitude larger than in the crystal discussed in Sec. V. As before we want to gauge this value by comparing it to the analytical estimate given in Eq. 24. The optimal coupling in the optimized crystal is |G F optimal | = 2π × 0.5 MHz. Again the magnetic mode volume bounds the coupling due to the smaller size of the magnetic mode compared to the optical mode which also extends to the Si 3 N 4 layers. This results in a overlap measure (see Eq. 25) of O = 0.04. Therefore the mode overlap is increased by 25% compared to the un-optimized crystal. Evaluating the directionality measure given in Eq. 26 gives D = 53% which is just slightly larger than in the un-optimized case. Taking both measures into account the analytical coupling estimate shrinks to |G F expected | = O · D · |G F optimal | ≈ 2π × 10 kHz which lies slightly above the numerically obtained value. Although the fine structure peaks of the optical spin density and the magnon mode still do not coincide (see Fig. 11), the coupling values are improved by "pulling" the optical and magnetic modes completely into the defect area by the insertion of elliptical holes, creating an overlap area with high density of both modes.
The Cotton-Moutton effect in this structure evaluates to |G C num | = 2π × 1 kHz and results in a total coupling |G num | = |G F num +G C num | = 2π ×2 kHz. We can conclude that in this case both effects interfere constructively for the total coupling.
The cooperativity per photon and magnon in this case is C op 0 ∼ 2 · 10 −11 , which can be enhanced to C op ∼ 0.5 · 10 −6 by the number of photons trapped in the cavity. Thus the cooperativity at maximum photon density is slightly lower as in the crystal presented above, a consequence of the reduced quality factor of the optical mode.

VII. CONCLUSION
We proposed an optomagnonic crystal consisting of a one-dimensional array with an abrupt defect. We showed that this structure acts as a Bragg mirror both for photon and magnon modes, leading to co-localization of the modes at the defect. By proper design and taking into account the required symmetries of the modes in order to optimize the coupling, we showed that coupling values in the kHz range are possible in these structures. This value is orders of magnitude larger than the experimental state of the art in the field, but still rather small compared to the theoretically predicted optimal value for micron sized structures, which is in the range of ∼ 10 −1 MHz [32].
We showed that the strength of the coupling in our proposed structure is still limited largely by the suboptimal mode overlap, < 5%. Further optimization in design is moreover needed in order to boost the cooperativity value, which is limited mainly by the optical losses. The simultaneous optimization is challenging, due to the complexity of the demagnetization fields in patterned geometries. Whereas it is well known that tapered defect can highly increase the optical quality factors, its effect on the magnetic modes is non-trivial and it is outside of the scope of this work. More complex geometries, including two-dimensional crystals, remain to be explored. The first results shown in this work however point to the promise of designing the collective excitations in optomagnonic systems via geometry, in order to boost the coupling strength and minimize losses, paving the way for applications in the quantum regime. As the eigenmodes should diagonalize the Hamiltonian to ω γ |β γ | 2 , we should havê and For µ = γ, this gives the normalization for magnons. For circularly polarized magnons with δm µ = δm(y + iz)/ √ 2 and m 0 = x, the normalization becomesˆd Appendix B: Numerical settings -optical simulations In the following we shortly discuss how the optical band structure and the optical modes can be evaluated numerically. In our work we used two different computational methods: for calculations done in the time domain we use the electromagnetic simulation tool MEEP [103], whereas for calculations done in the frequency domain we use the finite element solver Comsol [106]. We use two simulation tools since with MEEP it is much easier to obtain the band structure of the crystal, and with Comsol the exact mode shape.
a. MEEP MEEP in general solves for Maxwell's equations in the time domain within some finite composite volume. Therefore it essentially performs a kind of numerical experiment [103]. We use MEEP for simulating the band structure of a YIG crystal without defect in order to find its band gap and its corresponding mid-gap frequency. Furthermore we use a transmission spectrum simulation to optimize the defect size in order to get the least lossy localized mode. For finding the exact frequency of the localized mode we also simulate its spatial shape in the time domain. For simplicity we simulate the YIG crystal in a 2D model (see [40,105]). Its material parameters are set by the relative permittivity of ε = 5. In order to account for the leakage of the electromagnetic field the investigated crystal is surronded by a finite size air region large enough so that the leaking electric field decays before it reaches the boundaries, in order to avoid spurious reflection effects. This is achieved by choosing the distance between the surface of the crystal and the boundary of the air region as d air = 3 · λ 0 where λ 0 is the vacuum wavelength ( d air = 33.6 µm in our case). Furthermore, we need boundary conditions along the outside of the air region that are transparent to the leaking field such that the truncated air region represents a reasonable approximation of free space. Therefore we use a perfectly matched layer at the boundaries of the air region which absorbs all outgoing waves. The thickness of this layer should be at least a vacuum wave length [124]. The whole geometry is meshed by one single resolution parameter which discretizes the structure in time and space and gives the number of pixels per distance unit. For all band simulations we used a resolution of 40 pixels, whereas in case of the transmission spectrum we used a resolution of 20 pixels and a resolution of 50 pixels in case of the mode shapes [105].

Band structure simulations
For obtaining the band structure we use a YIG crystal without defect (d = a) and therefore we can simulate only one single unit cell with a side length of a containing one air hole, and apply an infinite repetition of this cell at each side inx-direction. Since we expect the mid-gap frequency of the crystal with defect to be around 240 THz, we excite the crystal with a gaussian pulse with a center-frequency of 225 THz and a width of 450 THz to cover all modes around the band gap. We center the pulse peak at an arbitrary postion (x = 0.00123, y = 0) in order to couple the pulse to an arbitrary mode. Since we want to simulate only TE-like modes, in order to save computational time the pulse only has a H z component. For decreasing the computation time even more, we apply an odd mirror symmetry plane for y = 0. The mirror symmetry for x = 0 is broken by a boundary condition for 0 < k x < π [105].

T ransmission spectrum simulations
For optimizing the defect size we simulate a transmission spectrum for frequencies at the band gap by measuring the flux at the end of the waveguide stemming from a source at the other end. The measured flux then is normalized to the flux of a waveguide without holes. We therefore simulate the transmission spectrum as a function of different defect sizes and use the defect size which gives the highest transmission. In order to consider only TE-like modes where the electric field lies in plane we need to excite the system with a J y -current source transverse to the propagation direction which is achieved by a gaussian pulse with only a E y -component. Its center frequency thereby is 222 THz (simulated mid-gap frequency) and its width is 90 THz (> band width). Also in this case we apply an odd mirror symmetry for y = 0 for decreasing the simulation time. We note that the mirror symmetry for x = 0 is broken by the source since it is located at the edge of the waveguide [105].

M ode shape simulations
For evaluating the mode frequency of the localized mode within the band gap we simulate the time evolution of this single mode by exciting it by a gaussian pulse with a center frequency of 203 THz (frequency of the peak in the transmission spectrum) and a width of 15 THz. Since in this simulation no symmetry is broken we also apply an odd mirror symmetry for x = 0 and y = 0 for obtaining only a TE-like mode [105].

b. Comsol
We use Comsol to find the spatial mode shape. Therefore we use the "Electromagnetic waves, Frequency domain" package of COMSOL's "RF module" which solves for the Helmholtz equation of the form where k 0 indicates the vacuum wave number, ω the angular frequency, µ r the relative permeability and ε 0 the vacuum permittivity. Contrary to the MEEP simulations above, we simulate the full geometry composite of a YIG layer sandwiched by two Si 3 N 4 layers. The used material parameters thereby are ε yig = 5, ε si = 4, µ yig = µ si = 0, and σ yig = σ si = 0 with µ the relative permeabilty and σ the conductivity. Again we also need to simulate an truncated air region around the crystal which is able to absorb the outgoing radiation. The corresponding material parameters are ε air = µ air = 1 and σ air = 0. Besides perfectly matched layers we also can use second order scattering boundary conditions at the air surfaces given by the expression [124] with n the normal vector to the considered plane. For large enough air regions both approaches are almost equivalent as long as the leaking field is propagating normal to the air surfaces. In order to account for a large enough air region we choose the distance between the surfaces of the crystal and the air boundaries as 4.5 µm. For reducing the simulation time we use the symmetry requirements of a TE-like mode. Therefore, we cut the geometry into an eighth of the whole structure and apply perfect electric conductor boundary conditions (n × E = 0) at the cut surfaces along x = 0 and y = 0 and a perfect magnetic conductor boundary condition (n × H = 0) at the cut surface along z = 0. The full solution is then obtained by using the symmetry requirements of a TE-like mode. The whole geometry is meshed by a physics-controlled tetrahedral mesh with a maximum element size of λ 0 /5 ≈ 0.3 µm [125]. We note that in case of a physics-controlled mesh Comsol automatically meshes the material areas of interest with a finer mesh and uses a coarser mesh e.g. for the air regions.
Appendix C: Numerical settings -magnetic simulations in this section we briefly discuss how the magnetic band structure and magnetic mode shape is obtained numerically. For evaluating the magnetization dynamics we use the finite difference tool MuMax3 [117] which solves for the Landau-Lifshitz-Gilbert equation of the form with m = M/M s the local reduced magnetization of one simulation cell, g the gyromagnetic ratio, α the damping parameter, and B eff an effective field which contributions can be found in [117]. As material parameters we used the parameters for YIG, M s = 140 kA/m (saturation magnetization), A ex = 2 pJ/m (exchange constant), K c1 = −610 J/m 3 (anisotropy constant) with the anisotropy axis alongẑ. In order to accelerate the simulations, we used an increased Gilbert Damping parameter α = 0.008 (compare to α ≈ 10 − 5 for YIG) [126]. The used meshgrid had (1024, 50, 5) cells in the (x,ŷ,ẑ)-direction what guarantees to take the exchange interaction into account (l ex ≈ 13 nm for YIG).
In general, in all our simulations the spin wave dynamics is excited via an external pulse field and the time evolution is recorded for all three magnetization components. For post processing the output of the form m i (x, y, z) with i = (x, y, z) saved for all simulated time steps separately we create for each magnetization component i = (x, y, z) a 4D-array of the form m i (t, x, y, z).

Band structure simulations
In order to obtain the band structure along a specific direction j with j = (x, y, z) e.g. chosen to be thex-direction, we reduce the four dimensional array to a two dimensional array of the form m i (t, x) = ny m nz n δm i (t, x, y m , z n ) and perform a 2D Fourier transform on this array δm i (f, k x ) = FT 2D [m i (t, x)] resulting in the band diagram along the chosen direction. For increasing the resolution in the band diagram we plot the quantity |δm i (t, x)|/max(|δm i (t, x)|).

M ode shape simulations
In order to obtain the mode shape we perform a space-dependent Fourier transform in time on each array entry separately δm i (f, x, y, z) = FT 1D [m i (t, x, y, z)].