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

Optics Express, Vol. 14, Issue 9, pp. 4151-4168 (2006)

http://dx.doi.org/10.1364/OE.14.004151

Acrobat PDF (647 KB)

### 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

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]

*λ*/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.

## 2. Theory

### 2.1. Advection equations from Maxwell’s equations

**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*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

*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)

*ψ*

^{+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

*c*

_{0}, and represented by

### 2.2. Numerical method

*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}=

*ξ*(

*x*

_{i}) 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:

*n*and

*x́*/

*n*as shown below:

*ξ*

_{i}-

*c*

_{0}Δ

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

*ψ*

^{-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

*ξ*

_{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.

*x*

_{i}+ Δ

*l*, the values

*ξ*(

*x*

_{i}- Δ

*x*) and

*ξ*(

*x*

_{i}+ Δ

*x*) are represented by

*ξ*(

*x*

_{i}+ Δ

*x*) to

*ξ*

_{i}, and the value Δ

*ξ*

^{-x}is defined as the distance

*c*

_{0}Δ

*t*- Δ

*x*

_{i}- Δ

*l*as shown in Fig. 1. Therefore, the factor of the reflectance derived from the forward wave is multiplied by

*ψ*

^{+x}(

*t*- Δ

*ξ*

_{i}- Δ

*l*) and the value obtained is added to the value of the transmitted backward wave at

*t*- Δ

*x*

_{i}- Δ

*ψ*

^{+x}(

*t*- Δ

*ξ*

_{i}- Δ

*l*) is unknown, it can be obtained using a cubic polynomial derived from known values near the position

*ξ*

_{i}- Δ

*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

*x*

_{i}respectively, derived from the backward wave and expressed by

*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

**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

*∂E*

_{x}/

*∂x*are known. The unknown value

*∂E*

_{x}(

*x,y*)/

*∂x*is obtained by central difference as follows:

## 3. One-dimensional calculation

*w*

_{t}is the Gaussian time width. For CIP calculation, the derivative of this equation is used and represented by

*w*

_{t}= 1 fs,

*t*

_{0}= 10

*w*

_{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

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]

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

*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

*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

*y*-direction are expressed by

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]

### 4.2. Gaussian pulse propagation including a dielectric disc

*G*propagates as expressed by

*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*

_{0}= 5

*w*

_{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.

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

### 5. Conclusion

## 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. , |

2. | A. Taflove and S. C. Hagness, |

3. | P. Rochon, E. Batalla, and A. Natansohn, “Optically induced surface gratings on azoaromatic polymer films,” Appl. Phys. Lett. |

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

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

6. | P. Lefin, C. Fiorini, and J. M. Nunzi, “Anisotoropy of the photoinduced translation diffusion of azo-dyes,” Opt. Mater. |

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

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

9. | D. Barada, M. Itoh, and T. Yatagai, “Computer simulation of photoinduced mass transport on azobenzene polymer films by particle method,” J. Appl. Phys. |

10. | D. Barada, T. Fukuda, M. Itoh, and T. Yatagai, “Numerical analysis of photoinduced surface relief grating formation by particle method,” Opt. Rev. |

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

12. | H. Takewaki, A. Nishiguchi, and T. Yabe, “The cubic-interpolated pseudo-particle (CIP) method for solving hyperbolic-type equations,” J. Comput. Phys. |

13. | H. Takewaki and T. Yabe, “The cubic-interpolated pseudo particle (CIP) method: application to nonlinear and multi-dimensional hyperbolic equations,” J. Comput. Phys. |

14. | T. Yabe and E. Takei, “A new higher-order Godunov method for general hyperbolic equations,” J. Phys. Soc. Jpn. |

15. | T. Yabe and T. Aoki, “A universal solver for hyperbolic-equations by cubic-polynomial interpolation. I. One-dimensional solver,” Comput. Phys. Commun. |

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

17. | T. Yabe and P. Y. Wang, “Unified numerical procedure for compressible and incompressible fluid,” J. Phys. Soc. Jpn. |

18. | T. Yabe, F. Xiao, and T. Utsumi, “Constrained interpolation profile method for multiphase analysis,” J. Comput. Phys. |

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

20. | Y. Ogata, T. Yabe, and K. Odagaki, “An accurate numerical scheme for Maxwell equation with CIP-method of characteristics,” Commun. Copmut. Phys. |

21. | G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations,” IEEE Trans. Electromagn. Compat. |

22. | J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. |

**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: Year | Journal | Reset

### References

- 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]
- A. Taflove and S. C. Hagness, Computational electrodynamics: the finite-difference time-domain method, 3rd edition, (Artech House, Norwood, MA, 2005)
- P. Rochon and E. Batalla and A. Natansohn, "Optically induced surface gratings on azoaromatic polymer films," Appl. Phys. Lett. 66,136-138 (1995) [CrossRef]
- 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]
- 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]
- P. Lefin, C. Fiorini and J. M. Nunzi, "Anisotoropy of the photoinduced translation diffusion of azo-dyes," Opt. Mater. 9,323-328 (1998) [CrossRef]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- T. Yabe and E. Takei, "A new higher-order Godunov method for general hyperbolic equations," J. Phys. Soc. Jpn. 57,2598-2601 (1988) [CrossRef]
- 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]
- 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]
- T. Yabe and P. Y. Wang, "Unified numerical procedure for compressible and incompressible fluid, " J. Phys. Soc. Jpn. 60,2105-2108 (1991) [CrossRef]
- T. Yabe, F. Xiao and T. Utsumi, "Constrained interpolation profile method for multiphase analysis," J. Comput. Phys. 169,556593 (2001) [CrossRef]
- 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]
- 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)
- 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]
- 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.