論文 / 著書情報
Article / Book Information
Title

Weakly guiding approximation of a three dimensional waveguide modelfor extreme ultraviolet lithography simulation

Authors

Hiroyoshi Tanabe, Akira Jinguji, Atsushi Takahashi

Citation

Journal of the Optical Society of America A, Vol. 41, Issue 8, pp. 14911499

Pub. date

2024, 7

Copyright(c) 2024 Optica Publishing Group. Users may use, reuse, and build
upon the article, or use the article for text or data mining, so long as
such uses are for non-commercial purposes and appropriate attribution
is maintained. All other rights are reserved.

DOI

http://dx.doi.org/10.1364/JOSAA.516610

Powered by T2R2 (Tokyo Institute Research Repository)

Research Article

Vol. 41, No. 8 / August 2024 / Journal of the Optical Society of America A

1491

Weakly guiding approximation of a threedimensional waveguide model for extreme
ultraviolet lithography simulation
Hiroyoshi Tanabe,* Akira Jinguji, AND Atsushi Takahashi
Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8550, Japan*tanabe.h.af@m.titech.ac.jp
Received 20 December 2023; revised 27 May 2024; accepted 11 June 2024; posted 12 June 2024; published 8 July 2024

A three-dimensional (3D) waveguide model is applied in extreme ultraviolet (EUV) lithography simulations. The
3D waveguide model is equivalent to rigorous coupled-wave analysis, but fewer field components are used to solve
Maxwell’s equations. The 3D waveguide model uses two components of vector potential, A x and A y , corresponding to the two polarizations. The electric field of the A x polarization is approximately parallel to the x axis, and theelectric field of the A y polarization is approximately parallel to the y axis. The 3D waveguide model solves a coupled vector wave equation for two polarizations. The refractive index of conventional EUV absorbers is close to thatof vacuum. The weakly guiding approximation in optical fiber theory is applied to the 3D waveguide model. The
coupled vector wave equations for the two polarizations are decoupled into two independent scalar wave equations.
Maxwell’s equations are simplified to a set of scalar wave equations. The weakly guiding approximation reduces thecomputation time to solve the equations. The computation time required to solve the weakly guiding approximation is about 1/5 of the time to solve the original 3D waveguide model. © 2024 Optica Publishing Group under the termsof the Optica Open Access Publishing Agreement
https://doi.org/10.1364/JOSAA.516610

1. INTRODUCTION
Extreme ultraviolet (EUV) lithography has begun to be used in
the mass production of logic and memory devices. The typicalresolution for EUV lithography is 16 nm assuming a numericalaperture (NA) of 0.33 and a moderate k1 factor of 0.4. In thiscase, the minimum size of the mask pattern will be 64 nm. The
typical thickness of the Ta absorber is about 60 nm, and the
aspect ratio of the mask pattern is nearly 1. High-aspect absorbers give rise to several mask three-dimensional (3D) effects,
such as critical dimension differences between horizontal andvertical lines and focus-dependent pattern shift [1,2]. Accurate
electromagnetic (EM) simulation is required to reproduce thesemask 3D effects.
There are two types in EM simulation methods. The finitedifference time-domain (FDTD) method [3] solves Maxwell’sequations in coordinate space, and the solution is the near-field
diffraction amplitude. Since the near-field diffraction amplitudedepends on the incident angle, the calculation must be repeated
for all source points (usually 100 or more).
On the other hand, rigorous coupled-wave analysis (RCWA)[4–7] and the 3D waveguide model [8,9] solve Maxwell’s equation in momentum (or frequency) space. Assuming a periodic
boundary condition, the momentum is discretized proportional to the inverse of the period length. These models solve a1084-7529/24/081491-09 Journal © 2024 Optica Publishing Group

coupled-wave equation that connects the incoming and outgoing waves. All relationships between the incoming and outgoingwaves are calculated simultaneously. The momentum of theincoming wave corresponds to the source position. If the period
length is large enough, the size of the momentum grid is small,and these models do not need to repeat calculations for different
source points.
The formulations of RCWA and the 3D waveguide model areequivalent as shown in Appendix A of this paper. The difference
between the two approaches is the number of the field components used in solving Maxwell’s equations. The RCWA method
uses four components of electric and magnetic fields: E x , E y ,Hx and Hy . On the other hand, the 3D waveguide model uses
two field components: vector potentials A x and A y .
Maxwell’s equations contain six field components, the electric field E , and the magnetic field H. However, photons only
have two polarizations. The 3D waveguide model uses two fieldcomponents A x and A y corresponding to the two polarizations.
The first part of this report provides a step-by-step explanationof how six field components of Maxwell’s equations are reduced
to two field components. First, the electric field and magnetic
field are replaced by four field components: the scalar potential
and vector potential. Gauge transformation eliminates the two
field components: the scalar potential and the z component of1492

Research Article

Vol. 41, No. 8 / August 2024 / Journal of the Optical Society of America A

the vector potential. Finally, Maxwell’s equations are reduced toa coupled-wave equation for A x and A y .
For EUV lithography simulations, the 3D waveguide model
can be further simplified [10–13]. In optical fiber theory, theweakly guiding approximation [14] decomposes all polarization
modes into uncoupled linearly polarized modes. This approximation is possible because the refractive index of the core isclose to that of the cladding. Similarly, the refractive index of the
Ta absorber used in EUV masks is close to that of the vacuum.
In the second part of this report, we apply the weakly guiding
approximation to the 3D waveguide model. We decompose
the coupled vector wave equation of two polarizations into two
independent scalar wave equations. Maxwell’s equations are
simplified to a set of scalar wave equations.
In Section 2 we explain the 3D waveguide model. In Section 3
we apply the 3D waveguide model to EUV lithography simulation. In Section 4 we verify the validity of the weakly guiding
approximation. Section 5 is the summary.
2. 3D WAVEGUIDE MODEL
This section explains the basics of the 3D waveguide model.

∇ × E = ik H,(7)∇ × H = −ik ε̂ E ,(8)
∇ · ε̂ E = 0,

(9)

∇ · H = 0,(10)where ε̂ is the complex dielectric constant defined byε̂ = ε + i

4π σ
.
ω

(11)

Maxwell’s equations have six field components E and H. We
reduce the number to four by using the vector potential A andthe scalar potential ϕ. For an explanation of vector and scalarpotentials and their gauge transformation, see Secs. 6.4-5 inRef. [15].
From Eq. (10) we can define the magnetic field in terms of avector potential:
H = ∇ × A.

(12)

Substituting this equation into Eq. (7), we obtain

A. Wave Equation of Vector Potential in DielectricMaterials

In this subsection we derive wave equations in dielectric materials. The dielectric constant ε varies with location. We assume
that the materials are nonmagnetic and the magnetic permeability µ = 1. This assumption applies to the materials used in
EUV masks. Then, Maxwell’s equations are written as follows(Eq. (6.28) in Ref. [15]):∇ × (E − ik A) = 0.

(13)

Then, the electric field E is written by using the scalar potential ϕ as follows:
E = ik A − ∇ϕ.

(14)

(1)

Equations (8) and (9) are not independent from each other.
If E is the solution of Eq. (8), it automatically satisfies Eq. (9).
Substituting Eqs. (12) and (14) into Eq. (8), we obtain thefollowing equation:(2)

1 A + k 2 ε̂ A − ∇ (∇ · A) + ik ε̂∇ϕ = 0.

∇ · (ε E ) = 4πρ,(3)

It is well known that vector potential A and scalar potential
ϕ have gauge transformation freedom. If we change the vector
potential A and scalar potential ϕ simultaneously as follows,∇ · H = 0,

(4)

A → A − ∇3,(16)

where c is the speed of the light, and ρ and J are the charge density and the current density, respectively. The charge density andthe current density obey the following continuity equation:

ϕ → ϕ − ik3,

(17)

1 ∂H
,
∇ × E =−c ∂t∇×H=

1 ∂(ε E )
4π
J +
,
c
c ∂t

∂ρ
+ ∇ · J = 0.
∂t

(5)

The Ohm’s law gives the relationship between the electricfield and the current density:
J = σ E,

(6)

where σ is the conductivity.
We define the plane wave by exp(i(k · x − ωt)), where k is
the wave vector and ω is the angular velocity. The wave number
k is related to the angular velocity by k = ω/c . In steady state forincoming and outgoing waves, E , H, J , and ρ have the sametime dependence exp(−iωt). Then, Maxwell’s equations are
simplified as follows:(15)

these changes do not affect the observables E and H, where 3
is an arbitrary gauge function. This freedom allows applying an
additional constraint on A and ϕ,∇ · A − ik ε̂ϕ = 0,(18)which is called the Lorenz condition.
Inserting the above equation into Eq. (15), we obtain theHelmholtz type wave equation [16] of the vector potential A as
follows:
1 A + k 2 ε̂ A − ∇(logε̂)∇ · A = 0.

(19)

The Lorenz condition eliminates the scalar potential fromthe equation. The number of the field components is reduced
to three. When the complex dielectric constant ε̂ is not spatiallyResearch Article

Vol. 41, No. 8 / August 2024 / Journal of the Optical Society of America A

A z = 0.

Fig. 1.

uniform, the third term in this equation does not allow solving
A x , A y , and A z independently.
B. Coupled-Wave Equations Inside a Layer

In this subsection, we simplify the wave equations, Eq. (17),inside a layer where the complex dielectric constant ε̂ is uniform
in the z direction. We still assume that ε̂ depends on x and y .
The structure of EUV masks satisfies these assumptions (Fig. 1).
For simplicity we assume that the absorber is a single layer. Bothinside the absorber layer and Mo/Si multilayers, the complex
dielectric constant ε̂ is uniform in the z direction. Inside eachlayer the z component of the wave equation is written as follows:(20)

The Lorenz condition, Eq. (18), partially restricts gaugetransformation freedom, but not all of them. If the gauge
function 3 satisfies the following equation,
13 + k 2 ε̂3 = 0,(27)

Maxwell’s equations are reduced to coupled-wave equationsof the two components of the vector potential A x and A y . The
electric field and the magnetic field are calculated from thevector potential as follows:i1
E = ik A + ∇ ∇ · A ,
(28)
k
ε̂

Cross section of an EUV mask.

1A z + k 2 ε̂ A z = 0.

1493

(21)

then Lorenz condition is automatically fulfilled. Therefore, anarbitrary function 3 that satisfies Eq. (21) can be added to A
and ϕ. We use this gauge transformation freedom to eliminate
A z in Eq. (19). We define the gauge function 3 as follows [8]:Z z
A z (x , y , z 0 )dz 0 + f (x , y ),(22)3(x , y , z) =H = ∇ × A.

(29)

Equations (25) and (26) are further simplified in the case of
one-dimensional mask patterns. These equations are decoupled
to transverse electric (TE) and transverse magnetic (TM) waves
as shown in Appendix B.
Inside the Mo/Si multilayers, ε̂ is uniform in the x and y
directions, so Eqs. (25) and (26) can be solved analytically.
The computation time is dominated by the time to solve thecoupled-wave equation in the absorber layer. In Fig. 1 we assume
that the absorber is a single layer. If the absorber has multiple
layers, the computation time becomes longer. If the sidewallsof the absorber are sloped, the staircase approximation can beapplied, but the computation time will be very long.
C. Boundary Conditions

At the boundary between the j th and ( j + 1)th layers(z = z j +1 ), the following boundary conditions must be
imposed on the electric and magnetic fields in the x − y plane:
E xj +1 x , y , z j +1 = E xj x , y , z j +1 ,(30)
E yj +1 x , y , z j +1 = E yj x , y , z j +1 ,

(31)



Hxj +1 x , y , z j +1 = Hxj x , y , z j +1 ,

(32)



Hyj +1 x , y , z j +1 = Hyj x , y , z j +1 .

(33)

z0

where z0 is an arbitrary fixed point inside the layer, and f (x , y )is the solution of the following differential equation:
 2

∂2
∂ Az
∂
+
f (x , y ) + k 2 ε̂(x , y ) f (x , y ) = −.
∂x2 ∂ y2∂z z=z0(23)
It can be easily proved that the gauge function 3 satisfies
Eq. (21). With this gauge transformation, the z component of
the vector potential can be eliminated because, from Eq. (22),∂3= Az .
∂z(24)

After the gauge transformation, Eq. (19) is simplified asfollows:∂ Ay∂logε̂ ∂ A x
2
1A x + k ε̂ A x −
+
= 0,
(25)
∂x
∂x
∂y


∂ Ay
∂logε̂ ∂ A x
1A y + k 2 ε̂ A y −
+
= 0,
(26)
∂y
∂x
∂y

We calculate the electric field and the magnetic field withineach layer by solving the coupled-wave equations, Eqs. (25) and(26). These equations contain two field components, A x andA y . Four boundary conditions, Eqs. (30)–(33), are required to+
fix the waves going forward A +x and A y , and the waves going−−backward A x and A y .
D. Incoming and Outgoing Waves in Vacuum

EUV masks are placed inside vacuum where the complex dielectric constant ε̂ is 1. The solution of the wave equation, Eq. (19),is a plane wave A(x ) = A 0 exp(ik · x ). From Eqs. (28) and(29), the electric field E and the magnetic field H of the incoming and outgoing waves in vacuum can be described using thevector potential A as follows:
E = ik A −i(k · A) k,k

H = ik × A.

(34)
(35)

1494

The wave vector, electric field, and magnetic field are perpendicular to each other:
E · k = H · k = E · H = 0.

(36)

The magnitudes of the electric field and magnetic field are thesame:|E |2 = |H|2 = k 2 | A|2 − |k · A|2 .

(37)

Gauge transformation freedom allows us to fix A z to be zero
if the media is uniform in the z direction. This gauge condition
can also be applied in vacuum. Therefore, two components ofthe vector potential, A x and A y , determine the electromagnetic
field of the light. We call these two components A x and A ypolarizations. The vector potential, electric field, and magnetic
field of the A x and A y polarizations are explicitly written as
follows.
A x polarization:
A = (A x , 0, 0) ,
E=

(38)


iA x 2k − k x2 , −k y k x , −kz k x ,
k

(39)


H = iA x 0, kz , −k y .

(40)A = 0, A y , 0 ,(41)

A y polarization:

E=


iA y
−k x k y , k 2 − k 2y , −kz k y ,k(42)

H = iA y (−kz , 0, k x ) .

(43)

Figure 2 shows schematic views of A x and A y polarizations.
When the wave vector is in the x − z plane, k y = 0. In thiscase, A x polarization corresponds to TM polarization and A ypolarization corresponds to TE polarization. In the same way,
when the wave vector is in the y − z plane, A x polarization corresponds to TE polarization and A y polarization corresponds to
TM polarization.
When k x , k y  k, E ∼ (ik A x , 0, 0) for A x polarizationand E ∼ (0, ik A y , 0) for A y polarization. The electric field forA x polarization is approximately parallel to the x axis, and theelectric field for A y polarization is approximately parallel to they axis.

Fig. 2.

Research Article

Vol. 41, No. 8 / August 2024 / Journal of the Optical Society of America A

Schematic views of A x and A y polarizations.

3. APPLICATION TO EUV LITHOGRAPHYSIMULATIONIn this section we apply 3D waveguide model to EUV lithography simulation.
A. Image Intensity on a Wafer

EUV scanners use Koehler illumination and partial-coherent
optics [17] like microscopes. EUV wavelength λ = 13.5 nm.
The magnification of the current EUV scanner is 1/4, and theNA is 0.33. According to Abbe’s theory [17], the total image
intensity on a wafer, I , is calculated by the incoherent sum of theimage intensities illuminated by all points of the effective source
S as follows:
ZZ
ZZ
E (k out ; k in )P ( p x , p y )I (x , y ) =S(s x , s y )× e i( p x x + p y y ) d p x d p y

2

ds x ds y ,(44)

where P is the pupil function of the projection optics. The
electric field E diffracted by an EUV mask is calculated bysolving the 3D waveguide model in Section 2. It depends on
the incoming wavevector k in and the outgoing wavevector k out .
These are calculated using the source position of the incoming
wave (s x , s y ) and the pupil position of the outgoing wave ( p x ,p y ) as follows:qk in = s x , −k sin θ + s y , k 2 − s x2 − (−k sin θ + s y )2 ,
(45)


q
k out = p x , −k sin θ + p y , − k 2 − p x2 − (−k sin θ + p y )2 .
(46)
The optics of EUV scanners is reflective, and the chief ray istilted θ = 6 deg in the negative y direction.
We calculate the intensity of the electric field because thephotoresist chemically reacts with the electric field. The pupil
position is limited by the radius of the pupil defined by NA. The
pupil function may include aberrations such as defocus.
B. Polarization Dependence of the Image Intensity

The image intensity depends on the polarization of the incoming wave. We calculate the image intensity of a mask (Fig. 3)that mimics the imec N3 metal layer [18]. The line widthis 56 nm on the mask and 14 nm on the wafer. The simulation uses the following optical settings: λ 13.5 nm, NA 0.33,
and dipole illumination σout /σin = 0.9/0.55 with the opening angle 90 deg. We assume a conventional Ta absorber of
60 nm thickness and a periodic boundary condition withthe period length L = 2048 nm. The size of the mask area is2048 nm × 2048 nm. The computational time to solve the3D waveguide model in this area is 2850 s using Intel Core
i9-10920X CPU and Nvidia RTX3090 GPU. We use the
MAGMA library [19] to speed up the eigenvalue calculation in
Eq. (A16). The calculation time highly depends on the cutoff

Research Article

Fig. 3.

Vol. 41, No. 8 / August 2024 / Journal of the Optical Society of America A

1495

Mask layout used in the simulation.

Fig. 5. Influence of the polarization change between the incoming
wave and outgoing wave.

4. WEAKLY GUIDING APPROXIMATION

Fig. 4. Image intensities of incoming waves with A x and A ypolarizations.

momentum. According to Ref. [11], we include the momentum(k x , k y ), which satisfies!
ky|k x |+1+ 1 ≤ 2,(47)
k xmaxk max
ywhere k xmax = k max= 6 · NA/4 · 2π/λ. This number is six
ytimes larger than the size of the pupil NA/4 · 2π/λ. When the
momentum is discretized by 2π/L, there are 8705 (k x , k y )pairs that satisfy Eq. (47). The size of the matrix solving thecoupled-wave equations is 17,410 × 17,410 because there are
two polarizations.
The region in Eq. (47) is a quasi-hyperbola that mimics thediffraction spectrum of a mask pattern consisting of verticaland horizontal lines or holes. Mask patterns are conventionallydesigned using X–Y coordinates. The minimum pattern pitch in
the X or Y direction is small compared to the minimum pattern
pitch in the diagonal direction. Therefore, in momentum space,the diffraction amplitude in the diagonal direction decreases
rapidly compared to the amplitude in the X or Y direction.
Figure 4 compares the image intensities of incoming waves
with A x and A y polarizations. There is a small difference
between the two intensities, at most 1.6%. The number is anabsolute value relative to the ML reflectance of 0.64. This isbecause the incident angle is tilted 6 deg in the y direction.
In the case of the unpolarized light, the image intensity canbe calculated averaging the image intensities of two polarizedincoming waves.
Even when the incoming wave has A x (or A y ) polarization,
the outgoing wave has both A x and A y polarizations. Figure 5
shows the influence of the polarization change between theincoming wave and outgoing wave. At the middle of Fig. 5,
we artificially remove the outgoing wave, which has a different
polarization from that of the incoming wave. The contributionof the polarization change is very small, less than 0.01%.

In optical fiber theory, the weakly guiding approximation [14]is often used because the refractive index of the core is close tothat of the cladding. The weakly guiding approximation ignores
the small cross term that appears in the coupled-wave equationof the transverse electric field [20]. The cross term includes thederivative of the dielectric constant, which is proportional tothe difference in refractive index between the core and cladding.
In this way, the weakly guiding approximation decomposes thecoupled vector wave equation of the electromagnetic field into
uncoupled scalar wave equations.
The complex
refractive index of the Ta absorber is√(n, k) = ε̂ = (0.9567, 0.0343) [21], which is close to therefractive index of vacuum, 1. The weakly guiding approximation can be applied to EUV masks when the Ta absorber is used.
This is the reason that the polarization change from the incoming wave to the outgoing wave is small. According to the weakly
guiding approximation, the coupled vector wave equations of
two polarizations, Eqs. (25) and (26), are decomposed into thetwo scalar wave equations as follows:
1A x + k 2 ε̂ A x −∂logε̂ ∂ A x
= 0,
∂x ∂x

(48)

1A y + k 2 ε̂ A y −∂logε̂ ∂ A y
= 0.
∂y ∂y

(49)

The third terms in these equations are small and could be
removed, but we retain these terms because they do not affect tothe speed of the computation.
Each equation requires two boundary conditions for thewaves going forward and backward. Since the electric fieldof the A x polarization is approximately parallel to the x axis,
we impose the following two boundary conditions on the A xpolarization:
E xj +1 x , y , z j +1 = E xj x , y , z j +1 ,(50)
Hyj +1 x , y , z j +1 = Hyj x , y , z j +1 .

(51)

In the same way, we impose the following two boundary conditions on the A y polarization:

1496

Fig. 6. Comparison of 3D waveguide model and weakly guiding
approximation. Ta absorber with 60 nm thickness is used.

E yj +1 x , y , z j +1 = E yj x , y , z j +1 ,

(52)



Hxj +1 x , y , z j +1 = Hxj x , y , z j +1 .

(53)

Fig. 7. Comparison of the 3D waveguide model and the weakly
guiding approximation. A low-n absorber, TP1, with 45 nm thickness
is used.

5. SUMMARY

Each equation can be solved independently, and it takes
289 s for a mask area of 2048 nm × 2048 nm. It takes 578 s tosolve two equations, which is about 1/5 of the time to solve theoriginal 3D waveguide model.
Figure 6 compares the image intensities calculated by the 3D
waveguide model and the weakly guiding approximation. The
difference is very small, less than 0.1%. This result is consistent
with the results in Fig. 5.
Table 1 summarizes the maximum intensity errors for several
test mask patterns, including vertical and horizontal L/S andhole patterns. The illumination used for L/S patterns is dipole
illumination, σout /σin = 0.9/0.55, with the opening angle
90 deg. The direction of the dipole is perpendicular to the L/Spattern. The illumination used for the hole pattern is annular
illumination, σout /σin = 0.9/0.55. In all cases the maximum
intensity error is less than 0.1%. The accuracy of the weakly
guiding approximation is very high when the Ta absorber isused.
Figure 7 shows the results using the low-n absorber TP1 in
Ref. [22]. The complex refractive index of the TP1 absorber is(0.91, 0.032). As shown in Fig. 6, the difference between the
3D waveguide model and the weakly guiding approximation
increases to 1.4%. The accuracy of the weakly guiding approximation is deteriorated when the low-n absorber is used. Notethat the low-n absorber is still under development and the mask
process has not yet been established, as discussed in Ref. [22].
Table 1.

Research Article

Vol. 41, No. 8 / August 2024 / Journal of the Optical Society of America A

We apply the 3D waveguide model to EUV lithography simulation. The 3D waveguide model is equivalent to RCWA. The 3D
waveguide model uses two field components of vector potential
A x and A y . The two field components A x and A y correspond to
the two polarizations. The electric field of the A x polarization isapproximately parallel to the x axis, and the electric field of theA y polarization is approximately parallel to the y axis.
The 3D waveguide model solves a coupled vector wave
equation for two polarizations. We apply the weakly guiding
approximation to the 3D waveguide model. We decompose
the vector wave equation of the two polarizations into two
independent scalar wave equations. Maxwell’s equations are
simplified to a set of scalar wave equations. The weakly guiding approximation reduces the computation time to solve theequations. The computation time to solve the weakly guiding
approximation is about 1/5 of the time to solve the original 3D
waveguide model.
APPENDIX A: EQUIVALENCE BETWEEN 3D
WAVEGUIDE MODEL AND RCWA
In this appendix we verify that the 3D waveguide model is equivalent to RCWA.
1. RCWA Formulation

According to the RCWA formulation in Ref. [7], Maxwell’sequations are reduced to the following coupled first orderdifferential equations of E x , E y , Hx and Hy :Maximum Image Intensity Errors between the 3D Waveguide Model and Weakly Guiding Approximation

Incoming
Wave

V L/S
L 14 nmP 28 nm

V L/S
L 14 nm
P 42 nm

V L/S
L 20 nm
P 40 nm

H L/S
L 14 nmP 28 nm

H L/S
L 14 nm
P 42 nm

H L/S
L 20 nm
P 40 nmHole
H 20 nm
P 40 nmHole
H 20 nm
P 80 nm

Ax
Ay

0.05%
0.09%

0.01%
0.05%

0.02%
0.05%

0.004%
0.07%

0.01%
0.1%

0.008%
0.09%

0.01%
0.03%

0.008%
0.02%

Research Article

Vol. 41, No. 8 / August 2024 / Journal of the Optical Society of America A 
0∂ S y /∂(z 0 )00 ∂ S x /∂(z )   ∂U /∂(z 0 )  = K̂ x K̂ yy∂Ux /∂(z 0 )K̂ 2 − Êx00Ê − K̂ y2− K̂ x K̂ ywhere S and U are the Fourier components of the electric andmagnetic fields defined byX
E x (x , y , z) =S xmn (z) exp i k xm x + k ny y ,
(A2)
m,n

X

E y (x , y , z) =S ymn (z) exp i k xm x + k ny y ,

(A3)



Uxmn (z) exp i k xm x + k ny y ,

(A4)



U ymn (z) exp i k xm x + k ny y .

(A5)

m,n

Hx (x , y , z) = iXK̂ y Ê −1 K̂ x
Iˆ − K̂ y Ê −1 K̂ y  S y K̂ x Ê −1 K̂ x − Iˆ − K̂ x Ê −1 K̂ y 
Sx   Uy  ,
0 0
Ux
0 0

1497

(A1)

The above equations can be solved in terms of the eigenvalues
and eigenvectors of the matrices. We do not need to solve both
equations because S and U are not independent fields. For
example, if we solve Eq. (A13) for the electric field S, then themagnetic field U can be calculated by using Eq. (A1) as follows:
−1 

  
∂ S y /∂(z 0 )K̂ y Ê −1 K̂ x Iˆ − K̂ y Ê −1 K̂ yUy.
=
∂ S x /∂(z 0 )Ux
K̂ x Ê −1 K̂ x − Iˆ − K̂ x Ê −1 K̂ y
(A15)

m,n

Hy (x , y , z) = i

Xm,nThe wave vectors for the mth component in the x direction
and the nth component in the y direction are written as follows:
2π
,
3x

(A6)

2π
,
3y

(A7)

k xm = k x + m
k ny = k y + n

where 3x and 3 y are the periodic lengths in the x and y directions, respectively.
In Eq. (A1) the normalized z coordinate z 0 = kz, and thematrices Iˆ, K̂ x , K̂ y , and Ê are dimensionless values as follows:( Iˆ)mn,m 0 n 0 = δmm 0 δnn 0 ,( K̂ x )mn,m 0 n 0 = δmm 0 δnn 0( K̂ y )mn,m 0 n 0 = δmm 0 δnn 0

(A8)

k xm
,
k
k ny
k

( Ê )mn,m 0 n 0 = ε̂m−m 0 n−n 0 ,

,

(A9)

(A10)
(A11)

where the Fourier component of the dielectric constant isdefined by
 

X
2π x
2π y
+n
.
(A12)
ε̂(x , y ) =ε̂mn exp i m
3x
3x
m,n
From Eq. (A1) we can derive two sets of the second order differential equations as follows:

2. 3D Waveguide ModelIn the 3D waveguide model, Eqs. (25) and (26) can be rewritten
in the form of an eigenvalue equation as follows: 2∂ F x /∂(z 0 )2
∂ 2 F y /∂(z 0 )2"
# K̂ x2 + K̂ y2 − Ê − Ẑ x K̂ x
− Ẑ x K̂ y
Fx
=
,
Fy
− Ẑ y K̂ x
K̂ x2 + K̂ y2 − Ê − Ẑ y K̂ y
(A16)where F is the Fourier component of the vector potential A:
A x (x , y , z) =
i X mnF x (z) exp i k xm x + k ny y , (A17)k m,nA y (x , y , z) =
i X mnF (z) exp i k xm x + k ny y . (A18)k m,n yIn Eq. (A16) two new matrices Ẑ x and Ẑ y appear. They aredefined by
2π(m − m 0 )(ζ )m−m 0 n−n 0 ,k3x2π n − n 0( Ẑ y )mn,m 0 n 0 =(ζ )m−m 0 n−n 0 ,k3 y( Ẑ x )mn,m 0 n 0 =

(A19)

(A20)

where ζ is the Fourier component of logε̂ as follows:
 

X
2π x
2π y
log ε̂(x , y ) =(ζ )mn exp i m
+n
.
3x
3x
m,n
(A21)


 # 
 2
 " 2
K̂ x + K̂ y Ê −1 K̂ y − Iˆ Ê K̂ y Ê −1 K̂ x Ê − K̂ x
∂ S y /∂(z 0 )2
Sy


=
,
Sx
∂ 2 S x /∂(z 0 )2
K̂ x Ê −1 K̂ y Ê − K̂ y
K̂ y2 + K̂ x Ê −1 K̂ x − Iˆ Ê

(A13)

# 


 2
 " 2
K̂ y + Ê K̂ x Ê −1 K̂ x − IˆK̂ x − Ê K̂ x Ê −1 K̂ y
∂ U y /∂(z 0 )2
Uy


=
.
Ux
∂ 2 Ux /∂(z 0 )2K̂ y − Ê K̂ y Ê −1 K̂ x K̂ x2 + Ê K̂ y Ê −1 K̂ y − Iˆ

(A14)

1498

Research Article

Vol. 41, No. 8 / August 2024 / Journal of the Optical Society of America A

The electric field is calculated by Eq. (28), and this equation isexpressed in the matrix form as
  
 
Sx
Iˆ − K̂ x Ê −1 K̂ x − K̂ x Ê −1 K̂ y
Fx=
. (A22)Sy
Fy− K̂ y Ê −1 K̂ x Iˆ − K̂ y Ê −1 K̂ yEquations (B2)–(B4) are rewritten as follows:∂2 E y∂2 E y++ k 2 ε̂E y = 0,
2
∂x
∂z 2

(B6)


E = 0, E y , 0 ,(B7)

The magnetic field can be simply calculated from Eq. (29).
Ux = ∂ F y /∂z 0 ,(A23)

U y = −∂ F x /∂z 0 .

(A24)

We can verify that the 3D waveguide model is equivalent to
RCWA. Differentiating Eq. (A16) by z 0 and using Eqs. (A23)and (A24) we get 2∂ U y /∂(z 0 )2
∂ 2 Ux /∂(z 0 )2"
# K̂ x2 + K̂ y2 − Ê − Ẑ x K̂ x
Ẑ x K̂ y
Uy=.
Ux
Ẑ y K̂ x
K̂ x2 + K̂ y2 − Ê − Ẑ y K̂ y

i ∂Ey
i ∂Ey
H=
, 0, −k ∂z
k ∂x
Replacing the field from A x to Hy ,

(A26)

(A27)

∂ Ay
∂ Ay
, 0,
H= −
∂z
∂x


(B2)

(B3)

.

(B4)

Replacing the field from A y to E y ,E y = ik A y .

∂ 2 Hy∂ 2 Hy∂logε̂ ∂ Hy+
+ k 2 ε̂ Hy −
= 0,
2
2
∂x
∂z
∂x ∂x


i ∂ Hyi ∂ HyE= −, 0,,k ε̂ ∂zk ε̂ ∂ x

(B12)

(B13)
(B14)

Funding. Ministry of Education, Culture, Sports, Science and Technology(MEXT) Initiative to Establish Next-generation Novel Integrated Circuits
Centers (X-NICS) (Grant JPJ011438).
Disclosures. The authors declare no conflicts of interest.

(B1)

The above two equations can be solved independently. The
general solution is given by a combination of the following two
special solutions:I) A x = 0, A y 6 = 0, A z = 0.
From Eqs. (28) and (29) the electric and magnetic fields are
E = (0, ik A y , 0),(B11)

The solution is known to be a TM wave.

If the mask pattern is uniform in the y direction and the incoming wave is in the x − z plane, ε̂, A x , and A y do not depend on
the y coordinate. In this case, Eqs. (25) and (26) are further simplified as follows:

∂2 Ay
∂2 Ay
+
+ k 2 ε̂ A y = 0.
2
∂x
∂z 2

∂ Ax
.
∂z


H = 0, Hy , 0 .

APPENDIX B: ONE-DIMENSIONAL MASKPATTERN∂ 2 Ax∂ 2 Ax∂logε̂ ∂ A x++ k 2 ε̂ A x −
= 0,
2
∂x
∂z 2
∂x ∂x

Hy =

Equations (B1), (B9), and (B10) are rewritten as follows:

and its matrix representation isÊ K̂ x Ê −1 K̂ x = K̂ x2 − Ẑ x K̂ x .

(B8)

.

The solution is known to be a TE wave:II) A x 6 = 0, A y = 0, A z = 0.
From Eqs. (28) and (29) the electric and magnetic fields arei ∂ 2 Axi ∂ 2 Ax
E= −
, 0,
,
(B9)
k ε̂ ∂z 2k ε̂ ∂z∂ x


∂ Ax
,0 .
(B10)
H = 0,∂z(A25)
This equation is identical to Eq. (A14) because


∂ 1 ∂
∂logε̂ ∂
∂2
ε̂
,
= 2−
∂ x ε̂ ∂ x
∂x
∂x ∂x



(B5)

Data availability. The data supporting the findings of this study are available within the paper.

REFERENCES
1. V. Philipsen, “Mask is key to unlock full EUV potential,” Proc. SPIE
11609, 1160904 (2021).
2. A. Erdmann, P. Evanschitzky, G. Bottiglieri, et al., “3D mask effects inhigh NA EUV imaging,” Proc. SPIE 10957, 109570Z (2019).
3. A. Wong and A. Neureuther, “Mask topography effects in projection
printing of phase-shifting masks,” IEEE Trans. Electron Devices 41,
895–902 (1994).
4. M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysisof planar-grating diffraction,” J. Opt. Soc. Am. 71, 811–818 (1981).
5. M. G. Moharam, E. B. Grann, and D. A. Pommet, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings and optics,” J. Opt. Soc. Am. A 12, 1068–1076(1995).
6. M. G. Moharam, D. A. Pommet, and E. B. Grann, “Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings:
enhanced transmittance matrix approach,” J. Opt. Soc. Am. A 12,
1077–1086 (1995).

Research Article

Vol. 41, No. 8 / August 2024 / Journal of the Optical Society of America A

7. P. C. Logofatu, “Rigorous coupled-wave analysis for twodimensional gratings,” Proc. SPIE 5972, 59720Q (2005).
8. H. Tanabe, “Modeling of optical images in resists by vector
potentials,” Proc. SPIE 1674, 637–649 (1992).
9. K. D. Lucas, H. Tanabe, and A. J. Strojwas, “Efficient and rigorous
three-dimensional model for optical lithography simulation,” J. Opt.
Soc. Am. A 13, 2187–2199 (1996).
10. H. Tanabe, S. Sato, and A. Takahashi, “Fast EUV lithography simulation using convolutional neural network,” J. Micro/Nanopattern.
Mater. Metrol. 20, 041202 (2021).
11. H. Tanabe and A. Takahashi, “Data augmentation in extreme ultraviolet lithography simulation using convolutional neural network,” J.
Micro/Nanopattern. Mater. Metrol. 21, 041602 (2022).
12. H. Tanabe, A. Jinguji, and A. Takahashi, “Evaluation of convolutional
neural network for fast extreme violet lithography simulation using
3nm node mask patterns,” J. Micro/Nanopattern. Mater. Metrol. 22,024201 (2023).
13. H. Tanabe, A. Jinguji, and A. Takahashi, “Accelerating extreme ultraviolet lithography simulation with weakly guiding approximation and1499

source position dependent transmission cross coefficient formula,”
J. Micro/Nanopattern. Mater. Metrol. 23, 014201 (2024).
14. D. Gloge, “Weakly guiding fibers,” Appl. Opt. 10, 2252–2258 (1971).
15. J. D. Jackson, Classical Electrodynamics, 2nd ed. (Wiley, 1975).
16. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968).
17. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge
University, 1999).
18. D. Xu, W. Gillijns, L. E. Tan, et al., “Investigation of low-n mask in 0.33
NA EUV single patterning at pitch 28 nm metal design,” Proc. SPIE
12051, 120510H (2022).
19. S. Tomov, J. Dongarra, and M. Baboulin, “Towards dense linear
algebra for hybrid GPU accelerated manycore systems,” Parallel
Comput. 36, 232–240 (2010).
20. A. W. Snyder and W. R. Young, “Modes of optical waveguide,” J. Opt.
Soc. Am. 68, 297–309 (1978).
21. E. Gullikson, “CXRO X-ray database,” https://henke.lbl.gov/optical_
constants/.
22. S. Lin, C. Lee, Y. Chen, et al., “EUV APSM mask prospects and challenges,” Proc. SPIE 12751, 127510N (2023).