Baets, Roel

Bienstman, Peter

Maes, Björn

Vandersteegen, Peter

English

Ghent University Academic Bibliography

DOI 10.1007/s11082-006-0021-xOptical and Quantum Electronics (2006) 38:35–44 © Springer 2006Using the complex Jacobi method to simulate Kerrnon-linear photonic componentsp e t e r vander s t e egen∗, b j o rn mae s, p e t e r b i en s tmanand roel baet sDepartment of Information Technology, Ghent University-IMEC, St.-Pietersnieuwstraat 41, 9000 Ghent,Belgium(*author for correspondence: E-mail: Peter.Vandersteegen@intec.UGent.be)Received 30 September 2005; accepted 31 January 2006Abstract. The iterative complex Jacobi technique has been extended to simulate the third order Kerreffect in wavelength scale dielectric structures. This method solves the Helmholtz equation in a dis-crete ﬁnite simulation space by an iterative process. An update equation reﬁnes the ﬁeld values at eachiteration step, until a desired accuracy is achieved. We have extended the iterative process with anextra calculation step which allows simulating materials with the non-linear third order Kerr effect. Ouradjustment of the discrete ﬁeld operators in the update equation also introduces PMLs as absorbingboundaries and the total ﬁeld/scattered ﬁeld formalism as ﬁeld source for this method.Key words: central differences, complex Jacobi iteration, Helmholtz equation, Kerr effect, PML, totalﬁeld/scattered ﬁeld1. IntroductionIn this paper several improvements to the recently introduced complexJacobi (CJ), (Hadley 2005) technique are discussed. The proposed exten-sions allow modeling of dielectric structures with the instantaneous non-linear Kerr effect. This Kerr effect is modelled as an intensity dependentrefractive index. High intensity is achievable by increased conﬁnement oflight in advanced dielectric structures, e.g. photonic crystals and photonicwires. The fabrication of these structures has been made possible by recentadvanced production techniques.Finite Difference Time Domain (FDTD), (Taﬂove 1995), Beam Propaga-tion (BPM), (Burzler et al. 1996) and eigenmode expansion (EME), (Maeset al. 2005) are able to simulate a wide variety of non-linear dielectric wave-length scale components. FDTD gives insight in the time-evolution of theﬁelds by numerically integrating Maxwell’s equations in the time domain.A steady state solution however requires all transient ﬁelds to be extinct.The extinction rate of these transient ﬁelds is completely dependent on thereﬂections in the structure. BPM only gives an accurate solution for thesteady state ﬁelds if reﬂections are negligible. The lack of bidirectionality is36 P. VANDERSTEEGEN ET AL.a limiting assumption. EME is a bidirectional frequency domain method.Its main application area are dielectric structures which consist of a few z-invariant sections.We propose several extensions of the complex Jacobi method, a generalfrequency domain method which numerically integrates the central differ-ence equivalent of the Helmholtz equation. Bidirectionality is automaticallyimplied by the Helmholtz equation. There is also no assumption ofz-invariance. The found ﬁeld-amplitudes are located on an equidistantmesh. These amplitudes are invariant in time, save a phase factor.Our improvements allow modeling of general non-linear components.Good absorption at the boundaries is achieved by extending the methodwith perfectly matched layers (PML). We inject a ﬁeld proﬁle in the simu-lation area by an adjusted total ﬁeld/scattered ﬁeld formalism (TFSF). ThePML and TFSF extensions are respectively based on (Chew and Weedon1994) and (Taﬂove 1995). Although the basic complex Jacobi method is afully vectorial solver in 3 dimensions, the proposed extensions have onlybeen tested for the scalar 2D Helmholtz equation.In Section 2 we will give an overview of the linear extensions for thecomplex Jacobi method, such as PML and TFSF. In Section 3 we willintroduce the non-linear extension. A quantative comparison of the pro-posed simulation technique has been made with non-linear EME in Sec-tion 4. Another example, a propagating soliton in homogeneous space, isdiscussed in Section 5.2. Linear complex Jacobi MethodThe basic complex Jacobi method (Hadley 2005) is an iterative methodwhich solves the Helmholtz equation on a mesh. Only the ﬁelds at equi-distant points are calculated. The derivatives in the Helmholtz equation arereplaced by central differences:xei,j ≡ ei+1,j − ei−1,j2x (1)The ﬁeld ei,j is located at the point (ix, jy), i and j are integers, xand y are the discretization steps in x and y directions. This gives us thefollowing scalar Helmholtz equation for a discrete mesh:∇2ei,j +k2ei,j =(2x2+ 2y2+ (k0ni,j )2)ei,j=0 (2)USING THE COMPLEX JACOBI METHOD 37ni,j =n(ix, jy)k0 = 2πλ0(3)Thewavelength in free space isλ0,ni,j is the position dependent refractive index.The two step iteration process reﬁnes the ﬁelds with a correction factoruntil a desired accuracy has been achieved.en+1i,j = eni,j +(2α1(1x2+ 1y2))−1((δ2x + δ2y +(2πλni,j)2)ei,j)(4)The superscript n is the iteration step. The iteration process consists ofrepeatedly updating the ﬁelds with this equation, ﬁrstly with α1, then withα2, until a desired accuracy is achieved.Using exp(−jωt) as convention for the time dependent phase factor ofthe invariant ﬁeld amplitudes results in the following optimal constants:α1 =√3− 1j and α2 =−α∗1 . Derivation of these optimal values for α1 andα2 can be found in the original paper (Hadley 2005).The implementation of the TFSF and PMLs results in the simulationarea of Fig. 1. The photonic structure under investigation is placed in thetotal ﬁeld region. These extensions are respectively described in subsection2.1 and 2.2.2.1. total ﬁeld/scattered ﬁeldThe TFSF formalism allows to inject a desired wave in the simulationspace. Figure 1 demonstrates how the scattered ﬁeld/total ﬁeld interfaceFig. 1. A simulation box with ‘total ﬁeld/scattered ﬁeld’ as ﬁeld source and PMLs as absorbing bound-ary conditions.38 P. VANDERSTEEGEN ET AL.Fig. 2. Interface between scattered ﬁeld–total ﬁeld.injects the desired ﬁeld proﬁle in the total ﬁeld region. In the scattered ﬁeldregion only the scattered ﬁeld, i.e. the total ﬁeld without the exciting ﬁeldis calculated.Equation 4 allows us to calculate en(i,j),tot in the total ﬁeld region,en(i,j),scatt in the scattered ﬁeld region. The use of the update equation in thescattered ﬁeld region implicitly assumes superposition of exciting ﬁeld andscattered ﬁeld. This also implies the exciting ﬁeld is a solution for the dis-cretized Helmholtz equation and the scattered ﬁeld region consists of linearmaterials. Non-linear materials have to be located in the total ﬁeld region.The calculation of an amplitude at the interface between the regions isillustrated by Fig. 2.en(i ′−1,j),scatt = en(i ′−1,j),tot − e(i ′−1,j),sourceen(i ′,j),tot = en(i ′,j),scatt + e(i ′,j),source, ∀n,∀j (5)The injected ﬁeld amplitudes e(i ′−1,j),source and e(i ′,j),source are the ampli-tudes of the exciting wave at respective positions (i ′ − 1, j ) and (i ′, j ).Substituting Equation (5) in the update Equation (4) in order to calculateen(i ′−1,j),scatt and en(i ′,j),tot gives the required update equations. The followingequations have omitted the j -index of the y-direction for simplicity.en+1i ′−1,scatt = eni ′−1,scatt+C(eni ′,tot − esource,i ′ + eni ′−2,scatt −2eni ′−1,scattx2−k2eni ′−1,scatt)en+1i ′,tot = eni ′,tot +C(eni ′+1,tot + eni ′−1,scatt + esource,i ′−1 −2eni ′,totx2−k2eni ′,tot)C =(2α(1x2+ 1y2))−1With these adjusted update formulae, it becomes possible to inject anydesired ﬁeld.2.2. absorbing boundaries based on PMLElimination of reﬂections at the boundaries of a ﬁnite simulation box isnecessary to achieve physical simulation results. Every wave incident onUSING THE COMPLEX JACOBI METHOD 39these boundaries should be absorbed for any polarization and incidentangle, without any spurious reﬂections. One of the most efﬁcient absorbingboundaries is the Perfectly Matched Layer (PML). In this paper we reporton the extension of the complex Jacobi method with a PML based on com-plex coordinate stretching (Chew and Weedon 1994). Complex Coordinatestretching results in the easiest implementation of PML for the complexJacobi method.The implementation of PMLs is based on adjusted difference operatorsin the absorbing boundary layer:(δ2x + δ2x)e=1sxx(1sxxe)+ 1syy(1syye)(6)For sx = 1.0 the previous formula reduces to the classical Helmholtzequation. A complex function sx results in absorption in the x-direction, acomplex function sy results in absorption in the y-direction.A ﬁnite difference equivalent of Equation (6) gives us:δ2xe=1sx(x1sx)xe+ 1s2x2x2e= 1sx,i( 1sx,i+1− 1sx,i−12x)(ei+1 − ei−12x)+ 1s2x,iei+1 + ei−1 −2eix2In our experience, optimal absorption was achieved with linearly increas-ing PML from sx = 1.0 to sx = 1.0+ 0.75j over 30 grid-points. These com-plex coordinates implicitly assume the time dependant phase factor of theamplitudes is exp(−jωt).3. Extending the iterative process to simulate Kerr non-linearitiesSimulating structures with the instantaneous Kerr effect can be modeledwith an extra extension besides TFSF and PML. The basic idea is an extraupdate-step of the refractive index after each iteration step. The instanta-neous Kerr effect is modeled by n=nlin +n2|E|2. This formula for the spa-tially dependent refractive index can be transformed in the following updateequation:ni,j =(ni,j,lin +ni,j,2|ei,j |2)2(7)This equation is used in combination with the iteration process, describedin Equation (4). After each iteration step, the ﬁeld at a certain locationresults in an adjusted refractive index.40 P. VANDERSTEEGEN ET AL.Fig. 3. Two Bragg mirrors encapsulate a cavity. (λresonance =1.56µm).Validation of this principle has been done by comparing with nonlineareigenmode expansion, as described in (Maes et al. 2005). The ﬁrst simula-tion test setup consists of a resonant cavity for both the linear and non-linear case. The simulation time needed by both methods is also discussed.The end of this section also addresses the inﬂuence of the discrete mesh.The second simulation test setup shows the calculation of a soliton in freespace.4. Comparison with non-linear eigenmode expansionComparison with eigenmode expansion (EME) has been done by simu-lating the 1-dimensional structure of Fig. 3. The refractive indices in thissetup are na =2.6 and nb =2.36. The distances da, db and dcavity are respecti-vily λresonance/4na, λresonance/4nb and λresonance/2nb. This structure is basicallya cavity encapsulated by two Bragg mirrors.The cavity causes a resonance peak in the middle of the bandgap, asillustrated in Fig. 4. A non-linear material in the central cavity causes ashift of the central resonance wavelength to a higher wavelength. Figure 5illustrates this behaviour. The incoming plane wave has an amplitude e =1Vm. The non-linear cavity uses n2 = 5.× 10−3 m2V 2 , n2 as deﬁned in formula7. Figure 5 also illustrates a good agreement between the results calculatedby EME and the results calculated by the complex Jacobi method (CJ).Another point under discussion is the simulation time. EME has anadvantage if large sections of the simulation space can be modeled aslinear. When using EME, these linear sections only need to be calculatedonce. EME uses an iterative process for the non-linear sections, whichresults in the re-calculation of the eigenmodes in these sections for everyUSING THE COMPLEX JACOBI METHOD 41Fig. 4. Power transmission and reﬂection, see structure in Fig. 3.Fig. 5. Power transmission for linear and non-linear cavity. (x =15 nm).iteration. This quickly becomes a bottleneck when increasing the amountof non-linear sections.Convergence of our method is guaranteed even for an extremely largenon-linearity of n2 = 0.025m2V 2 with the previously described simulationparameters. This results in an increase of refractive index of approximatelyn≈ 0.6. We did not notice a signiﬁcant decrease of convergence rate forthese large non-linearities in this setup.Therefore, our proposed extension to CJ is very well suited for struc-tures where the non-linearity is present in a large portion of the simulationdomain. In such a situation EME would not be very efﬁcient.Simulating the steady state solution of the previous problem with a timedomain simulation method – e.g. FDTD – would require extremely long42 P. VANDERSTEEGEN ET AL.Fig. 6. Resonance peak -cfr. Fig. 3- shifts in function of the discretization step. (The analytical disper-sion model uses Equation (8) to adjust na and nb in the calculated structure.).simulation times. This is especially the case for a structure with many reﬂec-tions, as is the case for this example. The bidirectionality of the componenteliminates BPM as a simulation tool.4.1. inﬂuence of discretization stepA coarse grid for simulating the linear structure of Fig. 3 results in a reso-nance peak shift to higher wavelengths, as can be seen in Fig. 6. The ‘ana-lytical dispersion model’ curve calculates the shift of the resonance peak byonly taking into account numerical dispersion. The effective refractive index−neff = ωck0 -, ‘felt’ by a plane wave in a discrete simulation space is differentthan the actual refractive index. Based on Fig. 6 we can conclude that theshift is mostly caused by numerical dispersion.The analytical dispersion model ignores the staircase phenomena whichis caused by only taking into account the refractive indices at an integer(x,y). This is also an intrinsic property of a discrete mesh.The analytical dispersion model curve calculates the resonance peak ofFig. 3 with the effective refractive indices neff,a and neff,b. These effectiverefractive indices are calculated with the dispersion relationship of planewave in a homogeneous wave.4 sin2(k0neff x2)x2= ω2c2neff =arcsin(k0nx2)k0x2(8)USING THE COMPLEX JACOBI METHOD 43Fig. 7. Injection of a gaussian ﬁeld proﬁle in a linear and non-linear medium. ((a): diffraction in linearmedium, (b): soliton in non-linear medium).For a small value formulae 8, neff reduces to the refractive indices incontinuous space: neff =n.5. Spatial soliton in a non-linear Kerr-materialWe have validated our extensions by simulating the self-focusing of a gauss-ian beam in non-linear space. In linear space, the injection of a Gaussianﬁeld proﬁle by the total ﬁeld/scattered ﬁeld formalism results in a diffrac-tion pattern. Both effects can be seen in Fig. 7. The Gaussian ﬁeld proﬁleis described by e=1.0 e x22∗0.252 Vµm. The material is described by the refractiveindex n=3.6 and n2=. 0.2m2V 2 .Eigenmode Expansion would require a large amount of eigenmodes tosimulate the same straight forward simulation setup. The lack of linearregions, which should only be calculated once, takes away the main advan-tage of non-linear EME.6. ConclusionOur extensions to the recently introduced complex Jacobi method allow thesimulation of 2D-components with Kerr-based materials in the frequencydomain. Our proposed extension to the complex Jacobi method is very wellsuited for structures where the non-linearity is present in a large portion ofthe simulation domain.AcknowledgementsBjorn Maes and Peter Bienstman acknowledge the Flemish Fund for Sci-entiﬁc Research (FWO-Vlaanderen) for their postdoctoral fellowships. Parts44 P. VANDERSTEEGEN ET AL.of this work were performed in the context of the Belgian DWTC projectIAP-Photon.ReferencesBurzler, J.M., S. Hughes and B.S. Wherrett. Appl. Phys. B 62(4) 389, 1996.Chew, W.C. and W.H. Weedon. Microwave Opt. Tech. Lett. 7 599, 1994.Hadley, R.G. J. Comput. Phys. 203(1) 358, 2005.Maes, B., P. Bienstman and R. Baets. Opt. Quan. Electron. 36(1–3) 15, 2005.Taﬂove, A. Computational Electrodynamics The Finite-Difference Time-Domain Method, Artech House,Boston-London, 111–117, 1995.

Using the complex Jacobi method to simulate Kerr non-linear photonic components.

https://biblio.ugent.be/publication/354296/file/1135954

Using the complex Jacobi method to simulate Kerr non-linear photonic components.

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography