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 Optical Society of America
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, 2] or control over the surface plasmon polariton propagation via an external magnetic field , were accompanied by frequency–domain simulations such as the scattering matrix method . 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.
The paper is organized as follows: In section 2, we discuss general properties of non–simple transition metals and model the isotropic part of the dielectric response as a generalized Drude model that takes retardation effects to lowest order into account. This is followed by section 3, where we motivate why most magneto–optic anisotropies can be expressed in terms of a Lorentz force term in the equations of motion for the polarization currents. We state the frequency–domain susceptibility tensor that arises from including the Lorentz force into an arbitrary dispersive material model. In section 4, we formulate the previously developed material models in terms of auxiliary differential equations that can be efficiently incorporated in time–domain codes such as the Discontinuous Galerkin time–domain (DGTD) method . We present parameter fits to experimental data for the isotropic material response of popular transition metals and for their magneto–optic response in the fifth section. Section 6 is devoted to the validation of our implementation by comparing numerical to analytical results and to experimental reference data.
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.
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 . 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
2.2. Isotropic dielectric susceptibility
The itinerant electrons in nickel are often modelled to first approximation as a Landau Fermi liquid (see e.g. , 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 susceptibilityFig. 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):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
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
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:
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.
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  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:
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 , it becomes considerably involved in 3D  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.
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  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 , 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 ). 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:
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: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
The addition of the Lorentz force to the differential equations of a polarization current leads to an anisotropic susceptibility tensor in frequency–domain. We can derive it without explicitly specifying the underlying isotropic material model. All previously presented equations of motion for the polarization current including the Lorentz force can be formulated in frequency–domain in this general form:Eq. (5). We reformulate this using the 3 × 3 identity matrix 1 and the Levy-Civita tensor ε̳: Eq. (3) and find
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 is4], García–Martin et al.):
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 ). 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 . 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 :
The equally popular Lorentz oscillator model, which is commonly used to model inter–band transitions, requires slightly more attention when introducing the Lorentz force. First, we once again state the ADE formulation of the basic, isotropic model as it is well known from the literature :
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):Eq. (3), the polarization current satisfies
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:
5. Parameter fits
We fitted the free parameters of our isotropic model to the experimental data of Johnson and Christy on transition metals . For the magneto–optic properties, we focused on nickel because this seems to be most popular for upcoming experimental work and because it appears to be better characterized than cobalt or iron. We kept all parameters of the isotropic model from the previous fit fixed and only regarded the three cyclotron frequencies (one for the retarded Drude and each Lorentz oscillator of the isotropic model) as parameters. They were chosen to optimally reproduce the polar MOKE measurements by Višňovský . It would be straightforward to fit to the off–diagonal parts of the susceptibility tensor provided in that publication. However, it is obvious from the figures, that these data strongly depend on which isotropic permittivity is assumed. In our case, this is a fit in itself and the approach outlined above would lead to excess fit errors. The proper way is to fit measured quantities; here, it is the complex Kerr angle, e.g. using the equations provided by Višňovský . The fitted parameters are shown in table 1 and a comparison of the experimental and fitted permittivities for nickel is shown in Fig. 2.
Overall, we find good agreement for the isotropic  and qualitative agreement for the magneto–optic response . 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 . 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.
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  of a metallic half-space under normal incidence with our numerical data obtained with the DGTD method . 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 . 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  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. . 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.
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 . 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.
This paper is definitely not the end of the story but provides some first and simple, yet very efficient magnetic transition metal models that should be sufficient for many simulations. To this end, we first analyzed the asymptotic behavior of the losses at high frequencies and found that they require the inclusion of a retardation term in the Drude model. Then, we collected arguments for modelling magneto–optic effects in general by adding the Lorentz force to the equations of motion of the polarization current. Finally, we derived auxiliary differential equation formulations for the generalized Drude model as well as for the magneto–optic material models, fitted the free parameters to literature values for three different transition metals and showed that the frequency–domain permittivity tensor and our implementation of the material model are consistent. We demonstrated the usefulness of our work by reproducing previously published experimental data. A next desirable step would be to focus on a semi–microscopic description of the electrons in metals with strong correlations. Upon linearization and spatial averaging, this should reproduce our retarded Drude model. Another question posed by the quality of the magneto–optics fits is whether the Lorentz–force model really is sufficient. It was derived directly from the quantum mechanics of a charged harmonic oscillator, but the deviations in the ultra–violet Kerr spectra of nickel spark some doubt about the validity of the harmonic oscillator approximation. Finally, spin–orbit coupling has been identified in the past to play a major role for the Kerr effect but is probably not fully contained in our models.
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.
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] .