The method of lines (MoL) has been developed to study coupling efficiency on hemispherical lens. In this paper, the physical shape of the lens is approximated by cascading a number of straight waveguide segments. The perfectly matched layer (PML) is applied as an absorber for the MoL to reduce numerical reflection in the simulation region. Analysis is done by calculating coupling efficiency at the plane of integration where the coupling efficiency is an overlap integral between laser diode field and fiber field. The result of coupling efficiency in this analysis is compared to the experiment and ABCD matrix. It is found that MoL gives good result accuracy.
©2009 Optical Society of America
Photonic device is a fundamental component in an optical communication network. Coupling of the laser between the fiber and the device plays an important role in achieving low insertion loss. A medium such as lens is often deployed in the device especially in a laser diode to obtain high coupling efficiency. To date, variety of lens designs have been produced which mostly focused on the fiber lens.
Design and analysis stage are the most important process in producing a workable fiber lens. It requires an accurate and reliable tool that can predict the optimum parameter of the fiber lens in order to achieve high coupling efficiency. At the moment, paraxial ray approximation technique based on ABCD matrix is the most widely used method in fiber lens analysis [1, 2]. This technique involves transformation of complex beam parameter by using A, B, C and D elements as a form of matrix which describes ray parameter of input and output planes . Its wide application is related to its simplicity and accuracy.
Apart from the paraxial ray approximation technique, method such as finite difference beam propagation method (FD-BPM), finite element method (FEM) and method of lines (MoL) are becoming popular due to an increasing computational power. The methods mentioned above use various discretization schemes, which some of them have an advantage over the other, depending on the application. The FD-BPM uses discretization scheme in the form of rectangular grid. Through finite difference approximation, the wave equation is discretized and propagated along the longitudinal axis by using Pade’ approximants operator or paraxial solution. The advantages of FD-BPM include a straightforward implementation of algorithm in various photonic structure and the ability to study polarization and nonlinearity . To date, this technique has been used in the simulation of lensed coreless fiber , tapered lens  and wedge-shaped lens . The second technique, the FEM has descretization scheme in the form of triangular shapes. Using variational formulation, each triangular element is analyzed individually in order to obtain full solution of the simulation region . Its flexible triangular shape yields more accurate discretization of the complex simulation structure. The successful implementation of the FEM can be found in the spot-size conversion application [9, 10]. The third technique, the MoL treats simulation regions as discrete lines in longitudinal direction. It solves wave equation by obtaining eigenvalue and eigenvector in each individual line analytically. Due to the nature of solving the wave equation, the MoL also can be considered as semianalytical . Various structures have been simulated by the MoL which include erbium doped waveguide amplifier , beam splitter  and multi-layer anti-resonant reflecting optical waveguide .
In this paper, we analyze the coupling efficiency of the hemispherical lens using the MoL technique. The MoL is chosen as a simulation tool due to its straightforward discretization and high accuracy calculation. Other than that, the backward propagation field can be simulated without compromising its computational stability. Furthermore, to our knowledge, this is the first time the MoL is adopted in calculating coupling efficiency on the hemispherical lens. For an absorber, our MoL is incorporated with the perfectly matched layer (PML) to suppress reflection at the edge of the simulation region. To verify its accuracy, the results produced by MoL are benchmarked with the results from the experiment and the ABCD matrix.
The paper is organized as follow; the basic of the MoL is reviewed in section 2. The implementation of PML in the simulation region is described in section 3. In section 4, the discretization scheme used in the lens is elaborated. Finally, results from the simulation are presented and analyzed in section 5.
2. Introduction to MoL
In this section, theory of MoL is explained using simple homogenous structure; a slab waveguide. Consider Helmholtz wave equation:
where k 0 is a wavenumber in freespace condition, n(x) is a refractive index distribution in transverse direction and E(x,z) is a transverse electric field.
In Fig. 1, the simulation region of the slab waveguide is discretized transversely into m number of lines where each lines represent the discrete E, n and x coordinate. The width between the two adjacent lines is represented by the stepsize, Δx. Using the central difference approximation, the second order derivative operator in Eq. (1) is discretized, so that the new discretized wave equation is:
where E is the column vector m × 1:
3. Perfectly matched layer
Scattering radiation always causes problem to MoL calculation. This radiation approaches the boundary of the simulation region and reflects back into the computed region. This phenomenon causes unwanted interference inside the simulation region. It will interfere the calculated results and might mislead the overall conclusion. To overcome this problem, few methods have been introduced in a form of artificial boundary to achieve minimum reflection. These boundary conditions are introduced where some absorb the incoming wave and act like a sponge . There is also method that behaves like a transparent boundary wall . Both boundary conditions mentioned above have disadvantages at which  fails to accommodate zero reflection at all angles  and  gives poor performance for highly diverging beam .
Our simulation involves propagation of laser field through the air medium. The laser field which is initially emitted from the laser diode is highly diverged. Without a proper absorbing boundary condition, a severe reflection will be created inside the simulation region. For this reason, the multiple layers of PML proposed by Jamid  is chosen as an absorber in our MoL due to its remarkable performance in minimizing reflection from highly diverging field.
The detail of PML structure in the simulation region is illustrated in Fig. 2. The MoL simulation region is split into two regions; real space and complex space. In the real space, the step size, h 0 and x 1,2…m-p coordinates are real number, hence represent the actual simulation region. In the case of complex space, this region consists of artificial layer where its step size, hi and x coordinate, x¯i can be defined by :
where ε is the PML strength, h 0 is the stepsize Δx in the real space, p is total number of PML layer, f(x¯i) is the loss profile function and i = 1,2....p. The function of the complex space is to handle the suppression of reflection by attenuating the amount of radiation field inside this region. To integrate the PML with MoL, the new Q 2 matrix is defined by changing the first and the last five rows of the first part of Eq. (4):
with ai 3b (b = 1,2...5) is defined by :
4. Fiber Lens Discretization
The discretization process used in section 2 is intended for homogeneous waveguide where its refractive index distribution is uniform along z axis. For the case of non-homogeneous waveguide, such as hemispherical lens, discretization at other direction is needed to approximate the varying refractive index distribution. To serve this requirement, the refractive index distribution is discretized into n number of segments along the z axis. Each discretized segment represents a straight waveguide structure, si whereby each segment is discretized transversely by following the scheme as explained in section 2. The step size in z axis for each segment is denoted by Δz. For better understanding, the discretization process performed on the hemispherical lens is illustrated in Fig. 3. Once the discretization is done, the field from the laser diode can be propagated through the lens by using the two equations  below:
where I is an identity matrix, Di = exp(UiΔz) and Ui = Qi. By referring to Fig. 4, the term Ai and Bi represent the transmission field and reflection field where both of these fields are related by this expression Bi = Γi DiAi. Both Eq. (11) and Eq. (12) satisfy the continuity of E and at the interface of each ith segment. Eq. (11) is solved recursively starting from the layer i = n till the layer i = 0 where it is assumed Γn+1 = 0. Once Γi are solved for each segments, the incident field, A 0 is propagated segment by segment by using Eq. (12).
The result of the calculation can be improved further if smaller Δx or Δz is chosen. The reason for this improvement can be explained by analyzing the relationship between Δx or Δz and the number of discretization lines. The smaller Δx or Δz value, the higher number of discretization line in x or z axis. The increased number of discretization lines leads to better approximation of the physical shape of hemispherical lens, thus increases simulation accuracy. However, this comes with long simulation time due to high demand for computational power and memory requirement.
5. Simulation result
The launched input field from the laser diode, E l0 is taken from the Gaussian field in Eq. (13) which has initial spotsize, ωox and ωoy.
Refer to Fig. 5 for the detail structure of hemispherical lens. In our entire analysis, the same simulation parameter is used consistently in MoL and PML algorithm which is listed in Table 1. The parameter of laser diode (LD) and single-mode fiber (SMF) are summarized in the Table 2.
The coupling efficiency, η between the LD and the lens is evaluated using overlap integral  defined as,
where E f0 is a fundamental mode of SMF and El is a transformed laser diode field at the interface of the lens and the SMF. Since the MoL simulation is done on the two separate axes, both fields in Eq. (14) can be decomposed into E f0(x,y) = E f0(x).E f0(y) and El(x,y)=El(x).El(y). Therefore, expression in Eq. (14) can be reduced to a much simpler form:
where ηx and ηy are single integral on the separate axes, x and y.
The first MoL simulation is done by varying the working distance, d between LD and lens. The coupling efficiency is then evaluated at the integral plane P1 for each iteration of d and the result is plotted in Fig. 6. For the comparison of our analysis, MoL result is compared with the experiment  and ABCD matrix analysis which includes the effect of Fresnel reflection and spherical aberration . In this graph, the maximum coupling efficiency predicted by MoL is 66% while the result from experiment produces 70%. Both maximum coupling efficiencies are recorded at an optimum working distance of 7 μm. By comparing both results, the experiment has higher coupling efficiency than MoL due to the effect of multiple reflection between the laser diode facet and the lens which is responsible to the increase of coupling efficiency. It should be noted that our MoL algorithm does not include the effect of multiple reflection. The other result produced by ABCD matrix has the same trend of coupling efficiency with MoL except the optimum working distance is shifted to 5.5 μm.
The factor that influences the coupling efficiency can be explained by considering two important parameters of El field which are spotsize, ωx,y and radius of curvature of the wave fronts, Rx,y. To analyze the contribution of the two parameters towards coupling efficiency, consider El field located at the position z=0 which has initial spotsize, ωox,oy and radius of curvature, Rox,oy = ∞. When the El field propagates through the air along the z axis, the parameter Rx,y and ωx,y are governed by:
Once the El field reaches the lens, Rx,y and ωx,y will be modified by the lens to match the Rof and ωof of E f0 field at plane P1. By matching together the two parameters of El field and E f0 field as close as possible, maximum coupling efficiency can be obtained. In MoL simulation, ωx,y and Rx,y  are obtained from :
Referring to Fig. 7, analysis is done in x axis to verify the influence of ωx and Rx towards coupling efficiency. By varying the working distance from d = 0 to d = 30 μm, it is found that Elx field has the highest value of Rx, approximately 200 μm at d = 7 μm. For the spotsize, its value increases gradually with d, thus the ωx at d = 7 μm is not the highest value. Therefore, Rx has greater influence on the coupling efficiency if compared to ωx. This can be explained by considering that Rx has a role in flattening the wave fronts of Elx field. The process of flattening the wave fronts is explained pictorially in Fig. 8(a). In this figure, we can see the wave fronts of Elx field in the air region initially has spherical form. However, after Elx field propagates through the lens, the spherical shape of the wave fronts changes to a nearly flat shape. This nearly flat shape corresponds to the large value of Rx as shown in Fig. 7. If the wave fronts of Elx field has the shape other than flat, Elx field will be out of phase with E f0 field. As a consequence of this, most of Elx field will not be coupled to E f0 field. For the spotsize conversion, Fig. 8(b) verifies the role of the lens in modifying ωx. By comparing Elx field at the plane P2 and P1, it clearly shows that ωx is reduced by the lens in order to have better match with ωof. The dicussion on the radius of curvature and the spotsize above is done on a single x-axis only. For ωy and Ry, the result differs little from ωx and Rx due to the nearly symmetrical shape of El field in both axes.
To prove MoL consistency on analyzing different parameter values of the lens, the radius of the hemispherical is varied from r = 3.5 to 6.0 μm. The coupling efficiency is derived from each radius at the optimum working distance and compared with the experimental result  and ABCD matrix. In Fig. 9, it is found that the result generated by MoL has good agreement with the experiment, although some variation occurs between them. The result deviation is contributed by the effect of multiple reflection between the laser diode and the lens which is not included in our MoL. Due to the multiple reflection, the rapid variation in the power for every half of λ change in working distance contributes difficulty in determining the real coupling efficiency value in the experiment. In the case of ABCD matrix, the method produces same pattern of coupling efficiency with experiment, but less accurate if compared to MoL. Despite the result variation, it is noticed that MoL simulation with smaller radius r = 4.0 to 3.5 μm has better accuracy than the large one. The improved accuracy at smaller radius is probably due to the reduced effect of multiple reflection between laser diode facet and hemispherical lens.
To further verify our MoL algorithm, the effective reflectivity, Re of the lens is simulated and analyzed. This parameter is derived based on a single axis  and is extended to both axes, x and y;
The calculation of Re is done at plane P2 where E ref and E inc are the reflection and incident laser fields at position z = d μm. In Fig. 10, the simulation result for all the radii show the decrease of Re along the working distance. This trend can be explained by analyzing laser propagation characteristic in the air medium. At the laser diode facet z = 0 μm, the emitted laser field, El starts to diverge during its propagation along the working distance. The divergence of the El field cause its field intensity to spread across the transverse dimension. When the El field or now known as E inc field reaches plane P2, the surface of the lens receives reduced intensity of the E inc field. The lower intensity means less E inc field is reflected from the lens surface. By increasing the working distance, the intensity of E ref field can be reduced accordingly, which leads to lower Re. This phenomenal has been similarly reported by other researcher  which confirms our result validity. The result also reveals another trend which shows the reduction of Re with the smaller lens radius. Smaller lens radius decreases the lens surfaces area, which in turn reduces the reflection by the lens. Our result on Re dependence on the lens radius is verified by paper  on tapered hemispherical lens.
The method of lines (MoL) algorithm which is integrated with the perfectly matched layer (PML) is constructed to analyze the hemispherical lens. Using overlap integral technique, the coupling efficiency is calculated in the MoL to evaluate the performance of the hemispherical lens. Further analysis by the MoL shows the effect of spotsize and radius of curvature of wave fronts on the coupling efficiency. Apart from studying coupling efficiency, the MoL simulation is extended to study backward reflection field by matching E and at waveguide discontinuities. By varying the gap distance and radius of the lens, the MoL is found to be reliable in producing accurate result. Due to the flexibility in discretization scheme, the MoL can be a suitable technique to study various fiber lens structures. In future, our MoL algorithm can be adapted to simulate more complex lens structures such as conical and tapered hyperbolic.
References and links
1. S. Gangopadhyay and S. Sarkar, “ABCD matrix for reflection and refraction of Gaussian light beams at surfaces of hyperboloid of revolution and efficiency computation for laser diode to single-mode fiber coupling by way of a hyperbolic lens on the fiber tip,” Appl. Opt. 36, 8582–8586 (1997). [CrossRef]
2. K. Sambanthan and F. A. Rahman, “Method to improve the coupling efficiency of a hemispherically lensed asymmetric tapered-core fiber,” Opt. Commun. 254, 112–118 (2005). [CrossRef]
3. W. L. Emkey and C. A. Jack, “Analysis and evaluation of graded-index fiber-lenses,” J. Lightwave Technol. 12, 1156–1164 (1987). [CrossRef]
4. R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical techniques for modeling guided-wave photonic devices,” IEEE J. Sel. Topics Quantum Electron. 6, 150–162 (2000). [CrossRef]
5. J. Yamauchi, K. Nishio, and H. Nakano, “Analysis of a lensed coreless fiber by a hybrid technique combining FD-BPM and FD-TDM,” J. Lightwave Technol. 16, 465–471 (1998). [CrossRef]
6. Z. Wang, B. Mikkelsen, B. Pedersen, K. E. Stubkjaer, and D. S. Olesen, “Coupling between angled-facet amplifiers and tapered lens-ended fibers,” J. Lightwave Technol. 9, 49–55 (1991). [CrossRef]
7. Y. He, S. K. Mondal, and F. G. Shi, “Design optimization of wedge-shaped lensed fibers for fiber-laser coupling: Fresnel reflection and non-Gaussian mode effects,” J. Lightwave Technol. 21, 2271–2275 (2003). [CrossRef]
8. K. Okamakoto, Fundamental of Optical Waveguides (Academic Press, California, 2000).
9. T. Wongcharoen, B. M. A. Rahman, M. Rajarajan, and K. T. V. Grattan, “Spot-size conversion using uniform waveguide sections for efficient laser-fiber coupling,” J. Lightwave Technol. 19, 708–716 (2001). [CrossRef]
10. M. Rajarajan, B. M. A. Rahman, and K. T. V. Grattan, “Numerical study of spot-size expanders for an efficient OEIC to SMF coupling,” IEEE Photonics Technol. Lett. 10, 1082–1084 (1998). [CrossRef]
11. M. N. O. Sadiku and C. N. Obiozor, “A simple introduction to the method of lines,” Int. J. Electr. Eng. Educ. 282–296 (2000).
12. W. Huang and R. R. A. Syms, “Analysis of folded erbium-doped planar waveguide amplifiers by the method of lines,” J. Lightwave Technol. 17, 2658–2664 (1999). [CrossRef]
13. H. A. Jamid, M. Z. M. Khan, and M. Ameeruddin, “A compact 90° three-branch beam splitter based on resonant coupling,” J. Lightwave Technol. 23, 3900–3906 (2005). [CrossRef]
14. A. Abdullah and M. A. Majid, “Analysis of multi-layer ARROW,” J. Microwaves and Optoelectronics 3, 1–8 (2003).
15. A. Dreher and R. Pregla, “Analysis of planar waveguides with the method of lines and absorbing boundary conditions,” IEEE Microwave Guided Wave Lett. 1, 138–140 (1991). [CrossRef]
16. B. Engquist and A. Majda, “Absorbing boundary conditions for the numerical simulation of waves,” Math. Comp. 31, 629–651 (1977). [CrossRef]
17. G. R. Hadley, “Transparent boundary condition for the beam propagation method,” IEEE J. Quantum Electron. 28, 363–370 (1992). [CrossRef]
18. D. Weiping and Z. Linchang, “An improvement algorithm of Mur’s First-Order absorbing boundary condition,” in IEEE 1997 International Symposium on Electromagnetic Compatibility, (Austin,USA1997), pp. 592–595.
19. C. Vassallo and F. Collino, “Highly efficient absorbing boundary conditions for the beam propagation method,” J. Lightwave Technol. 14, 1570–1577 (1996). [CrossRef]
20. H. A. Jamid, “Enhanced PML performance using higher order approximation,” IEEE Trans. Microwave Theory Tech. 52, 1166–1174 (2004). [CrossRef]
21. M. Z. M. Khan, “Analysis of one and two dimensional bandgap structures using automated method of lines with arbitrary longitudinal discontinuities,” Master dissertation, (King Fahd University of Petroleum and Minerals, Saudi Arabia, 2004).
22. A. A. Shittu, “Study of periodic waveguides by the finite-difference time domain method and the method of lines,” PhD dissertation, (King Fahd University of Petroleum and Minerals, Saudi Arabia, 1994).
23. J. John, T. S. M. Maclean, H. Ghafouri-Shiraz, and J. Niblett, “Matching of single-mode fibre to laser diode by microlenses at 1.5 μm wavelength,” IEE Proc.-Optoelectron. 141, 178–184 (1994). [CrossRef]
24. W. T. Chen and L. A. Wang, “Out-of-plane optical coupling between an elliptical Gaussian beam and an angled single-mode fiber,” J. Lightwave Technol. 16, 1589–1595 (1998). [CrossRef]
25. F. A. Rahman, K. Takahashi, and C. H. Teik, “A scheme to improve the coupling efficiency and working distance between laser diode and single mode fiber,” Opt. Commun. 208, 103–110 (2002). [CrossRef]
26. J. Alda, “Laser and Gaussian beam propagation and transformation,” in Encyclopedia of Optical Engineering, R. G. Driggers (Marcel Dekker, New York, 2003), pp. 999–1013.
27. T. Saitoh, T. Mukai, and O. Mikami, “Theoretical analysis and fabrication of antireflection coatings on laser-diode facets,” J. Lightwave Technol. LT-3, 288–293 (1985). [CrossRef]
29. H. Kuwahara, Y. Onoda, M. Goto, and T. Nakagami, “Reflected light in the coupling of semiconductor lasers with tapered hemispherical end fibers,” Appl. Opt. 22, 2732–2738 (1983). [CrossRef] [PubMed]