Finite-difference beam-propagation method for anisotropic waveguides with torsional birefringence

Xiuquan Ma; Shicheng Zhu; Li Li; Jinyan Li; Xinyu Shao; Almants Galvanauskas

doi:10.1364/OE.26.003995

1. Introduction

As one of the most powerful numerical modeling tools for optical waveguides, FD-BPM is widely used in designing new optical devices due to its excellent simulation speed with sufficient accuracy. So far, different kinds of FD-BPM algorithms are proposed to deal with different modeling needs. However, there are still some problems which demand FD-BPM’s speed but cannot be solved by existing FD-BPM algorithms, such as optical propagating in long fiber waveguides with torsional birefringence.

Recently, twisted optical waveguides have been receiving more attentions since they can be used in new types of fiber gratings [1], large mode area laser systems [2–5], fiber sensors [6], fiber dispersion controllers [7], spin-orbital couplers [8], and so on. However, when modeling optical propagation in twisted optical waveguides, dealing with torsional birefringence has always been difficult. There are only a few modeling tools available such as finite-difference time-domain (FDTD) method and finite-element (FE) BPM, which as we know are very time consuming and computationally expensive for simulating such long propagating problems.

The original FD-BPM algorithm based on the scalar wave equation ignores polarization dependence and cross-coupling effect [9]. Shortly afterwards, so-called “semi-vectorial” and “full-vectorial” FD-BPMs [10–13] are invented to solve optical propagation problems with polarization dependence and cross-coupling, which makes FD-BPM get more widely used. Based on above algorithmic frameworks, more advanced features of FD-BPM are invented, such as wide-angle formulation [14–17], bi-directional propagation [18–21] and so on. FD-BPM and its derivatives have become one of the most powerful modeling tools for simulating optical propagation in many different kinds of optical waveguides. However, for anisotropic waveguides with tensorial dielectric permittivities, FD-BPM algorithms are limited to deal with linear birefringence [22,23], while anisotropic waveguides with torsional birefringence still cannot be solved by FD-BPM algorithms. In this paper, we propose a new FD-BPM algorithm, in which both linear and torsional birefringences in weakly guiding anisotropic waveguides can be solved.

2. Algorithm derivation

We start from the wave equation for anisotropic waveguides

\nabla^{2} E + k_{0}^{2} {\tilde{ϵ}}_{r} E - \nabla [\nabla \cdot E] = 0 .

We separate the permittivity tensor

{\tilde{ϵ}}_{r}

into two parts: the scalar permittivity ϵ_s representing the isotropic refractive index of the waveguides and the perturbation tensor

\tilde{ϵ}

representing the linear and torsional birefringences. So we have

{\tilde{ϵ}}_{r} = ϵ_{s} + \tilde{ϵ},

where the perturbation

\tilde{ϵ}

is a full tensor as

\tilde{ϵ} = (\begin{matrix} Δ ϵ_{1} & ϵ_{6} & ϵ_{5} \\ ϵ_{6} & Δ ϵ_{2} & ϵ_{4} \\ ϵ_{5} & ϵ_{4} & 0 \end{matrix}) .

Now the wave equation can be written as

\nabla^{2} E + k_{0}^{2} {\tilde{ϵ}}_{r} E + \nabla [\frac{\nabla ϵ_{s}}{ϵ_{s}} \cdot E] + \frac{1}{ϵ_{s}} \cdot \nabla [\nabla (\tilde{ϵ} E)] = 0 .

Since ϵ_s is the scalar refractive index of the waveguides, the third term ∇ϵ_s/ϵ_s in Eq. (4) clearly describes the guiding effect of the waveguides, so we drop it since it should be a 2^nd-order small term for weakly guiding waveguides.

Then we focus on the fourth term $(1 / ϵ_{s}) \cdot \nabla [\nabla (\tilde{ϵ} E)]$ in Eq. (4). By dropping the 2^nd-order small terms, it can be approximated as

\frac{1}{ϵ_{s}} \cdot \nabla [\nabla (\tilde{ϵ} E)] \approx \frac{j β}{ϵ_{s}} \cdot \nabla [ϵ_{5} E_{x} + ϵ_{4} E_{y}] .

Meanwhile, for torsional birefringence induced by the shear strain, photoelastic theory gives [24]

ϵ_{4} = - σ τ y,

ϵ_{5} = σ τ x,

where σ is torsional photoelastic constant, and τ is twist rate. For convenience, we define

η = σ τ / ϵ_{s} .

Then, by taking Eqs. (5)–(8) into Eq. (4) and further dropping the 2^nd-order small terms, we finally have our core algorithmic equation for wave equation as

[\nabla^{2} + (\begin{matrix} k_{0}^{2} [ϵ_{s} + Δ ϵ_{1}] & k_{0}^{2} ϵ_{6} - j β η \\ k_{0}^{2} ϵ_{6} + j β η & k_{0}^{2} [ϵ_{s} + Δ ϵ_{2}] \end{matrix})] (\begin{matrix} E_{x} \\ E_{y} \end{matrix}) = 0 .

This is the form of wave equation we are going to use for weakly guiding anisotropic waveguides, where linear birefringence is included by Δϵ₁, Δϵ₂, and ϵ₆, and torsional birefringence is included by η. We need to point out that, this form of wave equation in Eq. (9) is energy conserved. The energy conservation is guaranteed by the 2 × 2 Hermitian matrix, which makes our BPM algorithm based on this equation working properly in terms of numerical simulation.

3. Algorithm verifications and a comparison between our algorithm and traditional FD-BPM algorithms

In the following discussion, we will examine a series of different scenarios in which our BPM algorithm can work properly. By adapting corresponding three-dimensional FD-BPM algorithm techniques such as the alternating-direction-implicit (ADI) method [25] and transparent boundary conditions [26], we obtain an efficient optical modeling tool whose calculation speed is comparable to commercial available FD-BPM-based software products. Comparisons between our algorithm with either existing algorithm formulations or analytical results will be presented, by doing which we can see that our BPM can deal with all kinds of tensorial dielectric permittivities very well including torsional birefringence.

First, we verify our algorithm with beam propagation problems in isotropic waveguides. As we can see, without linear and torsional birefringences, Eq. (9) goes

[\nabla^{2} + (\begin{matrix} k_{0}^{2} ϵ_{s} & 0 \\ 0 & k_{0}^{2} ϵ_{s} \end{matrix})] (\begin{matrix} E_{x} \\ E_{y} \end{matrix}) = 0 .

This is actually the scalar Helmholtz equation used in the scalar BPM for isotropic waveguides. It means our BPM goes back to the scalar BPM in this situation, so there is no need for more verification.

Second, we will benchmark our algorithm with rotational linear birefringence. In the previous work, researchers investigated the evolution of two orthogonal fields E_x and E_y when rotating linear birefringence presents. Assuming A(z) and B(z) are envelope functions of fields E_x and E_y in local coordinates, following results are achieved [27]

A (z) = [A (0) \cos (\sqrt{1 + X^{2}} τ z) + \frac{j X A (0) + B (0)}{\sqrt{1 + X^{2}}} \sin (\sqrt{1 + X^{2}} τ z)] e^{- j \bar{β} z},

B (z) = [B (0) \cos (\sqrt{1 + X^{2}} τ z) - \frac{j X B (0) + A (0)}{\sqrt{1 + X^{2}}} \sin (\sqrt{1 + X^{2}} τ z)] e^{- j \bar{β} z},

where we have A₀ = A(0) and B₀ = B(0) and average propagation constant

\bar{β}

and ratio X are defined as

\bar{β} = \frac{β_{x} + β_{y}}{2},

X = \frac{β_{x} - β_{y}}{2 τ} .

Then, the envelope functions A_c(z) and B_c(z) in the fixed laboratory Cartesian reference frame should take

A_{c} (z) = A (z) \cos τ z - B (z) \sin τ z,

B_{c} (z) = A (z) \sin τ z - B (z) \cos τ z .

Thus, we can make a comparison between our numerical modeling and equations above. For benchmarking, we take following parameters: wavelength at 1µm, helical pitch at 5mm, refractive index difference between the fast and slow axis of linear birefringence is dn_l = 2 × 10⁻⁴, simulation length at 2cm, and the initial value are A₀ = 1 and B₀ = 0. Then we plot both the simulation and analytical calculation in Fig. 1, in which red and blue curves are analytically calculated power evolutions of A_c(z) and B_c(z), blue and red circle points are simulated power evolutions of A_c(z) and B_c(z). We can see the simulation fits the analytical derivation very well. It means linear birefringence in our algorithm is indeed working properly.

Fig. 1 Simulation and analytical calculation for power evolutions of orthogonal fields E_x and E_y in waveguides with rotating linear birefringence. Red and blue curves are analytically calculated powers along the propagation direction. Blue and red circle points are simulated powers along the propagation direction.

Download Full Size | PDF

Third, we need to verify our algorithm in terms of torsional birefringence. From the previous work, we know that, given an optical fiber with only twisted torsional birefringence, amplitudes a_x and a_y of two orthogonal polarized modal fields should follow this coupled-mode equation [24]

\frac{d}{d z} [\begin{matrix} a_{x} \\ a_{y} \end{matrix}] = - j (\begin{matrix} 0 & - j η / 2 \\ j η / 2 & 0 \end{matrix}) [\begin{matrix} a_{x} \\ a_{y} \end{matrix}] .

Such a coupled-mode equation means that, given the torsional photoelastic constant σ and twist rate τ, torsional birefringence will induce an optical activity G whose value is half of the value of η as Eq. (8)

G = η / 2 .

We need to point out that, without linear birefringence, the coupled-mode equation of our algorithm derived from Eq. (9) shows the same formulation. Therefore, torsional birefringence is working properly as well.

The last and most critical question is whether our BPM algorithm works well for anisotropic waveguides with both linear and torsional birefringences. Actually, it has already been demonstrated: In our previous work, we presented a simulation for a Chirally-Coupled-Core (CCC) fiber [3], in which our FD-BPM algorithm successfully predicts all dips in the transmission spectrum. The CCC fiber structure consists of two waveguiding cores deposited in one glass cladding: a large central core and a small side core helically winding around the central core. Parameters of the CCC fiber are precisely selected so that high-order modes of central core can be coupled into side-core modes and then experience high bending loss of the side-core modes, while fundamental mode in the central core does not involve with any coupled effect, so those dips in the transmission spectrum of such a CCC fiber correspond to resonances among central-core fundamental mode and side-core lossy modes.

Here, to show that our BPM is the only FD-BPM algorithm that is capable of modeling waveguides with torsional birefringence, we compare our algorithm with two other FD-BPM algorithms. For CCC fibers, linear birefringence is caused by plane strain e_xx, e_yy, and e_xy during the high temperature drawing process [28, 29], and torsional birefringence is brought by shear strain e_xz, e_zx, e_yz, and e_zy caused by spinning the fiber during the high temperature drawing process [30,31]. The origins of linear and torsional birefringences are the differences of thermal-expansion coefficients and viscosity coefficients among the different doping areas such as central core, side core and cladding. The plane strain terms are well-known and easily calculable physical quantities, and we show the calculation results by “COMSOL Multi-Physics” in Figs. 2(a)–2(d). The strain energy distribution is also present in Fig. 2(e). We can see the plane strain is more heavily distributed in the region between two cores. The shear strain terms are not quantitatively calculable, however, we can legitimately assume these shear strain terms are also heavily distributed in the region between two cores because it is supposed to have the most dramatic change in viscosity coefficients as well. For simplicity, we assume that linear and torsional birefringences only exist in the elliptical black area between two cores and are evenly distributed, as shown in Fig. 2(f).

Fig. 2 COMSOL-calculated strain distributions in CCC fiber cross section and the simplified model for CCC fibers. (a)–(d), strain distributions of e_xx, e_yy, e_xy and e_xx e_yy. (e), strain energy distribution. (f), simpfied model of CCC fibers. The linear and−torsional birefringences are assumed to exist only in the elliptical black area and evenly distributed.

Download Full Size | PDF

Modeling results are shown in Fig. 3. Parameters of the fiber structure in this simulation example correspond to one of the fabricated fibers and are the following: central core diameter 35µm, side core diameter 13µm, central core numerical aperture (NA) 0.06, side core NA 0.1, core to core separation 3µm, and helical pitch 6.1mm. The discretizations in the simulations are set with Δx = Δy = 0.5µm and Δz = 20µm. The total transversal lengths of the computational windows are L_x = L_y = 100µm and propagation distance is L_z = 10cm. As we can see, only our algorithm can predict all resonances of the CCC fiber. The scalar BPM for isotropic waveguides can only predict scalar resonances marked with green lines in Fig. 3(a), while the vectorial BPM for anisotropic waveguides with linear birefringence can predict one linear birefringence resonance marked with red dash line and two scalar resonances in Fig. 3(b). Torsional birefringence resonances marked with blue dash lines can only be seen with our algorithm in Fig. 3(c). Also, it is worth noting that simulation time for such a 10cm fiber is about 28 seconds with our simulation program written in C++, which demonstrates the excellent computational performance of our algorithm in terms of the speed.

Fig. 3 Transmission spectrums in a CCC fiber simulated by different BPM algorithms: (a) Scalar BPM for isotropic waveguides shows only two scalar resonances marked with green dash lines. (b) Vectorial BPM for anisotropic waveguides with linear birefringence shows three resonances including one linear birefringence resonance marked with red dash line and two scalar resonances. (c) Our self-developed BPM for anisotropic waveguides with both linear and torsional birefringences shows six resonances including three torsional birefringence resonances marked with blue dash lines, one linear birefringence resonance and two scalar resonances.

Download Full Size | PDF

4. Eigenmodes calculation in helically-symmetric z-dependent waveguides with our algorithm

Moreover, we develop the technique to obtain eigenmode solutions in helically-symmetric z-dependent waveguides. Traditional BPM algorithms with correlation method [32] can be used in finding eigenmodes for z-independent optical waveguides since the definition for waveguide eigenmodes requires the profiles only depend on cross-section coordinates. Based on our BPM algorithm, we implement the capability of correlating the waveguide propagation in the helical coordinate system. Even though we have pointed it out in our previous work [3], it is still quite shocking to see that the directly solved eigenmodes with our BPM-based mode solving technique in such a CCC structure are the ones carrying optical angular momentums as shown in Fig. 4. Figures 4(a) and 4(b) are the amplitude and intensity profiles of BPM solved supermode in the CCC structure, while Fig. 4(c) is the phase distribution in the central core. Clearly, the ring-like profile of intensity and spiral wavefront show that the eigenmodes propagating in such a structure are carrying orbital angular momentums.

Fig. 4 Comparison of supermodes calculated by FEM and BPM in the CCC fiber. (a) and (b), amplitude and intensity solutions calculated by BPM. (c), phase distribution in the central core calculated by BPM. (d) and (e), amplitude and amplitude solutions calculated by FEM. (f), phase distribution in the central core calculated by FEM.

Download Full Size | PDF

For validity and comparison, we also performed Finite Element Method (FEM) simulation for solving eigenmodes in the static helicoidal coordinate system, which is discussed in more details in another publication [33]. Figures 4(d) and 4(e) are amplitude and intensity profiles of the solved eigenmodes by FEM, and Fig. 4(f) is phase distribution in the central core. As we can see, eigenmodes obtained by our BPM algorithm match our FEM calculation very well. This is, to the best of our knowledge, for the first time that both BPM propagation calculation and FEM modal-solving calculation are agreeing with each other to show the visualization of eigenmodes carrying optical angular momentums. Together with the theoretical framework carried out in [3] and [33], the agreeing simulation results by both BPM and FEM construct a self-consistent system of proving that eigenmodes propagating in the helically-symmetric structures such as CCC fibers are the ones carrying angular momentums.

5. Conclusion

In conclusion, the algorithm derivation, formulation, and implementation of a new FD-BPM algorithm capable of solving weakly guiding optical waveguides with torsional birefringence are presented in this paper. To verify the validity of this new algorithm, four different scenarios are discussed in details respectively. Particularly, by comparing our algorithm with existing FD-BPMs for simulating a CCC fiber, we show that only our BPM can predict all the resonances in the transmission spectrum. Moreover, this BPM can be extended to directly solve eigenmodes in such helically-symmetric z-dependent structures. This BPM mode solving technique, together with the FEM mode solving technique, are compared and combined to strongly prove that eigenmodes propagating in helically-symmetric structures are the ones carrying angular momentums.

Funding

National Basic Research Program of China (2014CB046703).

References and links

1. V.I. Kopp, V.M. Churikov, J. Singer, N. Chao, D. Neugroschl, and AZ. Genack, “Chiral Fiber Gratings,” Science 305, 74–75 (2004). [CrossRef] [PubMed]

2. X. Ma, I.-N. Hu, and A. Galvanauskas, “Propagation-length independent SRS threshold in chirally-coupled-core fibers,” Opt. Express 19(23), 22575–22581 (2011). [CrossRef] [PubMed]

3. X. Ma, C.-H. Liu, G. Chang, and A. Galvanauskas, “Angular-momentum coupled optical waves in chirally-coupled-core fibers,” Opt. Express 19(27), 26515–26528 (2011). [CrossRef]

4. X. Ma, C. Zhu, I.-N. Hu, A. Kaplan, and A. Galvanauskas, “Single-mode chirally-coupled-core fibers with larger than 50µm diameter cores,” Opt. Express 22(8), 9206–9219 (2014). [CrossRef] [PubMed]

5. S. Lefrancois, C.-H. Liu, M. L. Stock, T. S. Sosnowski, A. Galvanauskas, and F. W. Wise, “High-energy similariton fiber laser using chirally coupled core fiber,” Opt. Lett. 38(1), 43–45 (2013). [CrossRef] [PubMed]

6. X. M. Xi, G. K. L. Wong, T. Weiss, and P. St. J. Russell, “Measuring mechanical strain and twist using helical photonic crystal fiber,” Opt. Lett. 38(24), 5401–5404 (2013). [CrossRef] [PubMed]

7. G. K. L. Wong, M. S. Kang, H. W. Lee, F. Biancalana, C. Conti, T. Weiss, and P. S. J. Russell, “Excitation of orbital angular momentum resonances in helically twisted photonic crystal fiber,” Science 337(6093), 446–449 (2012). [CrossRef] [PubMed]

8. X. M. Xi, G. K. L. Wong, M. H. Frosz, F. Babic, G. Ahmed, X. Jiang, T. G. Euser, and P. St. J. Russell, “Orbital-angular-momentum-preserving helical Bloch modes in twisted photonic crystal fiber,” Optica 1(3), 165–169 (2014). [CrossRef]

9. S. T. Hendow and S. A. Shakir, “Recursive numerical solution for nonlinear wave propagation in fibers and cylindrically symmetrical systems,” Appl. Opt. 25, 1759–1764 (1986). [CrossRef]

10. W. P. Huang, C. L. Xu, S. T. Chu, and S. K. Chaudhuri, “A vector beam propagation method for guided-wave optics,” IEEE Photonics Technol. Lett. 3, 910–913 (1991). [CrossRef]

11. W. P. Huang, C. L. Xu, S. T. Chu, and S. K. Chaudhuri, “The finite-difference vector beam propagation method. Analysis and Assessment,” J. Lightwave Technol. 10, 295–305 (1992). [CrossRef]

12. P. L. Liu and B. J. Li, “Semivectorial beam propagation by Lanczos reduction,” IEEE J. Quantum Electron. 29, 2385–2389 (1993). [CrossRef]

13. W. P. Huang and C. L. Xu, “Simulation of three-dimensional optical waveguides by a full-vector beam propagation method,” IEEE J. Quantum Electron. 29, 2639–2649 (1993). [CrossRef]

14. D. Yevick and M. Glasner, “Analysis of forward wide-angle light propagation in semi-conductor rib waveguides and integrated-optic structures,” Electron. Lett. 25, 1611–1613 (1989). [CrossRef]

15. G. R. Hadley, “Wide-angle beam propagation using Padé approximant operators,” Opt. Lett. 17, 1426–1428 (1992). [CrossRef]

16. H. J. W. M. Hoekstra, G. J. M. Krijnen, and P. V. Lambeck, “New formulation of the beam propagation method based on the slowly varying envelope approximation,” Opt. Commun. 97, 301–303 (1993). [CrossRef]

17. Y. Chung and N. Dagli, “A wide-angle propagation technique using an explicit finite-difference scheme,” IEEE Photonics Technol. Lett. 6, 540–542 (1994). [CrossRef]

18. P. Kaczmarski and P. E. Lagasse, “Bidirectional beam propagation method,” Electron. Lett. 24, 675–676 (1988). [CrossRef]

19. Y. Chung and N. Dagli, “Modeling of guided-wave optical components with efficient finite-difference beam propagation methods,” IEEE Antennas Propag. Soc. Int. Symp. 1, 248–250 (1992).

20. Y. P. Chiou and H. C. Chang, “Analysis of optical waveguide discontinuities using the Padé approximants,” IEEE Photonics Technol. Lett. 9, 964–966 (1997). [CrossRef]

21. H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. 11, 830–832 (1999). [CrossRef]

22. C. L. Xu, W. P. Huang, J. Chrostowski, and S. K. Chaudhuri, “A full-vectorial beam propagation method for anisotropic waveguides,” IEEE J. Lightwave Technol. 12, 1926–1931 (1994). [CrossRef]

23. Q. Wang, G. Farrell, and Y. Semenova, “Modeling liquid-crystal devices with the three-dimensional full-vector beam propagation method,” J. Opt. Soc. Am. 23, 2014–2019 (2006). [CrossRef]

24. R. Ulrich and A. Simon, “Polarization optics of twisted single-mode fibers,” Appl. Opt. 18, 2241–2251 (1979). [CrossRef] [PubMed]

25. J. Yamauchi, T. Ando, and H. Nakano, “Beam propagation analysis of optical fibres by alternating direction implicit method,” Electron. Lett. 27(18), 1663–1666 (1991). [CrossRef]

26. G. R. Hadley, “Transparent boundary condition for beam propagation,” Opt. Lett. 16, 624–626 (1991). [CrossRef] [PubMed]

27. P. McIntyre and A. W. Snyder, “Light propagating in twisted anisotropic media: Application to photoreceptors,” J. Opt. Soc. Am. 68, 149–156 (1978). [CrossRef] [PubMed]

28. H. Tai and R. Rogowski, “Optical anisotropy induced by torsion and bending in an optical fiber,” Opt. Fiber Technol. 8, 162–169 (2002). [CrossRef]

29. A. D. Yablon, “Optical and mechanical effects of frozen-in stresses and strains in optical fibers,” IEEE J. Sel. Top. Quantum Electron. 10, 300–311 (2004). [CrossRef]

30. M. J. Li, X. Chen, and D. A. Nolan, “Effects of residual stress on polarization mode dispersion of fibers made with different types of spinning,” Opt. Lett. 29, 448–450 (2004). [CrossRef] [PubMed]

31. S. M. Pietralunga, M. Ferrario, M. Tacca, and M. Martinelli, “Local Birefringence in Unidirectionally Spun Fibers,” J. Lightwave Technol. 24, 4030–4038 (2006). [CrossRef]

32. M. D. Feit and J. A. Fleck, “Computation of mode properties in optical fiber waveguides by a propagating beam method,” Appl. Opt. 19, 1154–1164 (1980). [CrossRef] [PubMed]

33. L. Li, School of Mechanical Science and Engineering, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China, and S. Zhu, J. Li, X. Shao, A. Galvanauskas, and X. Ma are preparing a manuscript to be called “All-in-fiber method of generating orbital angular momentum with helically-symmetric fibers.”

Finite-difference beam-propagation method for anisotropic waveguides with torsional birefringence

Abstract

1. Introduction

2. Algorithm derivation

3. Algorithm verifications and a comparison between our algorithm and traditional FD-BPM algorithms

4. Eigenmodes calculation in helically-symmetric z-dependent waveguides with our algorithm

5. Conclusion

Funding

References and links

Cited By

Figures (4)

Equations (18)

Optics Express