JP5113572B2 - Simulation method and simulation program - Google Patents

Description

  The present invention relates to a simulation method and a simulation program using the FDTD method.

  As an example of a structure including a metal member, there is a metal diffraction grating configured by arranging a plurality of metal strips (fine metal wires) substantially in parallel. Such metal diffraction gratings are often utilized as wire grid polarizers or polarization controllers in display related optical devices. If the width or period of each metal strip is small compared to the wavelength of light incident on the metal diffraction grating, the diffraction grating characteristics depend on the width and period of the metal strip, the polarization of the incident beam, and the incident angle. Therefore, it is necessary to design the width and period of the metal strip according to the characteristics of the incident beam, etc. so that the desired diffraction grating characteristics can be obtained.

  A simulation using an analytical method is applied to the analysis of the electromagnetic field when electromagnetic waves are incident on a structure including a metal member. For example, analytical methods such as rigorous coupled wave analysis (RCWA) may be used for diffraction gratings assuming an infinite size. However, in the analytical method illustrated in RCWA, (1) it cannot be used in the case of an incident beam having a structure or finite width that should assume a finite size, (2) an irregular boundary surface (3) Although the dielectric constant of the medium is a complex number, it is assumed that the real part does not depend on the angular frequency of the incident beam. Therefore, for example, the above-described diffraction grating having a finite length in the arrangement direction, that is, a form closer to the actual diffraction grating cannot be accurately analyzed, and the metal member is made of a noble metal such as silver, aluminum, or gold. In this case, the difference between the electromagnetic field distribution obtained as a result of the simulation in the visible region and the electromagnetic field distribution in the real space becomes large.

On the other hand, as a simulation method that can solve the above problems, either a structure that has an irregular boundary surface and that has a finite size assumption or an infinite assumption can be used. There is known an FDTD method (Finite Differential Time Domain method) that can also be applied (see, for example, Patent Documents 1 and 2). The FDTD method solves the Maxwell equation directly by dividing the FDTD region into a plurality of cells and approximating the space and time derivatives in the Maxwell equation by finite differences, and is used to analyze various metal structures. Has been. And in the case of the FDTD method, for example, various analytical models for permittivity, such as Debye, Drude, or Lorentz, are incorporated into FDTD by using frequency domain convolution. Can do. Such a method using convolution is known as a recursive convolution method (hereinafter also simply referred to as RC method).
JP 2000-105259 A JP 2000-227450 A

  However, in the simulation using the conventional FDTD method, the simulation tends to become unstable depending on the condition setting such as the value of spatial discretization.

  Therefore, an object of the present invention is to provide a simulation method and a simulation program capable of more stably performing a simulation by the FDTD method.

において、プラズマ角周波数ω_ｐ，Ｆ及び衝突周波数ν_ｃ，Ｆに含まれる誘電率を、所定波長の電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、変数ｃ_１，Ｆ，ｃ_２，Ｆが、ＦＤＴＤ法における計算の安定化条件としての下記式（３）及び式（４）を満たすと共に、

解析用速度が０．７より小さくなるように、格子点の数及び解析用速度をコンピュータによって決定することを特徴とする。 The simulation method of the present invention is a simulation method for analyzing an electromagnetic field when electromagnetic waves are incident on a structure including a metal member according to the FDTD method, and is a space in an FDTD region in which an analysis model of the structure is arranged. The discretization condition is determined by a computer based on the dielectric constant of the metal member and the spatial discretization condition determined in the discretization condition determination process is used. And an analysis step of calculating an electromagnetic field by a computer according to an FDTD method using an inductive convolution method. In the discretization condition determination step, the dielectric constant is obtained by expressing the dielectric constant by a Drude dispersion formula in the inductive convolution method. a two variables c _1, c ₂ in a recursive relationship in which is the number of grid points for representing a predetermined wavelength in the FDTD region and Following equation using the plasma angular frequency omega _{p, F} and the collision frequency [nu _{c, F} including the dielectric constant of an electromagnetic wave analysis rate and metal member in the DTD region (1), the variable c ₁ represented by the formula (2) _{, F 2} , c _{2, F} :

, The dielectric constants included in the plasma angular frequency ω _{p, F} and the collision frequency ν _{c, F} are measured values for electromagnetic waves having a predetermined wavelength, or values determined based on the measured values, variables c _{1, F} , c _{2, F} satisfies the following formulas (3) and (4) as the stabilization conditions for the calculation in the FDTD method,

The number of grid points and the analysis speed are determined by a computer so that the analysis speed is less than 0.7.

において、プラズマ角周波数ω_ｐ，Ｆ及び衝突周波数ν_ｃ，Ｆに含まれる上記誘電率を、所定波長の電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、変数ｃ_１，Ｆ，ｃ_２，Ｆが、ＦＤＴＤ法における計算の安定化条件としての下記式（７）及び式（８）を満たすと共に、

解析用速度が０．７より小さくなるように、格子点の数及び解析用速度をコンピュータに決定せしめることを特徴とする。 The simulation program according to the present invention is a simulation program for causing a computer to execute a simulation for analyzing an electromagnetic field when electromagnetic waves are incident on a structure including a metal member according to the FDTD method. Based on the discretization condition determining step for determining the spatial discretization condition of the FDTD region in which the analysis model is arranged by the computer based on the dielectric constant of the metal member, and the spatial discretization condition determined in the discretization condition determining step Then, the computer executes an analysis step of calculating an electromagnetic field when electromagnetic waves are incident on the analysis model by a computer according to an FDTD method using an inductive convolution method. In the discretization condition determination step, an inductive convolution method is used. In the recursive relational expression obtained by expressing the above dielectric constant with the Drude dispersion formula Two variables c _{1, F,} a _c 2, _F, plasma angular frequency, including the dielectric constant of the number and the analysis speed of electromagnetic waves in FDTD regions and the metal member of the lattice points for representing a predetermined wavelength in the FDTD region The variables c _{1, F} , c _{2, F} represented by the following formulas (5) and (6) using ω _{p, F} and the collision frequency ν _{c , F} :

When the dielectric constant included in the plasma angular frequency ω _{p, F} and the collision frequency ν _{c, F} is an actual value for an electromagnetic wave having a predetermined wavelength or a value determined based on the actual value, the variables c _{1, F} , c _{2, F} satisfies the following formulas (7) and (8) as the stabilization conditions of the calculation in the FDTD method,

The computer is configured to determine the number of grid points and the analysis speed so that the analysis speed is less than 0.7.

  In the simulation method and the simulation program according to the present invention, the measured value of the dielectric constant of the metal member with respect to the electromagnetic wave having the predetermined wavelength λ is used, and the value determined based on the measured value is used to satisfy the calculation stabilization condition in the FDTD method. Thus, the spatial discretization of the FDTD region is attempted. Therefore, more accurate simulation can be stably performed.

としたとき、衝突周波数ν_ｃ，Ｆ及びプラズマ角周波数ω_ｐ，Ｆは、

として与えられることが好ましい。 In the simulation method according to the present invention, the predetermined wavelength of the electromagnetic wave in the real space is λ, the velocity of the electromagnetic wave is v, the angular frequency corresponding to the predetermined wavelength is ω _λ , the number of lattice points in the FDTD region is λ _F , and the analysis is performed. Suppose that the speed of operation is v _F , the values of the real part and the imaginary part of the dielectric constant with respect to the angular frequency ω _λ are ε ^r and ε ^i, and the analysis angular frequency ω _F in the FDTD region corresponding to the angular frequency ω _λ is

Where the collision frequency ν _{c, F} and the plasma angular frequency ω _{p, F} are

Is preferably given as

としたとき、衝突周波数ν_ｃ，Ｆ及びプラズマ角周波数ω_ｐ，Ｆは、

として与えられることが好ましい。 Similarly, in the simulation program according to the present invention, the predetermined wavelength of the electromagnetic wave in the real space is λ, the velocity of the electromagnetic wave is v and the angular frequency corresponding to the predetermined wavelength is ω _λ, and the number of lattice points in the FDTD region is λ _F , The analysis speed is v _F , the real part and imaginary part of the dielectric constant for the angular frequency ω _λ are ε ^r and ε ^i, and the analysis angular frequency ω _F in the FDTD region corresponding to the angular frequency ω _λ is

Where the collision frequency ν _{c, F} and the plasma angular frequency ω _{p, F} are

Is preferably given as

In this case, ν _{c, F} and ω _{P, F} are functions of λ _F and v _F. Therefore, the conditions satisfied by the variable _{c 1, F} and _{c 2, F, v F} is to take less than 0.7 value, it is possible to determine the lambda _F and _{v F.}

  In the simulation method according to the present invention, in the analysis step, for the TE mode in the electromagnetic wave, a first electromagnetic field analysis formula obtained by applying the inductive convolution method to the wave equation derived from the Maxwell equation Is preferably solved by a computer using an inductive relational expression. Similarly, in the simulation program according to the present invention, in the analysis step, the first electromagnetic field analysis obtained by applying the inductive convolution method to the wave equation derived from the Maxwell equation for the TE mode in the electromagnetic wave. It is preferable to solve the above formula by a computer using an inductive relational expression.

  In this case, the calculation amount can be reduced because the wave equation is used.

  Further, in the analysis step in the simulation method according to the present invention, for the TM mode in the electromagnetic wave, a second electromagnetic field analysis equation obtained by applying the inductive convolution method to the Maxwell equation is expressed as an inductive relational equation. It is also possible to adopt a mode of solving using Similarly, in the analysis process in the simulation program according to the present invention, the second electromagnetic field analysis formula obtained by applying the recursive convolution method to the Maxwell equation for the TM mode in the electromagnetic wave is represented by an inductive relationship. It is also possible to use an equation that is solved by a computer.

  Further, in the analysis step in the simulation method according to the present invention, when the electromagnetic wave includes both the TM mode and the TE mode, the above-described recursive convolution method is applied to the wave equation derived from the Maxwell equation for the TE mode. The second electromagnetic field obtained by solving the first electromagnetic field analysis formula obtained using the recursive relational expression and applying the recursive convolution method to the Maxwell equation for the TM mode. The analysis formula may be solved using an inductive relational expression, and the electromagnetic field in the FDTD region may be calculated based on the electromagnetic fields obtained for the TM mode and the TE mode. Similarly, in the analysis step in the simulation program according to the present invention, when the electromagnetic wave includes both the TM mode and the TE mode, the inductive convolution method is applied to the wave equation derived from the Maxwell equation for the TE mode. The first equation for electromagnetic field analysis obtained by application is solved by a computer using an inductive relational expression, and for TM mode, it is obtained by applying an inductive convolution method to the Maxwell equation. A mode in which the second electromagnetic field analysis equation is solved by a computer using an inductive relational expression, and the electromagnetic field in the FDTD region is calculated based on the electromagnetic field obtained for the TM mode and the TE mode. It can be.

  According to the present invention, simulation can be stably performed according to the FDTD method.

  Hereinafter, embodiments of a simulation method and a simulation program according to the present invention will be described with reference to the drawings. In the description of the drawings, the same elements are denoted by the same reference numerals, and redundant description is omitted. Further, the dimensional ratios and the like in the drawings do not necessarily match those described.

  A simulation method and a simulation program according to the present invention are for analyzing an electromagnetic field when electromagnetic waves are incident on a structure including a metal member according to an FDTD method (Finite Difference Time Domain Method). It is suitably used for design and the like. An example of a structure as an analysis target using the method and program according to the present invention is a metal diffraction grating. In this specification, unless otherwise specified, the optical element is described as a metal diffraction grating.

  FIG. 1 is a schematic diagram of a metal diffraction grating. The metal diffraction grating 10 has a plurality of metal thin wires (metal strips) 11 as a plurality of metal members extending in one direction (a direction perpendicular to the paper surface in FIG. 1) with respect to the extending direction of the metal thin wires 11. It is a diffraction grating arranged with a period Λ in a substantially orthogonal direction. For simplification of description, as shown in FIG. 1, the arrangement direction of the thin metal wires 11 is defined as the x-axis direction, the extending direction of the fine metal wires 11 is defined as the z-axis direction, and is orthogonal to the x-axis direction and the z-axis direction. The direction is the y-axis direction. Examples of the material for the fine metal wires 11 include silver, gold, and aluminum.

  The period Λ of the metal diffraction grating 10 is sufficiently small with respect to the predetermined wavelength λ of the incident beam 20 as an electromagnetic wave, for example, 5λ or less. When the incident beam 20 is incident on the metal diffraction grating 10, the incident beam 20 is diffracted by the metal diffraction grating 10. When the period Λ is sufficiently small with respect to the predetermined wavelength λ, when the incident beam 20 is incident on the metal diffraction grating 10, almost zero-order diffracted light is generated. In the following description, as shown in FIG. 1, the light diffracted with respect to the metal diffraction grating 10 to the side opposite to the incident side of the incident beam 10 is also referred to as a transmitted beam 21 and is the same side as the incident side of the incident beam 20. The light that is diffracted into the light is also referred to as a reflected beam 22. Note that the predetermined wavelength λ may be a wavelength in the visible region in the electromagnetic spectrum.

The characteristics of the metal diffraction grating 10 include the period Λ, the width (length in the x-axis direction) w and the thickness (length in the y-axis direction) h of the metal thin wire 11, the diffraction grating normal line N, and the incident beam 20. It depends on the incident angle θ _in which is formed, the polarization state of the incident beam 20, and the like. Therefore, the simulation method and program according to the present invention are effective for designing the metal diffraction grating 10 that realizes desired diffraction grating characteristics. In particular, the present invention is effective when the width W in the x-axis direction of the metal diffraction grating 10 is smaller than 30λ. This is because the simulation method and simulation program according to the present invention employs the FDTD method, which is effective when handling an object having a finite size, as its method.

  FIG. 2 is a schematic diagram of an analysis model in the FDTD region of the metal diffraction grating shown in FIG. In the FDTD method, as shown in FIG. 2, the Maxwell equation is spatially and temporally differentiated and directly numerically solved by dividing the FDTD region 100 into a plurality of cells. In the present embodiment, as shown in FIG. 2, the same reference numerals as those of the metal diffraction grating 10 shown in FIG. Moreover, the hatching in FIG. 2 is convenient for showing the metal thin wire 11 in the FDTD region 100 divided by cells.

  Usually, the metal diffraction grating 10 shown in FIGS. 1 and 2 extends in one direction (the z-axis direction in FIGS. 1 and 2), and the one in which the extension direction does not change is the analysis target. Can be treated as a two-dimensional problem. Therefore, a two-dimensional FDTD algorithm will be described.

μは透磁率である。また、Ｅは電界、Ｈは磁界である。更に、「∇」は、ｘ_ｕ及びｙ_ｕをｘ軸方向及びｙ軸方向の単位ベクトルとしたとき、ｘ_ｕ(∂／ｘ)＋ｙ_ｕ(∂／ｙ)で定義される微分演算子である。 In the FDTD method, the following Maxwell equations are used.

μ is the magnetic permeability. E is an electric field and H is a magnetic field. Furthermore, "∇" when the x _u and y _u and a unit vector in the x-axis and y-axis directions, is _{x u (∂ / x) +} y u (∂ / y) differential operators defined by .

式（１７）中のε（ｒ，ω）は金属細線１１の誘電率である。 D in the equation (16) represents the electric flux density and is given by the following equation.

In the equation (17), ε (r, ω) is the dielectric constant of the thin metal wire 11.

式（１８）において、ε（ｔ）＝Ｆ｛ε（ω）｝であり、Ｆはフーリエ変換を表す。 The electric flux density D in the time domain is expressed by the following equation by performing Fourier transform on both sides of Equation (17) and applying the convolution theorem.

In Expression (18), ε (t) = F {ε (ω)}, and F represents a Fourier transform.

  When the dielectric constant ε (ω) is expressed by a special function having an angular frequency ω, the convolution integral of the equation (18) can be calculated recursively. This is known as the RC method.

式（１９）において、ω_ｐはプラズマ周波数、ν_ｃは衝突周波数である。ε_∞は、周波数が無限大のときの誘電率であり１に等しい。χ（ω）は電気比感受率である。 In the present embodiment, a Drude primary approximation formula (Drude dispersion formula) is adopted as a function representing the dielectric constant ε (ω) of the thin metal wire 11. In this case, the dielectric constant ε (ω) is expressed as in Expression (19).

In Equation (19), ω _p is the plasma frequency and ν _c is the collision frequency. epsilon _∞ is equal to is 1 the dielectric constant when the frequency is infinite. χ (ω) is an electrical specific susceptibility.

Ｕ（ｔ）は、時間の単位ステップ関数であって、ｔ＝０の場合、Ｕ（ｔ）＝０であり、ｔ＞０の場合、Ｕ（ｔ）＝１である。 When χ (ω) is Fourier transformed, that is, F {χ (ω)} is expressed as χ (t), χ (t) is expressed by the following equation.

U (t) is a unit step function of time. When t = 0, U (t) = 0, and when t> 0, U (t) = 1.

が得られる。より具体的には、式（２１）は、区間［ｍΔｔ，（ｍ＋１）Δｔ］にわたってＥを定数としてとることによって、畳み込み積分から得られる。なお、ｍは０以上の整数である。 When Δt is a time step interval, when the convolution integral of Expression (18) is evaluated at time t = nΔt (n is an integer of 0 or more),

Is obtained. More specifically, equation (21) is obtained from convolution integral by taking E as a constant over the interval [mΔt, (m + 1) Δt]. Note that m is an integer of 0 or more.

ここで、Ψは、累積変数であり下記式（２３）が成り立つ。

Equation (21) can be written as:

Here, Ψ is a cumulative variable, and the following formula (23) is established.

と表すことができる。この式（２４）において、

である。式（２５）より、

が成立する。式（２６）において、

である。また、Ψ^０＝０、Ψ^１＝Ｅ^１χ^０、Ψ^２＝Ｅ^２χ^０＋ｃ_１Ｅ^１＋ｃ_２Ψ^１とする。 Equation (23) is

It can be expressed as. In this equation (24),

It is. From equation (25)

Is established. In equation (26),

It is. Further, ψ ⁰ = 0, ψ ¹ = E ¹ χ ⁰ , and ψ ² = E ² χ ⁰ + c ₁ E ¹ + c ₂ ψ ¹ are set.

式（２９）及び式（３０）では、ε_∞＝１、Δｔ＝１である。また、式（２９）及び式（３０）中、χ^０は式（２５）より、

で与えられる。更に、ｎ≦０であるすべてのｎについてΨの値は０とする。 From the expressions (24) to (28), the following inductive relational expression regarding the cumulative variable Ψ is obtained.

In the equations (29) and (30), ε _∞ = 1 and Δt = 1. Moreover, in Formula (29) and Formula (30), χ ⁰ is from Formula (25),

Given in. Furthermore, the value of Ψ is 0 for all n where n ≦ 0.

As will be described later, when the FDTD method is performed, c ₁ and c ₂ included in the equations (29) and (30) and expressed using ω _p and ν _c are ω _p and adopting _{c 1, F} and _{c 2, F} as a function of the plasma angular frequency omega _{p, F} and the collision frequency _{[nu c, F} in the FDTD region 100 corresponding to the [nu _c.

  In general, in two dimensions, an incident beam (electromagnetic wave) 20 in an arbitrary polarization state can be divided into two orthogonal polarization states, that is, a TM mode and a TE mode. Therefore, the algorithm of the FDTD method will be described below separately for the TM mode and the TE mode.

[TM mode]
In the TM mode, Equation (15) and Equation (16) are solved using the central difference.

この場合、式（１６）に等価な差分形式の方程式は、以下のように表される。

式（３３）中のＤ_１は∇に対応する２次の中心差分形式の演算子である。また、式（３３）中のΨ^ｎ＋１−Ψ^ｎは、式（２９）により与えられる。 Equation (32) is used as the differential form of the first derivative with respect to time in Equation (16).

In this case, the differential equation equivalent to equation (16) is expressed as follows.

D ₁ in the equation (33) is a second-order central difference form operator corresponding to ∇. In addition, Ψ ^{n + 1} −Ψ ⁿ in Expression (33) is given by Expression (29).

Here, an equation of a difference format equivalent to Equation (16) is shown, but for Equation (15), the following equation which is a difference format in the conventional FDTD method is used.

[TE mode]
Next, the case of the TE mode will be described.

本実施形態では∇・Ｅ＝０を仮定するが、高い伝導性を有する金属の場合、それらは自由電荷を長い間保持しないので、この仮定は妥当である。 In the TE mode, the following wave equation derived from the equations (15) and (16), which are Maxwell equations, is used.

In this embodiment, ∇ · E = 0 is assumed. However, in the case of metals having high conductivity, this assumption is valid because they do not retain free charge for a long time.

この場合、帰納的畳み込みを適用した式（３５）に等価な差分形式の方程式は、以下のように表される。

式（３７）中のＤ_１ ^２は∇^２に対応する２次の中心差分形式の演算子であり、ｖは位相速度である。また、式（３７）中のΨ^ｎ＋１＋Ψ^ｎ−１−２Ψ^ｎは式（３０）から与えられる。 Equation (36) is used as the differential form of the second derivative with respect to time of Equation (16).

In this case, the differential equation equivalent to the equation (35) to which the recursive convolution is applied is expressed as follows.

D ₁ ² in the formula (37) is an operator of second order central difference format corresponding to ∇ ^2, v is the phase velocity. In addition, Ψ ^{n + 1} + Ψ ⁿ⁻¹ −2Ψ ⁿ in Expression (37) is given from Expression (30).

  In the present embodiment, when the incident beam 20 is given in the TM mode, that is, when the incident beam 20 is a TM wave, an expression that is an electromagnetic field analysis expression (second electromagnetic field analysis expression). The electric field E and the magnetic field H are calculated using (33), Expression (34), and Expression (29). When the incident beam 20 is given in the TE mode, that is, when the incident beam 20 is a TE wave, the expression (37), which is an expression for electromagnetic field analysis (first electromagnetic field analysis expression), (30) is used to calculate the electric field E.

  Next, a method for determining a spatial discretization condition in the FDTD region in the present embodiment will be described.

  In the FDTD method, a grid composed of a plurality of cells is applied to the FDTD region 100 as shown in FIG. In the FDTD method solved by the simulation method of this embodiment, the discretization condition of the FDTD region 100 is determined based on the stability of the algorithm.

From Equation (29) and Equation (30), it is necessary for c ₁ and c ₂ to satisfy the following conditions for the stability of the algorithm.

c ₁ and c ₂ are expressed using ω _p and ν _c as shown in the equations (27) and (28).

が成り立つ。 On the other hand, from the equation (19), when the real part and the imaginary part of the dielectric constant ε (ω) are ε ^r (ω) and ε ⁱ (ω),

Holds.

Therefore, ω _p and ν _c are expressed as follows.

In the present embodiment, when the angular frequency corresponding to the predetermined wavelength λ of the incident beam 20 is ω _λ , the angles ε ^r (ω) and ε ⁱ (ω) in the equations (40) and (41) are The measured value for the frequency ω _λ is adopted.

When the measured value of the dielectric constant is not obtained for the predetermined wavelength λ, the dielectric constant for the predetermined wavelength λ, that is, ε ^r is obtained from the local fitting of the measured dielectric constant using Drude dispersion. Obtain (ω) and ε ⁱ (ω). The local fitting can be obtained by performing Drude's ε-ω plot so as to pass the measured value of the dielectric constant for the wavelength whose dielectric constant is measured closest to the predetermined wavelength λ.

For example, the measured value of the dielectric constant at a wavelength of 500 nm (angular frequency 3.77 × 10 ¹⁵ / sec) is not obtained, while the wavelength 496 nm (or angular frequency 3.8 × 10 ¹⁵ / sec) and the wavelength 551. It is assumed that dielectric constants are obtained for 111 nm (angular frequency 3.42 × 10 ¹⁵ / sec), respectively. In this case, the dielectric constant at a wavelength of 500 nm is calculated from Drude's ε-ω plot with dielectric constant ε = −9.564 + i0.309 for an angular frequency of 3.8 × 10 ₁₅ / sec (see FIG. 3).

  In order to obtain a more accurate dielectric constant, Drude's ε described above based on the dielectric constant for the wavelength whose dielectric constant is closest to the predetermined wavelength λ on the high wavelength side and the short wavelength side with respect to the predetermined wavelength λ. It is preferable to perform an -ω plot and calculate an average value of two estimated values obtained as dielectric constants of the predetermined wavelength λ to obtain an estimated value. When using, for example, light output from an LED as light incident on a metal diffraction grating, in many cases, the dielectric constant may not be actually measured for the wavelength of light output from the LED. The calculation method as described above is preferable. In the following description, the actual measurement value includes a value determined by performing the Drude ε-ω plot.

式（４４）及び式（４５）中のε^ｒ及びε^ｉは角周波数ω_λに対する実測値であり定数である。 Further, in the present embodiment, the angular frequency ω in the equations (42) and (43) is replaced with the angular frequency ω _F in the FDTD region, and the plasma angular frequency ω _p and the collision frequency ν _c are plasma angular frequencies in the FDTD region 100. ω _{p, F} and collision frequency ν _{c, F} are made to correspond. In this case, the plasma angular frequency ω _{p, F} and the collision frequency ν _{c, F} are expressed as follows.

In the equations (44) and (45), ε ^r and ε ⁱ are measured values for the angular frequency ω _λ and are constants.

式（４６）中のｖは入射ビーム２０の実空間での位相速度である。そして、λ_Ｆは、ＦＤＴＤ領域１００において波長λを表すのに必要とする格子点１１０の数であり、ｖ_ＦはＦＤＴＤ領域１００における入射ビーム２０の位相速度（解析用速度）である。 In the present embodiment, the angular frequency ω _λ in the real space is associated with the angular frequency (angular frequency for analysis) ω _F in the FDTD region 100 in the equations (44) and (45) by the following equation. Therefore, ω _{p, F} and ν _{c, F} are free parameters.

In Expression (46), v is the phase velocity of the incident beam 20 in real space. Λ _F is the number of lattice points 110 necessary to represent the wavelength λ in the FDTD region 100, and v _F is the phase velocity (analysis velocity) of the incident beam 20 in the FDTD region 100.

と表されることになる。 When ω _{p, F} and ν _{c, F} in the FDTD region 100 are defined corresponding to ω _p and ν _c in the real space as described above, the stabilization conditions shown in the equations (38) and (39) are as follows:

Will be expressed.

また、アルゴリズムの安定性として、式（４９）及び式（５０）と共に、ＦＤＴＤ法で通常知られている以下の関係式を満たす必要がある。

C _{1, F} and c _{2, F} are obtained by replacing ω _p and ν _c in c ₁ and c ₂ with ω _{p, F} and ν _{c, F} , and are represented by the following equations.

Further, as the stability of the algorithm, it is necessary to satisfy the following relational expression normally known by the FDTD method together with the expressions (49) and (50).

As described above, when a predetermined value for the angular frequency ω _λ is used as ε ^r and ε ⁱ and ω _{P, F} and ν _{c, F} are set as free parameters by using the equation (46), the equation (49) and to determine the lambda _F and _{v F} so as to satisfy the equation (50) and equation (51). Since λ _F is the number of lattice points 110 for representing the predetermined wavelength λ, the spatial discretization condition is determined.

  Next, an analysis apparatus, a simulation method, and a simulation program for executing the above algorithm will be described.

  FIG. 4 is a block diagram of an analysis apparatus for carrying out an example of the simulation method according to the present invention.

  The analysis device 30 shown in FIG. 4 is a computer device such as a personal computer device or a workstation, and includes an input unit 31, a memory 32, a central processing unit 33, and an output unit 34 as main components.

The input means 31 is a keyboard or the like, and accepts input data including parameters for performing a simulation. The above parameters include various conditions regarding the incident beam 20 including the wavelength λ of the incident beam 20 and the angular frequency ω _λ corresponding thereto, measured values of ε ^r and ε ^{i with} respect to the angular frequency ω _λ , conditions of the thin metal wire 11, and boundary conditions And the maximum time t _max (t _max = nΔt) in the simulation. The various conditions of the fine metal wire 11 are the width w and thickness h of the fine metal wire 11 and the period Λ. In addition to the above, the conditions regarding the incident beam 20 include the polarization state of the incident beam 20 (for example, TE mode or TM mode), the incident angle θ _{in of} the incident beam 20 to the metal diffraction grating 10 (metal diffraction grating normal) The angle between N and the incident beam 20).

The memory 32 is a FDTD program (simulation program) 40 and an FDTD program for performing FDTD analysis in accordance with various programs for operating the analysis device 30, input data from the input means 31, and a simulation method described later. The calculation result during or after execution of 40 is stored. The FDTD program 40 is based on a spatial discretization condition based on the equations (44) to (51), that is, a discretization condition determination module 41 that determines λ _F , and an algorithm described separately for the TM mode and the TE mode. And an analysis module 42 for executing the FDTD method.

  The central processing unit 33 reads the FDTD program 40 stored in the memory 32 and executes various calculations for simulation. The central processing unit 33 also has a function of reading various programs stored in the memory 32 and controlling components of the analysis device 30.

  The output unit 34 is, for example, a display and displays a simulation result obtained by being calculated by the central processing unit 33.

  Next, a simulation method according to the FDTD method will be described. FIG. 5 is a flowchart showing a simulation method. The simulation is performed by the central processing unit 33 reading and executing the FDTD program 40 stored in the memory 32 of the analysis device 30.

  As shown in FIG. 5, in step S10, the operator inputs simulation parameters via the input means 31 (input process). When the input means 31 receives a parameter, the memory 32 stores the parameter.

Next, in step S20, the discretization condition is determined (discretization condition determination step). In this discretization condition determination step, the central processing unit 33 uses the values of ε ^r and ε ⁱ input in step S10 and the equations (44) to (51) to obtain equations (47) and (48). ) And equation (51) satisfying λ _F and v _F. At this time, when there are a plurality of λ _F and ν _F satisfying the equations (47), (48), and (51), for example, the central processing unit 33 compares them and λ _F is minimized. Is preferably selected. This is because stable simulation is possible with a small number of grid points 110.

  In subsequent step S30, the metal diffraction grating 10 satisfying the various conditions input in step S10 is set in the FDTD region 100 defined by the spatial discretization condition determined in step S20, and the electromagnetic field in the FDTD region 100 is set. Is calculated for each time step interval Δt (analysis step). In this analysis step (step S30), the boundary condition input in step S10 is used as necessary. In this embodiment, Δt = 1.

  Since algorithms used in the TE mode and the TM mode are different, in step S30, (i) when the incident beam 20 is given as the TE mode, (ii) when the incident beam 20 is given as the TM mode, The analysis is performed in two steps.

  FIG. 6 is a flowchart of the analysis process of FIG. In the analysis process, first, in step S31, the polarization state of the incident beam 20 is determined. Hereinafter, the case where the incident beam 20 is given in the TE mode in step S10 and the case where the incident beam 20 is given in the TM mode will be described.

  (I) If the incident beam 20 is given in the TE mode, in step S32A, the electromagnetic field of the cell is initialized. Specifically, the components in the axial directions of the electric field E and the magnetic field H in each cell are set to zero.

Then, it progresses to step S33A, calculates the electric field E using Formula (37) and Formula (30), updates the electric field E, and updates the time t. As the value of _c 1, _{c 2} in equation (30), may be employed values of _{c 1, F} and _{c 2, F} defined by lambda _F and _{v F} determined in step S20.

When the electric fields E in all cells are calculated in step S33A, the process proceeds to step S34A, and it is determined whether or not the updated time t is the maximum time _tmax . If the time t is smaller than the maximum time _tmax (Y in step S34A), the process returns to step S33A.

Next, (ii) the case where the incident beam 20 is given in the TM mode will be described. In this case, the cell is initialized in step S32B. Specifically, the components in the axial directions of the electric field E and the magnetic field H in each cell are set to zero. Next, the process proceeds to step S33B, the electric field E and the magnetic field H of each cell are calculated using the equations (33), (34), and (29), the electric field E and the magnetic field H are updated, and the time t is set. Update. The values of c ₁ and c ₂ in equation (29) are the same as in the TE mode.

When the electric field E and the magnetic field H in all the cells are calculated in step S33B, the process proceeds to step S34B, and it is determined whether or not the updated time t is the maximum time _tmax . If the time t is smaller than the maximum time t _max (“Y” in step S34B), the process returns to step S33B. In FIG. 6, in order to explain separately for the TE mode and the TM mode, the determination step of the magnitude relationship between the time t and the maximum time t _max is divided into two, but step S34A and step S34B is substantially the same step.

In step S34A or step S34B, when the time t is equal to or greater than the maximum time t _max (in the case of “N” in steps S34A and S34B), the process proceeds to step S40 shown in FIG. Output (output process).

式（５２）及び式（５３）中において、ａ＝(σ_ｄΔｔ)／（２ε_ｄ）であり、σ_ｄ及びε_ｄは隣接する金属細線１１間の媒質の電導度及び誘電率である。また、μは隣接する金属細線１１間の媒質の透磁率である。 The implementation of the FDTD method in the present embodiment is the same as the simulation method of the conventional FDTD method except that the above-described equations are used, and thus detailed description thereof is omitted. As for the cells other than the thin metal wires 11, the following equation (52) is obtained for the TE mode, similarly to the FDTD method in the case of dealing with the Maxwell equation or the wave equation for electromagnetic waves propagating in normal air or vacuum. Using the TM mode, the equation (53) and the following equation (54) may be used.

In the equations (52) and (53), a = (σ _d Δt) / (2ε _d ), and σ _d and ε _d are the electric conductivity and dielectric constant of the medium between the adjacent metal thin wires 11. Further, μ is the magnetic permeability of the medium between the adjacent fine metal wires 11.

  In the simulation method and the FDTD program of the present embodiment described above, the algorithm is formulated by applying the Drude dispersion formula to the dielectric constant of the fine metal wire 11 of the metal diffraction grating 10 as the analysis target. In addition, the FDTD region 100 is spatially discretized by adopting measured values for the dielectric constant of the fine metal wires 11. And the stability of the algorithm is also guaranteed by spatial discretization. As a result, the electromagnetic field when the incident beam 20 is incident on the metal diffraction grating 10 can be more accurately and stably simulated.

  Further, in the simulation for the TE mode, the wave equation derived from the Maxwell equation is used, so that the simulation algorithm can be simplified. As a result, the time required for simulation is shortened. Furthermore, since the simulation method and the FDTD program of the present embodiment employ the RC method, the amount of data to be stored in the memory 32 can be reduced.

  Next, simulation results of the metal diffraction grating 10 using the simulation method and the FDTD program described in this embodiment will be described. Here, simulations were performed for two metal diffraction gratings 10 as Example 1 and Example 2, respectively.

The simulation conditions of Example 1 are as follows.
-Material of the fine metal wire 11: Period of the silver-metal diffraction grating 10 Λ: 176.125 nm
Fill factor (w / Λ) of metal diffraction grating 10: 0.5
-Thickness h of the fine metal wire 11: 100 nm
The wavelength λ of the incident beam 20: 704.5 nm (angular frequency ω _λ : 2.676 × 10 ¹⁵ / sec)
-Phase velocity v of incident beam 20: 3 × 10 ⁸ m / sec
- angular value of the frequency _ω λ for ε ^r: -23.405
The value of ε ⁱ for the angular frequency ω _λ : 0.387

The simulation conditions of Example 2 are as follows.
-Material of the fine metal wire 11: Period of the silver-metal diffraction grating 10 Λ: 496 nm
Fill factor (w / Λ) of metal diffraction grating 10: 0.5
・ Thickness h of the fine metal wire 11: 496 nm
The wavelength λ of the incident beam 20: 496 nm (angular frequency ω _λ : 3.8 × 10 ¹⁵ / sec)
Phase velocity v of incident beam 20: 3 × 10 ⁸ m / sec
The value of ε ^r with respect to the angular frequency ω _λ : −9.564
The value of ε ⁱ with respect to the angular frequency ω _λ : 0.309

In Examples 1 and 2, simulation was performed based on the flowcharts shown in FIGS. 5 and 6 while changing the incident angle θ _in when the incident beam 20 was given as the TE mode and the TM mode, respectively. In Examples 1 and 2, it is assumed that the size of the metal diffraction grating 10 is infinite in the x-axis direction.

In the case of the conditions of Example 1, in step S20, when λ is 704.5 nm, ε ^r is −23.405, and ε ⁱ is 0.387, λ _F is 35.227, and ν _F was 0.538. In this case, c1 _{, F} is 0.225, and c2 _{, F} is 1.002. In the case of the condition of the second embodiment, in step S20, when λ is 496 nm, ε ^r is −9.564, ε ⁱ is 0.309, λ _F is 35.227, and ν _F was 0.538. In this case, c1 _{, F} is 0.098, and c2 _{, F} is 1.003.

FIG. 7 shows a transmission spectrum and a reflection spectrum based on the simulation result in the case of the first embodiment. The horizontal axis of FIG. 7 indicates the incident angle θ _in (°), and the vertical axis indicates the transmittance (%) and the reflectance (%). The transmittance and the reflectance are defined as the ratio of the light intensity of the transmitted beam and the reflected beam to the light intensity of the incident beam 20. FIG. 8 shows a transmission spectrum and a reflection spectrum based on the simulation result of Example 2. The horizontal and vertical axes in FIG. 8 are the same as in FIG. The simulation results in FIGS. 7 and 8 are based on the far field transmittance and reflectance of an infinitely large metal diffraction grating 10 made of silver. The result in the far fill is calculated by integrating the complex electromagnetic field obtained by executing the FDTD method along a closed curve surrounding the metal diffraction grating 10.

As shown in FIGS. 7 and 8, the transmission characteristics and the reflection characteristics can be obtained for the TE mode and the TM mode by the simulation. From the above simulation, it can be seen that a difference in characteristics due to a difference in the configuration of the metal diffraction grating 10 appears. For example, as shown in FIG. 7, in the metal diffraction grating 10 of Example 1, most of the TE mode is reflected and transmission is suppressed by the simulation. Moreover, in the metal diffraction grating 10 of Example 1, while many TM modes can be permeate | transmitted, it is possible to reflect one part. Thus, it can be seen that polarization separation is achieved in the configuration of the metal diffraction grating 10. Further, as shown in FIG. 8, in the metal diffraction grating 10 of Example 2, both the TM mode and the TE mode can be transmitted, but the TE mode can transmit more, and the TM mode can be transmitted. It can be seen that the transmittance can be suppressed to 10% to 40% when the incident angle θ _in is 0 to 60 °.

  Therefore, if the simulation method and the FDTD program are used, the metal diffraction grating 10 having desired transmission characteristics and reflection characteristics can be designed by changing the parameters related to the metal diffraction grating 10.

  FIG. 9 is a side view schematically showing a configuration of a light guide plate unit to which the metal diffraction grating shown in FIG. 1 is applied. The light guide plate unit 50 includes a light guide plate 51 having a substantially rectangular parallelepiped shape and a metal diffraction grating 10 provided on a surface 51 a of the light guide plate 51. The light guide plate unit 50 can be used as a surface light source device by arranging a light source (not shown) on the side of the side surface 51 b of the light guide plate 51.

  Light (beam) output from the light source and incident into the light guide plate 51 from the side surface 51b propagates in the light guide plate 51 toward the side surface 51c. When the light propagating in the light guide plate 51 enters the metal diffraction grating 10, a part of the incident light is emitted from the light guide plate unit 50 by the diffraction of the metal diffraction grating 10, and the other part is mainly the light guide plate. 51 is returned.

  As described above, when the metal diffraction grating 10 is used as a part of the light guide plate unit 50, the metal diffraction grating 10 is preferably designed as a polarization separation element. Specifically, the metal diffraction grating 10 is designed so that a predetermined polarization component (for example, a TM mode as a P polarization component) is selectively transmitted and other light is reflected back into the light guide plate 51. It is preferable that In this case, since a beam (transmitted beam 21) in which a predetermined polarization component is dominant is emitted from the light guide plate unit 50, the emitted beam from the light guide plate unit 50 is used as a backlight in a liquid crystal display device, for example. It can be preferably used. In a liquid crystal display device using such a backlight, an optical element that removes unnecessary polarized light is not required, and thus it is possible to reduce the size and thickness.

  In order to design the metal diffraction grating 10 having the above-described characteristics, the above-described simulation method and FDTD program are effective. For example, when the metal diffraction grating 10 is used in a surface light source device included in a liquid crystal display device, the simulation may be performed with the predetermined wavelength λ as a wavelength in the visible region.

  When the metal diffraction grating 10 is provided on the light guide plate 51, as shown in FIGS. 1 and 2, a simulation may be performed with only the metal diffraction grating 10 as an analysis target. For example, as shown in FIG. As described above, the light guide plate unit 50 including the plate-like light guide plate 51 may be the analysis target.

式（５５）及び式（５６）中においては、ａ＝(σ_ｄΔｔ)／（２ε_ｄ）であり、σ_ｄ及びε_ｄは導光板５１の電導度及び誘電率である。また、μは導光板５１の透磁率である。 In this case, in the cell corresponding to the light guide plate 51, the following equation may be used as an algorithm of the FDTD method. That is, Equation (55) is used in the TE mode, and Equation (56) and Equation (57) are used in the TM mode.

In the equations (55) and (56), a = (σ _d Δt) / (2ε _d ), and σ _d and ε _d are the conductivity and dielectric constant of the light guide plate 51. Further, μ is the magnetic permeability of the light guide plate 51.

  Moreover, in the said embodiment, although the incident beam 20 was given as TE mode or TM mode, you may give the incident beam 20 of the non-polarization state containing TE mode and TM mode, for example. In this case, in the analysis step (step S30), for the TE mode in the incident beam 20, the electric field E is calculated in the same manner as described in step S33A, and the TM mode in the incident beam 20 is calculated. Thus, the electric field E and the magnetic field H may be calculated in the same manner as described in step S33B, and the electric field E and the magnetic field H of each cell may be combined.

As a method of combining the electric field E and the magnetic field H calculated for the TM mode and the TE mode, a weighted average of each of the electric field E and the magnetic field H calculated for each cell can be taken. For example, when the electric field calculated for the TM mode is E _TM and the electric field calculated for the TE mode is E _TE corresponding to one cell, the electric field in that cell is (E _TM + E _TE ) / 2 can be used. The same applies to the magnetic field. When taking the weighted average in this way, for example, the weighted average is calculated as described above for the electric field E and the magnetic field H calculated on the assumption that the incident beam 20 is light composed of the TE mode and the TM mode. Thus, the electric field E and the magnetic field H when the incident beam 20 in a non-polarized state is incident can be obtained.

  The metal diffraction grating 10 can be used for the light guide plate unit 50 as described with reference to FIG. 9. For example, the metal diffraction grating 10 can be used as a polarization separation element arranged at a predetermined position in the optical system. You can also. In this case, for example, it is conceivable to design the metal diffraction grating 10 so as to completely separate the TE mode and the TM mode.

  Moreover, although the simulation program according to the present invention has been described as a program that can handle both the TE mode and the TM mode, for example, the program for the TE mode and the program for the TM mode can also be performed. In the analysis step (step S30) in the simulation method, the determination step in step S32 shown in FIG. 6 may be omitted, and the simulation may be performed along the TE mode and TM mode steps shown in FIG. Is possible. Furthermore, although the two-dimensional FDTD method has been described in the above embodiment, the method can be extended to the three-dimensional FDTD method. Furthermore, although the metal diffraction grating 10 was illustrated and demonstrated in the said embodiment as a structure containing a metal member, the said structure is not limited to a metal diffraction grating, The simulation method and simulation program which concern on this invention are metal members. It can be applied to any structure including

It is a schematic diagram of the optical element for designing by applying one embodiment of the simulation method according to the present invention. FIG. 2 is a schematic diagram of the optical element shown in FIG. 1 set in an FDTD region. It is drawing which shows an example of the measured value of dielectric constant with respect to an angular frequency, and a fitting curve. It is a block diagram of the analysis device for applying one embodiment of the simulation method concerning the present invention. It is a flowchart of one Embodiment of the simulation method which concerns on this invention. It is a flowchart of the analysis process of FIG. It is drawing which shows the transmission spectrum and reflection spectrum based on the simulation result of Example 1. FIG. It is drawing which shows the transmission spectrum and reflection spectrum based on the simulation result of Example 2. FIG. It is a side view which shows typically the structure of the light-guide plate unit containing the optical element shown in FIG.

Explanation of symbols

  DESCRIPTION OF SYMBOLS 10 ... Metal diffraction grating (structure), 11 ... Metal fine wire (metal member), 20 ... Incident beam (electromagnetic wave), 100 ... FDTD area | region, 110 ... Lattice point.

Claims (10)

A simulation method for analyzing an electromagnetic field when electromagnetic waves are incident on a structure including a metal member according to the FDTD method,
A discretization condition determining step of determining a spatial discretization condition of the FDTD region in which the analysis model of the structure is arranged by a computer based on a dielectric constant of the metal member;
Based on the spatial discretization condition determined in the discretization condition determination step, an electromagnetic field when the electromagnetic wave is incident on the analysis model is calculated by the computer according to the FDTD method using an inductive convolution method. Analysis process,
With
In the discretization condition determination step, two variables c _{1, F 2,} c _{2, F} in an inductive relational expression obtained by expressing the dielectric constant by a Drude dispersion formula in the inductive convolution method, Plasma angular frequencies ωp _{, F} and collision frequency ν including the number of lattice points for representing the predetermined wavelength of the electromagnetic wave in the FDTD region, the speed for analyzing the electromagnetic wave in the FDTD region, and the dielectric constant of the metal member _The variables c _{1, F} , c _{2, F} represented by the following formulas (1) and (2) using c _{and F} :

In
When the dielectric constant included in the plasma angular frequency ω _{p, F} and the collision frequency ν _{c, F} is an actual value for the electromagnetic wave having the predetermined wavelength or a value determined based on the actual value, the variable c _{1 , F 1} , c _{2 and F} satisfy the following formulas (3) and (4) included in the stabilization conditions in the FDTD method,

The simulation method, wherein the number of grid points and the analysis speed are determined by the computer so that the analysis speed is less than 0.7. In the real space, the predetermined wavelength is λ, the electromagnetic wave velocity is v, the angular frequency corresponding to the predetermined wavelength is ω _λ , the number of lattice points in the FDTD region is λ _F , and the analysis speed is v _F. The values of the real part and the imaginary part of the dielectric constant with respect to the angular frequency ω _λ are ε ^r and ε ⁱ , and the angular frequency for analysis ω _F in the FDTD region corresponding to the angular frequency ω _λ is

Where the collision frequency ν _{c, F} and the plasma angular frequency ω _{p, F} are

The simulation method according to claim 1, wherein the simulation method is given as: In the analysis step,
For the TE mode in the electromagnetic wave, the first electromagnetic field analysis formula obtained by applying the recursive convolution method to the wave equation derived from the Maxwell equation is obtained using the recursive relational expression. The simulation method according to claim 1, wherein the simulation method is solved by a computer. In the analysis step,
For the TM mode in the electromagnetic wave, the second electromagnetic field analysis formula obtained by applying the recursive convolution method to the Maxwell equation is solved by the computer using the recursive relational expression. The simulation method according to any one of claims 1 to 3, wherein the simulation method is characterized. In the analysis step,
When the electromagnetic wave includes both the TM mode and the TE mode, the first electromagnetic field analysis for the TE mode is obtained by applying the recursive convolution method to the wave equation derived from the Maxwell equation. An equation is solved using the recursive relational expression, and for the TM mode, a second electromagnetic field analysis equation obtained by applying the recursive convolution method to the Maxwell equation is expressed as the recursive equation. 3. The simulation according to claim 1, wherein an electromagnetic field in the FDTD region is calculated based on an electromagnetic field obtained by solving using a relational expression and obtained for the TM mode and the TE mode. Method. A simulation program for causing a computer to execute a simulation for analyzing an electromagnetic field when electromagnetic waves are incident on a structure including a metal member according to the FDTD method,
A discretization condition determining step of determining a spatial discretization condition of the FDTD region in which the analysis model of the structure is arranged by a computer based on a dielectric constant of the metal member;
Based on the spatial discretization condition determined in the discretization condition determination step, an electromagnetic field when the electromagnetic wave is incident on the analysis model is calculated by the computer according to the FDTD method using an inductive convolution method. Analysis process,
To the computer,
In the discretization condition determination step, two variables c _{1, F 2,} c _{2, F} in an inductive relational expression obtained by expressing the dielectric constant by a Drude dispersion formula in the inductive convolution method, The number of lattice points for representing the predetermined wavelength of the electromagnetic wave in the FDTD region, the analysis speed of the electromagnetic wave in the FDTD region, and the plasma angular frequency ω _p and the collision frequency ν _{c, F} including the dielectric constant of the metal member The variables c _{1, F} , c _{2, F} represented by the following formulas (8) and (9) are used:

In
When the dielectric constant included in the plasma angular frequency ω _{p, F} and the collision frequency ν _{c, F} is an actual value for the electromagnetic wave having the predetermined wavelength or a value determined based on the actual value, the variable c _{1 , F 1} , c _{2 and F} satisfy the following formulas (10) and (11) included in the stabilization conditions in the FDTD method,

A simulation program that causes the computer to determine the number of grid points and the analysis speed so that the analysis speed is less than 0.7. The predetermined wavelength of the electromagnetic wave in the real space is λ, the speed of the electromagnetic wave is v, the angular frequency corresponding to the predetermined wavelength is ω _λ , the number of lattice points in the FDTD region is λ _F , and the analysis speed is v _F, and real and imaginary values of the permittivity for the angular frequency ω _λ are ε ^r and ε ^i, and the analysis angular frequency ω _F in the FDTD region corresponding to the angular frequency ω _λ is

Where the collision frequency ν _{c, F} and the plasma angular frequency ω _{p, F} are

The simulation program according to claim 6, which is given as: In the analysis step,
For the TE mode in the electromagnetic wave, the first electromagnetic field analysis formula obtained by applying the recursive convolution method to the wave equation derived from the Maxwell equation is obtained using the recursive relational expression. The simulation program according to claim 6, wherein the simulation program is solved by a computer. In the analysis step,
For the TM mode in the electromagnetic wave, the second equation for electromagnetic field analysis obtained by applying the recursive convolution method to the Maxwell equation is solved by the computer using the recursive relational expression. The simulation program according to any one of claims 6 to 8, wherein: In the analysis step,
When the electromagnetic wave includes both the TM mode and the TE mode, the first electromagnetic field analysis for the TE mode is obtained by applying the recursive convolution method to the wave equation derived from the Maxwell equation. A second equation for electromagnetic field analysis obtained by applying the recursive convolution method to the Maxwell equation for the TM mode by solving the equation with the computer using the recursive relational expression. The electromagnetic field in the FDTD region is calculated on the basis of the electromagnetic fields obtained for the TM mode and the TE mode by solving the problem with the computer using the inductive relational expression. The simulation program according to claim 6 or 7.

Images