Magnonic crystals — prospective structures for shaping spin waves in nanoscale

We have investigated theoretically band structure of spin waves in magnonic crystals with periodicity in one(1D), two- (2D) and three-dimensions (3D). We have solved Landau–Lifshitz equation with the use of plane wave method, finite element method in frequency domain and micromagnetic simulations...

Full description

Saved in:
Bibliographic Details
Published in:Физика низких температур
Date:2015
Main Authors: Rychły, J., Gruszecki, P., Mruczkiewicz, M., Kłos, J.W., Mamica, S., Krawczyk, M.
Format: Article
Language:English
Published: Фізико-технічний інститут низьких температур ім. Б.І. Вєркіна НАН України 2015
Subjects:
Online Access:https://nasplib.isofts.kiev.ua/handle/123456789/128078
Tags: Add Tag
No Tags, Be the first to tag this record!
Journal Title:Digital Library of Periodicals of National Academy of Sciences of Ukraine
Cite this:Magnonic crystals — prospective structures for shaping spin waves in nanoscale / J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, M. Krawczyk // Физика низких температур. — 2015. — Т. 41, № 10. — С. 959–975. — Бібліогр.: 65 назв. — англ.

Institution

Digital Library of Periodicals of National Academy of Sciences of Ukraine
id nasplib_isofts_kiev_ua-123456789-128078
record_format dspace
spelling Rychły, J.
Gruszecki, P.
Mruczkiewicz, M.
Kłos, J.W.
Mamica, S.
Krawczyk, M.
2018-01-05T17:39:42Z
2018-01-05T17:39:42Z
2015
Magnonic crystals — prospective structures for shaping spin waves in nanoscale / J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, M. Krawczyk // Физика низких температур. — 2015. — Т. 41, № 10. — С. 959–975. — Бібліогр.: 65 назв. — англ.
0132-6414
PACS: 75.30.Ds, 75.70.Cn, 75.75.–c
https://nasplib.isofts.kiev.ua/handle/123456789/128078
We have investigated theoretically band structure of spin waves in magnonic crystals with periodicity in one(1D), two- (2D) and three-dimensions (3D). We have solved Landau–Lifshitz equation with the use of plane wave method, finite element method in frequency domain and micromagnetic simulations in time domain to find the dynamics of spin waves and spectrum of their eigenmodes. The spin wave spectra were calculated in linear approximation. In this paper we show usefulness of these methods in calculations of various types of spin waves. We demonstrate the surface character of the Damon–Eshbach spin wave in 1D magnonic crystals and change of its surface localization with the band number and wavenumber in the first Brillouin zone. The surface property of the spin wave excitation is further exploited by covering plate of the magnonic crystal with conductor. The band structure in 2D magnonic crystals is complex due to additional spatial inhomogeneity introduced by the demagnetizing field. This modifies spin wave dispersion, makes the band structure of magnonic crystals strongly dependent on shape of the inclusions and type of the lattice. The inhomogeneity of the internal magnetic field becomes unimportant for magnonic crystals with small lattice constant, where exchange interactions dominate. For 3D magnonic crystals, characterized by small lattice constant, wide magnonic band gap is found. We show that the spatial distribution of different materials in magnonic crystals can be explored for tailored effective damping of spin waves
The research leading to these results has received funding from Polish National Science Centre project DEC-2- 12/07/E/ST3/00538 and from the EUs Horizon2020 research and innovation programme under the Marie Sklodowska-Curie GA No644348. The numerical calculation were performed at Poznan Supercomputing and Networking Center (grant No. 209).
en
Фізико-технічний інститут низьких температур ім. Б.І. Вєркіна НАН України
Физика низких температур
Специальный выпуск К 80-летию уравнения Ландау–Лифшица
Magnonic crystals — prospective structures for shaping spin waves in nanoscale
Article
published earlier
institution Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection DSpace DC
title Magnonic crystals — prospective structures for shaping spin waves in nanoscale
spellingShingle Magnonic crystals — prospective structures for shaping spin waves in nanoscale
Rychły, J.
Gruszecki, P.
Mruczkiewicz, M.
Kłos, J.W.
Mamica, S.
Krawczyk, M.
Специальный выпуск К 80-летию уравнения Ландау–Лифшица
title_short Magnonic crystals — prospective structures for shaping spin waves in nanoscale
title_full Magnonic crystals — prospective structures for shaping spin waves in nanoscale
title_fullStr Magnonic crystals — prospective structures for shaping spin waves in nanoscale
title_full_unstemmed Magnonic crystals — prospective structures for shaping spin waves in nanoscale
title_sort magnonic crystals — prospective structures for shaping spin waves in nanoscale
author Rychły, J.
Gruszecki, P.
Mruczkiewicz, M.
Kłos, J.W.
Mamica, S.
Krawczyk, M.
author_facet Rychły, J.
Gruszecki, P.
Mruczkiewicz, M.
Kłos, J.W.
Mamica, S.
Krawczyk, M.
topic Специальный выпуск К 80-летию уравнения Ландау–Лифшица
topic_facet Специальный выпуск К 80-летию уравнения Ландау–Лифшица
publishDate 2015
language English
container_title Физика низких температур
publisher Фізико-технічний інститут низьких температур ім. Б.І. Вєркіна НАН України
format Article
description We have investigated theoretically band structure of spin waves in magnonic crystals with periodicity in one(1D), two- (2D) and three-dimensions (3D). We have solved Landau–Lifshitz equation with the use of plane wave method, finite element method in frequency domain and micromagnetic simulations in time domain to find the dynamics of spin waves and spectrum of their eigenmodes. The spin wave spectra were calculated in linear approximation. In this paper we show usefulness of these methods in calculations of various types of spin waves. We demonstrate the surface character of the Damon–Eshbach spin wave in 1D magnonic crystals and change of its surface localization with the band number and wavenumber in the first Brillouin zone. The surface property of the spin wave excitation is further exploited by covering plate of the magnonic crystal with conductor. The band structure in 2D magnonic crystals is complex due to additional spatial inhomogeneity introduced by the demagnetizing field. This modifies spin wave dispersion, makes the band structure of magnonic crystals strongly dependent on shape of the inclusions and type of the lattice. The inhomogeneity of the internal magnetic field becomes unimportant for magnonic crystals with small lattice constant, where exchange interactions dominate. For 3D magnonic crystals, characterized by small lattice constant, wide magnonic band gap is found. We show that the spatial distribution of different materials in magnonic crystals can be explored for tailored effective damping of spin waves
issn 0132-6414
url https://nasplib.isofts.kiev.ua/handle/123456789/128078
citation_txt Magnonic crystals — prospective structures for shaping spin waves in nanoscale / J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, M. Krawczyk // Физика низких температур. — 2015. — Т. 41, № 10. — С. 959–975. — Бібліогр.: 65 назв. — англ.
work_keys_str_mv AT rychłyj magnoniccrystalsprospectivestructuresforshapingspinwavesinnanoscale
AT gruszeckip magnoniccrystalsprospectivestructuresforshapingspinwavesinnanoscale
AT mruczkiewiczm magnoniccrystalsprospectivestructuresforshapingspinwavesinnanoscale
AT kłosjw magnoniccrystalsprospectivestructuresforshapingspinwavesinnanoscale
AT mamicas magnoniccrystalsprospectivestructuresforshapingspinwavesinnanoscale
AT krawczykm magnoniccrystalsprospectivestructuresforshapingspinwavesinnanoscale
first_indexed 2025-11-26T13:15:49Z
last_indexed 2025-11-26T13:15:49Z
_version_ 1850622333728849920
fulltext Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10, pp. 959–975 Magnonic crystals — prospective structures for shaping spin waves in nanoscale J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk Faculty of Physics, Adam Mickiewicz University in Poznan, 85 Umultowska, Poznań 61-614, Poland E-mail: krawczyk@amu.edu.pl Received April 1, 2015, published online August 25, 2015 We have investigated theoretically band structure of spin waves in magnonic crystals with periodicity in one- (1D), two- (2D) and three-dimensions (3D). We have solved Landau–Lifshitz equation with the use of plane wave method, finite element method in frequency domain and micromagnetic simulations in time domain to find the dynamics of spin waves and spectrum of their eigenmodes. The spin wave spectra were calculated in linear approximation. In this paper we show usefulness of these methods in calculations of various types of spin waves. We demonstrate the surface character of the Damon–Eshbach spin wave in 1D magnonic crystals and change of its surface localization with the band number and wavenumber in the first Brillouin zone. The surface property of the spin wave excitation is further exploited by covering plate of the magnonic crystal with conductor. The band structure in 2D magnonic crystals is complex due to additional spatial inhomogeneity introduced by the demagnetizing field. This modifies spin wave dispersion, makes the band structure of magnonic crystals strongly dependent on shape of the inclusions and type of the lattice. The inhomogeneity of the internal magnetic field becomes unimportant for magnonic crystals with small lattice constant, where exchange interactions dominate. For 3D magnonic crystals, characterized by small lattice constant, wide magnonic band gap is found. We show that the spatial distribution of different materials in magnonic crystals can be explored for tailored effective damping of spin waves. PACS: 75.30.Ds Spin waves; 75.70.Cn Magnetic properties of interfaces (multilayers, superlattices, heterostructures); 75.75.–c Magnetic properties of nanostructures. Keywords: magnonic crystals, spin waves, eigenmodes. Introduction Material patterning is one of the major factors used for controlling propagation of waves with short wavelengths in solids. The use of periodic patterning for controlling the wave dynamics has increased significantly since the dis- covery of photonic crystals in 1987 [1,2]. Since that time periodicity has been extensively used for molding the flow of electromagnetic waves in the range from microwaves to optical wavelengths [3], which led to the discovery of new materials with properties unheard of in nature [4]. Ideas developed in photonics were transferred to other types of excitations, such as phonons [5], plasmons [6] and also spin waves (SWs) [7]. Periodicity was used for controlling SWs excitations in ferromagnetic materials already in 1970s [8,9]. The trans- mission of magnetostatic waves was tailored by periodic distribution of the saturation magnetization, realized by ion implantation [10], a regular lattice of etched grooves in a magnetic dielectric, periodic modulation by metallic stripes or dots on top of a ferromagnetic film [11], or peri- odic perturbation of the magnetic field [12]. Due to li- mitations in fabrication and technology the investigations were limited to large structures, in which the dipolar inter- action prevailed over the exchange interaction. The dis- covery of photonic crystals renewed the interest in magnet- ic periodic structures, inspiring many new ideas, providing abundant new physics, and pushing magnetization dynamics studies in unexplored directions. The concept of magnonic crystals (MCs) was proposed as a SW counterpart of pho- tonic crystals [13–16]. MCs are magnetic structures with periodic distribution of the constituent materials or period- ic modulation of some magnetic parameters (e.g. saturation magnetization, exchange interactions or magnetocrystalline © J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk, 2015 J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk anisotropy), or other parameters relevant to the propaga- tion of SWs, such as external magnetic field, film thick- ness, stress or a surrounding of the ferromagnetic film. In any periodic medium the eigensolutions of the wave equation in the linear regime fulfill Bloch’s theorem and, regardless of the type of excitation, form a band structure in the frequency–wave vector space. The same applies to the magnonic (i.e., spin-wave) band structure. Neverthe- less, the details of the band structure can be derived mostly from numerical calculations. Usually the magnonic band structure is studied in terms of its sensitivity to: structural changes of the MC, modifications of material parameters or application of external fields. SWs have a complex dispersion relation even in homo- geneous thin film [16] which is substantially different from the dispersion of SWs in a bulk material. The dispersion depends on the direction and magnitude of the wave vec- tor, the relative strength of short-range exchange interac- tion with respect to the strength of long-range dipolar in- teraction, the external shape of the sample, the magnitude of the applied static magnetic field and its orientation in relation to the direction of SWs propagation. Also the mag- netocrystalline anisotropy and magnetoelastic effects con- tribute to SWs dynamics. This means that the SW band structure of MCs will be influenced by many additional factors, both intrinsic and external ones, apart from those which are typical for the other types of artificial crystals. This makes MCs an intriguing system for scientific studies and the main subject of research in the field of magnonics [7]. This paper is dedicated to present various spin wave band structures which can be realized in MCs to demonstrate inter- esting properties of SWs resulting from the periodicity. We present the results of calculations for one-dimensional (1D), two-dimensional (2D) and tree-dimensional (3D) MCs ob- tained with different numerical methods. All spectra are ob- tained by solving Landau–Lifshitz (LL) equation. The plane wave method, finite element method in the frequency do- main and micromagnetic simulations based on finite dif- ference method in the time domain are demonstrated to be complementary for calculations of the SW spectra in MCs. Model and wave equation The dynamics of SWs can be analyzed using two main theoretical approaches [18]. One uses the discrete lattice model based on Heisenberg Hamiltonian with atomic structure of the ferromagnetic material directly taken into account. The other is based on solutions of the LL equation defined in continuous medium which describes precession of classical magnetic moments in effective magnetic field. The later approach is more suitable for systems with com- plex geometry in nano- and larger scales, where local pa- rameters (exchange integral, spin) can be expressed by macroscopic parameters (exchange length and magnetiza- tion saturation) [19]. The typical spatial resolution of nano- litography techniques is order of magnitude larger than interatomic distances. This justify use of an approach based on LL equation to these structures, which are widely fabri- cated nowadays [20]. We solve the LL equation, i.e., the equation of motion for the magnetization vector ( , ):tM r 0 eff ( , ) ( , ) ( , )t t t t ∂ = −γµ × + ∂ M r M r H r [ ]0 eff( , ) ( , ) ( , ) , S t t t M αγµ + × ×M r M r H r (1) where γ is gyromagnetic ratio, eff ( , )tH r denotes effec- tive magnetic field, 0µ is permeability of vacuum, r is position vector, t is time and SM is saturation magnetiza- tion. Relaxation processes are described by the last term on the right-hand side of the Eq. (1). We assume that MC is saturated, i.e., a collinear static magnetization in all the investigated structures is assumed. In this paper we use a coordinate system with the z axis being the direction of the external magnetic field and also the direction of the static magnetization in the saturated state. In the case of small disturbance of the magnetization from its equilibrium orientation, linear SWs are generated and the calculations can be conducted in linear approxima- tion. In this case the component of the magnetization vec- tor along equilibrium direction (z axis) is constant in time ˆ[ ( , ) ( ) ( , )]zt M z t= +M r r m r and can be approximated by spontaneous magnetization ( ) ( )z SM M≈r r . This assump- tion requires much larger magnitude of zM than those of the perpendicular components: ( , ) ( ),St Mm r r where ( , )tm r is a two-dimensional dynamic vector lying in the (x, y) plane: ˆ ˆ[ ( , ) ( , ) ( , ) ].x yt m t x m t y= +m r r r The effective magnetic field is assumed to be the sum of three terms: eff 0 ms ex ,= + +H H H H where 0H is the external static magnetic field; msH is the magnetostatic field with two components: static demagnetizing field dem ( )H r and dynamic components ( , )th r that are perpen- dicular to 0:H ms dem[ , , ];x yh h H=H exH is the exchange field which can be formally defined in saturation using linear space dependent exchange operator ex ˆ :H ex ex ˆ( ) ( ) ( ).=H r H r m r We neglect the contribution of the magnetic anisotropy. In magnetostatic approximation the magnetostatic field can be expressed as a gradient of the scalar magnetostatic potential: ms ,= −∇ϕH thus the curl of magnetostatic field (the one of the Maxwell equations) is always equal to zero: ms 0.∇× =H In order to find the dynamic components of the demagnetizing field ( , ),x yh h we need to solve the Gauss equation: ms[ ( ) ( )] 0.∇⋅ + =H r M r Putting to the Gauss equation magnetostatic field written as the gradient of ϕ we receive equation: 2 ( )( ) ( ) ( ) . yx Smm M x y z ∂∂ ∂ ∇ ϕ = + + ∂ ∂ ∂ rr r r (2) 960 Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 Magnonic crystals — prospective structures for shaping spin waves in nanoscale This equation defines the relation between magnetostat- ic potential (and magnetostatic field) and the magnetiza- tion. This shall be solved together with Eq. (1). In Eq. (1) with linear approximation we can separate space and time variables. Then, the time dependent equation has solutions in the form of monochromatic spin waves: ( , ) ( , ) e i t x yt m m − ω= ∝m r (similar for dynamic compo- nents of the magnetostatic field ( , ) ( , ) e i t x yt h h − ω= ∝h r ), where ω is the angular frequency of SW. Thus, the mono- chromatic SW is a solution of the following equations: 0 dem( ) ( ( ) )Ω x y yi m H m H m= + −r r r S ex( ) ( ) ( ) ( ) ( ), y S yM h M m− −r r r H r r 0 dem( ) ( ) (Ω  )y x xi m H m H m− = + −r r r ex( ) ( ) ( ) ,( ) ( )S x S xM h M m− −r r r H r r (3) where 0 Ω ω = γµ is reduced angular frequency of the SW. Planar magnonic crystals The simplest MC is a 1D periodic structure which has a form of magnetic multilayers [21]. These systems were intensively studied in the past and their fundamental fea- tures are well-understood using both discrete and continu- ous models [22]. More complex spin wave dynamics can be observed in planar structures with in-plane period- icity [23,24]. The planar geometry is the most common for the nano- and microstructures fabricated using top-down techniques (e.g. nanolithography) [20]. The obtained struc- tures can possess long-range order of high precision, which is crucial to explore the effects of the periodicity in SW dynamics. On the other hand, the bottom-up methods [25], are effective by utilizing processes of self-organization, however they often fail in fabrication of ideal 2D and 3D periodic systems because of incontrollable defects and dis- location. Because of the mentioned above reasons, both high quality 1D and 2D MCs are fabricated mainly in the form of planar systems with in-plane dimensions much larger than thickness of the structure [26]. MCs, as any other kind of periodic medium, has aniso- tropic dynamics (according with symmetry of the structure) but the geometry is not the only source of the anisotropy in SW dynamics. An inhomogeneity of the static magnetiza- tion distribution induces demagnetization and exchange fields, the additional two factors which influence the SW dynamics. In MCs these magnetic fields have also periodic distribution, however they introduce additional inhomoge- neity and anisotropy which can reduce symmetry of the lattice. The other factor which is inevitable for magnonic systems is an external magnetic field which can be arbitrar- ily oriented with respect to the crystallographic lattice axes or MC plane. The magnetic ground state (spatial distribu- tion of the static magnetization) can be very complex in MCs at low external magnetic field, below the saturation field. The investigation of SW dynamics in such system is complicated, because of complexity of the ground state and presence of remagnetization processes. To design MCs with desired SW dispersion, the stable magnetic configura- tion is usually required. Strong enough magnetic field can be always used to enter system in saturation state. Under this condition, the external magnetic field determines the direction of the static magnetization, which is almost ho- mogeneous and constant in different pieces of the same material. The magnetization dynamics can be then mostly attributed to SWs. Plane wave method In this section we present an adaptation of the plane wave method (PWM) to the calculations of the SW disper- sion relation in planar MCs with 2D periodicity. The gen- eral algorithm of this method is presented schematically in Fig. 1. If surface magnetic anisotropy on the surfaces of the MC is neglected and the thickness of the planar MC is small then the magnetocrystalline and dipolar pinning on its top and bottom surfaces is week. Therefore, we can as- sume that the magnetization is free on film’s surfaces. Moreover, because the ratio of the period to thickness for the MC is small, the lowest magnonic bands will represent the modes which are not quantized across thickness of the MC. These two assumptions allow to consider the SW amplitude as constant across the slab thickness in thin planar MCs. Using this approach we will proceed with PWM, Fig. 1. (Color online) Algorithm of the PWM to calculate the dispersion relation of SWs and their amplitude distribution in magnonic crystal. It is split into two parts, analytical derivation of the algebraic eigenproblem and numerical solutions. Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 961 J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk which applies 2D Fourier transform to linearized LL Eq. (3) supplemented with Gauss Eq. (2) and converts this equa- tions into the algebraic eigenproblem. This consist first part of the algorithm presented in Fig. 1. To illustrate PWM, we will discuss the solution for the system presented in Fig. 2. This is a square lattice (with lattice constant a) of the square dots (of side size s). Dots (ferromagnetic material A) and matrix (ferromagnetic ma- terial B) are saturated by the in-plane external magnetic field. The solutions of Eq. (3), which is differential equa- tion with coefficients being periodic function of the posi- tion vector, can be written in the form of the Bloch waves: ( ) ( )e ,i ⋅= k r km r m r (4) where (0, , )x yk k=k is the wavevector and ( )km r is an SWs complex amplitude with the same periodicity as the lattice of the MC, ( ) ( )+ =k km r a m r  , (0, , )y za a=a and ya and za are components of a two dimesional lattice vec- tor. For a given wave vector k, we expand the periodic factor of the Bloch function ( )km r and material parame- ters (i.e., saturation magnetization and exchange constant) in the basis of the plane waves: ei ⋅G r ( ) ( ) ( )e and ( )ei i S SM M⋅ ⋅= =∑ ∑G r G r k k G G r m G r Gm (5) where G is a reciprocal lattice vector, which for the square lattice takes the values: 2(0, , ) (0, , ),y z y zG G n n a π = =G with ,y zn n being integers. Also, in Eq. (2) the magnetostatic potential can be rep- resented as a Bloch function, and its periodic part can be expanded in Fourier series, together with SM or Bloch functions [ ( ), ( )]y xm mr r included in this equation. The so- lution of Eq. (2) for static demagnetizing field in the pla- nar structure with piecewise constant magnetization of ar- bitrary orientation was derived by Kaczer et al. in Ref. 27. These formulas were extended to dynamical components of the magnetization in the form of the Bloch function (4) in Ref. 28. Finally, dm, ( )zH r and ( ), ( )x yh hr r can be ex- pressed as a function of the Fourier coefficients ( ), SM G ( ) xm G and ( )ym G in the following form: ( ) ( ) 2 dm 2 ( ) 1 , e ,iS zM G H S x ⋅= − −  ∑ G r G G r G G (6) ( ) ( ) ( ), ,y y x y k G h im S x + = + − + ∑ k G r G k G k G ( ) ( ) ( ) , , e ,i xm C x + ⋅  − +   k G r k G k G (7) ( ) ( ) ( ), , y y y x k G h im S x  + = + −  + ∑ k G r G k G k G ( ) 2 , 2 ( )( ) [1 ( , )] e ,y y y im k G S x + ⋅ + − − + +  k k G rG k G k G (8) where 2( , ) sinh ( )e , d S x x − = κ κ κ 2( , ) cosh ( )e d C x x − = κ κ κ and d is the MC’s thickness. In PWM we use exchange field in the form 2 ex ex2 0 2 ,( ) ( ) ( ) S A l M     = ∇ ⋅ ∇ =∇ ⋅ ∇   µ   H m r r m r where ( )ex 2 0 2 S Al M = µ r is the exchange length. The ex- pression in the square brackets is the exchange operator. According to Bloch theorem (4) and using the expan- sions (5)–(8), the Eq. (3) can be transformed into the alge- braic eigenproblem, which finalize the first part of the al- gorithm shown in Fig. 1: Fig. 2. (Color online) (a) The schema of the structure of bi-component planar magnonic crystal. Square lattice (lattice constant, a) of square shaped dots (with size s) immersed in ferromagnetic material B. The MC is saturated by external magnetic field along the z-axis. (b) Reciprocal lattice of the square lattice with the first Brillouin zone marked with red dotted line. The irreducible part of the Brillouin zone is marked with bold solid line. 962 Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 Magnonic crystals — prospective structures for shaping spin waves in nanoscale ˆ .i= Ωk kMm m (9) The eigenvector , ,[ , ]=k k km m mx y consists of the two sub-vectors being the finite sets (of N elements) of Fourier coefficients for periodic factor of the Bloch functions: , , 1 , 2 ,[ ( ), ( ), , ( )]x x x x Nm m m=k k k km G G G and , , 1 , 2 ,[ ( ), ( ), , ( )].y y y y Nm m m=k k k km G G G The matrix M is a block matrix ˆ . xx xy yx yy M M M M    =     M (10) The blocks of the matrix (10) are defined as: _____________________________________________________ ( ) ( ), , , + y y j yyxx ij S i j j ij j k G M i M S x M + = − − + = −G G k G k G ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 ,2 0 ex 2 1 , y y jxy ij j l l j S i l S i j jij l j k G M H l M M S x +  = δ + + ⋅ + − − + − − + −  + ∑ k G k G G G G G G G k G k G ( ) ( ) ( ) 2 , , 2 1 , , z i z j S i j i j i j G G M S x +  − − − +  − G G G G G G ( ) ( ) ( ) ( ) ( ) ( )2 0 ex ,yx ij j l l j S i l S i j jij l M H l M M C x= − δ − + ⋅ + − − − − + +∑ k G G G G G G G G k Gk ( ) ( ) ( ) 2 , , 2 1 , . z i z j S i j i j i j G G M S x +  + − − +  − G G G G G G ________________________________________________ i and j index reciprocal lattice vectors used in the Fourier expansions. For the numerical solution of the Eq. (9), we need to limit a number of the reciprocal lattice vectors in all expansions and also calculate the coefficients of the Fourier expansion (5) of the magnetization saturation and exchange length, ( )SM G and 2 ex ( ),l G respectively. The ge- neral formula for the coefficients of Fourier expansion reads: ( ) ( ) 1 e ,i S F F dS S − ⋅= ∫ G rG r where ( )F r and ( )F G are periodic functions in real space describing the spatial distribution of material parameter (this is ( )SM r or 2 ex ( )l r in Eq. (3)) and its Fourier expan- sion coefficient for reciprocal lattice vector G (this is ( )SM G or 2 ex ( )l G in eigenvalue problem Eq. (9)), respec- tively. The S denotes the area of the unit cell. The Fourier coefficients for MC with inclusions of the regular shapes can be analytically calculated. For instance for 2D MC with inclusions of circular and square shapes, the coeffi- cients cir ( )F G and sq ( )F G take the following form: ( ) ( ) ( ) ( ) 2 2cir 1 for 0, 1 2 for 0, B A B A B F F F R F RS F F J RG G  + − π = =  π − ≠  G G G ( ) ( ) ( ) ( ) ( ) 2 2 sq 2 2 for 0, sinc for 0, 0, 2 1 sinc for 0, 0, 2 sinc sinc for , 2 2 0 B A B x A B x y y A B x y yx A B F F F s G s F F s G G G sF S F F s G G G sG s F F s  + − =    − ≠ =    =   − = =          − ≠         G G G where R is the radius of circular inclusions, G stands for the length of the reciprocal lattice vector G and 1J denotes the Bessel function of the first kind. AF and BF stands for the respective values of saturation magnetization (or ex- change length) in the inclusion material and host material, respectively. For the square lattice 2S a= (Fig. 2(a)). The algebraic eigenproblem (9) is solved numerically using standard routines to find eigenvalues and eigenvec- tors, and this constitute the second part of the algorithm. The solutions are eigenvalues Ω, from which frequency f of the SW modes for given k-vector are obtained, and ei- genvectors ,km which are used to obtain distribution of the SW amplitude in real space, ( ).km r The dispersion relation is found by solving eigenproblem (9) repetitively for successive values of the k along the high symmetry path in the first Brillouin zone (BZ). For square lattice this path is marked with red dotted line in Fig. 2(b). Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 963 J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk Magnonic band structure in thin 2D MCs The spectrum of MCs can be tailored by adjustment of structural and material parameters. The following parame- ters can be considered for system presented in Fig. 2: mate- rial composition of the matrix and inclusions, lattice con- stant a, size s and shape of the inclusions, and also the thickness of the slab d. It is worth to notice that due to in- terplay between dipolar and exchange interactions, the spectrum of magnonic crystals do not scale with the size of the system, as it is in photonic systems [3]. To demonstrate basic changes in the magnonic band structure due to struc- tural parameters, we investigate the variation of the abso- lute values of structural parameters but only their relative changes, even if we are interested in exploring qualitative features of the magnonic spectrum. Let’s consider MCs in two different limits of sizes [29]: exchange dominated re- gime — for small lattice constant and dipole dominated regime — for larger patterns. In Fig. 3 we present results of PWM calculations for MC with small lattice constant (a = 50 nm) for various shapes of inclusions: square, hex- agonal and circular, with fixed filling fraction 0 55,( .ff = the filling fraction is defined as a ratio of the area occupied by inclusion to the area of the unit cell, for structure from Fig. 2 this is 2 2 / ).ff s a= The first magnonic band does not change significantly with the changes of the shape of in- clusions. Also the corresponding profiles of dynamical magnetizations, which don’t reproduce all details of the inclusions, are similar. In this range of sizes, the exchange interactions, which are isotropic in terms of external field direction, dominate. This property makes the MC a close counterpart of the photonic crystal and allows to deduce a lot of its features in analogy to photonics. However, the impact of dipole interactions is still noticeable in these MCs. We can observe differences in the dispersions for the equivalent crystallographic directions — cf, the dispersion of the first band along Γ–X and Γ–X′, parallel and perpen- dicular to the external magnetic field. Figure 4 presents the dispersions for planar MCs of larger sizes (lattice constant a = 400 nm), with important influence of the dipole interactions. We observe the dra- matic change of the magnonic spectra for systems differing in the shape of inclusions (note that we kept the filling fraction constant for the three considered shapes). The one of the most important factors responsible for such differ- ences in the SW spectra is the static demagnetizing filed (see the right insets in Fig. 4). This field has a form of sharp wells/peeks located at the interfaces of materials differing in magnetization saturation. The amplitude of demagnetizing field depends both on the magnetization contrast (on both sides of interface) and the orientation of the interface with respect to the direction of the external magnetic field. The later feature is shape-dependent. The strongest impact of the demagnetizing field is noticed for the system with square inclusions where long sides of squares are oriented perpendicularly to the direction of H0. This induces deep and long wells of demagnetizing fields which can capture the modes and localize them. These modes (in Fig. 4(a) modes 1 and 2) are called edge modes and have frequency below the frequency of the fundamen- tal mode. We discussed the PWM and its adaptation to calcula- tion of the SW dispersion relation in the planar MC. One of the assumption we have made was about homogeneity of the sample across its thickness. However, there are pla- nar MCs which are nonuniform in out-of-plane direction, e.g. magnetic slab with an array of grooves or with mag- netic dots on its surface. In these cases, this assumption can be avoided, if the slab is with tiny periodic modulation on its top surface. Then we can apply the method presented in Ref. 31 where the Walker equation [18] with modulated Fig. 3. (Color online) (a) Dispersion relation for planar MC formed by Ni inclusions embedded in Fe matrix in the small lat- tice constant regime (a = 50 nm). The filling fraction was fixed to the value 0.55 for different shapes of inclusions. Both, the disper- sions (a) and (b) the profiles of the out-of-plane component of dynamical magnetization do not show significant difference for the systems with different shapes of inclusions. We presented only the profiles for the first (bottom row) and the second (top row) band in the center of the Brillouin zone (Γ point). The thickness of the slab is 5 nm. The external field magnitude is 50 mT. The fol- lowing values of the material parameters were assumed: MS,Fe = = 1.752⋅10–6 A/m, lex,Fe = 3.30 nm, and MS,Ni = 0.484⋅10–6 A/m, lex,Ni = 7.64 nm. The gyromagnetic ratio is assumed to have the same value: γ = 176 GHz/T for both materials [30]. 964 Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 Magnonic crystals — prospective structures for shaping spin waves in nanoscale boundary conditions was solved in the plane wave basis to find magnonic spectrum in the magnetostatic approxima- tion. The planar MCs of more sophisticated geometry, with competing magnetostatic and exchange interactions can be investigated using other methods. In the next two sections we demonstrate finite element method in frequency do- main and micromagnetic simulations with time domain fi- nite difference used to calculate magnonic band structure in planar 1D MC with inhomogeneity across the thickness taken into account. Finite element method in frequency domain In calculations which take into account the inhomoge- neity of the magnetization across the MC thickness we considered planar 1D MCs. The schematic structure con- sidered here is shown in Fig. 5(a). It is composed of infi- nitely long cobalt (Co) and permalloy (Py) stripes. Co and Py stripes are placed side by side. Their dimensions are the same: thickness d = 30 nm and width aCo = aPy = 100 nm, the lattice constant is a = 200 nm. The external magnetic field is directed along the stripes and has the magnitude 0 0 0.1Hµ = T. The material parameters of Co and Py are as follows: Co saturation magnetization MS,Co = 1.45⋅106 A/m, Py saturation magnetization: MS,Py = 0.7⋅106 A/m; Co ex- change constant: ACo = 3⋅10–11 J/m and Py exchange con- stant: ACo = 1.1⋅10–11 J/m. Studied configuration, in which external magnetic field is directed along the stripes (the z- axis), and SW propagation is perpendicular to the stripes (along the y-axis) in the film plane, is called Damon–Eshbach (DE) geometry [18]. For this geometry the static demagnetiz- ing field equals to zero: dem 0.H = Structure is finite in the x direction, which is a thickness of the structure. The finite element method (FEM) is a numerical tech- nique for finding approximate solutions to boundary value problem for partial differential equations. It uses subdivi- Fig. 4. (Color online) Dispersion relation for planar 2D MC formed by Fe inclusions of (a) square, (b) hexagonal, (c) circular shape em- bedded in Ni matrix in square lattice with lattice constant a = 400 nm. The filling fraction was fixed to the 0.55 for different shapes of inclu- sions. The profiles of the out-of-plane component of dynamical mag- netization for five lowest modes in the center of the Brillouin zone are presented in the bottom row in each figure. The spectra and pro- files depends significantly on the shape of inclusion due to the impact of the demagnetizing field dem zH (the insets on the right of each figure). The thickness of the MC is 20 nm. The external field H0 along the z-axis has magnitude 100 mT. The material parameters were assumed the same as in Fig. 3 [30]. Fig. 5. (Color online) (a) 1D MC composed of parallel alternately ordered, connected witch each other Co and Py stripes (the stripes have thickness d and width aCo and aPy; lattice constant is a = aPy + + aCo). The stripes are magnetically saturated along the z-axis by the external magnetic field H0. (b) The unit cell of the 1D MC used in FEM computations. Periodic Bloch boundary conditions (PBC) are assumed along the y axis. (c) The computational area used in FEM. This area extends to large distance along the x axis above and below of the MC. At the borders of this area the magnetostatic potential takes the value 0. Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 965 J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk sion of a computational domain into smaller parts, called finite elements (the computational domain used in our computations is shown in (Fig. 5(c)). FEM encompasses simpler methods for connecting many elementary equa- tions over small subdomains (finite elements), to approxi- mate solution of the complex equation over a large domain. This procedure can be easily realized using commercial software COMSOL Multiphysics. Derivation of the proper equations for COMSOL Multiphysics solver which will be used to solve LL equation in linear approach complement- ed with Maxwell equations is presented below. To take into account appropriate electromagnetic boundary conditions for magnetostatic potential we con- sidered MC which is surrounded by the nonmagnetic die- lectric above and underneath the structure (Fig. 5 (c)), and assumed 0ϕ = at the borders of the computational area, which are far from the MC. The magnetostatic potential shall fulfill Gauss Eq. (2) and because 0/SM z∂ ∂ = we get 2 ,yx mm x y ∂∂ ∇ ϕ = + ∂ ∂ (11) and from Eq. (3) we get ( ) ( ) 0 ex, ex, 0 Ω .  y S y yx y S x x x H m M h Hm i m M h H H m  − +   =   + −    (12) We take here advantage of the Bloch theorem and assumed periodic boundary conditions (PBC) along the y direction. PBC are applied on the boundaries of the unit cell com- posed of Co and Py stripes, which is shown in Fig. 5(b). Solutions of Eqs. (11) and (12), according to Bloch’s theo- rem for dynamic components of the magnetization, are in the form of Eq. (4) with km  being periodic functions of y and dependent on x (i.e. across the thickness of the struc- ture). Similar form for magnetostatic potential is taken: ( ) ( ), , e ,yik yx y x yϕ = ϕ (13) where ϕ is a periodic function of y and depends also on x. yk is a Bloch wavenumber, which can be limited to the first Brillouin zone, i.e., to the range from π/a− to π/ .a Due to symmetry of the structure, this range can be further re- duced to (0,π/ ).k a∈ Substituting Bloch functions of m and ϕ (Eqs. (4) and (13)) to the system of Eqs. (11)–(12) we obtain the eigenvalue problem which is solved with the use of COMSOL 4.3a. According to literature, we used the exchange field in the form appropriate for FEM calculations [46,49]: ( )ex 0 1 2 . S S A M M    = ∇ ⋅ ∇  µ    H m r This formulation of the exchange field superimposed on the interfaces between magnetic materials introduces con- tinuity of the dynamic magnetizations im and i s mA M x ∂ ∂ . The solutions of Eqs. (11) and (12) satisfy also electro- magnetic boundary conditions on the interfaces between magnetic materials and dielectric imposed by the Maxwell equations (i.e., tangential H component and normal B com- ponents are continuous). Magnonic band structure calculated with FEM is shown in Fig. 6 with solid diamonds. The significant dispersion is visible for the first band, the second band is folded-back from the second BZ and has opposite slope. The third band has again positive slope. Between bands are magnonic band gaps. Width of the gap between 1st and 2nd band amounts to 3 GHz. This magnonic band structure is similar to the one measured by Brillouin light scattering experi- ment for 1D MC composed of Co and Py stripes but with lattice constant equal to 500 nm [22]. The identification of the SW excitations related to the bands can be made by the analysis of the SWs’ amplitudes. These are obtained direct- ly from FEM calculations as eigenvectors. The spatial distribution of the selected modes on the ( , )x y plane is shown in Fig. 7. The 1st mode is a funda- mental mode without any nodal points in the BZ center. For this mode, the SW amplitude is higher in Py than in Co. It results from the fact that the frequency range Fig. 6. (Color online) Dispersion relation of SWs in 1D MC. Blue lines correspond to the results obtained in micromagnetic simula- tions. Red diamonds mark the results obtained using FEM calcu- lations in frequency domain. The vertical dashed line points the BZ boundary. Fig. 7. (Color online) Distribution of the SW amplitude’s modu- lus in the unit cell of the 1D MC composed of Co and Py stripes calculated with FEM. The amplitude for the 1st, 2nd, 3rd and 4th SW mode in the BZ center (a), for 0.5 /yk a= π (b) and in the BZ border (c) of the dispersion relation from Fig. 6 is shown. 966 Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 Magnonic crystals — prospective structures for shaping spin waves in nanoscale (around 9 GHz) of 1st band is well below FMR frequency of Co. It means that oscillations in Co stripes can be re- garded as forced by the magnetization oscillations in Py. The modes with increased frequency have increased num- ber of nodal points. For instance the 3rd mode which has two nodal points in the Py stripes, still don’t have oscilla- tions of the amplitude in Co. Interesting is also the 4th band. This mode has the nodal point across the thickness of the Py stripe (see, Fig. 7). This type of the excitations are named perpendicular standing spin waves (PSSW). Due to oscillations of the magnetization on short distance, ex- change interactions determine their properties. The oscilla- tions in neighbor Py stripes are not coupled via Co and the band does not change frequency with increasing .yk How- ever, for modes 1, 2 and 3 the change of the amplitude with increasing wavenumber is found. We can distinguish two kinds of changes with increasing Bloch wave vector ky. First is associated with the change of the phase in the plane wave envelope of the Bloch function Eq. (4). At the BZ boundary (Fig. 7(c)) /yk a= π and the phase of the oscillations changes sign in the neighbor unit cells. Be- cause harmonic oscillations appear only in Py, the zero of the SW amplitude in Co stripes for 1st and 3rd modes, and nonzero for the 2nd mode is observed at the BZ boundary. The second kind of change is caused by surface charac- ter of the DE wave [32]. In the homogeneous film the strength of the DE wave localization at the surface is pro- portional to the wavenumber. Thus, increase of the yk shall result in the surface localization of the SW also in MCs. Indeed, we observe asymmetric distribution of the SW amplitude across the thickness of the MC in the middle of the first BZ (Fig. 7(b)). The amplitude is larger at the top surface for the 1st and 3rd mode, and at the bottom surface for the 2nd mode. This distribution of the amplitude links results of the periodicity and non-reciprocity of the DE wave. The DE wave propagating in the positive y direction (for H0 field pointed at positive z axis) is localized at the top surface of the film, whereas the wave propagating in opposite direction is localized at the bottom surface. 2nd mode in the first BZ is obtained from the mode with 0.75 /yk a=− π (in the second BZ along negative direction of the wave vector) translated with reciprocal lattice vector 2 / G a= π to the first BZ, thus this mode is localized at the bottom surface of the MC. 3rd mode is obtained from wavenumber 1 .25 /yk a= π (3rd BZ) by translation with 2 / G a= − π to the first BZ. Thus, it is localized at the top surface and has larger localization than the 1st mode. The asymmetric distribution of the SW amplitude is lost at the BZ border, because at this wavenumber the superposition of the two counter-propagating waves is formed. The results of the frequency domain FEM calculations will be verified by micromagnetic simulations in the fol- lowing section. Micromagnetic simulations in time and real space domain Micromagnetic simulations (MMS) are very effective computational techniques extensively used for calculations of the magnetization dynamics in nano- and micro-sized ferromagnetic structures [33–39]. MMS are designed for solving numerically full LL Eq. (1) in the time domain and real space. There are many different MMS computational environments. They usually use one from the following nu- merical methods: Finite Difference Method (FDM) or FEM. Here, we use GPU accelerated MuMax3 software which is based on FDM with cuboidal space discretization [35]. Analysis of the SWs dynamics in time domain is very complex and time consuming task. The ideal situation is a single run of simulation which gives the full information about SWs dynamics and can be used to generate disper- sion relation. However, in many cases, the results of several MMS are needed to obtain the dispersion relation of suffi- ciently fine resolution in frequency and wave vector do- main, containing as many as possible branches of SW modes. In this section, the procedure for generation disper- sion relation, based on MMS, will be described. That pro- cess will be exemplified by generation of the dispersion in 1D MC (Fig. 5(a)). The dispersion relation for this system, evaluated using FEM in frequency domain, was already discussed in the previous subsection. MMS are performed in two stages. During the first stage, which can be called static simulations, equilibrium magnetic configuration is obtained. Then, this configura- tion is used as a starting point of the second stage, called dynamic simulations. During this step the static magnetic configuration is disturbed by adding relatively small (to stay in the linear regime and do not destroy equilibrium configuration) dynamic magnetic field. The dynamic field is usually directed orthogonally to the effective magnetic field. This perturbation of the magnetic field in turn induc- es coherent magnetization precession around the equilibri- um direction. The form of dynamic magnetic field deter- mines character of the excited SWs. It is important to use dynamic field which will excite as many SWs modes of different symmetries as possible, to collect information about full spectra of the SWs excitations. There are many possibilities of SWs excitations [40]. Very useful for that purpose are excitations in form of the sinc function, be- cause sinc function is transformed in Fourier space to a window function. Here, we use dynamic magnetic field in form of product of two sinc functions (in space and time domain): ( )dyn , , ,t x y z =b ( ) ( )0 cut 0 cut 0 ˆsinc sinc 2 ,b k y y f t t n   = − π −    (14) where 0b is the maximal value of the magnitude appearing at time 0t and at point 0y along the periodicity direction, n̂ Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 967 J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk is the unit vector perpendicular to the static magnetization. The 0b should be small enough to stay in linear regime. Space discretization of the MC and time steps for mag- netization dynamics, as well as total size of the structure and total time of simulations, determine the correctness of MMS and the resolution of magnonic band structure. How- ever, the requirement of the high resolution shall be com- promised with the size of the computational problem to handle it with available computer resourcers. Maximal value of the frequency taken into account in MMS is max samp1/(2Δ ) f t= and maximal value of the wave vector is max samp/Δ ,k y= π where sampΔt is an interval in time sampling and sampΔy is an interval in space sampling. Re- solution of the dispersion relation (in the frequency and wave vector space) is determined by number of sampl- ings taken into account: max samp simΔ 2 Δ /f f t t= and max sampΔ 2 Δ ,/ yk k y L= where Δf and Δk are frequency and wavenumber steps in dispersion relation, simt is total time of the simulation and yL is total length of the sample. It shows that in many cases sampling intervals can be larg- er than space and time steps used by solver to space dis- cretization and in time integration, respectively. Especially it is useful in the case of Runge–Kutta (RK) solvers [41]. For example, if we investigate SWs propagation up to 25 GHz, then optimal sampling frequency is samplΔt = = 2⋅10–11 s. However, such time steps can be too long in many cases and it may be better to use shorter ones. It is convenient to use maximal step equal to samplΔt and mini- mal step a few orders of magnitude shorter. Similarly is in the case of space discretization, but here, due to technical reasons, resampling is usually not convenient. Another important parameter in simulations is a number of periods in studied structure. This number of repetitions defines resolution of the dispersion relation in the wave vector domain (number of periods is the exact number of different k vectors inside 1st BZ). Summarizing, to obtain fine reso- lution of dispersion relation, it is necessary to simulate the system of many periods through long time. The evolution of the magnetization in the whole system, obtained using MMS, can be written in the form of 4 di- mensional matrix ( , , , )t x y zM . To reduce size of the prob- lem, it is convenient to choose one component of the mag- netization which is perpendicular to static magnetization — ( , , , )cM t x y z where c denotes chosen component of the magnetization, here we use y component. To obtain dis- persion relation for waves propagating along y axis for certain values of ix and iz we calculate two-dimensional Fourier transform from y and t coordinates to yk and .f This is effectively realized with Fast Fourier Transform (FFT) algorithm , ( ) { , , , } ( , , , ).c c t y yM t x y z M f x k z=  Then, if we plot module of the obtained result ( , , , )c i y iM f x k z we will get the dispersion relation. If we aren’t interested in modes visualization we can simply accelerate that process by choosing values ix and iz before calculating FFT: ,( , , , ) ( , ,{ }., )c c i y i t y i iM f x k z M t x y z=  Results obtained in such way depend strongly on the particular form of the excitation. Note that, the different lines in dispersion relation have unlike intensities. Due to that, the disper- sion relation could be unclear. Very helpful in extracting results from simulations are methods used for signal and image processing. More specific information about tech- nical details can be found in Ref. 40. Modes related with dispersion relation can be quite eas- ily visualized. Firstly, the values of frequency 0( )f and wave number ( ),yk corresponding to particular mode, should be selected. In next step we can reduce the size of obtained matrix ( , , , )c yM f x k z by leaving only elements corresponding to the selected frequency, 0( , , , )c yM f x k z = 0 ( , , ).c yfM x k z=  Subsequently, the matrix is filtered, i.e. all values corresponding to different wave numbers than 0 ,k nG+ where n∈ and G is a lattice vector, are replaced by zeros: 0 0 0 ( , , ) ( , , ).c c k nG y yf fM x k z M x k z+ ′δ =  In the further step one- dimensional inverse FFT is performed, separately for every value of x and :z 0 0 0 1 , ( , , ) { ( , , )}. y c f k yk fM x y z M x k z− ′=  Then, the real (imaginary) part of that matrix, 0 0,Re [ ( , , )]f kM x y z ( 0 0,Im [ ( , , )]f kM x y z ), can be plotted as profile of the particular mode corresponding to the points 0 0( , )f k in ( )f k space belonging to the dispersion branch. Algorithm of calculating the dispersion relation and visual- ization of the particular modes is presented in Fig. 8. To obtain dispersion relation already calculated with FEM (Fig. 6) we used 128 alternately ordered wires (Fig. 5(a)), which gives the total length 12.8yL = µm. We discretized this structure to 2600×1×4 cuboidal cells of sizes 5×20×7.5 nm. Due to symmetry along wires axis, it is possible to reduce number of cells along the z-axis during MMS by applying periodic boundary conditions. We’ve intended to study SWs propagation for frequen- cies up to 25 GHz, therefore time of sampling intervals Fig. 8. (Color online) The post-processing algorithm of the MMS results in order to obtain dispersion relation of SWs in MC and to visualize amplitude of the SWs excitations. 968 Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 Magnonic crystals — prospective structures for shaping spin waves in nanoscale was set on samplΔt = 2⋅10–11 s. For SWs excitation we used small, orthogonal to static magnetic field, external dynamic magnetic field defined in Eq. (14), where cutoff frequency was cut 30f = GHz and cutoff wave vector 8 cut 10k = m–1. The parameters which shift value of the maximal excita- tion field in space and time was taken as 0 /2 yy L= and 0t = 1⋅10–9 s. To avoid nonlinear effects, the maximal am- plitude was set as 0 0.05b = T. During dynamic simula- tions LL equation was solved using RK45 (Dormand– Prince) method [35] with maximal time step samplΔt and minimal time step equal to 10–18 s. Total simulation time was set as 100 ns. The dispersion relation obtained in MMS is shown in Fig. 6 with the blue color map. Good agreement with results of FEM in frequency domain is obtained. However, the intensity of the different bands is different. It is related to the shape of the excitation field and its relation to the SWs amplitude [42,43]. The intensity of the 1st band is high, as this band is connected to the fundamental oscillations (Fig. 9) according also with FEM results (Fig. 7(a)). The 2nd mode is invisible in the dispersion relation at 0yk = obtained with MMS. It is because, for the symmet- ric excitation field, used in simulations (Eq. (14)), the asymmetric SWs in Py (see Fig. 7(a)) are not excited. However, with increasing wave number, the intensity of the second band in the dispersion relation increases. This is, because at the BZ boundary the oscillations are in oppo- site phase in neighboring unit cells. Nonreciprocity in magnonic crystals In the previous two sections dispersion relation and amplitude of SW in 1D MC in the DE geometry were in- vestigated. SWs propagating in thin homogeneous films, perpendicular to the external in-plane magnetic field (DE geometry), possess nonreciprocal properties discussed al- ready in the previous section. This means that the SWs propagating in opposite directions (in yk or – yk direction) have amplitude of the dynamic magnetic fields distributed non-symmetrically across the film thickness and they have maximum near one of the surfaces of the film (maxima appear at opposite surfaces for )yk± [18]. The localization of the amplitude has already been observed in the results of FEM calculations for 1D MCs (Fig. 7(b)). However, the frequencies of the oppositely propagating waves were the same, ( ) ( )y yf k f k= − (see Fig. 9). Placing a metal overlayer atop of the film breaks spatial symmetry and also results in breaking symmetry of the dispersion relation with the respect to the change of the wavenumber sign [44]. The dispersion relation reveals significant change of frequency only for the spin waves propagating in one direction, resulting in nonreciprocal dispersion in the film ( ) ( ).y yf k f k≠ − In other words, placing the metal overlayer causes spin waves with equal frequency propagating in converse direction to possess different wave vector magnitudes. The effect of non-reciprocity has impact on the disper- sion of MCs and the magnonic band gaps. Placing a metal plate or perfect electric conductor (PEC) on the top of mag- nonic crystal might leads to destruction of the Bragg condi- tion, since the wavelengths of incident and reflected spin waves at fixed frequency are different [45]. However, the magnonic band gap can still form in this structure, as has been recently demonstrated theoretically and experimental- ly [46,47]. The spin waves propagating in opposite direc- tions and possessing different values of wave vector mag- nitudes can interact and fulfill the general Bragg condition required for band gap opening. FEM in the frequency domain has been implemented here to solve the LL and Maxwell equations in the magne- tostatic approximation (i.e., neglecting dynamical coupling of the magnetostatic field with the electric field), to find the dynamical components of the magnetization vector ( )m r and to obtain the dispersion relation in 1D MC shown in Fig. 10(a). The unit cell used in the calculation is Fig. 9. (Color online) Distribution of the SW amplitude in the unit cell of the 1D MC calculated with MMS. The amplitude for the 1st, and 3rd mode in the BZ center is shown. Good match with amplitudes obtained with FEM is shown [Fig. 7(a)]. Fig. 10. (Color online) (a) A structure of 1D MC with a layer of a perfect electric conductor (PEC) on the top surface. The MC is composed of alternating, infinitely long stripes of Py and Co. The external static magnetic field H0 is applied in the plane of the film, parallel to stripes. The SW propagate along y-axis. (b) The rectangular unit cell used in numerical calculations with periodic boundary conditions (PBC) assumed along the y axis. The effect of a PEC being in the direct contact with the MC is implemented via boundary condition at the top surface. The bot- tom border of the unit cell is far from the MC. Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 969 J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk shown in Fig. 10(b). The equations are defined by Eqs. (11) and (12). Because of PEC on the top surface we need to modify electromagnetic boundary condition at this surface. Here, the continuity of the normal component of the magnetic induction Bx is required. We set the Bx to zero on the top boundary of MC with PEC: 0. x xH m x ∂ϕ = − + = ∂ The dispersion calculated for MC with Co and Py stripes of width 250 nm, with dielectric surroundings is shown in Fig. 11(a) with dashed lines. Material parameters were taken from Refs. 48, 49 and the external magnetic field of 0.2 T was assumed. The dispersion is in good agree- ment with the experimental data from Ref. 23. The edges of the magnonic band gaps are indicated by points 1 and 2 at the BZ border in Fig. 11(a). This SW dispersion is re- ciprocal. The nonreciprocity is introduced with the asymmetric boundary conditions between the top and the bottom sur- face of the MC. The magnonic band structure calculated for the same MC but with PEC at the top surface is shown by solid lines in Fig. 11(a). The dispersion has now nonre- ciprocal character and the edges of the band gap are shifted towards the center of BZ (points 3 and 4 in Fig. 11(a)). The band structure possess an indirect band gap. The Fig. 11(b) shows the distribution of )e[ ]R (ym r along the y-axis. The profiles of the dynamic magnetization components are lat- erally quantized for higher bands, showing that effective wave vector parameter increases with the number of mode. Magnonic band gap in three-dimensional magnonic crystals In the last part we present the theory of 3D MCs. Con- cerning bi-component 3D MCs the big challenge is their fabrication, especially if the lattice constant is in the na- nometer range. Generally, there are two approaches to the fabrication of such structures: top-down and bottom-up [50]. In top-down methods holes are drilled in the bulk material (matrix) and filled with another magnetic material (inclusions). Such methods are common rather in 2D case. The idea of bottom-up techniques is to prepare the lattice of inclusions and then to fill empty spaces with the matrix. This idea gives possibility to make 3D MCs using self- assembling magnetic nanoparticles (NPs) as a template. One of the most extensively studied example of magnetic NPs is magnetoferritin (mFT), a biomimetic NP based on a ferritin, a protein used in living organisms to store an iron in a nontoxic form [51,52]. The usage of cage-like proteins to grow magnetic NPs has a number of advantages [53] (extremely high level of homogeneity, variety of the sizes and properties of protein cages, diversity in physical or chemical functionality of pro- tein shells) which allows to control the self-assembly pro- cess without modifying the NPs obtained inside the protein cages. Especially, mFT NPs can be filled with numerous magnetic materials resulting in different magnetic proper- ties of the NPs [54,55]. The protein crystallization tech- nique used to crystallize mFT NPs, allows to produce high- ly ordered 3D structures up to about 0.4 mm in size [56]. Obtained mFT crystals have high quality fcc structure and the lattice constant about 18.5 nm. An interesting effect is a reduction of the lattice constant to ca. 14 nm as a result of dehydration [57]. Moreover, it was shown theoretically that dried mFT crystals have the crystallographic structure and the lattice constant almost optimized for the occur- rence of a complete magnonic band gap [58]. The object of this part of our study is a bi-component MC based on mFT crystal, i.e., a crystal consisting of mFT NPs (inclusions) arranged periodically in a ferromagnetic host material (B). The geometry of such MC is limited to fcc structure, in which mFT NPs crystallize, and the diame- ter of the inclusions is fixed at 8 nm (the diameter of the magnetic core of fully loaded mFT). The minimal lattice constant of such MC is 11.314 nm (mFT cores are touching each other). Constituent materials are characterized by magnetic parameters: the saturation magnetization SM and the exchange length ex .l The contrast for each parameter is defined as the ratio of its value in mFT to its value in the matrix. An external magnetic field strong enough 0 0( 0.1 Hµ = T) to saturate the system is applied along one of the main crystal axes. The studied system is shown schematically in Fig. 12(c). To determine the spin-wave spectra of 3D MCs we use PWM already introduced in Sec. 3 for 2D MCs. The mag- netostatic field can be expanded in Fourier series using Eq. (2), as it was done for 2D MCs: Fig. 11. (Color online) (a) The dispersion relation of the 1D MC composed of alternating Co and Py stripes with one side metal- ized (solid lines) and with both dielectric surroundings (dashed lines) for the bias field of µ0H0 = 0.2 T. (b) The absolute value of the dynamic magnetization |Re [my(y)]| of the first three modes: I, II and III of the Co/Py metalized structure for the wavenumbers marked in (a) with circles. 970 Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 Magnonic crystals — prospective structures for shaping spin waves in nanoscale ( ) ( ) ( ) ( )( ) ( ) ( ) 2 , , 2 , e x x x x x y y y i x k G m k G k G m h + ⋅+ + + + =− + ∑ k k k G r G G G r k G ( ) ( ) ( ) ( )( ) ( ) ( ) 2 , , 2 , e x x y x x y y x i y k G m k G k G m h + ⋅+ + + + =− + ∑ k k k G r G G G r k G ( ) ( ) ( ) 2 dm 2 e ,S z i z M G H H ⋅≡ = −∑ G r G G r r G _______________________________________________ where position vector, the wave vector and reciprocal lat- tice vectors have now three components ( , , ),x y z=r ( , , )x y zk k k=k and ( , , )x y zG G G=G and the whole sys- tem is periodic in all three dimensions. In a consequence also elements of the matrix (10) in the eigenvalue problem (9) are modified: ,yyxx ij ijM M= − (15) , , 2 ( ) ( , ) ( )x x j y y jxx ij S i j j k G k G M M + + = − + G G k G (16) 2 0 ex( ) ( ) ( ) ( )xy ij j l l j S i jij l M H l M= δ + ⋅ + − ++ −∑ k G k G G G G G 2 2 , , , 2 2 ( ) ( ) ( ) ( ),y y j z i z j S i j S i j j i j k G G G M M + + + − − − + + G G G G k G G G (17) 2 0 ex( ) ( ) ( ) ( )yx ij j l l j S i jij l M H l M= − δ − + ⋅ + − − −∑ k G k G G G G G 2 2 , , , 2 2 ( ) ( ) ( ) ( ).x x j z i z j S i j S i j j i j k G G G M M + + − − + − + + G G G G k G G G (18) By diagonalization of the matrix (10) with sub-matrices defined in Eqs. (15)–(18), we obtained the reduced fre- quencies Ω (eigenvalues) and the eigenvectors — coeffi- cients of the Bloch expansion of the dynamic part of the magnetization (5), as in 2D MC. In Fig. 12 we present two spin-wave spectra plotted along the line connecting high-symmetry points in the first BZ. The lattice constant is assumed to be 18.5 nm (this value correspond to mFT crystal as prepared) and the mag- netic parameters of the mFT NPs are: SM = 0.346⋅106 A/m and A = 10–11 J/m [59]. Two matrix materials are used: Co in (a) or Ni in (b). For Co matrix a wide (ca. 100 GHz), Fig. 12. Spin-wave spectra (up to 500 GHz) of mFT MC with (a) cobalt and (b) nickel matrix for lattice constant 18.5 nm (mFT crystal as prepared). The spectra are plotted along the high-symmetry path in the 1st Brillouin zone shown in (d). Shaded area represents: in (a) the complete magnonic band gap and in (b) the overlapping of the first and second band. (c) Schematic depiction of the mFT MC struc- ture with coordinating system used. Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 971 J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk complete magnonic band gap appears. This gap is located between the first and the second band. In the case of Ni matrix the band gap failed to open which is the conse- quence of too weak contrast of saturation magnetization SM [60]. The existence of a critical SM contrast below which the complete gap does not open stems from the fact that in a spectrum for 3D structure a complete bandgap is indirect gap. The bottom and the top of this gap is related to different propagation directions (in contrast to the direc- tional gap, i.e., the gap for particular direction of the prop- agation). Roughly speaking, the width of the directional gap depends on the contrast of magnetic parameters, and its central frequency is determined by the value of the wave vector at the boundary of the first BZ. If directional gaps for different directions of propagation are too narrow, then the overlapping of neighboring bands appears. This very effect underlies the lack of complete magnonic gap for the Ni matrix (Fig. 12(b)). As we have already mentioned, the mFT NPs can be loaded with different magnetic materials which means that the contrast of SM can be modified in quite broad range. In Fig. 13(a) we show the evolution of the complete mag- nonic gap vs. the total magnetic moment of mFT NPs in MCs, which are based on four matrix materials: Fe, Co, Py and Ni. The lattice constant is fixed to 14 nm (dried mFT crystal). Presented dependences range from the magnetic moment twice the typical value of mFT 4(10 )Bµ [57] which is gradually reduced to half of this value. The gap widens quickly with increasing SM contrast (decreasing magnetic moment of the mFT NPs), consequently the complete magnonic gap opens even for Ni matrix at 9.9⋅103 (a value only 1% less than in typical mFT NPs). Another way to tailor the magnonic gap is adjusting to the MC lattice constant. As we already mentioned, the mFT crystal reduces its lattice constant from 18.5 to 14 nm while drying. Additionally, functionalization of the exter- nal surface of protein cage is also promising feature in terms of controlling the lattice constant of the MC. We plotted the evolution of edges of the complete magnonic gap vs. lattice constant for different matrix materials (Co, Fe, and Py) to explore the ranges in which the complete magnonic gap exists (see Fig. 13(b)). The magnetic mo- ment of inclusions is fixed to 104 µB. In this case there is no gap for Ni matrix for any value of the lattice constant. For all cases we studied, the complete magnonic gap changes a lot with the lattice constant and its maximal width is observed at approximately 13 nm (the precise val- ue depends on the matrix material). The existence of this maximum is related to the concur- rence of the exchange and dipolar interactions. Coexist- ence of long- and short-range interactions leads to very interesting phenomena in variety of systems [61–64]. In the MCs concerned in this study it strongly influences the spin wave profiles of two lowest modes (Fig. 14). For small lattice constant the modes of lowest frequency have a bulk character due to the strong exchange interactions be- tween excitations in neighboring areas of the MC. As the lattice constant grows, the significance of exchange inter- actions fades while dipolar interactions gain in importance. As a result for the larger lattice constant, stronger spin wave profiles are concentrated in mFT NPs (first mode) or Fig. 13. (Color online) The edges of the complete magnonic gap vs. (a) magnetic moment of inclusions and (b) lattice constant of MC for a Fe, Co, Py and Ni matrix. In (a) the lattice constant of MC is fixed at 14 nm (dried mFT crystal). In (b) magnetic moment of inclu- sions is fixed at the typical value for mFT NPs (104 µB) and there is no gap for Ni matrix in this case. Fig. 14. (Color online) Profiles of dynamic magnetization for the two lowest spin-wave modes in an mFT/Co MC in a plane perpendicular to the external field and passing through the centers of mFT NPs (the contours of which are represented by circles). 972 Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 Magnonic crystals — prospective structures for shaping spin waves in nanoscale in the matrix (second mode). For in-between region, where both types of interactions have comparable weightiness, maximum of the gap width appears. The profiles of the spin waves have great influence on the effective damping as well. Using PWM based on LL equation with damping included [65], we found the con- centration of the mode in one of the constituent material changes effective damping towards its value in this ma- terial [65]. This issue can lead to strong anisotropy of the ef- fective damping. As an example in Fig. 15 we show results for the lowest modes obtained in MC containing spherical scattering centers (material B, MS,B = 0.194⋅106 A/m, AB = 3.996⋅10–12 J/m — these values are close to yttrium iron garnet) disposed in the matrix (material A, MS,A = =1.752⋅106 A/m, AA = 2.1⋅10–11 J/m — the values close to Fe). The lattice is assumed to be simple cubic with the lattice constant a = 10 nm and the sphere radius equal to 3.628 nm. The Gilbert damping parameter was chosen as αA = 0.0019 and αB = 0.064. The eigenmodes being solutions of LL are characterized by complex eigenvalues i′ ′′Ω =Ω+ Ω . In Fig. 15(a) and (b) we plotted the real and imaginary part of Ω, respectively. The imaginary part  ′Ω is a measure of the life time of the particular SW excitation. However, the value ′′Ω has to be referred to ′Ω for direct comparison of the damping of different modes. In Fig. 15(c) we show so called “figure of merit” (FOM) which is the real part of the frequency divided by its imag- inary part. As we can see the damping of the lowest mode is much higher (lower FOM) than for the second mode. To explain this feature we plotted spin-wave profiles (the dis- tribution of the dynamical part of magnetization) for two lowest bands (see Fig. 15(d)). It is clearly seen that the lowest mode is strongly concentrated in scattering centers (with higher damping) while the second one — in the ma- trix (with lower damping). Therefore, the effective damp- ing for these two modes is much different. For the same reason the effective damping depends on the propagation direction as well. For example: comparing the profiles of the lowest mode for two different wave vectors, the con- centration of the dynamical magnetization in the spheres is a bit stronger for direction R than at point Γ. This results in the difference of effective damping: 0.059 for R and 0.045 for Γ. Conclusions We have demonstrated the application of different nu- merical methods: PWM, FEM and MMS for calculations of the magnonic band structure in 1D, 2D and 3D MCs. The PWM is an useful tool for calculation of SW proper- ties in periodic structures, independent on its dimensionali- ty, shape of the elements and crystallographic lattice. However, PWM is limited to fully saturated materials and it works only in linear approximation. In the case of PWM applied to planar MCs, the considered MCs have to be ho- mogeneous across the thickness and relatively thin. In the proposed FEM the last assumption is avoided. We showed usefulness of this method in the study of inhomogeneous excitations across the thickness of the homogeneous film of MC. We showed that FEM can be easily extended to include the boundary effect induced by conductive materi- als surrounding MC film. This method is suitable for in- vestigation of the nonreciprocal properties in SW dynamics in 1D and 2D MCs. MMS are the most general tool, which avoid assumptions used in PWM and FEM. However, this method requires the complex post-processing to extract information obtained directly in PWM or FEM. The other drawback of MMS is a long time of the calculations, how- ever, this is cut down by implantation of the GPU units for calculations. We have demonstrated usefulness of MMS to study magnonic band structure in planar 1D MC. Its exten- sion to the planar 2D MC can be easily done. We have demonstrated, that the surface character of the Damon–Eshbach wave is also preserved in magnonic crys- tals. Surface localization increases with the band number, however the localization exists only inside the Brillouin zone, this is excluding Brillouin zone center and border. Fig. 15. (Color online) The real and imaginary parts of the frequency in the first BZ for sc MC in (a) and (b), respectively. (c) The figure of merit (FOM) for the same structure. (d) The distribution of the amplitude of the dynamical components of the magnetization across the planes perpendicular to the external field. Planes (001) and (002) are between and across the spheres, respectively. The profiles from the first and the second band in Γ and R point of the first BZ are shown. (Figures taken from Ref. 65). Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 973 J. Rychły, P. Gruszecki, M. Mruczkiewicz, J.W. Kłos, S. Mamica, and M. Krawczyk The surface property of the spin wave excitation is further exploited by covering plate of the magnonic crystal with a conductor. Due to that the nonreciprocal dispersion rela- tion is introduced. The band structure in 2D magnonic crystals is very complex. It is caused not only by structurization, but also by the demagnetizing field. This field introduces additional spatial inhomogeneity for spin waves and modifies spin wave dispersion. It also makes the band structure strongly dependent on shape of inclusion and type of lattice for magnonic crystal. Pronounced effect is the localization of low frequency spin waves in the areas of the lowered in- ternal magnetic field. However, the inhomogeneity of the internal magnetic field becomes unimportant for magnonic crystals with small lattice constant where demagnetizing effect can be neglected. In this limit, the band structure is only slightly dependent on the shape of inclusions and type of the lattice. For 3D magnonic crystals, characterized by small lattice constant, for the structures with characteristic sizes comparable to diameter of magnetoferritin crystals, the wide band gap was found. Finally, we have pointed out that the spatial distribution of different materials could be explored for tailoring effective damping of spin waves. Acknowledgments The research leading to these results has received fund- ing from Polish National Science Centre project DEC-2- 12/07/E/ST3/00538 and from the EUs Horizon2020 re- search and innovation programme under the Marie Sklo- dowska-Curie GA No644348. The numerical calculation were performed at Poznan Supercomputing and Network- ing Center (grant No. 209). 1. E. Yablonovitch, Phys. Rev. Lett. 58, 2059 (1987). 2. S. John, Phys. Rev. Lett. 58, 2486 (1987). 3. J.D. Joannopoulos, S.G. Johnson, J.N. Winn, and R.D. Meade, Photonic Crystals—Moulding the Flow of Light, Princeton Univ. Press, Princeton (2008). 4. S. Guenneau, R.C. McPhedran, S. Enoch, A.B. Movchan, M. Farhat, and N.-A.P. Nicorovici, J. Opt. 13, 024014 (2011). 5. M. Maldovan, Nature 503, 209 (2013). 6. Y. Wang, E.W. Plummer, and K. Kempa, Adv. Phys. 60, 799 (2011). 7. M. Krawczyk and G. Grundler, J. Phys.: Condens. Matter 26, 123202 (2014). 8. C.G. Sykes, J.D. Adam, and J.H. Collins, Appl. Phys. Lett. 29, 388 (1976). 9. J.P. Parekh and H.S. Tuan, Appl. Phys. Lett. 30, 667 (1977). 10. P. Hartemann, J. Appl. Phys. 62, 2111 (1987). 11. J.M. Owens, J.H. Collins, C.V. Smith, and I.I. Chiang, Appl. Phys. Lett. 31, 781 (1977). 12. A.V. Voronenko, S.V. Gerus, and V.D. Haritonov, Sov. Phys. J. 31, 76 (1988). 13. J.O. Vasseur, L. Dobrzynski, B. Djafari-Rouhani, and H. Puszkarski, Phys. Rev. B 54, 1043 (1996). 14. M. Krawczyk and H. Puszkarski, Acta Physica Polonica A 93, 805 (1998). 15. Y.V. Gulyaev and S.A. Nikitov, Dokl. Phys. 46, 687 (2001). 16. H. Puszkarski and M. Krawczyk, Solid State Phenom. 94, 125 (2003). 17. B. Kalinikos and A. Slavin, J. Phys. C 19, 7013 (1986). 18. D.D. Stancil and A. Prabhakar, Spin Waves Theory and Applications, Springer (2009). 19. M. Krawczyk, M. Sokolovskyy, J.W. Kłos, and S. Mamica, Adv. Condens. Matter Phys. 2012, 1 (2012). 20. J.W. Lau and J.M. Shaw, J. Phys. D 44, 303001 (2011). 21. R.E. Camley, T.S. Rahman, and D.L. Mills, Phys. Rev. B 27, 261 (1983). 22. Linear and Nonlinear Spin Waves in Magnetic Films and Superlattices, M.G. Cottam (ed.), World Scientific Publish- ing Co., Singapore (1994). 23. Z.K. Wang, V.L. Zhang, H.S. Lim, S.C. Ng, M.H. Kuok, S. Jain, and A.O. Adeyeye, Appl. Phys. Lett. 94, 083112 (2009). 24. G. Gubbiotti, S. Tacchi, M. Madami, G. Carlotti, A.O. Adeyeye, and M. Kostylev, J. Phys. D 43, 264003 (2010). 25. W. Lu and C.M. Lieber, Nature Mater. 6, 841 (2007). 26. A.O. Adeyeye, S. Jain, and Y. Ren, IEEE Transact. Magnet. 47, 1639 (2011). 27. J. Kaczer and L. Murtinova, Phys. Status Solidi A 23, 79 (1974). 28. M.L. Sokolovskyy and M. Krawczyk, J. Nanoparticle Research 13, 6085 (2011). 29. J.W. Kłos, M.L. Sokolovskyy, S. Mamica, and M. Krawczyk, J. Appl. Phys. 111, 123910 (2012). 30. J.W.  Kłos  et al., J. Appl. Phys. 111, 123910 (2012). Copy- rights 2012 American Institute of Physics. 31. I. Lisenkov, D. Kalyabin, S. Osokin, J.W. Klos, M. Krawczyk, and S. Nikitov, J. Magn. Magn. Mater. 378, 313 (2015). 32. R.W. Damon and J. Eshbach, J. Phys. Chem. Solids 19, 308 (1961). 33. R. Hertel, W. Wulfhekel, and J. Kirschner, Phys. Rev. Lett. 93, 257202 (2004). 34. K.M. Lebecki, M.J. Donahue, and M.W. Gutowski, J. Phys. D 41, 175005 (2008). 35. A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, AIP Adv. 4, 107133 (2014). 36. S.-K.Kim, J. Phys. D 43, 264004 (2010). 37. W. Kim, S.-W. Lee, and K.-J. Lee, J. Phys. D 44, 384001 (2011). 38. S.-K. Kim, S. Choi, K.-S. Lee, D.-S. Han, D.-E. Jung, and Y.-S. Choi, Appl. Phys. Lett. 92, 212501 (2008). 39. G. Venkat, D. Kumar, M. Franchin, O. Dmytriiev, M. Mruczkiewicz, H. Fangohr, A. Barman, M. Krawczyk, and A. Prabhakar, IEEE Trans. Magn. 49, 524 (2013) 40. D. Kumar, Development of Nanoscale Magnetic Systems for Spin Wave Propagation, PhD Diss., Center for Research in 974 Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 Magnonic crystals — prospective structures for shaping spin waves in nanoscale Nanoscience and Nanotechnology, University of Calcutta, (2014). 41. J.C. Butcher, Runge–Kutta Methods in Numerical Methods for Ordinary Differential Equations, John Wiley & Sons (2008). 42. K.-S. Lee, D.-S. Han, and S.-K. Kim, Phys. Rev. Lett. 102, 127202 (2009). 43. K. Di, H.S. Lim, V.L. Zhang, M.H. Kuok, S.C. Ng, M.G. Cottam, and H.T. Nguyen, Phys. Rev. Lett. 111, 149701 (2009). 44. A.G. Gurevich and G.A. Melkov, Magnetization Oscillations and Waves, CRC Press (1996). 45. E.N. Beginin, Y.A. Filimonov, E.S. Pavlov, S.L. Vysotskii, and S.A. Nikitov, Appl. Phys. Lett. 100, 252412 (2012). 46. M. Mruczkiewicz, M. Krawczyk, G. Gubbiotti, S. Tacchi, Y.A. Filimonov, D.V. Kalyabin, I.V. Lisenkov, and S.A. Nikitov, New J. Phys. 15, 113023 (2013). 47. M. Mruczkiewicz, E.S. Pavlov, S.L. Vysotsky, M. Krawczyk, Yu.A. Filimonov, and S.A. Nikitov, Phys. Rev. B 90, 174416 (2014). 48. M. Sokolovskyy, J. Kłos, S. Mamica, and M. Krawczyk, J. Appl. Phys. 111, 07C515 (2012). 49. Z.K. Wang, V.L. Zhang, H.S. Lim, C. Ng, M.H. Kuok, S. Jain, and A.O. Adeyeye, ACS Nano 4, 643 (2010). 50. V.V. Kruglyak et al., Magnonic metamaterials, X.-Y. Jiang (ed.), Metamaterial (InTech, 2012), available from: http://www.intechopen.com/books/metamaterial/magnonic- metamaterials. 51. P.M. Harrison and P. Arosio, Biochim. Biophys. Acta 1275, 161 (1996). 52. S.C. Andrews, Biochim. Biophys. Acta 1800, 691 (2010). 53. M.L. Flenniken, M. Uchida, L.O. Liepold, S. Kang, M.J. Young, and T. Douglas, Curr. Top. Microbiol. Immunol. 327, 71 (2009). 54. I. Yamashita, K. Iwahori, and S. Kumagai, Biochim. Bio- phys. Acta 1800, 846 (2010). 55. H. Yoshimura, Colloids Surfaces A 282–283, 464 (2006). 56. M. Okuda, J.-C. Eloi, S.E.W. Jones, A. Sarua, R.M. Richardson, and W. Schwarzacher, Nanotechnology 23, 415601 (2012). 57. O. Kasyutich, D. Tatchev, A. Hoell, F. Ogrin, C. Dewhurst, and W. Schwarzacher, J. Appl. Phys. 105, 07B528 (2009). 58. S. Mamica, M. Krawczyk, M.L. Sokolovskyy, and J. Romero-Vivas, Phys. Rev. B 86, 144402 (2012). 59. O. Kasyutich, A. Sarua, and W. Schwarzacher, J. Phys. D 41, 134022 (2008). 60. S. Mamica, J. Appl. Phys. 114, 043912 (2013). 61. S. Mamica, R. Józefowicz, and H. Puszkarski, Acta Phys. Pol. A 94, 79 (1998). 62. S. Mamica, J.C.S. Lévy, P. Depondt, and M. Krawczyk, J. Nanopart. Res. 13, 6075 (2011). 63. S. Mamica, J.C.S. Lévy, and M. Krawczyk, J. Phys. D 47, 015003 (2014). 64. M. Krawczyk, H. Puszkarski, J.-C. S. Levy, S. Mamica, and D. Mercier, J. Magn. Magn. Mater. 246, 93 (2002). 65. J. Romero-Vivas, S. Mamica, M. Krawczyk, and V.V. Kruglyak, Phys. Rev. B 86, 144417 (2012). Low Temperature Physics/Fizika Nizkikh Temperatur, 2015, v. 41, No. 10 975 Introduction Model and wave equation Planar magnonic crystals Plane wave method Magnonic band structure in thin 2D MCs Finite element method in frequency domain Micromagnetic simulations in time and real space domain Nonreciprocity in magnonic crystals Magnonic band gap in three-dimensional magnonic crystals Conclusions Acknowledgments