Theory of electromagnetic energy transfer in three-dimensional structures

ultramicroscopy ELSEVIER

Ul~amicroscopy61 (1995) 21-27

Theory of electromagnetic energy transfer in three-dimensional structures Jean-Pol Vigneron *, Fatiha Forati, Damien Andr6, Annick Castiaux, Isabelle Derycke, Alain Dereux b~stitute for Studies in Interface Sciences, Facult~s Universitaires Notre-Dame de la Paix, Rue de Bruxelles 61, B-5000 Namur, Belgium

Received 9 May 1995; accepted 12 June 1995

Abstract

The practical computation of electromagnetic energy transfer through an inhomogeneous dielectric or conducting film is considered at the level of vector waves multiple scattering theory. Maxwell's equations are formulated in the Laue representation and the scattering boundary conditions are implemented using a transfer-matrix technique. This formulation leads to a fast algorithm if all convolution products involved are handled by fast Fourier transform techniques.

1. I n t r o d u c t i o n

The theoretical simulation [1] of devices such as near-field optical microscopes [2] or integrated optic gates requires the assessment of the electromagnetic energy flux through complex, three-dimensional structures, with scattering centers whose sizes can be as small as a fraction of a wavelength. The computation of the transmission, reflection or absorption coefficients of a three-dimensional aggregate of small particles require the explicit knowledge of a multiply-scattered field. Quite often, models for :the description of electromagnetic energy transfer in these structures can be formulated in terms of the radiation transmission of a thick, corrugated junction which connects, with parallel planar boundaries, homogeneous, isotropic half spaces. The present paper describes a fast algorithm which allows to compute this energy transfer with little limitations arising from the complexity of the three-dimensional coupling layer. This algorithm, partly borrowed from LEED dynamical theory [3] and recent scanning-tunneling microscopy computational procedures [4], uses a transfer-matrix technique to compute the needed electromagnetic transmission coefficients. In the following section, we give a detailed account of the new algorithm, with enough detail to facilitate its implementation. Existing applications and those currently under development will be presented elsewhere.

* Corresponding author. Fax: + 32 81 724707; E-mail: [email protected]. 0304-3991/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0304-3991(95)00097-6

22

J.-P. Vigneron et al. / Ultramicroscopy 61 (1995) 21-27

2. Algorithm for the computation of electromagnetic energy transfer through an inhomogeneous film We specifically consider the following simple situation: a finite-size inhomogeneous film makes the interface connecting two homogeneous, semi-infinite dielectric media. We will refer to the first of these media as "the wave emitter", and the second as "the wave collector". The incident wave and the series of reflected wave propagate in the emitter, while transmitted waves emerge into the collector. The (not necessarily thin) film separating the emitter from the collector is assumed to contain no current source, so that relevant Maxwell's equations reduce to the simple linear relations

V × E = iW/~oH,

(1)

VXH=

(2)

-iw~0s(r)E.

In this expression, e ( r ) is a space-varying dielectric function, which reduces to the constant values e~ in the emitter (for z < z]) and em in the collector (for z > zm). In the junction region, it assumes arbitrary spatial variations. Our objective is to evaluate the fraction of the incident energy which is transmitted through the junction layer, after the infinite number of scattering events encountered within the junction layer. In order to represent more effectively the scattered field, we make the assumption - with little loss of generality - that the junction layer exhibits a two-dimensional periodic symmetry. This may be the case when the system under study is actually periodic, but otherwise, we enforce a two-dimensional periodic repetition of the mainscattering centers, using a large two-dimensional supercell. In this latter case, a convergence study must be carried out in order to determine whether the size of the supercell is large enough to overcome intercell interactions. 2.1. Electromagnetic wave equation in the Laue representation The supercell periodicity implies that the stationary states of the electromagnetic field (and our scattered wave) propagates in the structure as a lateral Bloch wave, which we write as a Fourier series in the t coordinates O, keeping a real space representation for the normal coordinate z: E ( r ) = E E g , ( z ) e i(k''+g')'p,

(3)

g'

H(r)

= EHg,(z) e i(~''+#')'p.

(4)

g'

The constant two-dimensional Bloch vector kll belongs to the first supercell Brillouin zone and is set by the frequency and orientation of the incident wave. As will be seen later, a strong formal simplification can be achieved if we separate the electric field components Egll(Z) parallel to the layer interface from the component E g ± ( Z ) normal to that layer. The same separation is applied to the magnetic field. In terms of these components, the wave equation takes the form of a simple one-dimensional system of combined differential and algebraic equations, a standard problem in numerical analysis:

dHgl'( z) = dz

z ) E . . , ( z) - (k,, + g ) . g,

dEglt(z) dz

= -~°tx°tryHglI(z) + i ( k l l + g ) E g ± ( z ) '

WSo E e g - g ' ( z ) Eg,± ( z ) = - i ( k l l + g ) O'yHgll ( z ) . g'

00tl,0

z)

(k, +g).

(5) (6) (7)

J.-P. Vigneron et al. / Ultramicroscopy 61 (1995) 21-27

23

Here, gg(z) are Fourier coefficients of the periodic function e(r), and % is the two-by-two matrix (which happens to be its own inverse):

Eqs. (5)-(7), combined with an adequate ordinary differential equation integrator, provide of a new algorithm which can be used to propagate the electromagnetic fields from the emitter region into the three-dimensionally structured layer. The possibility of using fast Fourier transform techniques to evaluate the convolution product in (5) and to solve the system (7) turns this procedure into an extremely fast Maxwell equation solver, while its unusually modest memory requirement (one or two z-layer retained at any time) facilitates its practical implementation.

2.2. Poynting vector and energy transfer The expressions found for the lateral electric and magnetic fields allow to determine the Poynting vector, which can readily be shown to be periodic, assuming the supercell periodicity. The relevant quantity for the evaluation of the energy transfer is here the surface integral J of the Poynting vector over the area of a two-dimensional superlattice unit cell. This amounts to O"

S = ~ }Re ~] Egll ( z ) × Hg]l ( z ) " e z .

(9)

g

Such an expression only involves the field projections actually used in Eqs. (5)-(7) and is a sound basis for a widely portable computational scheme for assessing the transmission and reflection coefficients in very general three-dimensional film structures.

2.3. Scattering boundary conditions and transfer matrix On a plane of constant z, the tangential components of the electric field Elf( p, z) is continuous, even if the dielectric function g( p, z) experiences a discontinuity. Since the normal direction across such a plane is the z direction, the tangential component of the electric field is a continuous function Of z.' B y using the uniqueness of the Fourier expansion, we can then show that all coefficients Ell(Z) are, individually continuous functions of z. For non-magnetic (/x = 1) materials, the tangential component of the magnetic field is also continuous across the interface and all coefficients Hgll(Z) are continuous functions of z. The continuity of the two basic field components used by the present algorithm at the entrance and exit of the scattering layer is an important simplification. In the emitter region I and in the collector region III, the dielectric function, el and gill are constant over space. In such a region of homogeneous response, the electric field can be represented as a superposition of traveling (or evanescen0 plane waves, with polarizations vectors normal to the propagation direction. For a wave with lateral wavevector kll + g, in region I, we define the polarization krl + g

Ik,,+ l

Xez

(lO)

as that of a local transverse electric wave (referring to the normal vector, or z direction) while the complementary polarization takes one of the directions

Ik,, + -

¢?7, oo/ c ez + x , '

(11)

J.-P. Vigneron et al. / Ultramicroscopy 61 (1995) 21-27

24

(in the case of a forward propagation along

ez), and

[k,l + g[

I.t{g

f~r w/c ez - Xig'

(12)

(in the case of a backward propagation along vectors is given in terms of the vector

kIg =

ez). The

tangential part of these complementary polarization

kll + g

(13)

o # c Igll + g['

where the z-component of the wavevector is chosen to be

,7-1k,,+gl =, ~e[klg]>0.

(14)

We have similar expressions for region III. Having defined the polarization reference frames, we can write the general form of the fields in the emitter region as a superposition of three-dimensional plane-wave beams (labelled by g vectors). For the electric field, we write Eiig ( Z ) = N~-g'lqg e ik!g(z-zI) + Ni-g~g e -ikIg(z-zI) -It-X~g Xlg eikI~(z-zO -- Xlg XIg e-ikIg(Z-Zl), (15) while the magnetic field, for the same amplitudes, expresses itself as

nll ( z) = ~ii [ N~ XIg eiklg(z-zt) -- Nig Xlg e-iklg(z-zI) -- X~g~l~g eikIg(Z-ZI) -- Xlg~C~g e-iklg(z-zl)] (16)

On the

other hand, in the collector homogeneous region (lII), the Fourier coefficients take the form

Eiig (z) = Uii+ig~g e ikmg(z-z"0 + NHig Tqg e -i~m~(z-zm) + X~Ig

XliIg

eikIt'g(Z-Zm) -- XlIIg

XIIIg

e-ikmg(Z-Zm) (17)

and Hiig ( z ) = V~III [ N~g X lng eikHI~(z-zm)- NilIg X IIIg e-ikm"(z-zm)/Z0¢ -- Xliig ~flg e - iknlg(Z-ZIu)] ,

XI+IIgrlg eikmg(Z-ZlH) (18)

Applying scattering boundary conditions amounts to find the amplitude of the reflected waves so that, combined with a single incident wave, the field in region I provides initial conditions to propagate across the junction layer and leaves in region III only outgoing waves (N~ig and X~I e both zero). Obtaining the correct coefficients is merely a matter of choosing the appropriate superposition of plane waves in both regions, with the constraint that they connect according to the propagation machinery (5)-(7). This problem of finding the relevant amplitudes can be given a matrix form in the following way. The set of coefficients /VI+, X?, NI, X~- and Nii+i, KiWi, NIII, XIII (the overbar denotes an entire column of components, each labelled using reciprocal lattice vectors) are representations of the wave in, respectively, the incidence and emergence regions. If these coefficients are provided, the wave can be reconstructed, using the explicit field expressions, given above. In particular, if the coefficients NI~I, XIII, -+ NIII, XIII are known, the fields at z = ZliI are given by E l l g ( Z l I I ) = NiISgaC~g+ NHig~f~g + Xi+ig XIIIg -- XIIIg X i i i g ,

Htlg( ZIII) =

t.~0c

[ niI+rgXIIIg

--

N~,,

XIIIe -

X~Ig ~q, - Xnlg~Ig]"

(19) (20)

J.-P. Vigneronet al./ Ultramicroscopy61 (1995)21-27

25

Conversely, on the interface z = ZI we can reconstruct the coefficients from the field values by inverting the above formulas, relocalized at the edge of region I. We find

I NO. = XIg

J

Ey.(z,) 7---1

(zI)

-, ( l ' L 0 C / ,~ ] g- ," g -

(21)

,

( [..LoC/~)gyg (Z,)

where T-1

ky + gy 2lkll +gl

o~

k~ + gx

2lkll +g[

c

k x + gx 2kigikll + gl

~

O)

ky + gy 2kig[kll +gl

c

~II °9

kx + gx

~ I (.0

ky + gy

ky + gy

kx + gx

c

2k[gikll +gi

c

2klglkii +gl

21kli+gi

2lkll+gi

ky + gy 21kll + g[

kx + gx

k x + gx

2lkil + gl

c

2k[glkLi + gl

ky + gy c

2k[glktl + gl

f~i o9

k, + g x

(~i o)

ky + gy

ky + gy

kx + gx

c

2krglkli+gl

c

2klglkll+gl

21kll+gl

2lkll +gl (22)

Now, because of the linearity of Maxwell's equations, as given by Eqs. (5)-(7), there exists a transfer matrix which connects the set of coefficients in regions I and III. This relation will be written as

F] N~+

[ M+;

M~;

M~N

M~x

[ NI-

/ MN;

M;]

MNN

M;X

XI

[ MX;

Mx;

MXN

Mxx

x;I1

(23)

Column m of the transfer matrix can be computed by assigning the value zero to all coefficients NI~I, X~I , ~/HI, Xm, except the one appearing in position m, reconstructing the fields at z = Zm using Eqs. (19) and (20), integrating system (5)-(7) to z = zi and rebuilding the coefficients N+, X~-, N~-, X~- using Eqs. (21) and (22). Once the transfer matrix is known (actually, only half of it must be constructed), we can implement the scattering boundary conditions by specifying that the "incident wave coefficients" N~+, X~ are assumed to be known, by providing the proper incident wave, and that there is no back-traveling wave (or growing evanescent wave) in the emergence (collector) region: Nffi= 0, :~i = 0. Once ported into Eq. (23), this provides a direct relationship between the coefficients of the transmitted waves and the coefficients of the incident waves. This relation takes the form of a system of linear equations M++

MNX ] [NI~-I] = [/~I+ ]

m;~

M ; + ] [ "~II J

L "~;

(24) '

26

J.-P. Vigneron et al. / Ultramicroscopy 61 (1995) 21-27

whose solution only requires the knowledge of a fraction of the transfer matrix. The energy transferred can be obtained directly from the computed coefficients, using

J=

~(INI,+~=I2+Ix~,=I 2 2/~0w g

=,1i - I k

+gl2o = i i i T r - l k + g l

2

(25)

The Heaviside O-function (1 for positive arguments, 0 otherwise) insures that only radiative emerging waves will transfer energy.

3. Discussion Most methods used today to solve Maxwell's equations [5,6] are oriented towards a real-space representation of the scattered fields, and in the case of three-dimensional problems, require m handle an extremely large amount of data. For sake of efficiency, these should be kept in fast access storage. Memory size demand puts a stringent limit on the applicability of any numerical approach, as the grid size is not only related to the size of the scattering centers, but also to the wavelength of the scattered radiation. The representation developed in the present work is taken for one part in real space and, for the other part, in momentum space. This combination (the so-called Laue representation) allows to operate with an extremely fine (eventually adaptive) mesh in the transmission direction. The algorithm presented here has the advantage of treating all z-layers one at a time, and that it is never necessary to retain the information about a large number of these layers, as the computation " s w e e p s " through the scattering layer. This means that the memory requirement for this algorithm is drastically shortened, and that the limiting constraints on the scattering medium complexity is significantly relaxed. With this approach, the lateral coordinates are treated with Fourier series, and all costly operations can take advantage of the efficiency of fast Fourier transform techniques. The basic algorithm needed here is built on the fact that a convolution product can be carried out by first transforming both terms in direct space, make the direct local product, and return to Fourier space. This leads to a procedure which roughly scales as the number of reciprocal lattice vectors, to be executed a number of times equal to the number of points representing the scattered wave within one supercell in the scattering layer. A large amount of flexibility is needed to represent scattered waves well enough so that the energy flux is accurately computed. The fundamental reason for this is that the energy flux is not a variational quantity, i.e. a quantity which assumes an extremal value for the correct representation of the scattered field. The energy flux is then particularly exposed to errors made in the representation of the waves. Globally, the multiple scattered field can be described here with an extremely large number of coordinate points. It is not unreasonable, with this procedure, to analyze the field with more than a thousand reciprocal space vectors, and ten thousand mesh points in the transmission direction. This provides the scattered field with ten million mesh points in space, with the advantage that this total amount of data is never required to be present instantaneously in any near or far computer memory. A more traditional direct real-space approach which relaxes the three-dimensional fields with those ten million grid points would require to treat a matrix system with 1014 elements. Even if the matrix is expected to be substantially sparse, this is out of reach of most traditional computing systems available today.

Acknowledgements This work was funded partly by the Interuniversity Research Project on "Sciences of interfacial and mesoscopic structures" ( P A I / I U A P N. P3-49) initiated by the Belgian State Policy Programming. A. Castiaux is supported by the National Science Foundation (FNRS). This work was, for another part, supported by the Human Capital and Mobility Program initiated by the Commission of the European Community (Contract No. ERBCHRXCT930130).

J.-P. Vigneron et al. / Ultramicroscopy 61 (1995) 21-27

27

References [1] [2] [3] [4]

C. Girard, Theoretical analysis of scanning near-field optical microscopy, Scanning 16 (1994) 333, and references therein. D.W. Pohl, Scanning near-field optical microscopy, Adv. Opt. Electron Microsc. 12 (1991) 243, and references therein. J.B. Pendry, Low Energy Electron Diffraction (Academic Press, London, 1974), and references therein. J.P. Vigneron, I. Derycke, Th. Laloyaux, Ph. Lambin and A.A. Lucas, Computation of scanning tunneling microscope images of nanometer-sized objects physisorbed on metal surfaces, Scanning Microsc. Suppl. 7 (1993) 261. [5] A. Dereux and D. Pohl, The 90° prism as a model SNOM probe, in: Proc. NATO Adv. Res. Workshop on Near-Field Optics, Ed. D.W. Pohl (Kluwer, Dordrecht, 1993) p. 189. [6] O. Martin, A Numerical Green's Function Approach to Investigate Vectorial Field-Matter Interactions, Th~se No. 1199, Ecole Polytechnique F6d6rale de Lausanne, 1993.