OSA's Digital Library

Optics Express

Optics Express

  • Editor: Andrew M. Weiner
  • Vol. 21, Iss. 10 — May. 20, 2013
  • pp: 12022–12037
« Show journal navigation

Simple magneto–optic transition metal models for time–domain simulations

Christian Wolff, Rogelio Rodríguez-Oliveros, and Kurt Busch  »View Author Affiliations


Optics Express, Vol. 21, Issue 10, pp. 12022-12037 (2013)
http://dx.doi.org/10.1364/OE.21.012022


View Full Text Article

Acrobat PDF (898 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Efficient modelling of the magneto–optic effects of transition metals such as nickel, cobalt and iron is a topic of growing interest within the nano–optics community. In this paper, we present a general discussion of appropriate material models for the linear dielectric properties for such metals, provide parameter fits and formulate the anisotropic response in terms of auxiliary differential equations suitable for time–domain simulations. We validate both our material models and their implementation by comparing numerical results obtained with the Discontinuous Galerkin time–domain (DGTD) method to analytical results and previously published experimental data.

© 2013 OSA

1. Introduction

Fairly recently, nickel and other magnetic transition metals have gained some interest within the nano–photonic and plasmonic communities. Most effects that have been experimentally studied so far, such as plasmon–induced increase and control in the magneto–optic Kerr effect (MOKE) [1

1. E. Th. Papaioannou, V. Kapaklis, E. Melander, B. Hjörvarsson, S. D. Pappas, P. Patoka, M. Giersig, P. Fumagalli, A. García–Martín, and G. Ctistis, “Surface plasmons and magneto–optic activity in hexagonal Ni anti–dot arrays,” Opt. Express 19, 23867–23877 (2011) [CrossRef] [PubMed] .

, 2

2. E. Melander, E. Östman, J. Keller, J. Schmidt, E. Th. Papaioannou, V. Kapaklis, U. B. Arnalds, B. Caballero, A. García–Martín, J. C. Cuevas, and B. Hjörvarsson, “Influence of the magnetic field on the plasmonic properties of transparent Ni anti-dot arrays,” Appl. Phys. Lett. 101, 063107 (2012) [CrossRef] .

] or control over the surface plasmon polariton propagation via an external magnetic field [3

3. V. V. Temnov, G. Armelles, U. Woggon, D. Guzatov, A. Cebollada, A. Garcia–Martin, J.–M. García–Martín, T. Thomay, A. Leitenstorfer, and R. Bratschitsch, “Active magneto-plasmonics in hybrid metal-ferromagnet structures,” Nat. Phot. 4, 107–111 (2010) [CrossRef] .

], were accompanied by frequency–domain simulations such as the scattering matrix method [4

4. A. García–Martín, G. Armelles, and S. Pereira, “Light transport in photonic crystals composed of magneto–optically active materials,” Phys. Rev. B 71, 205116 (2005) [CrossRef] .

]. This is convenient as the material can be modelled directly by means of an experimentally determined dielectric tensor for each frequency independently. However, the non–perturbative numerical study of non–linear effects such as the generation of higher harmonics or the inclusion of active media via the Maxwell–Bloch equations demands time–domain simulations. It turns out that even in the absence of an external magnetic field, the losses in metals such as nickel cannot be satisfactorily described by the otherwise successful Drude model. Concerning magneto–optics, effects such as the Faraday effect and MOKE are sometimes accounted for by using a non–dispersive anisotropic dielectric tensor. This approach can be found in some studies using off-the-shelf finite difference time–domain codes, and may work for narrow–band excitations but provide at best rough estimates.

Thus, efficient time–domain formulations for the dispersive nature both of the isotropic and the magneto–optic properties of transition metals have been missing. The aim of this work is to provide a first yet powerful approach to these two issues. To this end, we generalize the Drude model to include retardation effects, and find that this ansatz is very well suited to model the metals in question. We then show that adding the Lorentz force to the equations of motion for the polarization currents provides a reasonable description of magneto–optic effects and can be incorporated easily and very efficiently into existing time–domain methods.

2. Ferromagnetic metals

We discuss the topic of ferromagnetic metals using the case of nickel as an example. From a microscopic point of view, nickel is quite different from silver or gold. Firstly, the most obvious difference—ferromagnetism—is but one symptom of the fact that magnetic metals represent highly correlated electron systems rather than quasi–free electron gases that are commonly associated with simple metals. Secondly, a ferromagnetic metal does not feature one but several electron species. In nickel, e.g., we find a hybridized 4s/4p–electron system and two systems of 3d–electrons; one with upward polarized spin and one with downward polarized spin. All species differ in many important properties such as effective mass or equilibrium density not only from free electrons but also from each other. Eventually, this fact is for example the origin of the spontaneous magnetization. The fundamental disagreement between optical measurements and the simple Drude model (see fig. 1) indicates that at least one of these properties is relevant for the dielectric response of magnetic metals in the visible range.

Fig. 1 Left panel: The asymptotic behavior of the measured permittivity’s imaginary parts of nickel, cobalt and iron (solid lines) [13] compared to the asymptotic behavior of the Drude model (dashed line, ω−3 scaling) and a curve with ω−1 scaling (dash–dotted).
Right panel: Comparison between the experimental permittivity of nickel (solid lines) [13] and a Drude model (dashed lines) that was fitted to the real part (parameters: ε = 1.78, h̄ωP = 7.66eV, h̄γD = 0.697eV). Apart from inter–band transitions, the real part fits acceptably. However, the gross mismatch in the imaginary part cannot be overcome by adding a few Lorentz terms, rendering the Drude model unsuitable.

2.1. Magnetic susceptibility

The magnetic susceptibilities of the ferromagnetic transition metals feature ferromagnetic resonance in the microwave regime. This effect is characterized by Lorentzian resonance lines. Due to the spectral position several orders of magnitude below the optical regime and the asymptotic behavior of this function (ω−2 in the real and ω−3 in the imaginary part), we choose to neglect them in our work. Another possible contribution to the dynamic magnetic susceptibility might be caused by reorientable magnetic dipoles. This can be expected to at least contribute to the static magnetic susceptibility and to have the frequency dependence of the Debye model [6

6. M. Y. Koledintseva, K. N. Rozanov, A. Orlandi, and J. L. Drewniak, “Extraction of Lorentzian and Debye parameters of dielectric and magnetic dispersive materials for FDTD modeling,” J. Electr. Eng.–Slovak 53, 97–100 (2002).

]. Its real part would vanish at high frequencies but the imaginary part would feature an ω−1–tail. Again, we assume that the relevant time scales are many orders of magnitude larger than the duration of an optical cycle; for the remainder of this manuscript, we thus assume
μ(ω)=1
(1)
at optical frequencies.

2.2. Isotropic dielectric susceptibility

The itinerant electrons in nickel are often modelled to first approximation as a Landau Fermi liquid (see e.g. [7

7. V. Korenman, J. L. Murray, and R. E. Prange, “Local-band theory of itinerant ferromagnetism. I. Fermi-liquid theory,” Phys. Rev. B 16, 4032–4047 (1977) [CrossRef] .

], Korenmann et al.). Thus, the metal could be regarded as a free gas of Landau–quasi–particles (QPs). As they essentially behave like the constituents of a free electron gas (albeit with different effective parameters), the dielectric response of the metal in the infrared may be expected to follow the ubiquitous Drude–Sommerfeld susceptibility
χDrude(ω)=ωP2ω(ω+iγD),
(2)
with the plasma frequency ωP/(2π) and the scattering rate γD. Using the well–known connection between a dynamic susceptibility χ(ω) and the corresponding dynamic conductivity
σ(ω)=iωχ(ω),
(3)
this can be equivalently formulated as a differential equation for the polarization current
tjD(t)+γDjD(t)=ωP2E(t).
(4)
Unfortunately, this model fails to reproduce the losses as can be seen in Fig. 1. The ω−3–asymptotics of the imaginary part of χDrude(ω) is in contradiction to the experimental data, which merely exhibit a ω−1–behavior. While the Drude model can overall very nicely reproduce the real part of the permittivity of nickel, the losses are grossly under–represented.

In a Drude metal, all electrons react immediately to any change of the driving E-field. The reason for this behavior is that the electrons are assumed to be non–interacting particles and the relaxation mechanism is a continuous scattering process that is not influenced by the external field. In a correlated system such as nickel, however, this is not the case. Instead, it takes some time to establish a new collective equilibrium state after the system has been perturbed from the outside, specifically by the optical field. We can effectively account for it in Eq. (4) by introducing a memory kernel Z(s):
tjD(t)+γDjD(t)=0Z(s)E(ts)ds.
(5)
This ansatz is very generic, and in principle, we need not make any assumptions about the origin of Z(s). In the case of a ferromagnet, however, we use it to effectively describe the process of relaxing electron correlation effects. We may expect this to happen on the time scale τ = lc/vF, where lc is the typical correlation length of a few nanometers and vF is the Fermi velocity, giving very roughly τ = 1fs. This first estimate provides a sanity check when fitting to experimental data as we do later on. More importantly, it allows us to make the model numerically manageable. Due to the convolution, Eq. (5) is ill–suited for time–domain simulations. However, the short time scale of Z(s) suggests to expand E into its Taylor series around t, which yields
tjD(t)+γDjD(t)=n=0ζ(n)(tnE)(t),
(6)
ζ(n)=(1)nn!0snZ(s)ds
(7)
where is the n-th moment of Z(s), apart from a prefactor. We truncate the series after the linear term and – upon substituting ωP2=ζ(0) and τ = ζ(1)/ζ(0) – end up with
tjD(t)+γDjD(t)=ωP2(1+τt)E(t).
(8)
The corresponding susceptibility for this simplest retarded Drude model is
χrD(ω)=ωP2(1iωτ)ω(ω+iγD),
(9)
and indeed has the desired ω−1 asymptotic of the imaginary part. Within our linear approximation, τ must be positive to ensure that χrD(ω) does not have a negative imaginary part at any frequency. This would lead to numerical instability even if the permittivity was correctly reproduced in the spectral range of interest.

As a side remark, we would like to mention that Eq. (9) can be decomposed into the sum of a conventional Drude model and a Debye model, where the latter arises due to the retardation effects. We conclude this paragraph discussing the inclusion of further terms. The second order term in the convolution leaves the model unchanged except for a modification of the values of τ and the background permittivity. A finite number of higher order terms leads to diverging losses at high frequencies. This is unphysical as any susceptibility should vanish at sufficiently short wavelength. Besides that, higher order terms are very likely to create negative losses (i.e. unphysical gain) in some frequency region, leading to numerical instabilities in time–domain simulations. In principle, we could also replace the simple Drude damping γD in Eq. (5) with another convolution. This would model a retarded dependence of a scattering mechanism on the movement of electrons (any dependence on the driving field is contained in Z(s)). In the case of QP–phonon scattering, this would model e.g. lattice heating and we assume that no such process is relevant at room temperature and on the time scale of an optical cycle. Even more, to first approximation, it would just renormalize ωP and γD without qualitatively changing anything.

After having settled for a treatment of correlation effects, we should consider the second main property of ferromagnetic metals: The presence of at least two spin–polarized electron species. In order to account for this, we assign individual retarded Drude models to either electronic subsystem. On closer inspection, we indeed find that the plasma frequency
ωP=e2nmeff
(10)
depends on the effective mass meff of the quasi–particles as well as on the partial carrier densities. This holds true for every electronic subsystem, such as spin–up electrons and spin–down electrons. Thus, we find individual plasma frequencies ωP, ωP for either subsystem. Likewise, we may expect that the corresponding retardation times τ, τ are quite different, too. At room temperature, however, the Drude damping γD is dominated by QP–phonon scattering because the QP life time is very large close to the Fermi surface and magnons are thermally less excited due to their higher energies. Because QP-phonon scattering is not sensitive to spin orientation, we may assume that both subsystems are subject to essentially the same scattering rate. As a result, the conduction electrons are described by two retarded Drude models with a common denominator. These terms can be replaced by a single retarded Drude term with effective plasma frequency and effective relaxation time
ωP2=[ωP]2+[ωP]2,
(11)
τ=([ωP]2τ+[ωP]2τ)ωP2.
(12)

Finally, we should comment on the treatment of inter–band transitions. In a simple metal, they are degenerate with respect to spin. This degeneracy is lifted in a ferromagnet. However, we do not expect this to have any fundamental impact such as retardation effects because the destination states are unoccupied. Thus, we describe them using the Lorentz oscillator model as it is good custom for simple metals, too:
χLorentz(ω)=ΔεLωL2ωL2iωγLω2.
(13)
The only difference is that we can expect Lorentz poles to appear in pairs because of the spin–related splitting of the destination state energies.

3. Magneto–optic materials

After having presented our model for the isotropic optical properties of nickel, we now turn towards the treatment of magneto–optical effects. First, we discuss the treatment of inter–band transitions because this is valuable for the description of dielectric Faraday–materials as well.

3.1. Insulators

Although insulators are not exactly the topic of this paper, we mention them here because they exhibit the Faraday effect if an inter–band transition is split by an external magnetic field. Recently, Berman showed in two independent quantum–mechanical calculations [8

8. P. R. Berman, “Optical Faraday rotation,” Am. J. Phys. 78, 270–276 (2009) [CrossRef] .

] that within the linear regime, the resulting optical properties are identical to those of a simple classical Lorentz oscillator with a Lorentz force added to the equation of motion. Thus, we describe magneto–optic inter–band transitions by a corresponding model, which in our notation reads as:
t2j+γtj+ωL2j=ΔεLωL2E+Ω×tj.
(14)
Here, we have introduced an auxiliary vector
Ω=emBext
(15)
that is aligned along the external magnetic field Bext and whose absolute value Ω = |Ω⃗| is the angular cyclotron frequency of the particle with charge −e and mass m.

Next, we would like to briefly comment on the value of Eq. (14) for simulations of the Faraday effect. Although the discussion of numerical technicalities is the topic of later parts of this paper, we would like to restrict the discussion on magneto–optic insulators to this section. Within a time–domain method and, in particular, within the the Discontinuous Galerkin time–domain method, Eq. (14) is a rather elegant way to handle Faraday materials rather than assuming a non–dispersive anisotropic dielectric. This is for two main reasons. First, it is the more realistic model because magneto–optic effects in insulators usually are dispersive with a resonant character. This argument holds for any time–domain method. Secondly, it has great technical advantages for our specific method: A non–dispersive anisotropic dielectric would directly appear in Maxwell’s curl equations and, thus, requires a modification of the numerical flux. Although this is straight forward in two dimensions [9

9. M. König, K. Busch, and J. Niegemann, “The Discontinuous Galerkin time–domain method for Maxwells equations with anisotropic materials,” Phot. Nano. Fund. Appl. 8, 303–309 (2010) [CrossRef] .

], it becomes considerably involved in 3D [10

10. J. Alvarez, L. D. Angulo, A. R. Bretones, and S. G. Garcia, “3–D Discontinuous Galerkin time–domain method for anisotropic materials,” IEEE Antenn. Wireless Propag. Lett. 11, 1182 (2012) [CrossRef] .

] and to the best of our knowledge, nobody has rigorously demonstrated convergence of such a 3D anisotropic electrodynamic DGTD method so far. However, a dispersive and local material is described by an ordinary differential equation without any spatial derivatives and does not interfere with the numerical flux. In section 6.1, we show that this is actually true.

3.2. Metals

In a metal, there are generally two sources of magneto–optic effects. The first contribution stems from inter–band transitions and can be treated as described in section 3.1. The second contribution comes from Drude–like currents. This has been previously modelled [4

4. A. García–Martín, G. Armelles, and S. Pereira, “Light transport in photonic crystals composed of magneto–optically active materials,” Phys. Rev. B 71, 205116 (2005) [CrossRef] .

] by adding the Lorentz force to Eq. (4), which arises naturally from the Drude model of charged point particles scattering from an ionic background. The applicability of this entirely classical picture to the degenerate Fermi system of a real metal might be questioned. Starting from the hydrodynamic (HD) model [11

11. J. E. Sipe, V. C. Y. So, M. Fukui, and G. I. Stegeman, “Analysis of second–harmonic generation at metal surfaces,” Phys. Rev. B 21, 4389–4402 (1980) [CrossRef] .

], we can show that – at least for simple metals – the classical result is in fact consistent with a more microscopic theory. The HD model itself can be derived from quantum mechanics, e.g. via a moment expansion of the Wigner function (see chapter 3 in the book by Wyatt [12

12. R. E. Wyatt, Quantum Dynamics with Trajectories (Springer, 2005).

]). Neglecting tunneling effects and after closure of the moment hierarchy by explicitly assuming a degenerate gas of Fermions with charge −e, it reads in our notation:
etndivj=0,
(16)
tj+γj+div[jv]emgradp(n)=em[enEμj×H].
(17)
Here, −e is the electron charge, n the local electron density, v the corresponding velocity field, j = −env is the current density and p(n) = αn5/3 is the density–dependent pressure of free fermions with some constant α. For the sake of readability, we suppressed the explicit space–and time–dependence of n, v, j, E and H. From this, the differential equation of the Drude model (which is a linear model) can be recovered by first linearizing the equations. To this end, we introduce equilibrium quantities n0, j0 = 0, E0 = 0, H0 = Bext/μ and small deviations δn, δj = j, δE = E, δH. We drop all products of the latter four quantities (i.e. the nonlinear terms) and find
etδndivj=0,
(18)
tj+γjemp(n0)gradδn=em[en0Eμj×H0],
(19)
where p′(n) = (∂p(n))/(∂n). By averaging over a sufficiently large volume, all non–local effects are smeared out, which is equivalent to dropping the spatial derivatives:
tδn=0,
(20)
tj+γj=n0e2mEemj×Bext
(21)
=ωP2E+Ω×j.
(22)
This is the Drude model with a Lorentz force added.

Unfortunately, our retarded Drude model is not a limit of the HD model or any other established model we are aware of. However, it is a generalization of the Drude model and, thus, should contain at least the Lorentz force term in its response to an external magnetic field:
tjD(t)+γDjD(t)=ωP2(1+τt)E(t)+Ω×jD(t).
(23)
This yields even quantitatively acceptable numerical results (see sections 5 and 6). Furthermore, the appearance of the Lorentz force in Eq. (23) is justified by the fact that the magneto–optic Drude model is regained in the limit of vanishing τ. In principle, there might be a τ–dependent prefactor to the Lorentz force but this is not of practical relevance as it can be absorbed in the cyclotron frequency Ω, which is eventually fitted to experimental data. Finally, higher order terms might be expected to arise from a more elaborate microscopic description.

3.3. Susceptibility tensor

For our parameter fits in section 5, we need the explicit forms of such tensors for a magnetic field along the z–axis. The general solution to this is
χ_(ω)=ζ(iω)iω(A2+B2)(AB0BA000A2+B2A)
(27)
with A = p(n)(−iω) and B = (−iω)n−1Ωk. Obviously, the tensor element χzz is just the isotropic model. The nontrivial entries of the Lorentz oscillator model with Lorentz force are:
χxx(ω)=ω2+iωγωL2(ω2+iωγωL2)2ω2Ω2,
(28)
χxy(ω)=iΔεLωL2ωΩ(ω2+iωγωL2)2ω2Ω2.
(29)
The nontrivial entries of the Drude model with Lorentz force are (in accordance with [4

4. A. García–Martín, G. Armelles, and S. Pereira, “Light transport in photonic crystals composed of magneto–optically active materials,” Phys. Rev. B 71, 205116 (2005) [CrossRef] .

], García–Martin et al.):
χxx(ω)=ω+iγω[(ω+iγ)2Ω2],
(30)
χxy(ω)=iωP2Ωω[(ω+iγ)2Ω2].
(31)
Finally, the nontrivial entries of our retarded Drude model with Lorentz force are:
χxx(ω)=(ω+iγ)(1iωτ)ω[(ω+iγ)2Ω2],
(32)
χxy(ω)=iωP2Ω(1iωτ)ω[(ω+iγ)2Ω2].
(33)
For the magnetic field along the x– or y–direction, the index subscripts have to be interchanged in the above result. As the material response is linear, an arbitrary magnetic field direction can be handled by superimposing individual magnetic fields along the three Cartesian directions.

4. Auxiliary differential equations

Dispersive media in time–domain simulations of Maxwell’s equations are commonly treated using the technique of auxiliary differential equations (ADE, see e.g. our DGTD review [5

5. K. Busch, M. König, and J. Niegemann, “Discontinuous Galerkin methods in nanophotonics,” Laser & Photon. Rev. 5, 773–809 (2011).

]). That is a set of ordinary differential equations with respect to time that describe the dynamics of one or several polarization currents. This approach remains well–suited for magneto–optic anisotropies that are caused by the Lorentz force. The straight forward approach would be to explicitly calculate the susceptibility tensor in frequency–domain (section 3.3), introduce individual polarization currents for each diagonal element and each pair of off–diagonal ones and to derive an ADE for each of them [9

9. M. König, K. Busch, and J. Niegemann, “The Discontinuous Galerkin time–domain method for Maxwells equations with anisotropic materials,” Phot. Nano. Fund. Appl. 8, 303–309 (2010) [CrossRef] .

]. The resulting inflation in the number of auxiliary fields can be avoided by simply adding the Lorentz force to the ADE of the isotropic model as stated in the model Eq. (14) and Eqs. (22) and (23). In this section, we present ADEs that arise in this fashion. It is noteworthy that for a specific model there may be several different ADE formulations; we tried to chose the ones that are least prone to numerical inaccuracies. So far, the implementation of our equations in our DGTD code has always appeared to be robust. We present a convergence analysis in section 6.

4.1. Drude and Lorentz oscillator models

We begin our discussion with the most basic metal model: the Drude model. Its ADE formulation is well known from the literature [5

5. K. Busch, M. König, and J. Niegemann, “Discontinuous Galerkin methods in nanophotonics,” Laser & Photon. Rev. 5, 773–809 (2011).

]:
tjD=γDjDωP2E.
(34)
The inclusion of the Lorentz force is trivial; we simply add it as another driving term:
tjD=γDjD+Ω×jDωP2E.
(35)
The parameter m in Ω⃗ = eB/m is regarded as an effective mass that must be fitted to experimental data.

4.2. Retarded Drude model

Next, we will derive an ADE formulation for the isotropic retarded Drude model presented in section 2.2. There, we already mentioned that this model can be decomposed into a conventional Drude and a Debye model. However, treating both contributions together reduces the number of required auxiliary fields and can be implemented with overall fewer CPU instructions, i.e. is computationally more efficient. We start with the susceptibility Eq. (9):
χrD(ω)=ωP2(1iωτ)ω(ω+iγD).
(44)
Using Eq. (3), the polarization current satisfies
(1+γDτ1τ1iω)jrD(ω)=ωP2τE(ω).
(45)
Once again, we make an ansatz for an auxiliary field
p(ω)=jrD(ω)τ1iω,
(46)
resulting in
jrD(ω)=(τ1γD)p(ω)+ωP2τE(ω),
(47)
iωp(ω)=γDp(ω)+ωP2τE(ω).
(48)
After applying a Fourier transform, we find our result:
tp(t)=γDp(t)+ωP2τE(t),
(49)
jrD(t)=(τ1γD)p(t)+ωP2τE(t).
(50)
Only the first equation has to be propagated in time; the expression for the current can be directly incorporated into the update equation for the electric field, so this model requires no more degrees of freedom than the conventional Drude model.

Finally, we derive an ADE formulation for the magneto–optic retarded Drude model. We start by adding the Lorentz force to the equation of motion:
tjrD(t)+γDjrD(t)Ω×jrD(t)=ωP2(1+τt)E(t).
(51)
After Fourier transformation, the polarization current satisfies
(1+γDτ1Ω×τ1iω)jrD(ω)=ωP2τE(ω).
(52)
Using the ansatz
p(ω)=(γDτ1)jrD(ω)Ω×jrD(ω)iωτ1,
(53)
we find the final ADE:
jrD(t)=p(t)+ωP2τE(t),
(54)
tp(t)=γDp(t)+(1τγD)ωP2E(t)+Ω×jrD(t).
(55)

5. Parameter fits

Table 1. Fit parameters of a retarded Drude model with two additional Lorentz oscillators for nickel, cobalt and iron. The cyclotron frequencies Ω are valid for magnetically saturated nickel films. The fits were performed for different windows of photon energies . The fit errors are upper bounds to the relative error throughout the respective energy windows.

table-icon
View This Table
Fig. 2 Panels (a) and (b): Literature [13] values for the complex permittivity of nickel and our fit using the retarded Drude model with two additional Lorentz oscillator terms.
Panel (c) and (d): Magneto–optic properties of nickel from the literature [14] and our fit. Starting from the isotropic model shown in the left part, the individual cyclotron frequencies Ω of the retarded Drude and the two Lorentz oscillator models (three free parameters all together) were used to fit the polar MOKE. For easier comparison with Fig. 4 in [14], the ordinate of panel (d) is scaled with the photon energy squared. In all panels, only the white energy ranges were considered for the fit; nonetheless, the model reproduces fairly the experimental data at photon energies above the fit range (red shaded areas). All model parameters are listed in Table 1.

Overall, we find good agreement for the isotropic [13

13. P. B. Johnson and R. W. Christy, “Optical constants of transition metals: Ti, V, Cr, Mn, Fe, Co, Ni, and Pd,” Phys. Rev. B 9, 5056–5070 (1974) [CrossRef] .

] and qualitative agreement for the magneto–optic response [14

14. Š. Višňovský, V. Pařízek, M. Nývlt, P. Kielar, V. Prosser, and R. Krishnan, “Magneto–optical Kerr spectra of nickel,” J. Magn. Magn. Mater. 127, 135–139 (1993) [CrossRef] .

]. Better agreement can be obtained if the target frequency range is restricted but we went with this moderate quality fit as we intend to compare to a broad band experiment in section 6.2. We attribute our low–frequency deviations to shortcomings in the isotropic fit. The only truly startling difference is the spectral offset of the upper magneto–optic resonance. This is in fact a problem of the off–diagonal elements, as can be seen by comparing to Fig. 4 of [14

14. Š. Višňovský, V. Pařízek, M. Nývlt, P. Kielar, V. Prosser, and R. Krishnan, “Magneto–optical Kerr spectra of nickel,” J. Magn. Magn. Mater. 127, 135–139 (1993) [CrossRef] .

]. We would like to point out the relatively large deviations between the curves from different experiments in that plot. These differences suggest that the preparation of the film plays a major role for the (magneto–)optical properties of nickel. We thus attribute the ostensibly poor quantitative accuracy of our fit to the fact that the isotropic and magneto–optic parameters were fitted to very different experimental works. We conclude that for quantitative simulations, both the isotropic and magneto–optic properties of nickel should be measured for films that were grown under the same conditions as those used for the intended device. Consistency between the literature fits and these calibration measurements should be checked.

Fig. 3 Panels (a) and (b): Literature [13] values for the complex permittivity of cobalt and our fit using the retarded Drude model with two additional Lorentz oscillator terms.
Panels (c) and (d): Literature [13] values for the complex permittivity of iron and our fit using the retarded Drude model with two additional Lorentz oscillator terms.
All model parameters are listed in table 1.
Fig. 4 Left panel: Vertical cut through the numerical domain to simulate the reflection spectrum of a nickel half space; the length unit is one nanometer. PEC indicates a perfect electric conductor boundary, PML stands for perfectly matched layers.
Central panel: Relative numerical error of the reflectivity in the absence of a magnetic field for orders p of the interpolation polynomials.
Right panel: Absolute numerical error of the polar MOKE for various polynomial orders p. The analytical curves to be numerically reproduced are shown in the third panel of Fig. 2. We plot the absolute error because of the change in sign.

Assuming that the results by Johnson and Christy and those by Višňovský in fact do not exhibit such inconsistencies, the high–frequency mismatch might indicate a breakdown of the Lorentz force model to magneto–optics and suggests a refined model.

6. Numerical results

To validate our implementation and judge its overall performance, we performed several numerical tests. For this, we first checked the isotropic nickel model by comparing the numerically calculated reflection spectrum of a half space under normal incidence with the analytical expectation. Then, we checked the off–diagonal elements of the dielectric tensor by comparing numerical and analytical polar MOKE for the same setup. With this, we tested the correctness of our implementation, the convergence and the numerical stability. Finally, we calculated a reflection spectrum of a recently published nanostructure to show that our fitting errors are within the accuracy margins of state–of–the–art experimental and numerical techniques. Our models do not require more degrees of freedom than a conventional Drude model with two Lorentz oscillator poles. In our implementation, the number of floating point operations per magneto–optic metallic tetrahedron is increased by 70% (11% without magneto–optic terms).

6.1. Halfspace problem

To test the material models, we compared the analytical reflection spectra [15

15. Š. Višňovský, “Magneto–optical Ellipsometry,” Czech. J. Phys. B 36, 625–650 (1986) [CrossRef] .

] of a metallic half-space under normal incidence with our numerical data obtained with the DGTD method [5

5. K. Busch, M. König, and J. Niegemann, “Discontinuous Galerkin methods in nanophotonics,” Laser & Photon. Rev. 5, 773–809 (2011).

]. The numerical reference system consisted of a slender domain with square cross section. A vertical cut is depicted in Fig. 4. For testing an isotropic material model, we would have used perfect electric and perfect magnetic conductor boundary conditions on the long domain faces and an appropriate polarization state of the incident pulse in order to mimic an infinitely extended film. This was not possible because the Kerr effect rotates the polarization state, so we imposed periodic boundaries, instead. The nickel film at the bottom had a thickness of 500nm, which is sufficient to suppress spurious reflections at the (perfect electric conductor) bottom domain boundary. We injected a broad band pulse at the top and started to record the Poynting flux through the injection plane when the pulse center reached the metal surface. Pulse width and total domain height were chosen such that the field in the injection plane was negligible in this moment. The air parts were divided into tetrahedra with a maximum edge length hmax = 100nm, the metallic parts into tetrahedra with hmax = 50nm; on each tetrahedron the field components were expanded into Lagrange polynomials [5

5. K. Busch, M. König, and J. Niegemann, “Discontinuous Galerkin methods in nanophotonics,” Laser & Photon. Rev. 5, 773–809 (2011).

]. We varied the polynomial order to show convergence.

Figure 4 shows the numerical error for some polynomial orders. Overall, we find that the numerical scheme is stable and converges to the analytical result based on our expressions for the susceptibility tensor. Each polynomial order reduces the error by about an order of magnitude, which we take as evidence for exponential convergence. This behavior is expected for DG methods [5

5. K. Busch, M. König, and J. Niegemann, “Discontinuous Galerkin methods in nanophotonics,” Laser & Photon. Rev. 5, 773–809 (2011).

] and demonstrates that an anisotropic, local and dispersive material model is compatible with the usual (isotropic) Upwind flux. This is in contrast to an anisotropic background permittivity, which would require a modification of the numerical flux.

In the reflection spectra, we find decreasing accuracy for higher photon energies. The reason for this is the increasingly unfavorable ratio between tetrahedron size and wavelength. Considering the Kerr spectra, it should be kept in mind that we plot the absolute error and that the correct result is of the order of 0.1. This means that the Kerr spectra seem to be much more prone to numerical inaccuracies. We attribute this to the rather small physical effect and the fact that the asymmetric numerical grid causes spurious polarization changes in lower polynomial orders. The results seem to be better in the center of the frequency window. This is due to the Gaussian shape of the excitation pulse and the correspondingly higher impact of numerical noise at extreme frequencies. In this context, it should be noted that the simulation covers nearly three octaves.

6.2. Experimental data

Finally, we calculated the polar MOKE spectrum of a structure similar to that investigated by Papaioannou et al. [1

1. E. Th. Papaioannou, V. Kapaklis, E. Melander, B. Hjörvarsson, S. D. Pappas, P. Patoka, M. Giersig, P. Fumagalli, A. García–Martín, and G. Ctistis, “Surface plasmons and magneto–optic activity in hexagonal Ni anti–dot arrays,” Opt. Express 19, 23867–23877 (2011) [CrossRef] [PubMed] .

]. To this end, we set up a numerical domain as sketched in Fig. 5. The basic shape was a hexagonal column with perfectly matched layers on the top and the bottom end and periodic boundaries on the six vertical faces. Close to the bottom end, we included one unit cell of a structured nickel film according to the experimental parameters on a silicon substrate. The grid consisted of tetrahedra with hmax = 100nm in air and hmax = 50nm in the silicon substrate and the metal. The expansion basis consisted of fifth order Lagrange polynomials.

Fig. 5 Left panel: Vertical cut through the numerical domain used to calculate the polar MOKE reported on in [1], Papaioannou et al.. The lateral cross section is a regular hexagon with an edge length of 271.4nm corresponding to a triangular lattice with a lattice constant of 470nm; a cylindrical hole with radius 137.5nm is cut out from the nickel film. The length unit is one nanometer.
Right panel: Polar MOKE of the hexagonal array of these nickel anti–dots. For better comparison, the plot is arranged in the same way as in the reference. Differing from the experiment, we studied normal incidence.

For the sake of simplicity, we studied normal incidence in our simulation with periodic boundaries, whereas the experiment was done at an angle of 8°. Furthermore, we neglected the very thin spacer layer of titanium and the very thin top gold layer and modelled the silicon substrate as a simple dielectric with constant permittivity ε = 11.9. We believe that all these modifications do not dramatically spoil our results (Fig. 5). This has to be compared to the experimental and numerical data in Fig. 3(b) of [1

1. E. Th. Papaioannou, V. Kapaklis, E. Melander, B. Hjörvarsson, S. D. Pappas, P. Patoka, M. Giersig, P. Fumagalli, A. García–Martín, and G. Ctistis, “Surface plasmons and magneto–optic activity in hexagonal Ni anti–dot arrays,” Opt. Express 19, 23867–23877 (2011) [CrossRef] [PubMed] .

]. Some features such as the strong feature above 3.5 eV are better reproduced by our method as compared to the frequency–domain results published therein, some other (e.g. the secondary resonance at 2.9 eV) not quite that well. We attribute the latter as well as the slight shift in the main resonance around 2.7 eV to the missing gold top layer as the surface plasmon polariton dispersion relation is very sensitive to the losses of the very surface. We find that (considering the errors introduced by our simplifications) the overall agreement between our method and the experiment is comparable to the published theoretical data, which proves the practical value of our models.

7. Summary

We demonstrated that our material models can be very efficiently incorporated into time–domain methods such as the DGTD method. We furthermore showed in sec. 5 that they reproduce the isotropic optical properties of a class of transition metals that had been previously barely accessible in time–domain simulations. Concerning the magneto–optical properties, we focused on nickel because of its practical relevance and its good experimental characterization.

Acknowledgments

We are deeply indebted to Prof. Wolfgang Nolting for inspiring discussions. This work was financially supported by the Deutsche Forschungsgemeinschaft (DFG) within project BU 1107/8-1.

References and links

1.

E. Th. Papaioannou, V. Kapaklis, E. Melander, B. Hjörvarsson, S. D. Pappas, P. Patoka, M. Giersig, P. Fumagalli, A. García–Martín, and G. Ctistis, “Surface plasmons and magneto–optic activity in hexagonal Ni anti–dot arrays,” Opt. Express 19, 23867–23877 (2011) [CrossRef] [PubMed] .

2.

E. Melander, E. Östman, J. Keller, J. Schmidt, E. Th. Papaioannou, V. Kapaklis, U. B. Arnalds, B. Caballero, A. García–Martín, J. C. Cuevas, and B. Hjörvarsson, “Influence of the magnetic field on the plasmonic properties of transparent Ni anti-dot arrays,” Appl. Phys. Lett. 101, 063107 (2012) [CrossRef] .

3.

V. V. Temnov, G. Armelles, U. Woggon, D. Guzatov, A. Cebollada, A. Garcia–Martin, J.–M. García–Martín, T. Thomay, A. Leitenstorfer, and R. Bratschitsch, “Active magneto-plasmonics in hybrid metal-ferromagnet structures,” Nat. Phot. 4, 107–111 (2010) [CrossRef] .

4.

A. García–Martín, G. Armelles, and S. Pereira, “Light transport in photonic crystals composed of magneto–optically active materials,” Phys. Rev. B 71, 205116 (2005) [CrossRef] .

5.

K. Busch, M. König, and J. Niegemann, “Discontinuous Galerkin methods in nanophotonics,” Laser & Photon. Rev. 5, 773–809 (2011).

6.

M. Y. Koledintseva, K. N. Rozanov, A. Orlandi, and J. L. Drewniak, “Extraction of Lorentzian and Debye parameters of dielectric and magnetic dispersive materials for FDTD modeling,” J. Electr. Eng.–Slovak 53, 97–100 (2002).

7.

V. Korenman, J. L. Murray, and R. E. Prange, “Local-band theory of itinerant ferromagnetism. I. Fermi-liquid theory,” Phys. Rev. B 16, 4032–4047 (1977) [CrossRef] .

8.

P. R. Berman, “Optical Faraday rotation,” Am. J. Phys. 78, 270–276 (2009) [CrossRef] .

9.

M. König, K. Busch, and J. Niegemann, “The Discontinuous Galerkin time–domain method for Maxwells equations with anisotropic materials,” Phot. Nano. Fund. Appl. 8, 303–309 (2010) [CrossRef] .

10.

J. Alvarez, L. D. Angulo, A. R. Bretones, and S. G. Garcia, “3–D Discontinuous Galerkin time–domain method for anisotropic materials,” IEEE Antenn. Wireless Propag. Lett. 11, 1182 (2012) [CrossRef] .

11.

J. E. Sipe, V. C. Y. So, M. Fukui, and G. I. Stegeman, “Analysis of second–harmonic generation at metal surfaces,” Phys. Rev. B 21, 4389–4402 (1980) [CrossRef] .

12.

R. E. Wyatt, Quantum Dynamics with Trajectories (Springer, 2005).

13.

P. B. Johnson and R. W. Christy, “Optical constants of transition metals: Ti, V, Cr, Mn, Fe, Co, Ni, and Pd,” Phys. Rev. B 9, 5056–5070 (1974) [CrossRef] .

14.

Š. Višňovský, V. Pařízek, M. Nývlt, P. Kielar, V. Prosser, and R. Krishnan, “Magneto–optical Kerr spectra of nickel,” J. Magn. Magn. Mater. 127, 135–139 (1993) [CrossRef] .

15.

Š. Višňovský, “Magneto–optical Ellipsometry,” Czech. J. Phys. B 36, 625–650 (1986) [CrossRef] .

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(160.3820) Materials : Magneto-optical materials
(260.3910) Physical optics : Metal optics

ToC Category:
Materials

History
Original Manuscript: February 19, 2013
Revised Manuscript: March 16, 2013
Manuscript Accepted: March 20, 2013
Published: May 10, 2013

Citation
Christian Wolff, Rogelio Rodríguez-Oliveros, and Kurt Busch, "Simple magneto–optic transition metal models for time–domain simulations," Opt. Express 21, 12022-12037 (2013)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-21-10-12022


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. E. Th. Papaioannou, V. Kapaklis, E. Melander, B. Hjörvarsson, S. D. Pappas, P. Patoka, M. Giersig, P. Fumagalli, A. García–Martín, G. Ctistis, “Surface plasmons and magneto–optic activity in hexagonal Ni anti–dot arrays,” Opt. Express 19, 23867–23877 (2011). [CrossRef] [PubMed]
  2. E. Melander, E. Östman, J. Keller, J. Schmidt, E. Th. Papaioannou, V. Kapaklis, U. B. Arnalds, B. Caballero, A. García–Martín, J. C. Cuevas, B. Hjörvarsson, “Influence of the magnetic field on the plasmonic properties of transparent Ni anti-dot arrays,” Appl. Phys. Lett. 101, 063107 (2012). [CrossRef]
  3. V. V. Temnov, G. Armelles, U. Woggon, D. Guzatov, A. Cebollada, A. Garcia–Martin, J.–M. García–Martín, T. Thomay, A. Leitenstorfer, R. Bratschitsch, “Active magneto-plasmonics in hybrid metal-ferromagnet structures,” Nat. Phot. 4, 107–111 (2010). [CrossRef]
  4. A. García–Martín, G. Armelles, S. Pereira, “Light transport in photonic crystals composed of magneto–optically active materials,” Phys. Rev. B 71, 205116 (2005). [CrossRef]
  5. K. Busch, M. König, J. Niegemann, “Discontinuous Galerkin methods in nanophotonics,” Laser & Photon. Rev. 5, 773–809 (2011).
  6. M. Y. Koledintseva, K. N. Rozanov, A. Orlandi, J. L. Drewniak, “Extraction of Lorentzian and Debye parameters of dielectric and magnetic dispersive materials for FDTD modeling,” J. Electr. Eng.–Slovak 53, 97–100 (2002).
  7. V. Korenman, J. L. Murray, R. E. Prange, “Local-band theory of itinerant ferromagnetism. I. Fermi-liquid theory,” Phys. Rev. B 16, 4032–4047 (1977). [CrossRef]
  8. P. R. Berman, “Optical Faraday rotation,” Am. J. Phys. 78, 270–276 (2009). [CrossRef]
  9. M. König, K. Busch, J. Niegemann, “The Discontinuous Galerkin time–domain method for Maxwells equations with anisotropic materials,” Phot. Nano. Fund. Appl. 8, 303–309 (2010). [CrossRef]
  10. J. Alvarez, L. D. Angulo, A. R. Bretones, S. G. Garcia, “3–D Discontinuous Galerkin time–domain method for anisotropic materials,” IEEE Antenn. Wireless Propag. Lett. 11, 1182 (2012). [CrossRef]
  11. J. E. Sipe, V. C. Y. So, M. Fukui, G. I. Stegeman, “Analysis of second–harmonic generation at metal surfaces,” Phys. Rev. B 21, 4389–4402 (1980). [CrossRef]
  12. R. E. Wyatt, Quantum Dynamics with Trajectories (Springer, 2005).
  13. P. B. Johnson, R. W. Christy, “Optical constants of transition metals: Ti, V, Cr, Mn, Fe, Co, Ni, and Pd,” Phys. Rev. B 9, 5056–5070 (1974). [CrossRef]
  14. Š. Višňovský, V. Pařízek, M. Nývlt, P. Kielar, V. Prosser, R. Krishnan, “Magneto–optical Kerr spectra of nickel,” J. Magn. Magn. Mater. 127, 135–139 (1993). [CrossRef]
  15. Š. Višňovský, “Magneto–optical Ellipsometry,” Czech. J. Phys. B 36, 625–650 (1986). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1 Fig. 2 Fig. 3
 
Fig. 4 Fig. 5
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited