1. Introduction

Recently, the analysis of electromagnetic fields in various nanostructures has been performed in order to predict unique optical functions or find an optimal structure. The most popular numerical method for the propagation of electromagnetic waves is the finite-difference time-domain (FDTD) method. The FDTD method is suitable for analyzing any optical situation because Maxwell’s equations, which govern the propagation of the electromagnetic wave, are solved directly in this numerical method. However, the time and space coordinate grids are shifted by a half step between the electric and magnetic fields in Yee’s algorithm[1, 2]. Therefore, time and space interpolations are required in order to calculate some physical quantities such as the Poynting vector, irradiance or optical force because these quantities need to be determined for the same time and position for both the electric and magnetic fields. Until now, many authors have investigated the interaction between light and materials and energy transformations from optical to mechanical energy such as photoinduced mass transport[3–11]. In order to analyze the energy transformation phenomena, physical quantities derived from light propagation are important. Therefore a multidisciplinary numerical method is required. In our previous paper, the photoinduced mass transport was analyzed numerically by proposing a novel physical model [9–11]. The numerical results were compared with various experimental results and the photoinduced mass transport was explained qualitatively. However, more accurate numerical methods are required for quantitative analysis.

The cubic interpolated propagation or constrained interpolation profile (CIP) method has been investigated as a numerical method for analyzing fluid dynamics [12–19]. The CIP method is a multiscale numerical method available for any advections such as mass transport, wave propagation and so on. Therefore, the CIP method is applicable for the propagation of electromagnetic waves[20]. In fact, it is easy to apply the CIP method to the propagation of electromagnetic waves in a medium with constant refractive index, for example in vacuum. Moreover, the CIP method has an advantage with regards to numerical diffusion, which occurs in finite differerence methods. However, no journal papers on the propagation of electromagnetic waves inside dielectric media have been published to our knowledge. In the FDTD method, numerical diffusion is caused by the short wavelengths in narrow pulse when the grid resolution, which is generally set to less than λ/10, is not sufficient for the short wavelength components. In order to investigate energy transformation on a molecular scale, it has been noted that the analysis of the force from such shorter pulses is important. In such case, the CIP method is useful because such numerical diffusion is negligible and the wave shape can be maintainted. Moreover, absorption or radiation boundary conditions at the edge of the calculation area are not necessary in the CIP method.

The CIP method has another advantage regarding the calculation area. In the CIP method, it is easy to introduce light sources in the far-field without absorbing boundary conditions. Therefore, the calculation area can be prepared around only a target object. In the FDTD method, the light sources in the far-field can also be introduced by using scattered-field formulation. However, it is necessary to choose a complex absorbing boundary condition in order to neglect the numerical reflection. Utilizing these advantages of the CIP method, more accurate calculation of electromagnetic fields may be performed. In addition, the CIP scheme is suitable for evaluating the optical force that acts on a dielectric medium by the electromagnetic field. In this paper, a novel technique for electromagnetic analysis by the CIP method is proposed and one-and two-dimensional calculations are demonstrated. The results are compared with the results obtained by total- or scattered-field FDTD method.

2. Theory

2.1. Advection equations from Maxwell’s equations

In CIP calculation, the form of Maxwell’s equations, which are the governing equations for electromagnetics, must be changed to advection equations. Maxwell’s equations are represented by

where E, H, ε and μ ₀ are the electric field vector, magnetic field vector, permittivity and permeability of vacuum, respectively. The square root of the permittivity is expressed by

where n and ε ₀ are the refractive index and permittivity of vacuum, respectively. In this paper, we assume that the permeability is constant and equal to μ ₀ because we are only considering optical conditions. By assuming that the permittivity is constant locally, we can rewrite Eqs. (1) and (2) can be rewritten as

where $c=1/\sqrt{\epsilon {\mu }_0}$ is the speed of light in a medium with refractive index of n. For simplification, we consider the one-dimensional case, E = (0,E_y ,0) and H = (0,0,H_z ), and Eqs. (4) and (5)

By adding Eq. (7) to (6) and subtracting Eq. (7) from (6), two advection equations are obtained as follows:

where

Equations (8) and (9) are the advection equations of the forward and backward wave, respectively. The quantities ψ ^+x and ψ ^-x move along the +x-direction and -x-direction, respectively. Here, we define the the optical length ξ as ξ = nx and the derivative with respect to ξ is given by

Using this equation, Eqs. (8) and (9) are rewritten as

where $c_0=1/\sqrt{{\epsilon }_0{\mu }_0}$ is the speed of light in vacuum. The derivatives of Eqs. (13) and (14) are also advection equations with the same fluid velocity c ₀, and represented by

These advection equations are evaluated by the CIP method.

2.2. Numerical method

Solving the advection equations numerically, Eqs. (13) and (14) are discretized as

and rewritten as

where Δt is the time step and the subscripts i denote the coordinates. The time step Δt must satisfy the condition Δt < Δx/c ₀, where Δx is the space interval in the x-axis. The amplitude of the forward wave at ξ_i - c ₀Δt is positioned between ξ_i = ξ(x_i ) and ξ _i-1 = ξ(x - Δx). Also, the backward wave at ξ_i + c ₀Δt is positioned between ξ_i and ξ _i+1 = ξ(x + Δx). The profile between the position ξ_i and neighboring point ξ _i∓1 is expressed by a cubic polynomial Ψ_i^±x(ξ) as follows:

Four coefficients of $a_{i,3}^{\pm x}$, $a_{i,2}^{\pm x}$, $a_{i,1}^{\pm x}$ and $a_{i,0}^{\pm x}$ are obtained by ${\psi }_i^{\pm x}$, ${\psi }_{i\mp 1}^{+x}$, ${\psi }_i^{+\stackrel{\mathit{´}}{x}}$/n and ${\psi }_{i\mp 1}^{+}$ x́/n as shown below:

where

When the refractive index at ξ_i - c ₀Δt is different from the index at ξ_i , the quantities ψ ^+x, ψ ^-x, ψ ^+x́ and ψ ^-x́ become discrete as shown in Eqs. (10) and (11) because the electric field, magnetic field and their derivatives are continuous. The continuous quantities can be obtained by considering the reflection.

Here, the forward wave is considered and the electric field $E_y^{+x}$ and magnetic field $H_z^{+x}$ are expressed by

From Eq. (29), the ratio of $E_y^{+x}$(ξ_i ) and $E_y^{+x}$(ξ_i - c ₀Δt) is obtained

Using Eq. (10), the value ψ ^+x(t - Δt, ξ_i - c ₀Δt) is represented by

Using Eq. (28), (29) and (30), Eq. (31) is rewritten as

where

As shown in Eq. (32), Eq. (20) is corrected by multiplying the right side of Eq. (20) by the coefficient $T_i^{+x}$ and the electromagnetic field becomes discontinuous at the interface between two media wiith different refractive indices. In addition, $\sqrt{\epsilon }{E}_y^{+x}$ is not equal to $\sqrt{{\mu }_0}{H}_z^{+x}$ when the refractive index is inhomogenous. Then, the value ψ ^-x becomes not equal to zero as shown in Eq. (11) and the reflected wave is generated. Using Eq. (11), the reflected wave ψ ^-x(ξ_i - c ₀Δt) is represented by

and using Eq. (30), ψ ^-x(t,ξ) is obtained as follows:

where

Equation (35) corresponds to Fresnel’s reflection formula. The theory mentioned above is accurate when the boundary of the refractive index is at only ξ_i , but invalid at other arbitrary positions. The reflection value at ξ(x_i + Δx) is needed in practical calculation because the backward wave is calculated by a cubic polynomial between ξ_i and ξ(x_i + Δx). In order to avoid these limitations, a reflection model is developed.

The schematic diagram of the reflection model is shown in Fig. 1. When the boundary is at x_i + Δl, the values ξ(x_i - Δx) and ξ(x_i + Δx) are represented by

Here, the value Δ$t_i^{-x}$ is defined as the time that it takes for the backward wave to move from ξ(x_i + Δx) to ξ_i , and the value Δξ ^-x is defined as the distance c ₀Δ$t_i^{-x}$. Then, the position of the forward wave at t - Δ$t_i^{-x}$ becomes x_i - Δ${\xi }_i^{-x}$ + 2Δl as shown in Fig. 1. Therefore, the factor of the reflectance derived from the forward wave is multiplied by ψ ^+x(t - Δ$t_i^{-x}$, ξ_i - Δ${\xi }_i^{-x}$ + 2Δl) and the value obtained is added to the value of the transmitted backward wave at t - Δ$t_i^{-x}$ and x_i - Δ${\xi }_i^{-x}$, which is a known value. Although the value of ψ ^+x(t - Δ$t_i^{-x}$, ξ_i - Δ${\xi }_i^{-x}$ + 2Δl) is unknown, it can be obtained using a cubic polynomial derived from known values near the position ξ_i - Δ${\xi }_i^{-x}$ + 2Δl. Regarding derivatives, the sign of reflectance is opposite because the direction of advection of the reflected wave is reversed in reflection. Finally, the next time step values of four advection equations are represented by

where $T_i^{-x}$ and $R_i^{-x}$ are the transmittance and reflectance at position x_i respectively, derived from the backward wave and expressed by

The right side of Eqs. (39)–(42) are obtained by interpolation in CIP scheme.

 

Fig. 1. Schematic diagram of the numerical model for reflection. Here, the position in optical length ξ_i is equal to x_i . The boundary between two media with different refractive indices is placed at x = x_i + Δl. The solid, dashed and dotted lines are ψ ^+x(t - Δt), ψ ^+x(t) and R _i+1 ψ ^+x, respectively. The value R ^+x ψ ^+x is obtained by multiplying the reflection of ψ ^+x(t) in which the axis of reflection is x_i + Δl with reflectance R _i+1.

Download Full Size | PDF

In the CIP method, more computer storage and arithmetric operations per one cell are required in comparison with the FDTD method or other methods because not only the electric and magnetic fields, but also their derivatives are required. For example, the values of E_y , H_z , ∂E_y /∂x and ∂H_z /∂x are used in one cell at a time in the case of one dimensional propagation. However, the CIP method has third-order accuracy in time and space because the CIP method is based on cubic polynomials, whereas the standard FDTD method has second-order accuracy. Thus, the CIP method has high accuracy although the CIP method has a disadvantage regarding the computer storage.

2.3. Optical force

As was described above, the electric field, magnetic field and their derivatives are solved using the CIP method. Then, these quantities are solved at a single location, while in the FDTD method with Yee’s algorithm the position and time of the electric field and magnetic field are shifted. Therefore, it is easy to obtain other physical quantities such as the Poynting vector, irradiance and optical force. The optical force per unit volume is represented by

where P = ε ₀ χ E is the electric polarization. The coefficient χ = n ² - 1 is the electric susceptibility. The x-component of the first term is represented by

The quantities except for ∂E_x /∂x are known. The unknown value ∂E_x (x,y)/∂x is obtained by central difference as follows:

Using Eq. (2), ∂ P/∂t in the second term of Eq. (45) is calculated by

It is not necessary to evaluate the finite difference to calculate the second term of eq. (45) because all quantities are known. Thus, the optical force is easily obtained.

3. One-dimensional calculation

Using this theory, the propagation of one-dimensional Gaussian pulses was calculated by the CIP method and the result was compared with the result obtained from the total-field FDTD method. In this discussion, a Gaussian pulse was used for the light source:

where w_t is the Gaussian time width. For CIP calculation, the derivative of this equation is used and represented by

Simulation parameters were set as follows: w_t = 1 fs, t ₀ = 10w_t , Δx = 100 nm and Δt = 0.2 fs. Figure 2 shows the numerical results obtained by the CIP and FDTD methods. In this case, the center position of the Gaussian pulse at t = 0 fs is x = -3 μm. The refractive index was set to 1.0 on the left side of the boundary at x = 5 μm and 1.5 on the right side of the boundary as shown in Fig. 2. The calculation range is from x = -5 μm to x = 10 µm in the case of the FDTD method. In the case of the CIP method, the light source can be placed beyond the calculation range, whereas the light source must be placed within the calculation range in the case of the total-field FDTD method. Therefore, the calculation range can be limited to the range x = 0 μm to x = 10 μm in the case of the CIP method although the position of the pulse at t = 0 is to the left of 0 μm. Then, the pulse does not appear at t = 0 as shown in Fig. 2(a). The pulse is generated at x = 0 μm and appears at t = 10 fs as shown in Fig. 2(b). Then, the positon of the pulse obtained by the CIP method coincides with the position of the pulse obtained by the FDTD method. In the FDTD method, the numerical diffusion increased with propagation, the high frequency components were gradually delayed and the numerical reflection at the right side of the calculation area was eliminated by the Mur’s absorbing boundary condition[21] as shown in Fig. 2(f). In the CIP method, numerical diffusion did not occur and the amplitude reflectance calculation agreed with theory (R = -0.2). No numerical reflection at the left and right edge of the calculation is observed even without absorbing boundary conditions. In this condition, the computational resources used by the CIP and FDTD methods are about 6.5 and 2.5 Kbytes, respectively. The running time of both methods at 100 fs on a 1.5 GHz Pentium M processor was about 0.06 s. No difference for the running time was observed between two methods in this condition.

In order to confirm the validity of our proposed numerical model in further, the result with modified boundary position was found and is shown in Fig. 3. The time difference of each pulse was found to be accurate in spite of a position shift of less than Δx. These results show that our proposed model is valid for analyzing media with complicated surfaces using Cartesian coordinates. However, when there are two or more interfaces between two grid points, our proposed technique cannot be used. In that case, the Soroban grid may be useful[19].

4. Two-dimensional propagation

4.1. Advection equations

For two-dimensional calculation, we consider optical waves propagating in the x-y plane. In the CIP scheme, the advection direction is decomposed into x and y components. Using Eqs. (1) and (2), the advection equations for the x-direction are expressed by

where

Also, the advection equations for the y-direction are expressed by

where

These equations and their derivatives are solved by the type-M scheme in this paper[13]. Each advection equation can be solved by the technique proposed in Sec. 2.

 

Fig. 2. Comparison of Gaussian pulse propagation. The solid red and dashed green lines are the pulse shapes obtained by the CIP and FDTD methods, respectively. (a), (b), (c) and (d) show the y-components of the electric field after 0 s, 10 s, 30 s and 50 s, respectively. (e) shows the result after 50 s in the range from 8 μm to 10 μm. [Media 1]

Download Full Size | PDF

 

Fig. 3. Time difference of Gaussian pulse propagation by various boundary position. The red, green, blue, magenta and cyan lines are the profile of the Gaussan pulse when the shift amounts of the boundary position are 0 nm, 10 nm, 20 nm, 50 nm and 100 nm, respectively.

Download Full Size | PDF

4.2. Gaussian pulse propagation including a dielectric disc

In the two-dimensional calculation, a Gaussian pulse G propagates as expressed by

and irradiates a dielectric disc as shown in Fig. 5. Here, s is the propagation vector and the angle between the vector and the y-axis was set to 30°. in this simulation. The polarization state of incident pulse is set to TM or TE polarization. The calculation area was set to 5 μm×5 μm and the dielectric disc with radius 2 μm was placed with its central point at (2.5, 2.5) μm. The refractive indices of the dielectric disc and the ambient medium were set to 1.5 and 1,0, respectively. The simulation parameters were as follows: w_t = 1 fs, t ₀ = 5w_t , Δx = Δy = 100 nm and Δt = 0.2 fs. The radiation layer, which is used to interface with the external electromagnetic field, was placed on the edge of the calculation area. The radiation layer values are determined by the electromagnetic field from the incident pulse propagating through the ambient medium. The incident pulse used here is infinitely broad along the direction perpendicular to the propagation direction and the infinite width of the pulse can be expressed by the radiation layer.

The calculation results obtained by our proposed technique was validated by comparing with the result by scattered-field FDTD method. In the case of FDTD method, Mur’s second order absorbing boundary condition was used. The result with TM polarization after 20 fs is shown in Fig. 5. In this condition, quite similar results were obtained between two methods. As shown in Fig. 5(a), the electric field distribution becomes continuous on the boundary and no discrete profile nor numerical diffusion were observed. The calculations of TE polarzation case were also demonstrated by our proposed technique and scattered-fied FDTD method and the magnetic field distribution is shown in Fig. 6. In the case of the FDTD method, slight numrical reflection was observed as shown in Fig. 6(b). In contrast, no numerical reflection was observed in the case of the CIP method as shown in Fig. 6(b). Figure 7, 8 and 9(a) show the propagation of the pulse with TM polarization and the force vectors acted on the dielectric disc in the case of the CIP method. The pulse was gradually focused by the dielectric disc. After 8 fs, the incident pulse appeared in the calculation area as shown in Fig. 7(a). After 12 fs, the reflected pulse from the dielectric disc was observed as shown in Fig. 7(b). After 20 fs, the reflected pulse observed at 12 fs moved out of the calculation area and no numerical reflections at the edge of the calculation area were observed as shown in Fig. 8(a). After 30 fs, a part of the pulse moved out of the calculation area and again no numerical reflections were observed as shown in Fig. 8(b). After 40 fs, the incident light is out of the calculation area and only the light reflected at the edge of the disc is remained. The result after 40 fs was compared with the result obtained by the FDTD method. As shown in Fig. 9, the numerical reflection was observed in the case of the FDTD method, whereas no numerical reflection was observed in the case of the CIP method. The direction of the optical force acted on the dielectric disc was from the high field region to the low field region. This result is consistent with the theoretical result predicted from Eq. (45). By the numerical result, the validity of our numerical model was confirmed.

In both TM and TE polarzation cases, the computational resourses of the CIP method and FDTD method were 27.5 Kbytes and 7.9 Kbytes, respectively. The running times of the CIP method and FDTD method at 60 fs on 1.5 GHz Pentium M processor were about 2 s and 1 s, respectively. The computational resources and running times were almost the same between the cases of TE and TE polarizations. In this simulation, an adequate result could be obtained irrespective of the limited size of the calculation area that was only slightly larger than the size of the dielectric disc. Of course, the numerical reflection can be reduced by more complex absorbing boundary conditions such as perfectly matched layer[22] in the FDTD method. Howerver, the optimal boundary condition must be chosen in case by case. In the CIP method, no absorbing boundary condition is necessary. This is an advantage of the CIP method. Reflected waves are not generated even if the radiation layer is placed at the edge of the calculation area. Therefore, this numerical method is suitable for analyzing the local field and force.

5. Conclusion

In this paper a novel technique for the calculation of light propagation including the case of dielectric media was proposed. The reflection of a Gaussian pulse at the interface between media with different refractive indices was simulated. The CIP method has an advantage in that no absorbing boundary condition is necessary at the edge of the calculation area. Although a lot of memory is required per grid point in this method, the calculation area can be reduced by calculating a local area because the reflection at the edge of calculation area does not occur in principle. Also, it is easy to obtain other physical quantities, for example optical force, because the values of the electric field, magnetic field and its derivatives are positioned in the same grid. Our proposed method can be easily extended to three-dimensional systems although only one- and two-dimensional calculations were demonstrated in this paper. We believe that the CIP method with our proposed numerical model might be useful in the future for special situations such as the propagation of super short pulses, local electromagnetic analysis and energy transformation.

 

Fig. 4. Schematic diagram of the two-dimensional light propagation case. A dielectric disc with two-micron radius is placed at the center of the calculation area. A Gaussian pulse originates from the bottom of the figure and the angle between the propagation direction and the y-axis is 30°. The amplitude of the electromagnetic field is constant along the direction perpendicular to the propagation direction. The values of the radiation layer placed at the edge of the calculation area are determined by the electromagnetic field from the Gaussian pulse propagating through the medium with n = 1.

Download Full Size | PDF

 

Fig. 5. Electric field distribution obtained by the numerical result of the propagation of the Gaussian pulse with TM polarization after 20 fs. (a) and (b) are the results obtained by the CIP and FDTD methods, respectively. [Media 2]

Download Full Size | PDF

 

Fig. 6. Magnetic field distribution obtained by the numerical result of the propagation of the Gaussian pulse with TE polarization after 20 fs. (a) and (b) are the results obtained by the CIP and FDTD methods, respectively. [Media 3]

Download Full Size | PDF

 

Fig. 7. Result of optical field and force calculation. (a) and (b) are results after 8 and 12 fs, respectively. Dotted circle shows dielectric disc and arrows are the directions of optical force. The length of the arrows is the strength of the optical force. [Media 4]

Download Full Size | PDF

 

Fig. 8. Result of optical field and force calculation. (a) and (b) are results after 20 and 30 fs, respectively.

Download Full Size | PDF

 

Fig. 9. Result of optical field and force calculation after 40 fs. (a) and (b) are the results obtained by the CIP and FDTD methods, respectively. The optical force calculation is performed by only CIP method.

Download Full Size | PDF