Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Computational methods in applying the angular spectrum algorithm to optical fibers

Open Access Open Access

Abstract

We have been investigating the angular spectrum method for analysis of single-mode and multimode optical fibers. The justification is simple because the approach requires less computational time than other methods because it is heavily dependent on the fast Fourier transform. This method should become much more important as new approaches are being studied for higher data rate transmission. We examine two basic components in the application of this technique to propagation of light in optical fibers – an absorbing layer at the outer radius and design of the interface between the cord and cladding of the optical fiber.

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

1. Introduction

In recent work, we have studied the angular spectrum-based wave propagation method [1] for analysis of optical fibers using a random array of cores [2], for single mode optical fibers and different optical coupler geometries [3], and for multimode fibers [4]. Our motivation is simple - this approach is computationally fast because it makes use of the fast Fourier transform (FFT) twice in each iteration. It is difficult to justify this claim because most papers do not give specific computational times, nor do they discuss their computer programs. Our only comparison was in [2] where others [5] reported times of 1000 hours using the finite difference beam propagation method (FD-BPM) [6] and multicore computers. By comparison, using an iMac with a 2.66 GHz Intel Core 2 Duo processor, a computation of 100,000 steps only took about 15 minutes!

This computational time is becoming much more relevant as optical fiber networks are reaching capacity limits [7]. A variety of new approaches are being investigated including multimode optical fibers, some with elliptical cores [8,9], and vortex or vector beams that carry orbital angular momentum and spiral phase [10,11]. More complex coupling schemes would be required including spatial mode multiplexers [12], broadband polarization splitters [13] and photonic lanterns [14,15]. However, study of these new approaches would benefit from accurate and fast computational techniques.

In this work, we show the effects of adding an absorbing boundary layer at the outer radius of the fiber and smoothing the interface between the core and the cladding on the propagation of the light in the optical fiber. These results have not been presented earlier to our knowledge. The reduction in computational time in our approach allows such investigations to be completed in a timely manner.

The paper is organized as follows: Section 2 reviews weakly guided modes in an optical fiber. In Section 3, we briefly review the Fourier transform-based Fresnel diffraction algorithm. Section 4 discusses the role of an absorbing mask in the outside of the grid and Section 5 shows results when we soften the sharp index boundary between the core and the cladding. In all cases, we show the effects of the propagation of the beam through the fiber. Section 6 presents the conclusions of this study.

2. Review of weakly guided modes in an optical fiber

The derivation of the waveguide modes of an optical fiber using Maxwell’s equations is well understood [16,17]. The difficulty in the exact approach involves matching the Bessel functions within the core with the modified Hankel functions for the cladding.

We use the notation in Yariv [16]. The optical fiber has a core of radius a with a higher index of refraction ${n_1}$. The cladding is considered to have an infinite diameter with a lower index of refraction ${n_2}$. The light has a vacuum wavelength of $\lambda$. The mode propagation is governed by a parameter V defined as

$$V = ak\sqrt {n_1^2 - n_2^2}.$$

For values of $V < 2.403$, only a single mode will propagate. As the value of V increases, more modes are allowed to propagate. These modes are usually shown as a plot of the normalized propagation constant as a function of V [16]. In each case, the propagation constant $\beta$ for each mode increases as V increases from a minimum value of ${n_2}k$ to a maximum of ${n_1}k$ as

$${n_1}k \ge \beta \ge {n_2}k.$$

Gloge [18] showed that the problem is much simplified in the weakly guided case where the difference between the core and cladding indices of refraction is small.

Next we discuss the Fresnel diffraction algorithm.

3. Brief review of the Fourier transform-based Fresnel diffraction algorithm

Light transmission in optical fibers can be examined using a number of computational techniques [19]. Often, the wave propagation is studied using the finite difference beam propagation method (FD-BPM) in [5,20], or the split-step Fourier method [21].

We have been using [2225] the Fast Fourier transform-beam propagation method (FFT-BPM). In this algorithm, the beam passes through a two-dimensional phase mask representing the effects of the fiber medium over a given propagation distance. Then this electric field is propagated over that distance using a Fresnel diffraction algorithm. This process is repeated over a large number of steps. Note that this phase mask can change with propagation to allow for such ideas as tapered waveguides.

We use the angular spectrum method [1] for computing the Fresnel diffraction. From our investigations, this appears to be the first reference for the Fourier transform-based approach. An explanation of this Fresnel diffraction algorithm in terms of the ray matrix formalism was provided in [26]. Finally a digital representation of this algorithm was provided in [3].

Next we briefly review this Fresnel diffraction algorithm using the notation introduced by Goodman [1]. We begin with a complex two-dimensional image representing the product of the electric field at the ${z_1}$ plane with the phase mask as $g({{x_1},{y_1},{z_1}} )$. The output at the ${z_2}$ plane is given according to the sequence in Eq. (3).

$$g({{x_2},{y_2},{z_2}} )= {\Im ^{ - 1}}\{{H({p,q,{z_2} - {z_1}} )\Im [{g({{x_1},{y_1},{z_1}} )} ]} \}.$$

Here $\Im$ and ${\Im ^{ - 1}}$ denote the Fourier and inverse Fourier transforms and $H({p,q,{z_2} - {z_1}} )$ is the angular spectrum transfer function. This angular transfer function can be written as

$$H(p,q,{z_2} - {z_1}) = \exp\left[ {\frac{{i\pi \lambda ({z_2} - {z_1})({p^2} + {q^2})}}{{{N^2}{\Delta ^4}}}} \right].$$

This angular transfer function can then be rewritten, in the paraxial approximation, as a lens function [26] as

$$H(p,q,{z_2} - {z_1}) = \exp\left[ {\frac{{i\pi ({p^2} + {q^2})}}{{\lambda F}}} \right].$$

Here the focal length F of the lens function is defined as

$$F = {\left( {\frac{{N{\Delta ^2}}}{\lambda }} \right)^2}\frac{1}{{({z_2} - {z_1})}}.$$

The array size is $NxN$ pixels, the pixel size is given by $\Delta $ and the wavelength in the medium is given by $\lambda$. Note that the focal length of this lens increases as the step distance decreases. We considered using the nonparaxial version of this function, but found that it did not produce any improvements.

So each step in the iteration in Eq. (3) consists of 4 steps. We begin with an initial electric field and transmit it through the phase mask. Then to perform the Fresnel diffraction, we form the complex-two dimensional Fourier transform of this input. Next we multiply it by the angular spectrum transform function. Finally, we form the inverse Fourier transform of this product.

This algorithm has several aspects that contribute to its speed. First, it uses two fast Fourier transform (FFT) operations. In addition, both the phase mask and the angular spectrum transfer function are only calculated once and then used in each subsequent operation as long as the waveguide geometry remains constant.

Next we discuss our model and how we overcome problems introduced by the FFT algorithm.

4. Importance of a boundary absorbing layer

For our studies, we use a grid of 1024 × 1024 pixels. The fiber has a core diameter of 8.2 microns formed on the center 64.3 pixels. Therefore each pixel has a size of 0.127 micron. The indices of refraction are 1.45 and 1.4432 for the core and cladding giving a numerical aperture of NA = 0.14 and corresponding to Corning SMF-028 single mode fiber. The vacuum wavelength was kept at 1.625 microns. These parameters give a value of V = 2.22 using Eq. (1) and we are in the single mode regime.

Our Macintosh computer program was written by Don Cottrell where various computer routines such as the Fourier transform or lens functions are represented by icons as shown in Fig. 1 and the flow is indicated by the black lines.

 figure: Fig. 1.

Fig. 1. Icon based program used in this computational study of optical fiber propagation.

Download Full Size | PDF

The patterns representing the phase mask and angular spectrum transform function are previously calculated and are stored in the dark blue icons. Each is copied into separate multiplying paste functions with the red arrows.

We then form the initial Gaussian input function. The initial beam is a Gaussian beam having a Gaussian width of 64 pixels and whose center is offset by 32 pixels. We use this as a way of forcing the fiber to make all of the allowed propagation modes.

Now we begin a loop of 50 iterations. In each iteration, we multiply the output of the previous iteration by the phase function using the red arrow to connect the two. The Fourier transform operation is performed followed by the multiplication of the angular transform function, again represented by the red arrow. The inverse Fourier transform is then formed.

After the 50 iterations, the output is saved using the light blue icon. We then enter the second loop that returns us to the 50-loop iteration. This second loop repeats 4096 times to give us the final output. We used steps of 1.12 micron covering a total distance of 23 cm for a total of 204,800 iterations. We saved 4096 images of the propagating beam at 0.056 mm steps. The total computational time including the saving of the images is roughly 1 hour.

Using the parameters of our system, the focal length of the lens in Eq. (6) is 193 microns and is shown in Fig. 2 where black represents a phase of $\pi$ radians and white represents a phase of $- \pi$ radians.

 figure: Fig. 2.

Fig. 2. Phase function for lens in Eq. (6) having a focal length of 193 micron with pixel sizes and wavelength for our computational system.

Download Full Size | PDF

In the application of this algorithm, the Fourier transform operation produces higher spatial frequencies that undergo aliasing. These higher frequency components move off of the computational grid and reappear returning from the opposite edge because the domain of the fast Fourier transform is toroidal. Note that the number of steps for aliases that go on the horizontal or vertical directions is different from those going in a diagonal direction.

Figure 3 shows 512 equally spaced one-dimensional electric field curves for the case where there is no boundary layer. The initial beam immediately decreases because the input beam is not equivalent to the single mode supported by the fiber. However as the beam propagates, we see noise fluctuations because of the aliasing.

 figure: Fig. 3.

Fig. 3. Amplitude of the electric field as a function of distance through the optical fiber where there is no absorbing layer at the boundary of the fiber.

Download Full Size | PDF

Next we added an absorbing layer around the edge of the fiber as shown in Fig. 4. Our absorbing layer has a width of about 20 pixels and has a modified hyperbolic tangent form.

 figure: Fig. 4.

Fig. 4. Amplitude of absorbing layer near the border of the 1024 × 1024 grid.

Download Full Size | PDF

Results showing the propagation similar to Fig. 3 are shown in Fig. 5. The initial beam, shown at the left, immediately decreases as before. The electric field profile quickly agrees with theory as shown earlier in Ref [3]. However the beam never reaches an equilibrium value and slowly decreases as a function of distance.

 figure: Fig. 5.

Fig. 5. Amplitude of the electric field as a function of distance through the optical fiber where we apply an absorbing layer. The beam retains the exact solution of the wave equation, but the amplitude slowly decays.

Download Full Size | PDF

The reason for this decrease is the sharp index change at the core- cladding edge. This sharp interface continues to produce scattering into the boundary layer, and draining energy from the center. Others have studied various alternative designs [27,28].

Next we describe how to reduce this effect.

5. Importance of smoothing the fiber core/cladding index boundary

To reduce the scattering caused by the sharp core-cladding index change, we convolve the core-clad interface with a Gaussian apodizing function. We start with a function having values of ${n_1}k$ inside and ${n_2}k$ outside. This was softened by convolving this function with a Gaussian function having a width of 0.4 propagation wavelengths (around 0.4 microns). The smoothing of the core-cladding phase interface is shown in Fig. 6. This gradient is very small compared to the 8.2-micron core diameter of the fiber.

 figure: Fig. 6.

Fig. 6. Smoothing of the phase step boundary between the core and cladding.

Download Full Size | PDF

This was converted into a phase function using the propagation distance. The final phase mask is shown in Fig. 7. The central core is graded and the absorbing layer is applied.

 figure: Fig. 7.

Fig. 7. Phase mask for the transmission through the propagation distance with the absorbing layer shown in black. Note that the phase representing the fiber core is graded.

Download Full Size | PDF

The effect of this smoothing is apparent in Fig. 8. As before, we see the sharp decline in the amplitude because the initial beam does not match the mode of the fiber. However, the amplitude of the beam quickly becomes stable and does not vary with propagation.

 figure: Fig. 8.

Fig. 8. The amplitude does not decay noticeably over the 23 cm length of fiber.

Download Full Size | PDF

Figure 9 compares the total energy in the beam as a function of propagation distance over a total of 23 cm. In both cases, the energy initially decreases rapidly because the input beam does not match the modes of the fiber. The subsequent drop in the unsoftened case is due to the introduction of short wavelength components produced by the sharp edges of the central core. These components are partly at too great an angle to be contained in the fiber system and make their way to the absorbing boundary. However, as suggested by Fig. 9, the energy remains constant as the beam propagates because of the softened boundary between the core and cladding.

 figure: Fig. 9.

Fig. 9. Energy of the beam as a function of distance through the optical fiber comparing the case where the core-cladding interface is a step function (in red) and where we apply a smoothing function to this interface (in black).

Download Full Size | PDF

Other groups have examined different approaches for softening this boundary. However, to our knowledge, these kinds of measurements have not been made.

6. Conclusions

In conclusion, we have continued our investigation of the angular spectrum method for examining the transmission of light in optical fibers. Extension to more complex configurations is very easy. In this work, we explain the evolution of this algorithm so that it can more easily implemented by others.

We continue to emphasize that this is an extremely simple algorithm where each repetitive segment consists of 4 components as discussed in Section 3. The output electric field from the previous iteration segment is multiplied by a phase mask shown in Fig. 7 that represents the phase difference between the light traveling through the core and cladding. Our major result is to show two ways in which this phase mask can be modified. First, we introduce an absorbing layer at the outer boundary of the cladding as shown in Fig. 4. Secondly we smooth the index of refraction between the core and the cladding as shown in Fig. 6. After this multiplication, we form the Fourier transform of this product. This Fourier transform is then multiplied by the angular transfer function of Eq. (6). However we find that this angular transfer function is equivalent to the multiplication by a lens function whose focal length is given by Eq. (7) and shown in Fig. 2. Finally the inverse Fourier transform of this product is obtained.

We note that the phase mask can be modified to allow for an input coupler or to allow more complicated cases where the fiber diameter is fusion-tapered or even given an elliptical format. The two dimensional FFT is then formed and multiplied by the angular spectrum transfer function. Finally the inverse FFT is formed. Then the segment is repeated.

We added some new features to the program and the total computational time now is roughly 60 minutes for the 1024 × 1024 array.

We find excellent agreement with the electric field results using the weakly guided model reported in [4]. However our results show the importance of the absorbing boundary layer and the softening of the index of refraction at the core-cladding boundary. To our knowledge, these effects have not been covered previously. However any investigation of these optical fiber characteristics would be severely limited by extensively long computational times. We are deeply disappointed in the lack of information regarding computational times in other works. We also believe that Ref [1]. shows the first work regarding these double Fourier algorithms and should be regarded as the main prior art.

The main conclusion of this study is that this angular spectrum approach using the Fresnel diffraction fast Fourier transform algorithm is extremely powerful and should be considered as a way of exploring new designs for optical fiber transmission and optical coupling. This could serve as a prelude to further explorations of these structures using more complex computational approaches.

We would welcome collaborations with other groups that wish to examine this approach.

Disclosures

The authors declare no conflicts of interest.

References

1. J. W. Goodman, Introduction to Fourier Optics, (McGraw Hill, 1969)

2. J. A. Davis and D. M. Cottrell, “Simulation of Anderson localization in a random fiber using a fast Fresnel diffraction algorithm,” Opt. Eng. 55(6), 066122 (2016). [CrossRef]  

3. D. M. Cottrell and J. A. Davis, “Simulation of optical fiber couplers using the angular spectrum algorithm,” Appl. Opt. 57(19), 5319–5327 (2018). [CrossRef]  

4. D. M. Cottrell and J. A. Davis, “Simulation of multimode optical fibers using the angular spectrum algorithm and a Fourier analysis,” Appl. Opt. 58(17), 4585–4591, (2019) [CrossRef]  

5. S. Karabasi, C. R. Mirr, R. J. Frazier, P. G. Yarandi, K. W. Koch, and A. Mafi, “Detailed investigation of the impact of the fiber design parameters on the transverse Anderson localization of light in disordered optical fibers,” Opt. Express 20(17), 18692–18706 (2012). [CrossRef]  

6. W. P. Huang and C. L. Xu, “Simulation of three-dimensional waveguides by full vector-beam propagation method,” J. Lightwave Technol. 29(10), 2639–2649 (1993). [CrossRef]  

7. R. Essiambre, G. Kramer, P. J. Winzer, G. J. Foschini, and B. Goebel, “Capacity limits of optical fiber networks,” J. Lightwave Technol. 28(4), 662–701 (2010). [CrossRef]  

8. J. Liang, Q. Mo, S. Fu, M. Tang, P. Shum, and D. Liu, “Design and fabrication of elliptical-core few-mode fiber for MIMO-less data transmission,” Opt. Lett. 41(13), 3058–3061 (2016). [CrossRef]  

9. G. Milione, E. Ip, M.-J. Li, J. Stone, G. Peng, and T. Wang, “Mode crosstalk matrix measurement of a 1 km elliptical core few-mode optical fiber,” Opt. Lett. 41(12), 2755–2758 (2016). [CrossRef]  

10. J. Richardson, J. Fini, M. Yao, and M. J. Padgett, “Orbital angular momentum: origins, behavior, and applications,” Adv. Opt. Photonics 3(2), 161–204 (2011). [CrossRef]  

11. G. Milione, M. P. J. Lavery, H. Huang, Y. Ren, G. Xie, T. A. Nguyen, E. Karimi, L. Marrucci, D. A. Nolan, R. R. Alfano, and A. E. Willner, “4×20 Gbit/s mode division multiplexing over free space using vector modes and a q-plate mode (de)multiplexer,” Opt. Lett. 40(9), 1980–1983 (2015). [CrossRef]  

12. Y. Wu and K. S. Chiang, “Ultra-broadband mode multiplexers based on three-dimensional asymmetric waveguide branches,” Opt. Lett. 42(3), 407–410 (2017). [CrossRef]  

13. T. Zhao, S. Lou, X. Wang, M. Zhou, and Z. Lian, “Ultrabroadband polarization splitter based on three-core photonic crystal fiber with a modulation core,” Appl. Opt. 55(23), 6428–6434 (2016). [CrossRef]  

14. S. G. Leon-Saval, N. K. Fontaine, J. R. Salazar-Gil, B. Ercan, R. Ryl, and R. J. Bland-Hawthon, “Mode-selective photonic lanterns for space-division multiplexing,” Opt. Express 22(1), 1036–1044 (2014). [CrossRef]  

15. L. Shen, L. Gan, L. Huo, C. Yang, W. Tong, S. Fu, M. Tang, and D. Liu, “Design of highly mode group selective photonic lanterns with geometric optimization,” Appl. Opt. 57(24), 7065–7069 (2018). [CrossRef]  

16. A. Yariv, Optical Electronics in Modern Communications, 5th ed. (Oxford University, 1997)

17. B.E.A. Saleh and M. C. Teich, Fundamentals of Photonics, 2nd Edition, (Wiley & Sons, 2007)

18. D. Gloge, “Weakly guiding fibers,” Appl. Opt. 10(10), 2252–2258 (1971). [CrossRef]  

19. J. Bures, Guided Waves, (Wiley-VCH, 2009)

20. W. P. Huang and C. L. Xu, “Simulation of three-dimensional waveguides by full vector-beam propagation method,” J. Lightwave Technol. 29(10), 2639–2649 (1993). [CrossRef]  

21. O. V. Sinkin, R. Holzlohner, J. Zweck, and C. R. Menyuk, “Optimization of the split-step Fourier method in modeling optical-fiber-communications systems,” J. Lightwave Technol. 21(1), 61–68 (2003). [CrossRef]  

22. J. A. Fleck, J. R. Morris, and M. D. Feit, “Time-dependent propagation of high-energy laser beams through the atmosphere,” Appl. Phys. 10(2), 129–160 (1976). [CrossRef]  

23. M. D. Feit and J. A. Fleck, “Light propagation in graded-index optical fibers,” Appl. Opt. 17(24), 3990–3998 (1978). [CrossRef]  

24. C. Yeh, W. P. Brown, and R. Szejn, “Multimode inhomogeneous fiber couplers,” Appl. Opt. 18(4), 489–495 (1979). [CrossRef]  

25. J. Van Roey, J. van der Donk, and P. E. Lagasse, “Beam propagation method; analysis and assessment,” J. Opt. Soc. Am. 71(7), 803–810 (1981). [CrossRef]  

26. J. A. Davis and D. M. Cottrell, “Ray matrix analysis of the fast Fresnel transform with applications towards liquid crystal displays,” Appl. Opt. 51(5), 644–650 (2012). [CrossRef]  

27. D. Yevick and J. Yu, “Optimal absorbing boundary conditions,” J. Opt. Soc. Am. A 12(1), 107–110 (1995). [CrossRef]  

28. R. Learn and E. Feigenbaum, “Adaptive step-size algorithm for Fourier beam-propagation method with absorbing boundary layer of auto-determined width,” Appl. Opt. 55(16), 4402–4407 (2016). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (9)

Fig. 1.
Fig. 1. Icon based program used in this computational study of optical fiber propagation.
Fig. 2.
Fig. 2. Phase function for lens in Eq. (6) having a focal length of 193 micron with pixel sizes and wavelength for our computational system.
Fig. 3.
Fig. 3. Amplitude of the electric field as a function of distance through the optical fiber where there is no absorbing layer at the boundary of the fiber.
Fig. 4.
Fig. 4. Amplitude of absorbing layer near the border of the 1024 × 1024 grid.
Fig. 5.
Fig. 5. Amplitude of the electric field as a function of distance through the optical fiber where we apply an absorbing layer. The beam retains the exact solution of the wave equation, but the amplitude slowly decays.
Fig. 6.
Fig. 6. Smoothing of the phase step boundary between the core and cladding.
Fig. 7.
Fig. 7. Phase mask for the transmission through the propagation distance with the absorbing layer shown in black. Note that the phase representing the fiber core is graded.
Fig. 8.
Fig. 8. The amplitude does not decay noticeably over the 23 cm length of fiber.
Fig. 9.
Fig. 9. Energy of the beam as a function of distance through the optical fiber comparing the case where the core-cladding interface is a step function (in red) and where we apply a smoothing function to this interface (in black).

Equations (6)

Equations on this page are rendered with MathJax. Learn more.

V = a k n 1 2 n 2 2 .
n 1 k β n 2 k .
g ( x 2 , y 2 , z 2 ) = 1 { H ( p , q , z 2 z 1 ) [ g ( x 1 , y 1 , z 1 ) ] } .
H ( p , q , z 2 z 1 ) = exp [ i π λ ( z 2 z 1 ) ( p 2 + q 2 ) N 2 Δ 4 ] .
H ( p , q , z 2 z 1 ) = exp [ i π ( p 2 + q 2 ) λ F ] .
F = ( N Δ 2 λ ) 2 1 ( z 2 z 1 ) .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.