# Nonlinear Coherent Directional Coupler: Coupled Mode Theory and BPM Simulation

Authors: Dharmadas Kumbhakar
Journal: International Journal of Optics
 (2012)

---

## Abstract

Finite difference beam propagation method is an accurate numerical procedure, used here to explore the switching dynamics of a nonlinear coherent directional coupler. The coupling lengths derived from this simulation are compared with coupled mode theories. BPM results for the critical power follow the trend of the coupled mode theories, but it lies in between two coupled mode theories. Coupled mode theory is sensitive to numerical approximations whereas BPM results practically do not depend on grid size and longitudinal step size. Effect of coupling-region-width and core-width variations on critical power and coupling length is studied using BPM to look at the aspects of optical power-switch design.


---

## Body

## 1. Introduction

Beam propagation method [1] is a versatile tool to investigate or model various optical phenomena in photonic devices. Earlier, the method was popular to find the guided modes in dielectric waveguide structures [2]; even in microstructure fibers like holey fibers it can be used to find the modes [3]. With suitable boundary conditions, it gave proper estimates of the propagation constant in leaky structures like arrow waveguides [4]. The method was used in linear directional couplers to simulate switching of a CW signal and filtering TE/TM mode signals [5]. With growing interest in nonlinear fibers, it was the only tool to model the pulse propagation [6] in conjugation with split-step Fourier transform. It has been also used to model various other nonlinear phenomena like second harmonic generation [7], leakage loss in buried silicon substrate [8], and so forth. Beam propagation method can be incorporated either in the finite element or finite difference framework. In cases where the interfaces are parallel to the axes, the two methods are almost the same in terms of accuracy, but the latter is sometimes preferred due to its simplicity in implementation.Nonlinear directional coupler (NLDC) is a very useful device in photonic circuits; it can be used as all optical modulator, switches, logic gates, and so forth. Jenson [9] first showed that total power exchange in a linear coherent coupler is lost in a nonlinear coherent coupler above a critical input power in one guide. This fact was used to design power-controlled all optical switches [10, 11]. The device can also be used as phase-controlled switches [12] and optical modulator [13]. So, the knowledge of the dependence of the critical power and phase variations on the NLDC parameters is important. There are numerous studies on this device; analytic methods include coupled mode theory [9, 14], variational method [15], and power flow portrait [16]. Numerical studies are mainly based on finite element [17, 18] or finite difference beam propagation methods [19–21]. There are very few experimental results reported so far (see [22] and references therein).In some works, explicit FDBPM is used to model some nonlinear phenomena [18, 20], we preferred the implicit formulation, and it has been proved to be robust and fast. In this paper, we have simulated 2D nonlinear coherent directional coupler with implicit finite difference beam propagation method (IFDBPM). The results are compared with coupled mode theory. Some aspects of power-switching design are also investigated.


## 2. Coupled Mode Theory

Figure1 shows the structure of a 2D symmetric or coherent directional coupler. The refractive index is given by(1)ndc2(x)=n12for|x-xc±(w2+d2)|≤d2=n22elsewhere.Figure 1
Structure of a coherent directional coupler.We assume that the coupling is small so that the power flow in the two waveguide components can be obtained from the individual waveguide modes. The individual waveguide structure is obtained removing the second core region. The index distribution in individual guides is given byguidea(2a)nw2(x)=n12for|x-xc+(w2+d2)|≤d2=n22elsewhere,
guide b(2b)nw2(x)=n12for|x-xc-(w2+d2)|≤d2=n22elsewhere.The field evolution along the length of the directional coupler can be expressed as(3)Ψ(x,z)=[A(z)ψa(x)+B(z)ψb(x)]e-iβz.ψa(x) and ψb(x) are the normalized mode fields of individual waveguides a and b, respectively. A(z) and B(z) are the amplitudes of the fields in waveguide a and b, respectively. Using coupled mode theory, one gets (4a)dAdz=i[K0A+K1B+(K2|A|2+2K3|B|2)A],(4b)dBdz=i[K0B+K1A+(K2|B|2+2K3|A|2)B],where(5a)K0=k022β∫bΔn2|ψa(x)|2dx,(5b)K1=k022β∫aψa(x)Δn2ψb(x)dx,(5c)K2=k022β∫α|ψa(x)|4dx,(5d)K3=k022β∫α|ψa(x)|2|ψb(x)|2dx,k0=2π/λ, λ is the free space wavelength and β is the propagation constant in the coupler. α=ε0c0n2n2′, ε0 is the free space electric permittivity, and c0 is  the velocity of light in free space; n2′ is the nonlinear coefficient of the nonlinear medium either core or cladding or both as the case may be and n is the linear refractive index of that medium. Δn2 is the perturbing refractive index which will be explained in Section 4.The indexa or b in the integrals indicates integration within the core region of waveguide a or b. K0 is the contribution of the evanescent field of one waveguide in the other; it just modifies the propagation constant β [23]. K1 is the overlap of the evanescent field with the field into the other; it has the pivotal role in the exchange of power in a linear directional coupler. K2 is the contribution due to self phase modulation and K3 is due to cross-phase modulation.For weak coupling, one can neglect the small contributionsK0 and K3; this is justified with the evaluation of these parameters in Section 4.IfP0 is the initial power fed into the waveguide a, P0=FA(0)2, where F=(1/2)ε0c0ne, ne is the effective index of the waveguide. Then the power P(z) in that guide at a distance z from the input end can be expressed as(6)(dPdz)2+4K1P(P0-P)-K2′P2(P0-P)2=0,
where K2′=K2/F. The above equation is obtained from (4a) and (4b) and using P0=F(|A|2+|B|2)=constant, assuming the coupler to be loss-free.Equation (6) can be easily converted to Jacobi elliptic equation(7)(dP′dz′)=(1-P′2)(1-m+mP′2),
where P′=(P-P0/2)/P0/2, z′=2K1z, and m=(P0K2′/4K1)2.Form≤1, the solution of (7) is(8)P=12P0[1+cn(2K1z∣m)].cn is the relevant Jacobi elliptic function. The value of the input power P0 for which m=1 is called the critical power. So, critical power Pc = 4K1/K2′.The coupling lengthL is the minimum value of z at which power in the waveguide falls from initial value P0 to zero. So,(9)cn(2K1L∣m)=-1.For linear directional coupler,m=0, the coupling length L0 is given by(10)cos⁡(2K1L0)=-1orL0=π2K1.Above the critical input power(m>1), (7) can be written as(11)(dP′d(mz′))=(P′2-1)(1-m-1-P′2).
Its solution is
(12)P=12P0[1+dn(2K1mz∣m-1)].


## 3. Beam Propagation Method

In the linear regime, the paraxial scalar wave equation in finite difference form results in a linear matrix equation; evolution of the field in thez-direction is then obtained successively solving the linear equation through small steps. With Kerr or non-Kerr-type nonlinearity, the equation is not linear. However, the nonlinear contribution in the paraxial wave equation, which contains the electric field value, can be used to modify the refractive index profile. The electric field can then be obtained as in the linear case. The modal field and propagation constant of a three-layer or multilayer waveguide can then be obtained in this process with finite difference inverse iteration method or shifted power method or imaginary distance beam propagation method. But when field evolution is required as in a directional coupler, to keep or develop the proper phase of the fields, the successive solution of the original paraxial wave equation with imaginary quantities is essential.Inz-invarinat or slow z-variant structures, the light propagation is given by the paraxial wave equation(13)2iβ∂Ψ∂z=∂2Ψ∂x2+k02(ndc2(x)+α|Ψ|2-β2)Ψ.Ifψn(x) is the field distribution in the directional coupler at the nth step, in the (n+1)th step, the index distribution can be taken to be effectively(14)n2(x)=ndc2(x)+α|ψn(x)|2.
The field distribution at the (n+1)th step can be expressed as(15)ψn+1=(1+GΔz2iβ)(ξψn+1+(1-ξ)ψn).
In finite difference form, x is discretised into N segments; ψn is an N-dimensional column matrix and G is an N×N matrix operator. n2 is revised at each step with the latest ψn using (14) and then ψn is updated through (15). ξ is a parameter required for numerical convergence. In Crank-Nicholson scheme, it is 0.5. ξ<0.5 is called explicit formulation and ξ≥0.5 is called implicit formulation. In fully explicit form (ξ=0), ψn+1 can be obtained directly with (15). This requires very small steps of Δz for convergence and consequently long runtime. In implicit formulation, a matrix inversion is necessary at each step; however, it is very stable and large Δz can be used [24]. Matrix inversion is very easy and fast using the sparse matrix techniques in Matlab; for that reason we have used implicit formulation for fast and stable computation.


## 4. Coupling Length Variations with Power

Coupled mode theory [9] considered only the linear perturbation contribution in Δn2(16)Δn2=(ndc2-nw2),
and the individual waveguide modes were also taken to be the solution of Helmholtz equation without the nonlinear term:(17)d2ψdx2+k02{nw2(x)-ne2}ψ=0.ne is the effective index of the waveguide.Power flow in the coherent nonlinear directional coupler was derived from coupled mode equations.Later improved coupled mode theory [14] derived the power flow equations from Lorentz reciprocity theorem and used the nonlinear perturbation in Δn2(18)Δn2=(ndc2-nw2)+α|ψb|actual2.
Individual waveguide modes were obtained using the nonlinear Helmholtz equation(19)d2ψdx2+k02{nw2(x)-ne2+α|ψ|2}ψ=0.
Actual field |ψ|actual and normalized ψ are related as(20)|ψ|actual=ψPε0c0ne.Two couplers with different nonlinearity configurations are used in this paper to show the versatility of our BPM simulation model.The mode fields are determined with highly accurate finite difference inverse iteration method which we have developed and coded in Matlab. The transverse grid size used in this evaluation is 0.001μm and it takes only a few seconds. These fields are used to find the coupling constants K1 and K2.K0,K3 are found to be two orders smaller in magnitude than K1,K2, respectively. One can solve the coupled equations (4a) and (4b) numerically; it shows that these terms practically have no effect on coupling length and critical power within the range of power level we consider.A smaller transverse grid size like 0.04μm is used to determine critical power with BPM simulation. A longitudinal step size of 1 μm is found to be enough. Refining transverse grid size to 0.01 and longitudinal step size to 0.1 μm does not produce any significant change in critical power estimate.


### 4.1. Nonlinear Coupling Region

The 2D directional coupler parameters are taken as in [14]: n1 = 1.57; n2 = 1.55; λ= 1.064 μm; d = 2 μm; w= 2.4 μm; α = 6.3786 m2/V. Central cladding layer is taken to be nonlinear.Figure2 shows the variation of coupling length in this coupler as a function of input power as obtained from two coupled mode theories and our BPM simulation. The results of two coupled mode theories are widely different; the BPM simulation results follow the trend of improved coupled mode theory [14]; the critical power is slightly higher in BPM.Figure 2
Comparison of coupling length variation with power: coupled mode theory 1 (Jenson): dashed curve, coupled mode theory 2 (Okamoto): dotted curve, our BPM simulation: solid curve.


### 4.2. Nonlinear Core Regions

We have also studied a directional coupler with nonlinear core regions. The coupler parameters aren1 = 1.56, n2 = 1.55, d = 3 μm, W = 5 μm, λ = 1.3 μm, and α = 0.15 m2/V. The variation of coupling length with power for this coupler is shown in Figure 3.Figure 3
Power versus coupling length. Line properties are same as in Figure2.The critical power in our BPM simulation results lies in between the results of two coupled mode theories. For high nonlinearity, BPM simulation curve for the variation of coupling length with power shifts towards improved coupled mode theory and, for low nonlinearity, it shifts towards the other coupled mode theory.Although BPM algorithm is quite different from coupled mode theories, the functional form of the power flow in the coupler as obtained with both the approaches is very similar, apart from some shift as input power is increased. To justify this, we compare BPM simulation and coupled mode theory with the coupler structure of Section4.2. Figure 4(a) shows the BPM simulation picture of power evolution in the coupler for low input power (0.1 W/m). In this low power regime, only one coupling constant K1 describes power flow:(21)K1=π2L0,
where L0 is the low power coupling length. Figure 4(b) shows the plot of the function P=P0[1+cos(2K1z)] and the computed total power in the guide a from the field evolution obtained with BPM simulation shown in Figure 4(a). As revealed in Figure 4(b), agreement is excellent.(a) Field evolution at 0.1 W/m input power: BPM. (b) Power flow pattern, (line curve) computed from BPM shown left, (crossed points) coupled mode theory. (c) Field evolution at 70 W/m power: BPM. (d) Power flow pattern (line curve), computed from BPM shown left, (crossed points) coupled mode theory fitted. (e) Field evolution at 100 W/m power: BPM. (f) Power flow pattern, (line curve) computed from BPM shown left, (crossed points) coupled mode theory fitted.
(a)(b)(c)(d)(e)(f)Figure4(c) shows the field evolution with BPM for an input power 70 W/m, near but below the critical power (75 W/m). The nature of the power flow in BPM is different from cosine function as the look of the field evolution (Figure 4(c)) tells. The power flow pattern according to coupled mode theory is (8), that is, P=(1/2)P0(1+cn(2K1z∣m)). We check whether power flow obeys this relation. We take a slightly different (lower) value of K1, than given by (21); find K2, m from K1 and critical power determined from BPM. Putting these values in (8), plot of P-z clearly matches the BPM power flow, shown in Figure 4(d). It clearly shows that power flow happens to be of cn nature below critical power.Figure4(e) shows the field evolution at a power 100 W/m which is above the critical power. The power is almost contained in the guide a; almost 90% on the average. Only about 10% power is exchanged between the two arms of the coupler. The match of the BPM power flow curve with the coupled mode theory equation (12) in Figure 4(f) verifies the fact that the undulations in the peak power are of dn nature according to (12).In the experiments [22] on the nonlinear switching of directional couplers, GaAs/AlGaAs multiple quantum well structures were used and operated at wavelengths near the band edge to tap sufficient Kerr nonlinearity due to exitonic resonances. However, these structures exhibit large losses. To reduce the losses, later, only the coupling region was made nonlinear. Even with such improved structures, the loss is of the order of 10 cm−1. With high initial powers, although straight-through condition is achieved, some power is coupled to the other guide due to power loss in the lossy nonlinear medium. The experimental data reveals that power flow is different from Jenson’s theory but it is insufficient to differentiate further.In coupled mode theory, the coupling constantsK0-K3 are determined with the solution of the Helmholtz equation in individual wave guide ψa and ψb at a fixed power. However, the power in each waveguide gets changing all along the length due to power transfer. Coupled mode theory did not incorporate the phase and field overlap changes continuously. But in the BPM, at any step, the field is computed with the instantaneous power. So, it incorporates the power and phase change continuously. Power transfer is likely to be more complicated than that predicted by existing theoretical frameworks. Some mechanism to get proper coupling constants in the coupled mode theory may be one step forward.


## 4.1. Nonlinear Coupling Region

The 2D directional coupler parameters are taken as in [14]: n1 = 1.57; n2 = 1.55; λ= 1.064 μm; d = 2 μm; w= 2.4 μm; α = 6.3786 m2/V. Central cladding layer is taken to be nonlinear.Figure2 shows the variation of coupling length in this coupler as a function of input power as obtained from two coupled mode theories and our BPM simulation. The results of two coupled mode theories are widely different; the BPM simulation results follow the trend of improved coupled mode theory [14]; the critical power is slightly higher in BPM.Figure 2
Comparison of coupling length variation with power: coupled mode theory 1 (Jenson): dashed curve, coupled mode theory 2 (Okamoto): dotted curve, our BPM simulation: solid curve.


## 4.2. Nonlinear Core Regions

We have also studied a directional coupler with nonlinear core regions. The coupler parameters aren1 = 1.56, n2 = 1.55, d = 3 μm, W = 5 μm, λ = 1.3 μm, and α = 0.15 m2/V. The variation of coupling length with power for this coupler is shown in Figure 3.Figure 3
Power versus coupling length. Line properties are same as in Figure2.The critical power in our BPM simulation results lies in between the results of two coupled mode theories. For high nonlinearity, BPM simulation curve for the variation of coupling length with power shifts towards improved coupled mode theory and, for low nonlinearity, it shifts towards the other coupled mode theory.Although BPM algorithm is quite different from coupled mode theories, the functional form of the power flow in the coupler as obtained with both the approaches is very similar, apart from some shift as input power is increased. To justify this, we compare BPM simulation and coupled mode theory with the coupler structure of Section4.2. Figure 4(a) shows the BPM simulation picture of power evolution in the coupler for low input power (0.1 W/m). In this low power regime, only one coupling constant K1 describes power flow:(21)K1=π2L0,
where L0 is the low power coupling length. Figure 4(b) shows the plot of the function P=P0[1+cos(2K1z)] and the computed total power in the guide a from the field evolution obtained with BPM simulation shown in Figure 4(a). As revealed in Figure 4(b), agreement is excellent.(a) Field evolution at 0.1 W/m input power: BPM. (b) Power flow pattern, (line curve) computed from BPM shown left, (crossed points) coupled mode theory. (c) Field evolution at 70 W/m power: BPM. (d) Power flow pattern (line curve), computed from BPM shown left, (crossed points) coupled mode theory fitted. (e) Field evolution at 100 W/m power: BPM. (f) Power flow pattern, (line curve) computed from BPM shown left, (crossed points) coupled mode theory fitted.
(a)(b)(c)(d)(e)(f)Figure4(c) shows the field evolution with BPM for an input power 70 W/m, near but below the critical power (75 W/m). The nature of the power flow in BPM is different from cosine function as the look of the field evolution (Figure 4(c)) tells. The power flow pattern according to coupled mode theory is (8), that is, P=(1/2)P0(1+cn(2K1z∣m)). We check whether power flow obeys this relation. We take a slightly different (lower) value of K1, than given by (21); find K2, m from K1 and critical power determined from BPM. Putting these values in (8), plot of P-z clearly matches the BPM power flow, shown in Figure 4(d). It clearly shows that power flow happens to be of cn nature below critical power.Figure4(e) shows the field evolution at a power 100 W/m which is above the critical power. The power is almost contained in the guide a; almost 90% on the average. Only about 10% power is exchanged between the two arms of the coupler. The match of the BPM power flow curve with the coupled mode theory equation (12) in Figure 4(f) verifies the fact that the undulations in the peak power are of dn nature according to (12).In the experiments [22] on the nonlinear switching of directional couplers, GaAs/AlGaAs multiple quantum well structures were used and operated at wavelengths near the band edge to tap sufficient Kerr nonlinearity due to exitonic resonances. However, these structures exhibit large losses. To reduce the losses, later, only the coupling region was made nonlinear. Even with such improved structures, the loss is of the order of 10 cm−1. With high initial powers, although straight-through condition is achieved, some power is coupled to the other guide due to power loss in the lossy nonlinear medium. The experimental data reveals that power flow is different from Jenson’s theory but it is insufficient to differentiate further.In coupled mode theory, the coupling constantsK0-K3 are determined with the solution of the Helmholtz equation in individual wave guide ψa and ψb at a fixed power. However, the power in each waveguide gets changing all along the length due to power transfer. Coupled mode theory did not incorporate the phase and field overlap changes continuously. But in the BPM, at any step, the field is computed with the instantaneous power. So, it incorporates the power and phase change continuously. Power transfer is likely to be more complicated than that predicted by existing theoretical frameworks. Some mechanism to get proper coupling constants in the coupled mode theory may be one step forward.


## 5. Effect of Core and Coupling Width Variations

A nonlinear directional coupler can find application as a power-controlled switch. If power is input at the waveguidea, with proper choice of the length of the directional coupler, at low powers the optical power is switched out through the other guide b; above a critical power, the optical power is retained and out through the guide a itself.Obviously, the nonlinear medium—either the coupling region or core regions—must have high Kerr-type nonlinearity to operate the coupler as an optical switch at moderately low critical power. The coupling length of a coupler is also an important factor to optimize the size of the device. So, one must look at the coupler dimensions which govern these two quantities—critical power and coupling length. We have studied the variation of these two quantities with the change of the width of the coupling region and the width of the core regions for both the types of nonlinear couplers studied in Section4. Figure 5(a) shows the variation of critical power with the change in w (width of the coupling region), core width d as a parameter—for the coupler with nonlinear coupling region (Section 4.1). Figure 5(b) shows those variations for the nonlinear core coupler (Section 4.2). As a general rule, the critical power falls with increase in the separation between the guides. If d is kept fixed, K2 remains unchanged as the mode profile ψa (or ψb) does not change (see (5c)); however, as w is increased, the value of K1 decreases due to reduced overlap of the fields (see (5b)). This is the reason that critical power  (Pc=4K1/K2′) decreases with increasing w.(a) Variation of critical power with core separation: nonlinear coupling region. (b) Variation of critical power with core separation: nonlinear core regions.
(a)(b)But the effect of core-widthd variation is different for two types of nonlinear couplers. With the decrease in d, the field inside the core decreases in the first type of coupler—where the coupling region is nonlinear. The decrease of field inside the core leads to field spreading that increases the overlap to some extent, but more is the increase in K2 as the field increases in the nonlinear region. For the other type, just the opposite happens—field inside the core increases with decrease in d. That shrinks the fields, increasing the overlap but keeping K2 almost the same.Figures6(a) and 6(b) show the variations in the low power coupling length. As the separation between the coupler arms increases, K1 decreases due to reduction in the field overlap. Coupling length (π/2K1) increases as a result whatever the value of d is.(a) Variation of coupling length with core separation: nonlinear coupling region. (b) Variation  of coupling length with core separation: nonlinear core regions.
(a)(b)To summarize, critical power decreases as the core separation is increased in both the two types of couplers we consider. But as the core separation is increased, coupling length also increases. So, to contain the size of the switching device, one has to optimize between the critical power and the coupling length.


## 6. Conclusion

We have developed an implicit BPM algorithm which gives the effective power flow in a nonlinear directional coupler, easily and correctly. The method is fairly insensitive to transverse and longitudinal step size. It appears to be more effective than the coupled mode theories which use some approximations. Hence IFDBPM should find wide applications in investigation of NLDC.



---

*Source: 173250-2012-04-11.xml*