- Email: [email protected]
Coordinate transformation and coordinate mapping for numerical simulation of Z-scan measurements
Optics Communications 244 (2005) 71–77 www.elsevier.com/locate/optcom
Coordinate transformation and coordinate mapping for numerical simulation of Z-scan measurements Wei-Ping Zang, Jian-Guo Tian *, Zhi-Bo Liu 1, Wen-Yuan Zhou, Feng Song, Chun-Ping Zhang Photonics Center, College of Physics, Nankai University, Tianjin 300071, China Received 1 April 2004; received in revised form 31 August 2004; accepted 9 September 2004
Abstract We apply a new method consisting of a coordinate transformation and a coordinate mapping technique, to simulate numerically Z-scan measurements of a nonlinear medium. A comparison of numerical results obtained by this new method with those by conventional numerical method is made. Results show that the new scheme yields a large, orders-of-magnitude range improvement in computational speed and accuracy for the analyses of closed-aperture and open-aperture Z-scan measurements, compared with conventional numerical methods. 2004 Elsevier B.V. All rights reserved. PACS: 42.65.k; 42.65.Hw; 42.25.Bs; 45.10.Db Keywords: Numerical simulation; Coordinate transformation; Coordinate mapping; Z-scan
1. Introduction The Z-scan technique has become a popular tool [1,2] for measurements of both amplitude and sign of nonlinear index and nonlinear absorption. The beam propagation in thin nonlinear media has been analyzed and modeled extensively. *
Corresponding author. E-mail addresses: [email protected] (J.-G. Tian), [email protected] (Z.-B. Liu). 1 Tel.: +862223507207; fax: +862223508379.
But the analyses of beam propagation in thick nonlinear media are usually more complicated and difficult. Various analytic approaches have been applied to such analyses [3–11], some of them involve an aberrationless approximation [3–7], and the results are correct to first order for irradiance; the others are approximate correct to higher orders for irradiance by using some approximate methods [8–11]. In higher-power regimes, numerical techniques are usually applied [12,13]. But for conventional numerical methods, a significant amount of computer CPU time is always required.
0030-4018/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2004.09.033
72
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
Therefore, the improvement of the numerical efficiency is one of important tasks for numerical analysis of Z-scan characteristics. In this paper, a new numerical method is applied to simulate numerically Z-scan measurements of nonlinear media. This method consists of a coordinate transformation and a coordinate mapping technique. Compared with conventional numerical methods, the method can enhance the computational speed and accuracy of Z-scan at least two orders of magnitude.
2. Methods 2.1. Coordinate transformation
Kerr coefficient and the two-photon absorption as follows: RefvNL ðr; tÞg ¼ 2n0 n2 jEj2 ; ImfvNL ðr; tÞg ¼
n0 b2 2 jE j ; k0
ð2aÞ ð2bÞ
where n2 is Kerr coefficient, b2 is two-photon absorption coefficient, n0 is linear refractive index of medium and k0 is wave vector in vacuum. Assuming a TEM00 Gaussian beam traveling in the +z-direction, we can write E as following: w0 Eðr; z; tÞ ¼ E0 ðtÞ wðzÞ r2 ikr2 þ itan1 ðz=z0 Þ ; exp 2 w ðzÞ 2RðzÞ ð3Þ
The coordinate transformation we selected is based on the vacuum behavior of a focused Gaussian beam [14]. Some other authors used similar systems [15,16] to enhance both the computational speed and accuracy in analysis of laser resonator cavity modes and beam propagating in atmosphere. The basic concept of the coordinate transformation is to calculate the amplitude changes induced in the vacuum rather than to calculate the complete amplitude. Therefore, this procedure removes the necessity for sampling the high frequency oscillation in the phase induced by focusing. The coordinate transformation alters the independent variables so that the dependent variables take a different functional form, but the dependent variables are numerically identical to the untransformed amplitude at equivalent points in the space of the two coordinate systems. The paraxial wave equation in the Fresnel approximation can be written as 1 o oE oE r ika0 E 2ik r or or oz 2 þ k 0 vNL ðr; z; tÞ E ¼ 0; ð1Þ where vNL is the nonlinear susceptibility of the medium, and a0 is linear absorption coefficient. For a nonlinear medium with Kerr effect and two-photon absorption, real and imaginary parts of the nonlinear susceptibility are related to the
2
¼ w20 ð1 þ z2 =z20 Þ is the z20 =z2 Þ is the radius of
beam radius, where w ðzÞ RðzÞ ¼ zð1 þ curvature of the wave-front at z, w0 is the radius of beam at the focus and z0 ¼ kw20 =2 is the Rayleigh length of the beam, k = 2p/k is the wave vector. E0(t) denotes the electric field at the focus and contains the temporal envelope of the laser pulse. The existence of the quadratic term in the imaginary part of the exponential will cause a high-frequency oscillation in the amplitude for r > w(z), which can exceed the sampling frequency of the computational mesh. Additionally, the amplitude will increase as z ! 0 due to the existence of w0/w(z). As the beam focuses and becomes smaller, the mathematical mesh should shrink at the same rate to model the transverse distribution accurately. These considerations result in the following sequences of coordinate transformations: Eðr; zÞ E0 ðtÞC ðr; zÞAð~r; ~zÞ; w0 iz 2 1 r þ itan ðzÞ ; exp C ðr; zÞ ¼ 1 þ z2 wðzÞ d~z 1 1 ¼ ; ¼ dz 1 þ z2 DðzÞ
ð4aÞ ð4bÞ
ð4cÞ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where z ¼ z=z0 ; r ¼ r=w0 and ~r ¼ r= ð1 þ z2 Þ. The phase of function Að~r; ~zÞ changes very slowly
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
with ~r. The quadratic term in the imaginary part of the exponential in Eq. (4b), which will cause a high-frequency oscillation in the amplitude for r > w(z), and the term w0/w(z) have been included in the function Cðr; zÞ. Function Cðr; zÞ can be solved analytically. If the function Að~r; ~zÞ is obtained, we can get the electric field through Eq. (4a). In addition, We can see clearly from Eq. (4c) that alteration of step size automatically takes smaller steps in the region of high beam intensity. Rewriting Eq. (1) for the function Að~r; ~zÞ in these new coordinates, the following equations will be obtained: 1 o oA oA ~r pðAÞA ¼ 0; þ 4 1 ~r2 A 4i ~r o~r o~r o~z ð5aÞ pðAÞ ¼ 2iDðzÞ a0 þ
4DU0 ðtÞjAj2 2iQ0 ðtÞjAj2 ; n0 ð5bÞ
where DU0(t) = kn2I0(t)z0 is the phase shift at focus for the nonlinear medium with one Rayleigh length, Q0(t) = b2I0(t)z0 is the nonlinear absorption at focus for the nonlinear medium with one Rayleigh length and a0 ¼ a0 z0 is the linear absorption within one Rayleigh length. For the case of p = 0, which corresponds to the beam propagation in vacuum, the function Að~r; ~zÞ has the following form: Að~r; ~zÞ ¼ expð~r2 Þ:
equ ¼ L
Z
L
z2 ¼0
d~zðz1 ; z2 Þ ¼
73
Z 0
L
dz2 1 þ ðz1 þ z2 Þ2
tan1 ðz1 Þ; ¼ tan1 ðz1 þ LÞ
ð7Þ
where z ¼ z=z0 ¼ z1 þ z2 (where z1 is the separation of the front face of the medium from the waist and z2 is the propagation distance from the front face), ¼ L=z0 and L equ is the dimensionless quantity. L It can be seen from Eq. (7) that the equivalent medium length in new coordinate system is the function of the front surface position of medium and the medium length. The equivalent medium length as a function of the position of medium is shown in Fig. 1. We can see from Fig. 1 that the equ is the function of equivalent medium length L the position of medium and medium length L. When the focal point of beam is located in middle equ obtain its maxof medium, for fixed length L; L imum. The equivalent length Lequ is approximately proportional to T1 1 for an open Z-scan [8]. equ increases slowly with medium length L and deL crease quickly as the distance between the middle of medium and focal point of beam increases, and if medium length L ! 1 and the focal point of beam is located in middle of medium, equ ¼ 2tan1 ðL=2Þ L is equal to p. So the length of medium in new coordinates is shrunken greatly. Eqs. (1) and (5) are parabolic partial differential
ð6Þ
It means that the form of the function Að~r; ~zÞ in vacuum keeps unchanged while the beam propagates in this coordinate system. Therefore, the use of the transformations in Eqs. (4) can significantly reduce the size of the data arrays and the computational time, since fewer points can express the radial field distribution and the propagation distance is shrunken greatly. The Z-scan experimental arrangement is the same as that in [2]. If the front surface of the medium is located at z and the medium length is L, one equ by can obtain the equivalent medium length L Eq. (4c), which is defined as the distance between the front surface and the exit surface of medium in this new coordinate system,
equ Fig. 1. The relation between the equivalent medium length L and the position of medium x. The real medium length is L = 1, 2, 4, respectively. The unit of medium length is Rayleigh length z0.
74
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
equations and can be solved numerically by Crank–Nicholson finite difference scheme [17]. Here we call the numerical simulation of Eq. (1) by Crank–Nicholson finite difference scheme as conventional method (CM), and the numerical simulation of Eq. (1) by Crank–Nicholson finite difference scheme as coordinate transformation method (CTM).
It is seen that when i = N, we can get S =2. However, for other value of i, there zN ¼ L are two values of zi . We can use Eqs. (4) to get the real electric q field Eðr; zi Þ ffiby Aði; ~rÞ; z ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S and r ¼ ~r 1 þ ðzi þ L S Þ2 . zi þ L
3. Results and discussion 2.2. Coordinate mapping Below we apply coordinate mapping method (CMM) to numerical simulations of Z-scan measurements. It can be seen that there is no explicit dependence on the medium position in the Eq. (5a). Therefore, the function Að~r; ~zÞ, solved from Eq. (5a), is not the explicit function of the medium position. The function Að~r; ~zÞ is only determined by the equivalent medium length that is the function of medium position and medium length. For a nonlinear medium with fixed length, different medium position in z-direction corresponds to different equivalent medium length. Using this rule, we can obtain the function Að~r; ~zÞ at the exit surface of the medium for different medium position by only computing beam propagation in the medium when the position of center point of the medium is located at origin point. The procedure is described below: S ¼ L=z0 , First, assuming the medium length is L when the position of front surface of medium is loS =2 the equivalent medium length is cated at L 1 Lequ ¼ 2tan ðLS =2Þ, which corresponds to the largest equivalent medium length. Second, slicing the medium into many thin slices, the beam propagates through each slice, and the function Aði; ~rÞ can be obtained, here Aði; ~rÞ denotes the function Að~r; ~zÞ after the beam propagation over ith slice, and is also Að~r; ~zÞ at the exit surface of medium when the medium is located at zi which can be determined by Eq. (7) as following: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 " #ffi9 u = 1< u L S 2S 4 1 2i 1 ; zi ¼ LS tL S =2Þ ; 2: tan N tan ðL ð8Þ
where N is the number of slices.
In order to check the range of validity as well as the accuracy of CM and CTM, we simulate Gaussian beam propagation in vacuum. The main parameters of the Gaussian beam are calculated numerically and then compared with the analytical values. The beam size which is evaluated as the second moment of intensity, can be written as: " R1 #1=2 2p 0 r2 Iðr; zÞr dr R1 wðzÞ ¼ : ð9Þ 2p 0 Iðr; zÞr dr The conditions of numerical computing are listed as following: the step size in r direction is Dr = w0/100, and the step sizes in z direction are equ ¼ 1= DL = z0/100 for CM and D~r ¼ 1=100; DL 100 for CTM, respectively. For Gaussian beams, the size of the numerical window has to be at least four times of the beam size w(z), in order that the reflection from the boundary can be neglected completely. Therefore, the size of the numerical window is selected to be four times of the beam size w(z) or w0 for CM or CTM. In our computations we used a PC with a 2.4-GHz Intel Pentium CPU and 512M DDR memory. Matlab was used for numerical computing. The calculated results of beam propagation over a distance of 20z0 are shown in Fig. 2. It is seen that CTM enhances the accuracy of computing at least one order of magnitude for the above computational conditions as CMs numerical computing. The comparison between computational time of two methods for different medium length is given in Fig. 3. Results show that CTM enhances the speed of numerical computing one or two orders of magnitude for the same conditions as CMs numerical computing. The reason of enhancement of accuracy and speed of numerical computing come from the size of the data arrays
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
75
Typical closed-aperture Z-scan curves without two-photon absorption and open-aperture Z-scan curves with two-photon absorption are given in Fig. 4, where the medium length LS = 10z0. The solutions from CTM and CMM agree very well. The computational time as a function of medium length for three methods is given in Fig. 5. The results show that the numerical efficiency of CTM enhances at least two orders of magnitude more than that of CM, and the numerical efficiency of
Fig. 2. The relative error ew of beam radius of the Gaussian beam propagation through the air for CM and CTM.
Fig. 4. Closed-aperture and open-aperture Z-scan curves. The medium length is 10z0. Solid line represents CTM, and circle represents CMM.
Fig. 3. Computational time taken by CM and CTM for different propagation distance. The propagation distance is scaled by Rayleigh length z0, and the temporal unit is second.
significantly reduced and fewer points can express the radial field distribution and the propagation distance shrunken greatly. Further, we compare the numerical simulations of Z-scan measurements obtained by CMM, CTM and CM, respectively. The parameters used in our calculations are listed as followings [2]: linear refractive index n0 = 1, Kerr coefficient n2r = 6.8 · 1014 cm2/GW, two-photon-absorption coefficient b2r = 5.8 cm/GW, beam waist w0r = 25 lm, and peak intensity I0 = 0.21 GW/cm2. All computations involve only a cw laser beam.
Fig. 5. Computational time taken by CM, CTM and CMM for numerical simulation of Z-scan measurements. The medium length is scaled by Rayleigh length z0, and the temporal unit is second. Circle, CM; square, CTM; triangle, CMM.
76
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
(a)
and each slice propagates as if it were a continuous wave (CW) laser beam. Second task is to accumulate all effects of time slices. The computational time is determined by the number of time slices, the total CPU time is obtained by multiplying the CPU time of cw laser beam and the number of time slices. For Gaussian-shape pulse and hyperbolic-sech shape pulse, the peak-to-valley normalized transmittance difference DTp–v of closed-aperture Z-scan curves without two-photon absorption as a function of medium length are given in Fig. 6(a) and (b). It can be seen that the ratio of DTp–v of Gaussian-shape pulse to that of pffiffiffi CW laser beam is 1= 2, and the ratio DTp–v of hyperbolic-shape pulse to that of CW laser beam is 2/3 in the case of ‘‘thin’’ medium [18]. The ratio decreases as medium length increases, finally reaches the value 0.676 for Gaussian-shape pulse and 0.64 for hyperbolic-shape pulse.
4. Conclusion
(b) Fig. 6. (a) The peak–valley normalized transmittance difference DTp–v vs. medium length L for CW laser beam, Gaussian-shape and hyperbolic-sech shape pulsed laser beam. (b) The peak– valley normalized transmittance difference DT p–v =DT CW p–v vs. medium length L for Gaussian-shape and hyperbolic-sech shape pulsed laser beam. DT CW p–v is the peak–valley normalized transmittance difference of CW laser beam.
CMM enhances at least one order of magnitude more than that of CTM. Therefore, compared with CM, the total numerical efficiency enhances over three orders magnitude. All the discussions above involve only the cw laser beam. For pulsed laser beam with the temporal envelope of Gaussian-shape or hyperbolic-sech shape, Eqs. (5) are hold correct if we ignore the effect of dispersion. The numerical simulation of laser pulse propagation in a thick nonlinear media can be split into two separate tasks. The first task is to split laser pulse into a number of time slices,
We apply CTM to improve the numerical efficiency of beam propagation and CMM to improve the numerical efficiency of numerical simulation of Z-scan measurements. The numerical analyses show that the numerical efficiency of the analyses on beam propagation enhances at least two orders of magnitude and total numerical efficiency of numerical simulation of Z-scan characteristics enhances at least three orders of magnitude, compared with that of CM. For CW laser beam and pulsed laser beam, we applied these new methods to analyze the peak-to-valley normalized transmittance difference for closed-aperture Z-scan as the function of medium length and nonlinear refraction of medium.
Acknowledgements This research is supported by Project 60025512 supported by the National Natural Science Foundation of China, Natural Science Foundation of Tianjin (Grant 0436012 11), a key project of the Ministry of Education (Grant 00026), the Foundation for University Key Teachers of the Ministry
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
of Education, and the Fok Ying Tung Education Foundation (Grant 71008). References [1] M. Sheik-Bahae, A.A. Said, E.W. Van Stryland, Opt. Lett. 14 (1989) 955. [2] M. Sheik-Bahae, A.A. Said, T.H. Wei et al., IEEE. J. Quantum Electron. 26 (1990) 760. [3] M. Sheik-Bahae, A.A. Said, T.H. Wei et al., Opt. Eng. 30 (1991) 1228. [4] J.A. Hermann, R.G. McDuff, J. Opt. Soc. Am. B 10 (1993) 2056. [5] W.-P. Zang, J.-G. Tian, C.-P. Zhang et al., Acta Phys. Sin. 43 (1994) 476 (in Chinese). [6] J.-G. Tian, W.-P. Zang, G.Y. Zhang et al., Acta Phys. Sin. 43 (1994) 1712 (in Chinese). [7] J.-G. Tian, W.-P. Zang, G.-Y. Zhang et al., Appl. Opt. 34 (1995) 4331.
77
[8] J.A. Hermann, J. Opt. Soc. Am. B 14 (1997) 814. [9] W.-P. Zang, J.-G. Tian, Z.-B. Liu et al., Opt. Lett. 28 (2003) 722. [10] W.-P. Zang, J.-G. Tian, Z.-B. Liu et al., J. Opt. Soc. Am. B 32 (2004) 63. [11] P.B. Chapple, J.A. Hermann, R.G. McDuff, Opt. Quantum Electron. 31 (1999) 555. [12] S. Hughes, J.M. Burzer, G. Spruce, B.S. Wherrett, J. Opt. Soc. Am. B 12 (1995) 1888. [13] D.I. Kvosh, S. Yang, D.J. Hagan et al., Appl. Opt. 38 (1999) 5168. [14] P.B. Ulrich, Naval Research Laboratory Report 7706, May 29, 1974. [15] E.A. Sziklas, A.E. Siegman, Appl. Opt. 14 (1975) 1874. [16] J.A. Fleck, J.R. Morris, M.D. Feit, Appl. Phys. 10 (1976) 129. [17] C.F. Gerald, P.O. Wheatley, Applied Numerical Analysis, Addison Wesley, Reading, MA, 1984 (Chapter 8). [18] R.L. Sutherland, Handbook of Nonlinear Optics, Marcel Dekker, New York, 1996 (Chapter 7).
Coordinate transformation and coordinate mapping for numerical simulation of Z-scan measurements Wei-Ping Zang, Jian-Guo Tian *, Zhi-Bo Liu 1, Wen-Yuan Zhou, Feng Song, Chun-Ping Zhang Photonics Center, College of Physics, Nankai University, Tianjin 300071, China Received 1 April 2004; received in revised form 31 August 2004; accepted 9 September 2004
Abstract We apply a new method consisting of a coordinate transformation and a coordinate mapping technique, to simulate numerically Z-scan measurements of a nonlinear medium. A comparison of numerical results obtained by this new method with those by conventional numerical method is made. Results show that the new scheme yields a large, orders-of-magnitude range improvement in computational speed and accuracy for the analyses of closed-aperture and open-aperture Z-scan measurements, compared with conventional numerical methods. 2004 Elsevier B.V. All rights reserved. PACS: 42.65.k; 42.65.Hw; 42.25.Bs; 45.10.Db Keywords: Numerical simulation; Coordinate transformation; Coordinate mapping; Z-scan
1. Introduction The Z-scan technique has become a popular tool [1,2] for measurements of both amplitude and sign of nonlinear index and nonlinear absorption. The beam propagation in thin nonlinear media has been analyzed and modeled extensively. *
Corresponding author. E-mail addresses: [email protected] (J.-G. Tian), [email protected] (Z.-B. Liu). 1 Tel.: +862223507207; fax: +862223508379.
But the analyses of beam propagation in thick nonlinear media are usually more complicated and difficult. Various analytic approaches have been applied to such analyses [3–11], some of them involve an aberrationless approximation [3–7], and the results are correct to first order for irradiance; the others are approximate correct to higher orders for irradiance by using some approximate methods [8–11]. In higher-power regimes, numerical techniques are usually applied [12,13]. But for conventional numerical methods, a significant amount of computer CPU time is always required.
0030-4018/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2004.09.033
72
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
Therefore, the improvement of the numerical efficiency is one of important tasks for numerical analysis of Z-scan characteristics. In this paper, a new numerical method is applied to simulate numerically Z-scan measurements of nonlinear media. This method consists of a coordinate transformation and a coordinate mapping technique. Compared with conventional numerical methods, the method can enhance the computational speed and accuracy of Z-scan at least two orders of magnitude.
2. Methods 2.1. Coordinate transformation
Kerr coefficient and the two-photon absorption as follows: RefvNL ðr; tÞg ¼ 2n0 n2 jEj2 ; ImfvNL ðr; tÞg ¼
n0 b2 2 jE j ; k0
ð2aÞ ð2bÞ
where n2 is Kerr coefficient, b2 is two-photon absorption coefficient, n0 is linear refractive index of medium and k0 is wave vector in vacuum. Assuming a TEM00 Gaussian beam traveling in the +z-direction, we can write E as following: w0 Eðr; z; tÞ ¼ E0 ðtÞ wðzÞ r2 ikr2 þ itan1 ðz=z0 Þ ; exp 2 w ðzÞ 2RðzÞ ð3Þ
The coordinate transformation we selected is based on the vacuum behavior of a focused Gaussian beam [14]. Some other authors used similar systems [15,16] to enhance both the computational speed and accuracy in analysis of laser resonator cavity modes and beam propagating in atmosphere. The basic concept of the coordinate transformation is to calculate the amplitude changes induced in the vacuum rather than to calculate the complete amplitude. Therefore, this procedure removes the necessity for sampling the high frequency oscillation in the phase induced by focusing. The coordinate transformation alters the independent variables so that the dependent variables take a different functional form, but the dependent variables are numerically identical to the untransformed amplitude at equivalent points in the space of the two coordinate systems. The paraxial wave equation in the Fresnel approximation can be written as 1 o oE oE r ika0 E 2ik r or or oz 2 þ k 0 vNL ðr; z; tÞ E ¼ 0; ð1Þ where vNL is the nonlinear susceptibility of the medium, and a0 is linear absorption coefficient. For a nonlinear medium with Kerr effect and two-photon absorption, real and imaginary parts of the nonlinear susceptibility are related to the
2
¼ w20 ð1 þ z2 =z20 Þ is the z20 =z2 Þ is the radius of
beam radius, where w ðzÞ RðzÞ ¼ zð1 þ curvature of the wave-front at z, w0 is the radius of beam at the focus and z0 ¼ kw20 =2 is the Rayleigh length of the beam, k = 2p/k is the wave vector. E0(t) denotes the electric field at the focus and contains the temporal envelope of the laser pulse. The existence of the quadratic term in the imaginary part of the exponential will cause a high-frequency oscillation in the amplitude for r > w(z), which can exceed the sampling frequency of the computational mesh. Additionally, the amplitude will increase as z ! 0 due to the existence of w0/w(z). As the beam focuses and becomes smaller, the mathematical mesh should shrink at the same rate to model the transverse distribution accurately. These considerations result in the following sequences of coordinate transformations: Eðr; zÞ E0 ðtÞC ðr; zÞAð~r; ~zÞ; w0 iz 2 1 r þ itan ðzÞ ; exp C ðr; zÞ ¼ 1 þ z2 wðzÞ d~z 1 1 ¼ ; ¼ dz 1 þ z2 DðzÞ
ð4aÞ ð4bÞ
ð4cÞ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where z ¼ z=z0 ; r ¼ r=w0 and ~r ¼ r= ð1 þ z2 Þ. The phase of function Að~r; ~zÞ changes very slowly
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
with ~r. The quadratic term in the imaginary part of the exponential in Eq. (4b), which will cause a high-frequency oscillation in the amplitude for r > w(z), and the term w0/w(z) have been included in the function Cðr; zÞ. Function Cðr; zÞ can be solved analytically. If the function Að~r; ~zÞ is obtained, we can get the electric field through Eq. (4a). In addition, We can see clearly from Eq. (4c) that alteration of step size automatically takes smaller steps in the region of high beam intensity. Rewriting Eq. (1) for the function Að~r; ~zÞ in these new coordinates, the following equations will be obtained: 1 o oA oA ~r pðAÞA ¼ 0; þ 4 1 ~r2 A 4i ~r o~r o~r o~z ð5aÞ pðAÞ ¼ 2iDðzÞ a0 þ
4DU0 ðtÞjAj2 2iQ0 ðtÞjAj2 ; n0 ð5bÞ
where DU0(t) = kn2I0(t)z0 is the phase shift at focus for the nonlinear medium with one Rayleigh length, Q0(t) = b2I0(t)z0 is the nonlinear absorption at focus for the nonlinear medium with one Rayleigh length and a0 ¼ a0 z0 is the linear absorption within one Rayleigh length. For the case of p = 0, which corresponds to the beam propagation in vacuum, the function Að~r; ~zÞ has the following form: Að~r; ~zÞ ¼ expð~r2 Þ:
equ ¼ L
Z
L
z2 ¼0
d~zðz1 ; z2 Þ ¼
73
Z 0
L
dz2 1 þ ðz1 þ z2 Þ2
tan1 ðz1 Þ; ¼ tan1 ðz1 þ LÞ
ð7Þ
where z ¼ z=z0 ¼ z1 þ z2 (where z1 is the separation of the front face of the medium from the waist and z2 is the propagation distance from the front face), ¼ L=z0 and L equ is the dimensionless quantity. L It can be seen from Eq. (7) that the equivalent medium length in new coordinate system is the function of the front surface position of medium and the medium length. The equivalent medium length as a function of the position of medium is shown in Fig. 1. We can see from Fig. 1 that the equ is the function of equivalent medium length L the position of medium and medium length L. When the focal point of beam is located in middle equ obtain its maxof medium, for fixed length L; L imum. The equivalent length Lequ is approximately proportional to T1 1 for an open Z-scan [8]. equ increases slowly with medium length L and deL crease quickly as the distance between the middle of medium and focal point of beam increases, and if medium length L ! 1 and the focal point of beam is located in middle of medium, equ ¼ 2tan1 ðL=2Þ L is equal to p. So the length of medium in new coordinates is shrunken greatly. Eqs. (1) and (5) are parabolic partial differential
ð6Þ
It means that the form of the function Að~r; ~zÞ in vacuum keeps unchanged while the beam propagates in this coordinate system. Therefore, the use of the transformations in Eqs. (4) can significantly reduce the size of the data arrays and the computational time, since fewer points can express the radial field distribution and the propagation distance is shrunken greatly. The Z-scan experimental arrangement is the same as that in [2]. If the front surface of the medium is located at z and the medium length is L, one equ by can obtain the equivalent medium length L Eq. (4c), which is defined as the distance between the front surface and the exit surface of medium in this new coordinate system,
equ Fig. 1. The relation between the equivalent medium length L and the position of medium x. The real medium length is L = 1, 2, 4, respectively. The unit of medium length is Rayleigh length z0.
74
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
equations and can be solved numerically by Crank–Nicholson finite difference scheme [17]. Here we call the numerical simulation of Eq. (1) by Crank–Nicholson finite difference scheme as conventional method (CM), and the numerical simulation of Eq. (1) by Crank–Nicholson finite difference scheme as coordinate transformation method (CTM).
It is seen that when i = N, we can get S =2. However, for other value of i, there zN ¼ L are two values of zi . We can use Eqs. (4) to get the real electric q field Eðr; zi Þ ffiby Aði; ~rÞ; z ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S and r ¼ ~r 1 þ ðzi þ L S Þ2 . zi þ L
3. Results and discussion 2.2. Coordinate mapping Below we apply coordinate mapping method (CMM) to numerical simulations of Z-scan measurements. It can be seen that there is no explicit dependence on the medium position in the Eq. (5a). Therefore, the function Að~r; ~zÞ, solved from Eq. (5a), is not the explicit function of the medium position. The function Að~r; ~zÞ is only determined by the equivalent medium length that is the function of medium position and medium length. For a nonlinear medium with fixed length, different medium position in z-direction corresponds to different equivalent medium length. Using this rule, we can obtain the function Að~r; ~zÞ at the exit surface of the medium for different medium position by only computing beam propagation in the medium when the position of center point of the medium is located at origin point. The procedure is described below: S ¼ L=z0 , First, assuming the medium length is L when the position of front surface of medium is loS =2 the equivalent medium length is cated at L 1 Lequ ¼ 2tan ðLS =2Þ, which corresponds to the largest equivalent medium length. Second, slicing the medium into many thin slices, the beam propagates through each slice, and the function Aði; ~rÞ can be obtained, here Aði; ~rÞ denotes the function Að~r; ~zÞ after the beam propagation over ith slice, and is also Að~r; ~zÞ at the exit surface of medium when the medium is located at zi which can be determined by Eq. (7) as following: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 " #ffi9 u = 1< u L S 2S 4 1 2i 1 ; zi ¼ LS tL S =2Þ ; 2: tan N tan ðL ð8Þ
where N is the number of slices.
In order to check the range of validity as well as the accuracy of CM and CTM, we simulate Gaussian beam propagation in vacuum. The main parameters of the Gaussian beam are calculated numerically and then compared with the analytical values. The beam size which is evaluated as the second moment of intensity, can be written as: " R1 #1=2 2p 0 r2 Iðr; zÞr dr R1 wðzÞ ¼ : ð9Þ 2p 0 Iðr; zÞr dr The conditions of numerical computing are listed as following: the step size in r direction is Dr = w0/100, and the step sizes in z direction are equ ¼ 1= DL = z0/100 for CM and D~r ¼ 1=100; DL 100 for CTM, respectively. For Gaussian beams, the size of the numerical window has to be at least four times of the beam size w(z), in order that the reflection from the boundary can be neglected completely. Therefore, the size of the numerical window is selected to be four times of the beam size w(z) or w0 for CM or CTM. In our computations we used a PC with a 2.4-GHz Intel Pentium CPU and 512M DDR memory. Matlab was used for numerical computing. The calculated results of beam propagation over a distance of 20z0 are shown in Fig. 2. It is seen that CTM enhances the accuracy of computing at least one order of magnitude for the above computational conditions as CMs numerical computing. The comparison between computational time of two methods for different medium length is given in Fig. 3. Results show that CTM enhances the speed of numerical computing one or two orders of magnitude for the same conditions as CMs numerical computing. The reason of enhancement of accuracy and speed of numerical computing come from the size of the data arrays
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
75
Typical closed-aperture Z-scan curves without two-photon absorption and open-aperture Z-scan curves with two-photon absorption are given in Fig. 4, where the medium length LS = 10z0. The solutions from CTM and CMM agree very well. The computational time as a function of medium length for three methods is given in Fig. 5. The results show that the numerical efficiency of CTM enhances at least two orders of magnitude more than that of CM, and the numerical efficiency of
Fig. 2. The relative error ew of beam radius of the Gaussian beam propagation through the air for CM and CTM.
Fig. 4. Closed-aperture and open-aperture Z-scan curves. The medium length is 10z0. Solid line represents CTM, and circle represents CMM.
Fig. 3. Computational time taken by CM and CTM for different propagation distance. The propagation distance is scaled by Rayleigh length z0, and the temporal unit is second.
significantly reduced and fewer points can express the radial field distribution and the propagation distance shrunken greatly. Further, we compare the numerical simulations of Z-scan measurements obtained by CMM, CTM and CM, respectively. The parameters used in our calculations are listed as followings [2]: linear refractive index n0 = 1, Kerr coefficient n2r = 6.8 · 1014 cm2/GW, two-photon-absorption coefficient b2r = 5.8 cm/GW, beam waist w0r = 25 lm, and peak intensity I0 = 0.21 GW/cm2. All computations involve only a cw laser beam.
Fig. 5. Computational time taken by CM, CTM and CMM for numerical simulation of Z-scan measurements. The medium length is scaled by Rayleigh length z0, and the temporal unit is second. Circle, CM; square, CTM; triangle, CMM.
76
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
(a)
and each slice propagates as if it were a continuous wave (CW) laser beam. Second task is to accumulate all effects of time slices. The computational time is determined by the number of time slices, the total CPU time is obtained by multiplying the CPU time of cw laser beam and the number of time slices. For Gaussian-shape pulse and hyperbolic-sech shape pulse, the peak-to-valley normalized transmittance difference DTp–v of closed-aperture Z-scan curves without two-photon absorption as a function of medium length are given in Fig. 6(a) and (b). It can be seen that the ratio of DTp–v of Gaussian-shape pulse to that of pffiffiffi CW laser beam is 1= 2, and the ratio DTp–v of hyperbolic-shape pulse to that of CW laser beam is 2/3 in the case of ‘‘thin’’ medium [18]. The ratio decreases as medium length increases, finally reaches the value 0.676 for Gaussian-shape pulse and 0.64 for hyperbolic-shape pulse.
4. Conclusion
(b) Fig. 6. (a) The peak–valley normalized transmittance difference DTp–v vs. medium length L for CW laser beam, Gaussian-shape and hyperbolic-sech shape pulsed laser beam. (b) The peak– valley normalized transmittance difference DT p–v =DT CW p–v vs. medium length L for Gaussian-shape and hyperbolic-sech shape pulsed laser beam. DT CW p–v is the peak–valley normalized transmittance difference of CW laser beam.
CMM enhances at least one order of magnitude more than that of CTM. Therefore, compared with CM, the total numerical efficiency enhances over three orders magnitude. All the discussions above involve only the cw laser beam. For pulsed laser beam with the temporal envelope of Gaussian-shape or hyperbolic-sech shape, Eqs. (5) are hold correct if we ignore the effect of dispersion. The numerical simulation of laser pulse propagation in a thick nonlinear media can be split into two separate tasks. The first task is to split laser pulse into a number of time slices,
We apply CTM to improve the numerical efficiency of beam propagation and CMM to improve the numerical efficiency of numerical simulation of Z-scan measurements. The numerical analyses show that the numerical efficiency of the analyses on beam propagation enhances at least two orders of magnitude and total numerical efficiency of numerical simulation of Z-scan characteristics enhances at least three orders of magnitude, compared with that of CM. For CW laser beam and pulsed laser beam, we applied these new methods to analyze the peak-to-valley normalized transmittance difference for closed-aperture Z-scan as the function of medium length and nonlinear refraction of medium.
Acknowledgements This research is supported by Project 60025512 supported by the National Natural Science Foundation of China, Natural Science Foundation of Tianjin (Grant 0436012 11), a key project of the Ministry of Education (Grant 00026), the Foundation for University Key Teachers of the Ministry
W.-P. Zang et al. / Optics Communications 244 (2005) 71–77
of Education, and the Fok Ying Tung Education Foundation (Grant 71008). References [1] M. Sheik-Bahae, A.A. Said, E.W. Van Stryland, Opt. Lett. 14 (1989) 955. [2] M. Sheik-Bahae, A.A. Said, T.H. Wei et al., IEEE. J. Quantum Electron. 26 (1990) 760. [3] M. Sheik-Bahae, A.A. Said, T.H. Wei et al., Opt. Eng. 30 (1991) 1228. [4] J.A. Hermann, R.G. McDuff, J. Opt. Soc. Am. B 10 (1993) 2056. [5] W.-P. Zang, J.-G. Tian, C.-P. Zhang et al., Acta Phys. Sin. 43 (1994) 476 (in Chinese). [6] J.-G. Tian, W.-P. Zang, G.Y. Zhang et al., Acta Phys. Sin. 43 (1994) 1712 (in Chinese). [7] J.-G. Tian, W.-P. Zang, G.-Y. Zhang et al., Appl. Opt. 34 (1995) 4331.
77
[8] J.A. Hermann, J. Opt. Soc. Am. B 14 (1997) 814. [9] W.-P. Zang, J.-G. Tian, Z.-B. Liu et al., Opt. Lett. 28 (2003) 722. [10] W.-P. Zang, J.-G. Tian, Z.-B. Liu et al., J. Opt. Soc. Am. B 32 (2004) 63. [11] P.B. Chapple, J.A. Hermann, R.G. McDuff, Opt. Quantum Electron. 31 (1999) 555. [12] S. Hughes, J.M. Burzer, G. Spruce, B.S. Wherrett, J. Opt. Soc. Am. B 12 (1995) 1888. [13] D.I. Kvosh, S. Yang, D.J. Hagan et al., Appl. Opt. 38 (1999) 5168. [14] P.B. Ulrich, Naval Research Laboratory Report 7706, May 29, 1974. [15] E.A. Sziklas, A.E. Siegman, Appl. Opt. 14 (1975) 1874. [16] J.A. Fleck, J.R. Morris, M.D. Feit, Appl. Phys. 10 (1976) 129. [17] C.F. Gerald, P.O. Wheatley, Applied Numerical Analysis, Addison Wesley, Reading, MA, 1984 (Chapter 8). [18] R.L. Sutherland, Handbook of Nonlinear Optics, Marcel Dekker, New York, 1996 (Chapter 7).