## Abstract

A new three-dimensional finite-difference (FD)-based beam propagation method (BPM) is proposed for simulating optical propagation in weakly guiding waveguides with torsional birefringence, which cannot be simulated by existing FD-BPM algorithms. We also demonstrate that this new BPM algorithm is capable of obtaining eigenmode solutions in a helically-symmetric z-dependent waveguide structure with torsional birefringence.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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

*ϵ*representing the isotropic refractive index of the waveguides and the perturbation tensor $\tilde{\mathit{\u03f5}}$ representing the linear and torsional birefringences. So we have where the perturbation $\tilde{\mathit{\u03f5}}$ is a full tensor as

_{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

_{s}*-order small term for weakly guiding waveguides.*

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

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

where*σ*is torsional photoelastic constant, and

*τ*is twist rate. For convenience, we define Then, by taking Eqs. (5)–(8) into Eq. (4) and further dropping the 2

*-order small terms, we finally have our core algorithmic equation for wave equation as*

^{nd}*ϵ*

_{1}, Δ

*ϵ*

_{2}, and

*ϵ*

_{6}, 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

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*when rotating linear birefringence presents. Assuming

_{y}*A*(

*z*) and

*B*(

*z*) are envelope functions of fields

*E*and

_{x}*E*in local coordinates, following results are achieved [27]

_{y}*A*

_{0}=

*A*(0) and

*B*

_{0}=

*B*(0) and average propagation constant $\overline{\beta}$ and ratio X are defined as

Then, the envelope functions *A _{c}*(

*z*) and

*B*(

_{c}*z*) in the fixed laboratory Cartesian reference frame should take

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

^{−4}, simulation length at 2cm, and the initial value are

*A*

_{0}= 1 and

*B*

_{0}= 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.

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*of two orthogonal polarized modal fields should follow this coupled-mode equation [24]

_{y}*σ*and twist rate

*τ*, torsional birefringence will induce an optical activity G whose value is half of the value of

*η*as Eq. (8)

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*, and

_{yy}*e*during the high temperature drawing process [28, 29], and torsional birefringence is brought by shear strain

_{xy}*e*,

_{xz}*e*,

_{zx}*e*, and

_{yz}*e*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).

_{zy}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*= 100

_{y}*µm*and propagation distance is

*L*= 10

_{z}*cm*. 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.

## 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.

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.”