Optics Communications 266 (2006) 505–511 www.elsevier.com/locate/optcom
A Modified full-vectorial finite-difference beam propagation method based on H-fields for optical waveguides with step-index profiles Jinbiao Xiao *, Xiaohan Sun Lab of Photonics and Optical Communications, Electronic Engineering Department, Southeast University, Nanjing 210096, China Received 20 February 2006; received in revised form 11 May 2006; accepted 11 May 2006
Abstract A modified full-vectorial finite-difference beam propagation method based on H-fields in solving the guided-modes for optical waveguides with step-index profiles is described. The propagation is split into two substeps. In the first substep, the field propagates in the absence of the cross-coupling terms, and then they are evaluated and double used in the second substep. An improved six-point finite-difference scheme is constructed to approximate the cross-coupling terms along the transverse directions. By using the imaginary-distance procedure, the field patterns and the normalized propagation constants of the fundamental modes for a buried rectangular waveguide and rib waveguide are presented, and the hybrid nature of the full-vectorial guided-modes is demonstrated. Solutions are in good agreement with the benchmark results from film mode matching method, which tests the validity and utility of the present method. 2006 Elsevier B.V. All rights reserved. PACS: 42.82.m; 42.79.Gn
1. Introduction The beam propagation method (BPM) has become one of the most popular and powerful numerical techniques for modeling electromagnetic wave propagation in photonic integrated circuits/planar lightwave circuits (PICs/ PLCs) and its imaginary-distance procedure has also been successfully used to solve the guided-modes for optical waveguides [1]. For complicated dielectric structures in which large refractive index contrasts and/or polarization effects are involved, full-vectorial (FV) BPM (FV-BPM) becomes necessary. So far, numerous versions of FVBPM have been proposed [1]. Among them, the FV-BPM based on finite-difference (FD) method (FV-FD-BPM) is attractive due to the simplicity of its implementation and the sparsity of its resulted matrix [2–8]. The discretization for FV-BPM formulation along the longitudinal direction should be dealt carefully because of
*
Corresponding author. Fax: +86 2583793686. E-mail address:
[email protected] (J. Xiao).
0030-4018/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2006.05.027
the cross-coupling terms (CCTs) which are indispensable in an FV-BPM. Improper treatment of the CCTs may cause instability and inaccuracy especially in strongly guiding structures. Huang et al. [2] made use of the iterative solver to perform this task, which has been applied to both the isotropic [2] and anisotropic [3] dielectric structures, but the iterative algorithm is time-consuming. Instead, the alternating direction implicit (ADI) method is a noniterative procedure for solving multidimensional partial differential equations, so the calculation is efficient and it does not suffer from nonconvergence [9]. The ADI algorithm with the Crank–Nicholson (CN) scheme was first proposed by Yamauchi et al. [10] to solve the scalar BPM formulation, and same approach was also used by Liu et al. [11] to solve the semi-vectorial BPM formulation. More recently, with operator decomposition technique the ADI algorithm with CN scheme has been applied to FV-BPM formulation [7,8]. Alternatively, Petra´cˇek et al. [6] proposed a modified FD scheme to discretize the FV-BPM formulation based on E-fields along the longitudinal direction. The propagation is also split into two substeps. In the first substep, only semi-vectorial formulation is used, whilst in the second
506
J. Xiao, X. Sun / Optics Communications 266 (2006) 505–511
substep, full-vectorial formulation is adopted in which the CCTs are double used. The order of two substeps is reversed for each polarized field component so that the CCTs are always expressed in implicit form. Therefore, the calculation is stable. Moreover, the resulted matrix is tridiagonal, thus the calculation is efficient. However, numerical results are not presented in [6]. On the other hand, the accuracy of an FV-BPM strongly depends on the FD scheme for the CCTs along the transverse directions. A direct discretization for the CCTs has been proved to result in poor convergence behavior [2,3]. To the best of our knowledge, there have been only two attempts so far to develop improved FD scheme for the CCTs based on E-fields [5,8]. In this paper, the modified FD scheme proposed by Petra´cˇek et al. in [6] is applied to the FV-BPM formulation based on H-fields for optical waveguides with step-index profiles. Although all field components can be calculated from either transverse electric (E-formulation) or magnetic formulation (H-formulation), the resulted field patterns do not always look exactly the same [12,13]. As a general rule, if H-fields are desired, it is best to choose the H-formulation. Moreover, an improved six-point FD scheme is also developed to approximate the CCTs along the transverse directions, which is independent of specific structures of waveguide, so the procedure is simple and versatile. In order to test the validity and utility of the present method, the fundamental modes for a buried rectangular waveguide and rib waveguide are analyzed by using the imaginary-distance procedure on its original implementation.
2 2k 0 n0 ohx o hx oz oz2 2 2k 0 n0 ohy o hy 2 oz oz Eq. (2) reduces to " x P xx þ P yxx o hx ¼ P yx oz hy
ð4aÞ ð4bÞ
P xy P xyy þ P yyy
#
hx
hy
ð5Þ
with P xxx hx P yxx hx P xyy hy P yyy hy P xy hy P yx hx
2 1 o hx 1 2 2 2 ¼ þ k n n0 hx 2jk 0 n0 ox2 2 0 1 1 ohx 1 2 2 2 o 2 ¼ n þ k 0 n n0 hx 2jk 0 n0 oy n2 oy 2 1 1 ohy 1 2 2 2 o 2 ¼ n þ k 0 n n0 hy 2jk 0 n0 ox n2 ox 2 2 1 o hy 1 2 2 2 ¼ þ k n n0 hy 2jk 0 n0 oy 2 2 0 2 1 o hy o 1 ohy n2 ¼ 2jk 0 n0 oyox oy n2 ox 2 1 o hx o 1 ohx n2 ¼ 2jk 0 n0 oxoy ox n2 oy
ð6aÞ ð6bÞ ð6cÞ ð6dÞ ð6eÞ ð6fÞ
where P xxx and P xyy are x-dependent operators, P yxx and P yxx are y-dependent operators, and Pxy and Pyx correspond to the CCTs. (5) together with (6) is the basic formulation for FV-BPM based on H-fields. 2.2. Discretization along the longitudinal direction
2. Mathematical formulation 2.1. Basic formulation for FV-BPM based on H-fields The vectorial wave equation based on H-fields is directly derived from Maxwell’s equations and it is given as below * * 1 r2 H þ 2 ðrn2 Þ ðr H Þ þ k 20 n2 H ¼ 0 ð1Þ n pffiffiffiffiffiffiffiffiffi where k 0 ¼ x e0 l0 is the wavenumber in vacuum and n = n(x, y, z) is the refractive index of the medium. If the refractive index varies slowly along the direction of the wave propagation (say, z-direction), the transverse components of (1) can be approximated by 1 on2 oH x oH y r2 H x þ k 20 n2 H x 2 ¼0 ð2aÞ n oy oy ox 1 on2 oH y oH x ¼0 ð2bÞ r2 H y þ k 20 n2 H y 2 n ox ox oy *
By letting H x ¼ hx expðjk 0 n0 zÞ
ð3aÞ
H y ¼ hy expðjk 0 n0 zÞ
ð3bÞ
where n0 is the reference index, and assuming the slowlyvarying envelope approximation, i.e., [2]
The modified FD scheme in [6] is used to discretize (5) along the longitudinal direction, in which the propagation step is also split into two substeps. The algorithm is summarized as follows: (1) First substep: z ¼ lDz ! l þ 12 Dz Dz Dz y 1 2 ¼ P 1 P xxx hp;q;lþ 1 þ hp;q;l x 2 2 xx x Dz Dz y 1 p;q;lþ12 2 DzP h P ¼ 1 þ 1 P xyy hp;q;lþ hp;q;l yx x y 2 2 yy y (2) Second substep: z ¼ 1 þ 12 Dz ! ðl þ 1ÞDz Dz Dz x p;q;lþ12 P 1 P yyy hp;q;lþ1 ¼ 1 þ y yy hy 2 2 Dz y 1 P xx hxp;q;lþ1 DzP xy hyp;q;lþ1 2 Dz x 1 2 ¼ 1 þ P xx hp;q;lþ x 2 1
ð7aÞ ð7bÞ
ð7cÞ
ð7dÞ
2 , and hp;q;lþ1 denote the component of h at where hxp;q;l , hp;q;lþ x x x plane z = lDz, z ¼ 1 þ 12 Dz, and z = (l + 1)Dz with x = pDx and y = qDy, respectively, in which p, q, and l are integers, Dx · Dy is the mesh size in transverse directions,
J. Xiao, X. Sun / Optics Communications 266 (2006) 505–511
and Dz is the propagation step size. The same denotations are also used for hy. In the present algorithm, the propagation step is split into two substeps z ¼ z = lDz ! (l + 1)Dz lDz ! l þ 12 Dz and z ¼ l þ 12 Dz ! ðl þ 1ÞDz. In the first substep, the field propagates in the absence of the CCTs, e.g. (7a) for hx, and then they are evaluated and double used in the second substep, e.g. (7d). In fact, the first and second substep corresponds to the semi- and full-vectorial formulation, respectively. The order of the two substeps is reversed for each polarized field component, e.g. (7b) and (7c) for hy, so that the CCTs are always expressed in implicit form. As a result, the calculation is stable. Moreover, the operator inversion for the CCTs is avoided through the operator splitting technique. In addition, the resulted matrix is tridiagonal, thus the calculation is efficient.
If the CCTs in (5) are neglected, semi-vectorial BPM formulation is obtained, so the discretizations for ðP xxx þ P yxx Þhx and ðP xyy þ P yyy Þhy are the same as those used in [2]. As mentioned above, the accuracy is severely restricted by the FD scheme for the CCTs which are indispensable in an FV-BPM. Here, we construct an improved FD scheme, which is derived and presented in Appendix A as (A5) and (A6), to approximate the CCTs along the transverse directions. In comparison with the direct discretization for the CCTs shown in [2], two more adjacent points corresponding to magnetic fields are added on the RHS in (A5) and in (A6). Therefore, the present formula is a six-point FD scheme and is independent of specific structures of waveguide since they are expressed in explicit forms. So, the procedure is simple and versatile. For boundaries of the computational window, the transparent boundary condition (TBC) [14] is used in this algorithm in order to eliminate the non-physical reflection. 3. Numerical results We first calculate the normalized propagation constants and the field patterns of the fundamental mode for a typical buried rectangular waveguide [15] whose geometry is shown in Fig. 1, in which w = h = 1.0 lm, n1 = 1.0, n2 = 1.5 with k = 1.5 lm. The fundamental mode can be calculated by the imaginary-distance procedure in which the coordinate z in the propagation direction is changed to js [16,17]. The effective index neff is evaluated by the growth in the field amplitude neff ðsÞ ¼ n0 þ
ln
R
R /ðx; y; s þ DsÞdxdy ln /ðx;y; sÞdxdy k 0 Ds ð8Þ
and the normalized propagation constant B is here defined as
Y0 n1 h
X0
y
n2
x Core w Claddings Fig. 1. Geometry of a buried rectangular waveguide.
B¼
n2eff n21 n22 n21
ð9Þ
As the injected field, we chose a step function as 1; for jxj 6 0:5 lm and jy j 6 0:5 lm /ðx; y; 0Þ ¼ 0; otherwise ð10Þ for the major component and zero for the minor component. The transverse mesh size is Dx = Dy = 0.025 lm, the propagation step size is Dz = 0.1 lm, and the reference index is n0 = (n1 + n2)/2. The computational window parameters are X0 = 1.5 lm and Y0 = 1.5 lm. Fig. 2 shows the normalized propagation constants B for the fundamental mode as a function of propagation distance s. It is noted that the convergence values of the normalized propagation constants for both Quasi-TE and Quasi-TM mode are obtained after propagation only 10 lm. The values for Quasi-TE mode are identical to those for Quasi-TM at each propagation step. Moreover, the solutions keep convergence after we propagate the injected fields till 50 lm. The resulted equation systems are solved by using MATLAB subroutines and implemented on a Pentium IV CPU PC with a CPU clock of 2.0 GHz. It is founded that the required CPU time for
0.5000
Quasi-TE Mode Quasi-TM Mode 0.4975
B
2.3. Discretization along the transverse directions
507
Δx=0.025 μm Δy=0.025 μm Δτ=0.1 μm
0.4950
0.4925 0
10
20
30
40
propagation distances τ (μm)
50
Fig. 2. Normalized propagation constants B of the fundamental mode for a buried rectangular waveguide as a function of propagation distance s.
508
J. Xiao, X. Sun / Optics Communications 266 (2006) 505–511
one propagation step is less than 1 s (the number of transverse mesh points is 80 · 80). The computational time for the same problem by using iterative solver in [2] is over 10 s. Therefore, the present method is efficient. Fig. 3 illustrates the typical field patterns for the fundamental mode in Quasi-TE and Quasi-TM modes. Because of the symmetric structure of the waveguide, only onequarter of the field patterns are displayed. It is seen that a single peak appears in the guiding region for the major component, and four peaks locate at the corners for the minor component. Moreover, the discontinuous feature of the derivatives of magnetic field component hx and hy at horizontal and vertical interfaces, respectively, are clearly illustrated. Therefore, the hybrid nature of the full-vectorial guided-modes is demonstrated. The second example is a rib waveguide [2] as shown in Fig. 4. The refractive indices of the superstrate, core, and substrate are n1 = 1.0, n2 = 3.44, and n3 = 3.40. The height and the width of the rib are h = 1.0 lm and w = 3.0 lm. The wavelength is k = 1.15 lm. The slab waveguide height t varies from 0.1 to 0.9 lm. The computational window parameters are Ya = 0.5 lm above the top of the rib, Yb = 2.5 lm below the guiding layer, and Xa = 2.5 lm (or 5.0 lm for t = 0.9 lm) aside from the rib lateral side. We also calculated the fundamental mode of this rib waveguide by using imaginary-distance procedure. The convergence of the numerical solutions was evaluated by
w
Ya
Xa
h n1
Superstrate
y t
n2 x
n3
Core Substrate
Yb
Fig. 4. Geometry of a rib waveguide.
reducing the transverse mesh size. Fig. 5 shows the percentage error of the normalized propagation constants as a function of the transverse mesh size (Dx = Dy) in which the slab waveguide height is t = 0.5 lm and the propagation step size is fixed at Ds = 0.1 lm. The percentage error is defined here as Bcalculated BFMMM errorð%Þ ¼ ð11Þ B FMMM
where BFMMM is the solutions from film mode matching method (FMMM) [18] which is so far regarded as the benchmark results for comparing the accuracy. It can be
1.0
y(μm)
y(μm)
1.0
0.5
(a)
0.5
0.5
1.0
x(μm)
0.5
(b)
1.0
1.0
x(μm)
y(μm)
y(μm)
1.0
0.5
(c)
0.5
0.5
x(μm)
1.0
(d)
0.5
x(μm)
1.0
Fig. 3. Field patterns of the fundamental mode for a buried rectangular waveguide: Quasi-TE mode (a) hx and (b) hy; Quasi-TM mode (c) hx and (d) hy.
J. Xiao, X. Sun / Optics Communications 266 (2006) 505–511
509
0.16
present method FV-FD-BPM in [2]
25
15
0.12
Δτ =0.1μm
error,%
error, %
20
Quasi-TM Mode
10
0.08
0.04
Δx=0.025μm
5
0.00 Quasi-TE Mode
0.0
0 0.05 0.10 0.15 transverse mesh size Δx=Δy (μm)
0.5
0.20
Fig. 5. Percentage error of the normalized propagation constants B of the fundamental mode for a rib waveguide as a function of transverse mesh size Dx = Dy.
observed that, for both Quasi-TE and Quasi-TM modes, the numerical results converge to the solutions from FMMM as the transverse mesh size decreases. Moreover, the percentage error of our results is smaller than that of the solutions obtained by the FV-FD-BPM in [2]. Fig. 6 plots the percentage error of the calculation as a function of the propagation step size Ds in which the transverse mesh size is Dx = Dy = 0.025 lm. It is seen that the variation of the propagation step size does not have much influence on the accuracy of the normalized propagation constants, so relatively large propagation step size Ds may be used. Tables 1 and 2, respectively, present the normalized propagation constants in Quasi-TE and Quasi-TM mode as a function of slab waveguide height after the propagation reaches the convergence. The transverse mesh size is Dx = Dy = 0.025 lm and the propagation step size is Ds = 0.1 lm. It is seen that the normalized propagation constants for both Quasi-TE and Quasi-TM modes increase
Quasi-TE Mode Quasi-TM Mode
Δy=0.025μm 1.0
1.5
Δτ (μm)
2.0
2.5
3.0
Fig. 6. Percentage error of the normalized propagation constants B of the fundamental mode for a rib waveguide as a function of propagation step size Ds.
with the increment of slab waveguide height. In order to examine the accuracy of present method, the results of the analysis under study together with those of FV-FDBPM [2], FMMM [18], modal matching method (MMM) [19], finite element method (FEM) [19], and mapped Galerkin method (MGM) [20] are also summarized in the stables. Our results are closer to the values from FMMM than those from FV-FD-BPM in [2], especially; the results in Quasi-TM obtained by the present method and the FMMM agree up to the third decimal digit. So the present FD scheme for CCTs improves the computational accuracy. Fig. 7 illustrates the field patterns of the fundamental mode in Quasi-TE and Quasi-TM modes with slab waveguide height t = 0.5 lm. Half of the patterns are displayed because of the symmetric structure with respect to x = 0. Here again the required discontinuities for hx and hy are clearly evident.
Table 1 Normalized propagation constants of Quasi-TE mode for a rib waveguide as a function of the slab waveguide height computed by various methods t (lm)
MMM [19]
FEM [19]
MGM [20]
FV-FD-BPM [2]
FMMM [18]
Our results
0.1 0.3 0.5 0.7 0.9
0.3019 0.3111 0.3270 0.3512 –
0.3019 0.3110 0.3270 0.3512 –
0.2959 0.3059 0.3231 0.3487 0.3871
0.3039 0.3144 0.3303 0.3533 0.3879
0.3019 0.3110 0.3270 0.3510 0.3883
0.3028 0.3118 0.3277 0.3515 0.3886
Table 2 Same as Table 1 except that the results are obtained for Quasi-TM mode t (lm)
MMM [19]
FEM [19]
MGM [20]
FV-FD-BPM [2]
FMMM [18]
Our results
0.1 0.3 0.5 0.7 0.9
0.2675 0.2751 0.2890 0.3107 –
0.2675 0.2751 0.2890 0.3106 –
0.2543 0.2626 0.2773 0.2999 0.3361
0.2690 0.2776 0.2915 0.3127 0.3451
0.2674 0.2751 0.2890 0.3107 0.3455
0.2678 0.2754 0.2893 0.3109 0.3456
J. Xiao, X. Sun / Optics Communications 266 (2006) 505–511
1
1
0
0
y(μm)
y(μm)
510
-1
-1
-2
-2 1
2
3
x(μm)
(a)
1
2 x(μm)
3
1
2 x(μm)
3
(b)
0
0
y(μm)
1
y(μm)
1
-1
-1
-2
-2
1 (c)
x(μm)
2
3 (d)
Fig. 7. Field patterns of the fundamental mode for a rib waveguide: Quasi-TE mode (a) hx and (b) hy; Quasi-TM mode (c) hx and (d) hy.
4. Conclusions The modified FD scheme is applied to discretize the FVBPM formulation based on H-fields along the longitudinal direction and an improved six-point FD scheme is also developed to approximate the cross-coupling terms along the transverse directions. The field patterns and the normalized propagation constants of the fundamental mode for a buried rectangular and a rib waveguide are presented by using the imaginary-distance procedure which is programmed by MATLAB subroutines. Solutions are in good agreement with those from film mode matching method, which shows that the present method improve the computational accuracy, although the traditional FD scheme is used to approximate the semi-vectorial operators along the transverse directions. Moreover, the results do not show up any instabilities and the calculation is efficient. Appendix A The improved FD scheme for the cross-coupling terms along the transverse directions is derived as follows. The first term on the RHS in (6e) is approximated by p;q;l
1 o2 hy 1 hpþ1;qþ1;l hyp1;qþ1;l 2jk 0 n0 oyox 8jk 0 n0 DxDy y hypþ1;q1;l þ hyp1;q1;l ðA1Þ
whilst the second term can be approximated by p;q;l 1 1 ohy 2 o n 2jk 0 n0 oy n2 ox " p;qþ1;l p;q1;l # n2p;q;l 1 ohy 2 1 ohy 2 2 n2 ox n ox 2jk 0 n0 Dy
ðA2Þ
with p;qþ1;l 1 ohy 2 2 1 pþ1;qþ12;l 1 hy 2 hyp1;qþ2;l 2 2 n ox np;q;l þ np;qþ1;l 2Dx 2 1 pþ1;qþ1;l h 2 þ hypþ1;q;l 2 np;q;l þ np;qþ1;l 4Dx y hyp1;qþ1;l hyp1;q;l ðA3Þ 1 p;q ;l 1 ohy 2 2 1 pþ1;q12;l p1;q12;l h h y n2 ox n2p;q;l þ n2p;q1;l 2Dx y
2 1 hpþ1;q1;l þ hypþ1;q;l 2 2 np;q;l þ np;q1;l 4Dx y hyp1;q1;l hyp1;q;l ðA4Þ where np,q,l stands for the index distribution at plane z = lDz with x = pDx and y = qDy. So the FD scheme for Pxyhyjp,q,l along the transverse directions in (6e) is expressed by
J. Xiao, X. Sun / Optics Communications 266 (2006) 505–511
p;q;l P xy hy
1 8jk 0 n0 DxDy
"
! 2n2p;q;l 1 2 hpþ1;qþ1;l np;q;l þ n2p;qþ1;l y !
2n2p;q;l hp1;qþ1;l n2p;q;l þ n2p;qþ1;l y ! 2n2p;q;l 1 2 hpþ1;q1;l np;q;l þ n2p;q1;l y ! 2n2p;q;l 1 2 hp1;q1;l np;q;l þ n2p;q1;l y
! 2n2p;q;l 2n2p;q;l hp;qþ1;l n2p;q;l þ n2pþ1;q;l n2p;q;l þ n2p1;q;l x ! # 2n2p;q;l 2n2p;q;l p;q1;l þ 2 h np;q;l þ n2pþ1;q;l n2p;q;l þ n2p1;q;l x
1 þ
! 2n2p;q;l 2n2p;q;l hpþ1;q;l n2p;q;l þ n2p;qþ1;l n2p;q;l þ n2p;q1;l y ! # 2n2p;q;l 2n2p;q;l p1;q;l þ 2 h np;q;l þ n2p;qþ1;l n2p;q;l þ n2p;q1;l y ðA5Þ Similarly, the FD scheme for Pyxhxjp,q,l in (6f) is expressed by " ! p;q;l 2n2p;q;l 1 1 2 P yx hx hpþ1;qþ1;l 8jk 0 n0 DxDy np;q;l þ n2pþ1;q;l x ! 2n2p;q;l hp1;qþ1;l 1 2 np;q;l þ n2p1;q;l x ! 2n2p;q;l 1 2 1 hxpþ1;q1;l np;q;l þ n2pþ1;q;l ! 2n2p;q;l hp1;q1;l þ 1 2 np;q;l þ n2p1;q;l x
511
ðA6Þ References [1] R. Scarmozzino et al., IEEE J. Sel. Top. Quantum Electron. 6 (2000) 150. [2] W.P. Huang et al., IEEE J. Quantum Electron. 29 (1993) 2639. [3] C.L. Xu et al., J. Lightwave Technol. 12 (1994) 1926. [4] I. Mansor et al., J. Lightwave Technol. 14 (1996) 908. [5] J. Yamauchi et al., J. Lightwave Technol. 16 (1998) 2458. [6] J. Petra´cˇek et al., Proc. SPIE 3320 (1998) 294. [7] Y.L. Hseuh et al., J. Lightwave Technol. 17 (1999) 2389. [8] Y.Z. He et al., IEEE Photon. Technol. Lett. 15 (2003) 1381. [9] W.H. Press et al., Numerical Recipes in C – The Art of Scientific Computing, second ed., Cambridge University Press, New York, 1992. [10] J. Yamauchi et al., Electron. Lett. 27 (1991) 1663. [11] P.L. Liu et al., IEEE J. Quantum Electron. 29 (1993) 1205. [12] D. Marcuse, IEEE J. Quantum Electron. 28 (1992) 459. [13] C. Vassallo, Optical Waveguide Concepts, Elsevier, Amsterdam, 1991. [14] G.R. Hadley, IEEE J. Quantum Electron. 28 (1992) 363. [15] S. Subdbø, J. Lightwave Technol. 10 (1992) 418. [16] S. Ju¨ngling et al., IEEE J. Quantum Electron. 30 (1994) 2098. [17] Y.Z. He et al., Opt. Commun. 225 (2003) 151. [18] S. Subdbø, Pure Appl. Opt. 3 (1994) 381. [19] S. Selleri et al., Opt. Quantum Electron. 33 (2001) 373. [20] J. Xiao et al., J. Opt. Soc. Am. B 21 (2004) 798.