OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 14, Iss. 9 — May. 1, 2006
  • pp: 4151–4168
« Show journal navigation

Cubic interpolated propagation scheme in numerical analysis of lightwave and optical force

Daisuke Barada, Takashi Fukuda, Masahide Itoh, and Toyohiko Yatagai  »View Author Affiliations


Optics Express, Vol. 14, Issue 9, pp. 4151-4168 (2006)
http://dx.doi.org/10.1364/OE.14.004151


View Full Text Article

Acrobat PDF (647 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

A novel technique using a cubic interpolated propagation or constrained interpolation profile (CIP) scheme for numerical analysis of light propagation in dielectric media is proposed. One- and two-dimensional calculations of the propagation of short Gaussian pulses are performed. The validity of the proposed technique is confirmed by applying it to the examination of the reflection from dielectric media. Using the CIP scheme, the optical force acted upon a dielectric disc is also calculated and it is shown that the direction of the calculated force is consistent with the direction predicted from theory.

© 2006 Optical Society of America

1. Introduction

The cubic interpolated propagation or constrained interpolation profile (CIP) method has been investigated as a numerical method for analyzing fluid dynamics [12–19

12. H. Takewaki, A. Nishiguchi, and T. Yabe, “The cubic-interpolated pseudo-particle (CIP) method for solving hyperbolic-type equations,” J. Comput. Phys. 61, 261–268 (1985) [CrossRef]

]. 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

20. Y. Ogata, T. Yabe, and K. Odagaki, “An accurate numerical scheme for Maxwell equation with CIP-method of characteristics,” Commun. Copmut. Phys. 1, 311–335 (2006)

]. 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

×E=μ0Ht,
(1)
×H=εEt,
(2)

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

ε=nε0,
(3)

where n and ε 0 are the refractive index and permittivity of vacuum, respectively. In this paper, we assume that the permeability is constant and equal to μ 0 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

μ0Ht=c×(εE),
(4)
εEt=c×(μ0H),
(5)

where c=1/εμ0 is the speed of light in a medium with refractive index of n. For simplification, we consider the one-dimensional case, E = (0,Ey ,0) and H = (0,0,Hz ), and Eqs. (4) and (5)

μ0Hzt=cεEyx,
(6)
εEyt=cμ0Hzx.
(7)

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

ψ+xt+cψ+xx=0,
(8)
ψxtcψxx=0,
(9)

where

ψ+x=μ0Hz+εEy,
(10)
ψx=μ0HzεEy,.
(11)

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

dxdξ=1n.
(12)

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

ψ+xt+c0ψ+xξ=0,
(13)
ψxtc0ψxξ=0,
(14)

where c0=1/ε0μ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 0, and represented by

ψ+x′/nt+c0ψ+x′/nξ=0,
(15)
ψx′/ntc0ψx′/nξ=0,
(16)
ψ+x′=ψ+xx,
(17)
ψx′=ψxx.
(18)

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

ψ±xtξiψ±x(tΔt,ξic0Δt)Δt=0,
(19)

and rewritten as

ψ±x(t,ξi)=ψ±x(tΔt,ξic0Δt),
(20)

where Δt is the time step and the subscripts i denote the coordinates. The time step Δt must satisfy the condition Δt < Δx/c 0, where Δx is the space interval in the x-axis. The amplitude of the forward wave at ξi - c 0Δt is positioned between ξi = ξ(xi ) and ξ i-1 = ξ(x - Δx). Also, the backward wave at ξi + c 0Δ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:

Ψi±x(ξ)=ai,3±x(ξξi)3+ai,2±x(ξξi)2+ai,1±x(ξξi)+ai,0±x.
(21)

Four coefficients of ai,3±x, ai,2±x, ai,1±x and ai,0±x are obtained by ψi±x, ψi1+x, ψi+x´/n and ψi1+ /n as shown below:

ai,0±x=ψi±x.
(22)
ai,1±x=ψi±x′,
(23)
ai,2±x=3(ψi1±xψi±x)(Δξi±x)22ψi±x′+ψi1±x′Δξi±x,
(24)
ai,3±x=ψi±x′ψi1±x′(Δξi±x)22(ψi1±xψi±x)(Δξi±x)3,
(25)

where

Δξi+x=ξi1ξi,
(26)
Δξix=ξi+1ξi.
(27)

When the refractive index at ξi - c 0Δt is different from the index at ξi , the quantities ψ +x, ψ -x, ψ + and ψ - 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 Ey+x and magnetic field Hz+x are expressed by

Hz+x=ψ+x2μ0,
(28)
Hy+x=ψ+x2ε.
(29)

From Eq. (29), the ratio of Ey+x(ξi ) and Ey+x(ξi - c 0Δt) is obtained

Ey+x(t,ξi)Ey+x(tΔt,ξic0Δt)=ε(ξi)ε(ξic0Δt)=n(ξi)n(ξic0Δt).
(30)

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

ψ+xtΔtξic0Δt=μ0Hz+xtΔtξic0Δt+ε(ξi)Ey+xtΔtξic0Δt.
(31)

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

ψ+xtξi=Ti+xψ+xtΔtξic0Δt,
(32)

where

Ti+x=2n(ξi)n(ξi)+n(ξic0Δt).
(33)

As shown in Eq. (32), Eq. (20) is corrected by multiplying the right side of Eq. (20) by the coefficient Ti+x and the electromagnetic field becomes discontinuous at the interface between two media wiith different refractive indices. In addition, εEy+x is not equal to μ0Hz+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 0Δt) is represented by

ψxtΔtξic0Δt=μ0Hz+xtΔtξic0Δt+ε(ξi)Ey+xtΔtξic0Δt,
(34)

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

ψxtξ=Ri+xψ+xtΔtξic0Δt,,
(35)

where

Ri+x=n(ξi)n(ξic0Δt)n(ξi)+n(ξic0Δt).
(36)

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 ξ(xi + Δx) is needed in practical calculation because the backward wave is calculated by a cubic polynomial between ξi and ξ(xi + Δ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 xi + Δl, the values ξ(xi - Δx) and ξ(xi + Δx) are represented by

ξ(xiΔx)={xi+n(xi)Δln(xi+Δx)(Δx+Δl)Δl<0xin(xi)ΔxΔl0,
(37)
ξ(xi+Δx)={xi+n(xi)Δl+n(xi+Δx)(ΔxΔl)Δl0xi+n(xi)ΔxΔl<0.
(38)

ψi+x(t)=Ti+xψ+xtΔtξic0Δt+RixψxtΔtξi+c0Δt2Δl,
(39)
ψix(t)=TixψxtΔtξi+c0Δt+Ri+xψ+xtΔtξic0Δt+2Δl,
(40)
ψi+x′(t)=Ti+xψ+xtΔtξic0ΔtRixψxtΔtξi+c0Δt2Δl,
(41)
ψix′(t)=TixψxtΔtξi+c0ΔtRi+xψ+xtΔtξic0Δt+2Δl,
(42)

where Tix and Rix are the transmittance and reflectance at position xi respectively, derived from the backward wave and expressed by

Tix=2n(ξi)n(ξi)+n(ξ+c0Δt),
(43)
Rix=n(ξi)n(ξi+c0Δt)n(ξi)+n(ξi+c0Δt).
(44)

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 xi . The boundary between two media with different refractive indices is placed at x = xi + Δ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 xi + Δl with reflectance R i+1.

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 Ey , Hz , ∂Ey /∂x and ∂Hz /∂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

F=PE+μ0Pt×H,
(45)

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

PEx=ε0χ[ExExx+EyExy+EzExz].
(46)

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

Exxyx=Exx+ΔxyExxΔxy2Δx.
(47)

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

Pt=χ×H.
(48)

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:

ε0Eyxt=μ0Hzxt
(49)
=ε0exp[(tt0x/c0)22wt2],
(50)

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

ε0Eyxtx=μ0Hzxtx
(51)
=ε0c0(tt0)xc02wt2exp[(tt0x/c0)22wt2].
(52)

Simulation parameters were set as follows: wt = 1 fs, t 0 = 10wt , Δ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

21. G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations,” IEEE Trans. Electromagn. Compat. 23, 377–382 (1981) [CrossRef]

] 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

19. T. Yabe, H. Mizoe, K. Takizawa, H. Morikia, H. N. Ima, and Y. Ogata, “Higher-order schemes with CIP method and adaptive Soroban grid towards mesh-free scheme,” J. Comput. Phys. 194, 55–77 (2004) [CrossRef]

].

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

ψ+xt+cψ+xx=0,
(53)
ψxtcψxx=0,
(54)

where

ψ+x=[μ0Hz+εEyμ0HyεEz],
(55)
ψx=[μ0HzεEyμ0Hy+εEz].
(56)

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

ψ+yt+cψ+yx=0,
(57)
ψytcψyx=0,
(58)

where

ψ+y=[μ0HzεExμ0Hx+εEz],
(59)
ψy=[μ0Hz+εExμ0HxεEz].
(60)

These equations and their derivatives are solved by the type-M scheme in this paper[13

13. H. Takewaki and T. Yabe, “The cubic-interpolated pseudo particle (CIP) method: application to nonlinear and multi-dimensional hyperbolic equations,” J. Comput. Phys. 70, 355–372 (1987) [CrossRef]

]. 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]
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.

4.2. Gaussian pulse propagation including a dielectric disc

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

Gxyt={exp[(tt0(ŝ·r)/c0)22wt2]tŝr/c00t<ŝr/c0,
(61)

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: wt = 1 fs, t 0 = 5wt , Δ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

22. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994) [CrossRef]

] 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.
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]
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]
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]
Fig. 8. Result of optical field and force calculation. (a) and (b) are results after 20 and 30 fs, respectively.
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.

References and links

1.

K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat. , 14, 302–307 (1966) [CrossRef]

2.

A. Taflove and S. C. Hagness, Computational electrodynamics: the finite-difference time-domain method, 3rd edition, (Artech House, Norwood, MA, 2005)

3.

P. Rochon, E. Batalla, and A. Natansohn, “Optically induced surface gratings on azoaromatic polymer films,” Appl. Phys. Lett. 66, 136–138 (1995) [CrossRef]

4.

D. Y. Kim, S. K. Tripathy, L. Li, and J. Kumar, “Laser-induced holographic surface relief gratings on nonlinear optical polymer films,” Appl. Phys. Lett. 66, 1166–1168 (1995) [CrossRef]

5.

C. J. Barrett, P. L. Rochon, and A. L. Natansohn, “Model of laser-driven mass transport in thin films of dye-functionalized polymers,” J. Chem. Phys. 109, 1505–1516 (1998) [CrossRef]

6.

P. Lefin, C. Fiorini, and J. M. Nunzi, “Anisotoropy of the photoinduced translation diffusion of azo-dyes,” Opt. Mater. 9, 323–328 (1998) [CrossRef]

7.

T. G. Pedersen, P. M. Johansen, N. C. R. Holme, and P. S. Ramanujam, “Mean-field theory of photoinduced formation of surface reliefs in side-chain azobenzene polymers,” Phys. Rev. Lett. 80, 89–92 (1998) [CrossRef]

8.

J. Kumar, L. Li, X. L. Jiang, D. Y. Kim, T. S. Lee, and S. Tripathy, “Gradient force: the mechanism for surface relief grating formation in azobenzene functionalized polymers,” Appl. Phys. Lett. 72, 2096–2098 (1998) [CrossRef]

9.

D. Barada, M. Itoh, and T. Yatagai, “Computer simulation of photoinduced mass transport on azobenzene polymer films by particle method,” J. Appl. Phys. 96, 4204–4210 (2004) [CrossRef]

10.

D. Barada, T. Fukuda, M. Itoh, and T. Yatagai, “Numerical analysis of photoinduced surface relief grating formation by particle method,” Opt. Rev. 12, 271–273 (2005) [CrossRef]

11.

D. Barada, T. Fukuda, M. Itoh, and T. Yatagai, “Proposal of novel model for photoinduced mass transport and numerical analysis by electromagnetic-induced particle transport method,” Jpn. J. Appl. Phys. 45, 465–469 (2006) [CrossRef]

12.

H. Takewaki, A. Nishiguchi, and T. Yabe, “The cubic-interpolated pseudo-particle (CIP) method for solving hyperbolic-type equations,” J. Comput. Phys. 61, 261–268 (1985) [CrossRef]

13.

H. Takewaki and T. Yabe, “The cubic-interpolated pseudo particle (CIP) method: application to nonlinear and multi-dimensional hyperbolic equations,” J. Comput. Phys. 70, 355–372 (1987) [CrossRef]

14.

T. Yabe and E. Takei, “A new higher-order Godunov method for general hyperbolic equations,” J. Phys. Soc. Jpn. 57, 2598–2601 (1988) [CrossRef]

15.

T. Yabe and T. Aoki, “A universal solver for hyperbolic-equations by cubic-polynomial interpolation. I. One-dimensional solver,” Comput. Phys. Commun. 66, 219–232 (1991) [CrossRef]

16.

T. Yabe, T. Ishikawa, P. Y. Wang, T. Aoki, Y. Kadota, and F. Ikeda, “A universal solver for hyperbolic-equations by cubic-polynomial interpolation. II. Two- and three-dimensional solvers,” Comput. Phys. Commun. 66, 233–242 (1991) [CrossRef]

17.

T. Yabe and P. Y. Wang, “Unified numerical procedure for compressible and incompressible fluid,” J. Phys. Soc. Jpn. 60, 2105–2108 (1991) [CrossRef]

18.

T. Yabe, F. Xiao, and T. Utsumi, “Constrained interpolation profile method for multiphase analysis,” J. Comput. Phys. 169, 556593 (2001) [CrossRef]

19.

T. Yabe, H. Mizoe, K. Takizawa, H. Morikia, H. N. Ima, and Y. Ogata, “Higher-order schemes with CIP method and adaptive Soroban grid towards mesh-free scheme,” J. Comput. Phys. 194, 55–77 (2004) [CrossRef]

20.

Y. Ogata, T. Yabe, and K. Odagaki, “An accurate numerical scheme for Maxwell equation with CIP-method of characteristics,” Commun. Copmut. Phys. 1, 311–335 (2006)

21.

G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations,” IEEE Trans. Electromagn. Compat. 23, 377–382 (1981) [CrossRef]

22.

J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994) [CrossRef]

OCIS Codes
(000.3860) General : Mathematical methods in physics
(000.4430) General : Numerical approximation and analysis
(260.2160) Physical optics : Energy transfer

ToC Category:
Physical Optics

History
Original Manuscript: February 9, 2006
Revised Manuscript: April 17, 2006
Manuscript Accepted: April 18, 2006
Published: May 1, 2006

Citation
Daisuke Barada, Takashi Fukuda, Masahide Itoh, and Toyohiko Yatagai, "Cubic interpolated propagation scheme in numerical analysis of lightwave and optical force," Opt. Express 14, 4151-4168 (2006)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-9-4151


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media," IEEE Trans. Antennas Propagat.,  14,302-307 (1966) [CrossRef]
  2. A. Taflove and S. C. Hagness, Computational electrodynamics: the finite-difference time-domain method, 3rd edition, (Artech House, Norwood, MA, 2005)
  3. P. Rochon and E. Batalla and A. Natansohn, "Optically induced surface gratings on azoaromatic polymer films," Appl. Phys. Lett. 66,136-138 (1995) [CrossRef]
  4. D. Y. Kim, S. K. Tripathy, L. Li and J. Kumar, "Laser-induced holographic surface relief gratings on nonlinear optical polymer films," Appl. Phys. Lett. 66,1166-1168 (1995) [CrossRef]
  5. C. J. Barrett, P. L. Rochon and A. L. Natansohn, "Model of laser-driven mass transport in thin films of dyefunctionalized polymers," J. Chem. Phys. 109,1505-1516 (1998) [CrossRef]
  6. P. Lefin, C. Fiorini and J. M. Nunzi, "Anisotoropy of the photoinduced translation diffusion of azo-dyes," Opt. Mater. 9,323-328 (1998) [CrossRef]
  7. T. G. Pedersen, P. M. Johansen, N. C. R. Holme and P. S. Ramanujam, "Mean-field theory of photoinduced formation of surface reliefs in side-chain azobenzene polymers," Phys. Rev. Lett. 80,89-92 (1998) [CrossRef]
  8. J. Kumar, L. Li, X. L. Jiang, D. Y. Kim, T. S. Lee and S. Tripathy, "Gradient force: the mechanism for surface relief grating formation in azobenzene functionalized polymers," Appl. Phys. Lett. 72,2096-2098 (1998) [CrossRef]
  9. D. Barada,M. Itoh and T. Yatagai, "Computer simulation of photoinduced mass transport on azobenzene polymer films by particle method," J. Appl. Phys. 96,4204-4210 (2004) [CrossRef]
  10. D. Barada, T. Fukuda, M. Itoh and T. Yatagai, "Numerical analysis of photoinduced surface relief grating formation by particle method," Opt. Rev. 12,271-273 (2005) [CrossRef]
  11. D. Barada, T. Fukuda, M. Itoh and T. Yatagai, "Proposal of novel model for photoinduced mass transport and numerical analysis by electromagnetic-induced particle transport method," Jpn. J. Appl. Phys. 45,465-469 (2006) [CrossRef]
  12. H. Takewaki, A. Nishiguchi and T. Yabe, "The cubic-interpolated pseudo-particle (CIP) method for solving hyperbolic-type equations," J. Comput. Phys. 61,261-268 (1985) [CrossRef]
  13. H. Takewaki and T. Yabe, "The cubic-interpolated pseudo particle (CIP) method: application to nonlinear and multi-dimensional hyperbolic equations," J. Comput. Phys. 70,355-372 (1987) [CrossRef]
  14. T. Yabe and E. Takei, "A new higher-order Godunov method for general hyperbolic equations," J. Phys. Soc. Jpn. 57,2598-2601 (1988) [CrossRef]
  15. T. Yabe and T. Aoki, "A universal solver for hyperbolic-equations by cubic-polynomial interpolation. I. Onedimensional solver," Comput. Phys. Commun. 66,219-232 (1991) [CrossRef]
  16. T. Yabe, T. Ishikawa, P. Y. Wang, T. Aoki, Y. Kadota and F. Ikeda, "A universal solver for hyperbolic-equations by cubic-polynomial interpolation. II. Two- and three-dimensional solvers," Comput. Phys. Commun. 66,233-242 (1991) [CrossRef]
  17. T. Yabe and P. Y. Wang, "Unified numerical procedure for compressible and incompressible fluid, " J. Phys. Soc. Jpn. 60,2105-2108 (1991) [CrossRef]
  18. T. Yabe, F. Xiao and T. Utsumi, "Constrained interpolation profile method for multiphase analysis," J. Comput. Phys. 169,556593 (2001) [CrossRef]
  19. T. Yabe, H. Mizoe, K. Takizawa, H. Morikia, H. N. Ima and Y. Ogata, "Higher-order schemes with CIP method and adaptive Soroban grid towards mesh-free scheme," J. Comput. Phys. 194,55-77 (2004) [CrossRef]
  20. Y. Ogata, T. Yabe and K. Odagaki, "An accurate numerical scheme for Maxwell equation with CIP-method of characteristics," Commun. Copmut. Phys. 1,311-335 (2006)
  21. G. Mur, "Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations," IEEE Trans. Electromagn. Compat. 23,377-382 (1981) [CrossRef]
  22. J. P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," J. Comput. Phys. 114,185-200 (1994) [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Supplementary Material


» Media 1: GIF (128 KB)     
» Media 2: GIF (631 KB)     
» Media 3: GIF (241 KB)     
» Media 4: GIF (690 KB)     

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited