Fully vectorial laser resonator modeling of continuous-wave solid-state lasers including rate equations , thermal lensing and stress-induced birefringence

The computer-aided design of high quality mono-mode, continuous-wave solid-state lasers requires fast, flexible and accurate simulation algorithms. Therefore in this work a model for the calculation of the transversal dominant mode structure is introduced. It is based on the generalization of the scalar Fox and Li algorithm to a fully-vectorial light representation. To provide a flexible modeling concept of different resonator geometries containing various optical elements, rigorous and approximative solutions of Maxwell’s equations are combined in different subdomains of the resonator. This approach allows the simulation of plenty of different passive intracavity components as well as active media. For the numerically efficient simulation of nonlinear gain, thermal lensing and stress-induced birefringence effects in solid-state active crystals a semi-analytical vectorial beam propagation method is discussed in detail. As a numerical example the beam quality and output power of a flash-lamp-pumped Nd:YAG laser are improved. To that end we compensate the influence of stress-induced birefringence and thermal lensing by an aspherical mirror and a 90◦ quartz polarization rotator. © 2015 Optical Society of America OCIS codes: (140.3410) Laser resonators; (140.3570) Lasers, single-mode; (140.4780) Optical resonators; (260.1440) Birefringence; (350.6830) Thermal lensing. References and links 1. J. R. Leger, D. Chen, and Z. Wang, “Diffractive optical element for mode shaping of a Nd:YAG laser,” Opt. Lett. 19(2), 108–110 (1994). 2. O. Svelto, Principles of Lasers (Springer, 2010). 3. W. W. Rigrod, “Saturation effects in high gain lasers,” Appl. Phys. 36(8), 2487–2490 (1965). 4. A. G. Fox and T. Li, “Effects of gain saturation on the oscillating modes of optical masers,” IEEE J. Quantum Electron. 2(12), 774–783 (1966). 5. W. Koechner, “Thermal lensing in a Nd:YAG laser rod,” Appl. Opt. 9(11), 2548–2553 (1970). 6. J. D. Foster and L. M. Osterink, “Thermal effects in a Nd: YAG laser,” Appl. Phys. 41(9), 3656–3663 (1970). 7. T. Taira, A. Mukai, Y. Nozawa, and T. Kobayashi, “Single-mode oscillation of laser-diode-pumped Nd:YVO4 microchip lasers,” Opt. Lett. 16(24), 1955–1957 (1991). 8. D. S. Kliger and J. W. Lewis, Polarized Light in Optics and Spectroscopy (Elsevier, 1990). 9. Y.-Z. Huang, W.-H. Guo, and Q.-M. Wang,“Analysis and numerical simulation of eigenmode characteristics for semiconductor lasers with an equilateral triangle micro-resonator,” IEEE J. Quantum Electron. 37(1), 100–107 (2001). #240663 Received 11 May 2015; revised 2 Jul 2015; accepted 3 Jul 2015; published 13 Jul 2015 © 2015 OSA 27 Jul 2015 | Vol. 23, No. 15 | DOI:10.1364/OE.23.018802 | OPTICS EXPRESS 18802 10. P. Nyakas, “Full-vectorial three-dimensional finite element optical simulation of vertical-cavity surface-emitting lasers,” J. Lightw. Technol. 25(9), 2427–2434 (2007). 11. A. Christ, N. Kuster, M. Streiff, A. Witzig, and W. Fichtner, “Correction of the numerical reflection coefficient of the finite-difference time-domain method for efficient simulation of vertical-cavity surface-emitting lasers,” J. Opt. Soc. Am. B 20(7), 1401–1408 (2003). 12. M. Streiff, A. Witzig, M. Pfeiffer, P. Royo, and W. Fichtner, “A comprehensive VCSEL device simulator,” IEEE J. Sel. Topics Quantum Electron. 9(3), 879–891 (2003). 13. J. Pomplun, S. Burger, F. Schmidt, A. Schliwa, D. Bimberg, A. Pietrzak, H. Wenzel, and G. Erbert, “Finite element simulation of the optical modes of semiconductor lasers,” Phys. Status Solidi B 247(4), 846–853 (2010). 14. H. Kogelnik and T. Li, “Laser beams and resonators,” Appl. Opt. 5(10), 1550–1567 (1966). 15. A. E. Siegman, Lasers (University Science Books, 1986). 16. C. N. Kurtz and W. Streifer, “Guided waves in inhomogeneous focusing media, part I: formulation, solution for quadratic inhomogeneity,” IEEE Trans. Microw. Theory Techn. 17(1), 11–15 (1969). 17. W. J. Firth, “Propagation of laser beams through inhomogeneous media,” Opt. Commun. 22(2), 226–230 (1977). 18. H. Kogelnik, “On the propagation of Gaussian beams of light through lenslike media including those with a loss or gain,” Appl. Opt. 4(12), 1562–1569 (1965). 19. T. Y. Fan and R. L. Byer, “Diode laser-pumped solid-state lasers,” IEEE J. Quantum Electron. 24(6), 895–912 (1988). 20. D. G. Hall, R. J. Smith., and R. R. Rice. “Pump-size effects in Nd:YAG lasers,” Appl. Opt. 19(18), 3041–3043 (1980). 21. L. W. Casperson, “Laser power calculations: sources of error,” Appl. Opt. 19(3), 422–431 (1980). 22. P. F. Moulton, “An investigation of the Co:MgF2 laser system,” IEEE J. Quantum Electron. 21(10), 1582–1588 (1985). 23. J. Junghans, M. Keller, and H. Weber, “Laser resonators with polarizing elements eigenstates and eigenvalues of polarization,” Appl. Opt. 13(12), 2793–2798 (1974). 24. T. Graupeter, R. Hartmann, and C. Pflaum, “Calculations of eigenpolarization in Nd:YAG laser rods due to thermally induced birefringence,” IEEE J. Quantum Electron. 50(12), 1035–1043 (2014). 25. A. G. Fox and T. Li, “Resonant modes in a maser interferometer,” Bell Syst. Tech. J. 40(2), 453–488 (1961). 26. A. G. Fox and T. Li, “Modes in a maser interferometer with curved and tilted mirrors,” Proc. IEEE 51(1), 80–89 (1963). 27. D. Asoubar, S. Zhang, F. Wyrowski, and M. Kuhn, “Laser resonator modeling by field tracing: a flexible approach for fully vectorial transversal eigenmode calculation,” J. Opt. Soc. Am. B 31(11), 2565–2573 (2014). 28. F. Wyrowski and M. Kuhn, “Introduction to field tracing,” J. Mod. Opt. 58(5-6), 449–466 (2011). 29. L. Thylen and D. Yevick, “Beam propagation method in anisotropic media,” Appl. Opt. 21(15), 2751–2754 (1982). 30. J. A. Fleck and M. D. Feit, “Beam propagation in uniaxial anisotropic media,” J. Opt. Soc. Am. 73(7), 920–926 (1983). 31. L. Xu, P. Huang, J. Chrostowski, and S. K. Chaudhuri, “Full-vectorial beam propagation method for anisotropic waveguides,” J. Lightw. Technol. 12(11), 1926–1931 (1994). 32. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 1995). 33. J. Saijonmaa and D. Yevick, “Beam-propagation analysis of loss in bent optical waveguides and fibers,” J. Opt. Soc. Am. 73(12), 1785–1791 (1983). 34. A. E. Siegman and E. A. Sziklas, “Mode calculations in unstable resonators with flowing saturable gain. 1: Hermite-Gaussian expansion,” Appl. Opt. 13(12), 2775–2792 (1974). 35. E. A. Sziklas and A. E. Siegman, “Mode calculations in unstable resonators with flowing saturable gain. 2: fast Fourier transform method,” Appl. Opt. 14(8), 1874–1889 (1975). 36. B. A. E. Saleh and M. C. Teich, Fundamentals of Photonics (Wiley Interscience, 2007). 37. Y. Sato and T. Taira, “Saturation factors of pump absorption in solid-state lasers,” IEEE J. Quantum Electron. 40(3), 270–280 (2004). 38. Meta Numerics, Library for advanced scientific computation in the .NET Framework, www.meta-numerics.net, accessed (2/18/2015). 39. R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey and D. E. Knuth, “On the LambertW function,” Adv. Comput. Math. 5(1), 329–359 (1996). 40. A. Hoorfar and M. Hassani, “Inequalities on the Lambert W function and hyperpower function,” J. Inequal. Pure and Appl. Math 9(2), 5–9 (2008). 41. N. Hodgson and H. Weber, Optical Resonators: Fundamentals, Advanced Concepts, Applications (Springer Science, 2005). 42. C. Pfistner, R. Weber, H. P. Weber, S . Merazzi, and R. Gruber, “Thermal beam distortions in end-pumped Nd: YAG, Nd: GSGG, and Nd : YLF rods,” IEEE J. Quantum Electron. 30(7), 1605–1615 (1994). 43. U. O. Farrukh, A. M. Buoncristiani, and C. E. Byvik, “An analysis of the temperature distribution in finite solidstate laser rods,” IEEE J. Quantum Electron. 24(11), 2253–2263 (1988). 44. S. C. Tidwell, J. F. Seamans, M. S. Bowers, and A. K. Cousins, “Scaling cw diode-end-pumped Nd: YAG lasers #240663 Received 11 May 2015; revised 2 Jul 2015; accepted 3 Jul 2015; published 13 Jul 2015 © 2015 OSA 27 Jul 2015 | Vol. 23, No. 15 | DOI:10.1364/OE.23.018802 | OPTICS EXPRESS 18803 to high average powers,” IEEE J. Quantum Electron. 28(4), 997–1009 (1992). 45. W. Koechner, Solid-State Laser Engineering (Springer-Verlag, 1996). 46. M. Schmid, T. Graf, and H. P. Weber, “Analytical model of the temperature distribution and the thermally induced birefringence in laser rods with cylindrically symmetric heating,” J. Opt. Soc. Am. B 17(8), 1398–1404 (2000). 47. S. Chenais, S. Forget, F. Druon, F. Balembois, and P. Georges, ”Direct and absolute temperature mapping and heat transfer measurements in diode-end-pumped Yb:YAG,” Appl. Phys. B 79(2), 221–224 (2004). 48. J. Didierjean, S. Forget, S. Chenais, F. Druon, F. Balembois, P. Georges, K. Altmann, and C. Pflaum, “Highresolution absolute temperature mapping of laser crystals in diode-end-pumped configuration,” Proc. SPIE 5707, 370-379 (2005). 49. J. D. Foster and L. M. Osterink, “Index of refraction and expansion thermal coefficients of Nd:YAG,” Appl. Opt. 7(12), 2428–2429 (1968). 50. J. L. Blows, J. M. Dawes, and T. Omatsu, “Thermal lensing measurements in line-focus end-pumped neodymium yttrium aluminium garnet using holographic lateral shearing interferometry,” Appl. Phys. 83(6), 2901–2906 (1998). 51. W. Koechner and D. K. Rice, “Effect of birefringence on the performance of linearly polarized YAG:Nd lasers,” IEEE J. Quantum Electron. 6(9), 557–566 (1970). 52. R. Weber, B. Neuenschwander, and H. P. Weber, “Thermal effects in solid-state laser materials,” Opt. Mater. 11(2), 245–254 (1999). 53. M. Schmid, R. Weber, Th. Graf, M. Roos, and H. P. Weber, “Numerical simulation and analytical description of the thermally induced birefringence in laser rods,” IEEE J. Quantum Electron. 36(5), 620–626 (2000). 54. R. W. Dixon, “Photoelastic properties of selected materials and their relevance for applications to acoustic light modulators and scanners,” Appl. Phys. 38(13), 5149–5153 (1967). 55. W. R. Cook, D. F. Nelson, and K. Vedam, “High frequency properties of dielectric crystals: piezooptic and e


Introduction and state-of-the-art simulation techniques
The computer-aided analysis and optimization of the beam quality and power of mono-mode, continuous-wave (cw) laser sources is of increasing importance to reduce development cycle times and costs of high performance lasers.Especially the increasing number of active media and passive optical components used for the resonator setup and the resulting increase of free parameters requires fast software-based optimization algorithms.The fundamental aim of any optimization approach is to obtain a flexible, fast and accurate model of the laser resonator including all the physical effects which have a significant influence on the transversal structure and the power of the dominant resonator mode.In principle the power and lateral structure of the dominant transversal mode are influenced by several physical effects introduced by passive or active intracavity components.Passive intracavity components like lenses, mirrors, prisms, diffractive optical elements, do not amplify light and mainly influence the transversal resonator mode by refraction, reflection, absorption or diffraction effects [1,2].On the other hand active optical components amplify light.The dominant effects of solid-state active components like Nd:YAG on the transversal resonator mode are mainly nonlinear and inhomogeneous gain [3,4], thermal lensing [5] and stress-induced [6] or natural birefringence [7,8].Consequently all of the above mentioned physical effects must be included in a simulation technique enabling the accurate and flexible simulation of beam quality and power of the dominant transversal resonator mode.
In literature there are mainly three simulation approaches available for the calculation of the dominant transversal laser resonator mode: rigorous Maxwell solvers, the Gaussian beam propagation method and the method of Fox and Li.In principle Maxwell solver based approaches [9][10][11] are rigorous simulation techniques, which calculate the dominant resonator mode by the discretization of Maxwell's equations directly or the wave equation and the usage of periodic boundary conditions.One important class of Maxwell solvers is e.g. based on the Finite Element Method (FEM) [12,13].Basically these rigorous techniques can simulate any type of resonator setup containing various active and passive components.However in practice rigorous Maxwell solvers suffer from high numerical effort, resulting in a very restricted computational volume of the resonator cavities.Typically these techniques are used for microcavity lasers only.
Methods based on the Gaussian beam propagation technique [14,15], which are also sometimes called ABCD-Matrix approaches, are restricted in their applicability to paraxial resonator setups containing passive intracavity components introducing quadratic phase terms, like paraxial, thin lenses.Furthermore these approximated techniques only include physical effects of simplified active components.So thermal lensing is restricted to quadratic temperature distributions, ending up with the concept of Gaussian ducts [16,17].Also gain effects on the beam profile can be modeled only by Gaussian ducts, meaning that a quadratic gain profile [18] must be assumed.All of the above mentioned approximations of the Gaussian beam propagation method are necessary to ensure that the transversal shape of the resonator mode is Gaussian.This Gaussian mode shape can be just found in a limited number of real resonators setups.For the calculation of the laser beam power the Gaussian beam propagation technique is combined with a semi-classical separation approach [2,[19][20][21][22], meaning that the laser beam power is separately calculated from the transversal mode structure by rate equations and an integration of the photon density over the resonator volume.This separation results in approximations which are not fulfilled in all laser resonators, e.g. the effects on the beam profile caused by nonlinear gain saturation as well as diffraction losses caused by intracavity apertures are neglected.For the simulation of natural or stress-induced birefringence effects on the beam quality, power and the polarization state of the resonator mode, the Gaussian beam propagation method can be extended by a Jones matrix analysis [23] as shown in [24].In this case the polarization state of the beam and transversal mode structure are calculated separately.Consequently still fundamental or higher order Gaussian mode shapes must be assumed.
A more general approach which is not restricted to laser modes with Gaussian shape is the method of Fox and Li [25,26].The scalar nature of this method was in the past its main limiting factor.Therefore in a previous publication we have presented a very flexible, fully vectorial generalization of the Fox and Li approach for the calculation of the shape of the dominant laser resonator eigenmode [27].So far we have shown that, in principle, any resonator setup containing all kinds of passive optical components can be simulated.In this work we will show how thermal lensing and stress-induced birefringence effects introduced by an active cw laser crystal can be included in the calculation of the transversal resonator mode.Furthermore we include gain in our concept to enable the calculation of the laser beam power.In section 2 we will review the generalized and fully-vectorial Fox and Li algorithm.In section 3 we will show how active cw solid-state components can be included in our model using the vectorial beam propagation method (vBPM).Finally we will apply the model in a numerical example, showing how the resonator geometry and the application of intracavity components can decrease the negative effects of stress-induced birefringence and thermal-lensing on the beam quality and output power of the dominant laser resonator mode.

Generalization of the Fox and Li algorithm
In our recent work [27] we have shown that the fundamental mode of any kind of laser resonator can be calculated by the fully-vectorial Fox and Li algorithm.The main idea is that an electric field given at an arbitrary plane inside the resonator must reproduce itself after propagating once through the hole resonator.Mathematically this idea can be written as Here V 1 is the E x and V 2 is the E y polarization component of the harmonic field T representing the transversal resonator mode.Please note that the remaining electric field component E z and the three magnetic components H x , H y , H z can be calculated from V 1 and V 2 on demand by Maxwell's equations [28].The propagation of the field through the resonator is described by the resonator round trip operator R.This round trip operator consists of a sequence of approximated or rigorous component operators C m and free space operators P m,m−1 for different subdomains m of the resonator [27]: Figure 1 illustrates the round trip operator concept given by Eq. ( 1) and Eq. ( 2): The fully- or a non-diagonal structure dependent on the simulation technique applied to propagate the light through the active or passive intracavity component m.The non-diagonal operator components C m,21 and C m,12 in Eq. ( 4) describe the polarization cross-talk of the light propagation through the component m.In our recent work [27] we have focused on the general structure of the component and free space operators and the resulting structure of Eq. (1).It was discussed in detail that the presence of one or more non-diagonal component operators result in a non-diagonal round trip operator in Eq. ( 1), ending up with a coupled eigenvalue problem, where γ 1 = γ 2 .Up to now we have mainly shown resonators consisting of passive optical components and simplified active optical components neglecting birefringence and thermal lensing.In the following section we will show a non-diagonal component operator for the simulation of active components including stress-induced birefringence, rate equations and refractive index inhomogeneities due to temperature distributions.The applied technique is based on the vectorial beam propagation method (vBPM).The combination of the vBPM with the Fox and Li algorithm enables the calculation of the beam power and the fully vectorial transversal mode shape of a resonator including thermal lensing, birefringence and nonlinear gain saturation.

Simulation of active components by the vectorial beam propagation method
In this section we will first introduce the vectorial Helmholtz equation describing the light propagation through a weakly inhomogeneous, nonlinear and birefringent medium.Then we will discuss a split-step beam propagation algorithm which is a fast Fourier Transformation (FFT) based numerical solver of this wave equation.Finally in this section we will show how the different material models describing the light amplification, stress-induced birefringence and thermal lensing are included in the wave equation.Please note that in principle most of the following material models in this section are known in literature separately.However in this work it is the fist time to our knowledge that these material models are combined and used for the simulation of light inside a laser cavity.Furthermore in section 3.2 we will add a novel semi-analytical approach for the numerically efficient inclusion of light amplification effects within the vBPM.

Fundamentals of the vectorial beam propagation method
The light propagation through a weakly birefringent and nonlinear medium with weak refractive index inhomogeneities given by the dielectric constant can be described by the Helmholtz equation in matrix form [29]: with the identity matrix I I I and the matrix Here ) is the transversal nabla operator, r r r = (x, y, z) T is the position vector and k 0 is the wavenumber in vacuum.Using the assumption of weak birefringence, nonlinearity and inhomogeneity effects in the active medium, we can approximate the tensor components of ε ε ε(r r r,V V V ) in Eq. ( 7) by [29] with the indices α = 1, 2 and β = 1, 2. n 0 is the complex linear, isotropic and homogeneous refractive index share of ε ε ε(r r r,V V V ).ñαβ (r r r,V V V ) represents the complex nonlinear, anisotropic and inhomogeneous refractive index shares and can be splitted in Here Δn therm αβ (r r r), Δn birefring αβ (r r r,V V V ) and Δn gain αβ (r r r,V V V ) separately describe the effect of thermal lensing, stress-induced birefringence and nonlinear gain on the refractive index.The actual structure of these variables will be discussed in detail in the following subsections.
The Helmholtz equation given by Eq. ( 6) can be solved by a split-step beam propagation method [29,30] or a finite-difference beam propagation method [31] to calculate the light propagation through the active component.In the following we will give a numerically efficient vectorial split-step beam propagation method which is based on the work of Thylen and Yevick [29].The complex amplitude V (x, y, z 0 + L) of a harmonic field in a plane perpendicular to the optical axis, located at the axial position z 0 + L at the end of the active component, can be calculated from the field V (x, y, z 0 ) in front of the active component by a split-step vBPM in symmetrized form: ) with the number of vBPM steps S = L/Δz.Here L represents the total propagation distance and Δz a single split-step distance.Please note that we have skipped the (x, y) variable dependency of the field only for better readability.P SPW s is the angular spectrum of plane waves operator [32], given by the diagonal operator: with Here F is the Fourier transformation with the conjugate variables ρ ρ ρ = (x, y) T and κ κ κ with the wavenumber k = 2π λ n 0 and the vacuum wavelength λ .Figure 2 illustrates the split-step vBPM given by Eq. (10).The operator Δz on the right hand side of Eq. ( 10) describes the influence of the nonlinear, inhomogeneous and birefringence terms of the active component and is modeled by the non-diagonal operator: with for α = 1, 2 and β = 1, 2. ñαβ (z ,V V V ) can be interpreted locally as the nonlinear, anisotropic and inhomogeneous contribution of the active component's refractive index and is equal to the expression given by Eq. ( 9).For better readability the (x,y) dependency was skipped in the notation.Consequently the integral on the right hand side of Eq. ( 14) describes the optical path length introduced to rays propagating parallel to the optical axis due to this residual refractive index.The contribution of gain can be separated by substituting Eq. ( 9) into Eq.( 14) resulting in: with and with For small slab sizes Δz the integral in Eq. ( 17) can be approximately solved numerically, ending up with the nonlinear operator In principle also the integral in Eq. ( 16) can be approximated in the same way as Eq. ( 18).However especially for large small-signal gain or strong nonlinear gain saturation inside the active medium, this approximation would require a very small vBPM step size Δz which would increase the computational effort.That is why in the following subsection we will introduce a novel semi-analytical solution of the integral in Eq. ( 16), ending up with a larger acceptable step size Δz and consequently a numerically more efficient vBPM.
To avoid aliasing effects within the vBPM due to propagation of the electric field to the edge of the computational window in a transverse plane (x, y) an absorbing boundary condition must be introduced.The absorber function, which is applied on both field components after each nonlinear operation C NL s , is given in Eq. ( 3) of Jaijonmaa's work [33].For the realistic simulation of light propagation through the active laser medium by the vBPM, the effects of thermal lensing, birefringence and gain saturation have to be included by the corresponding terms in Eqs. ( 16) and (18).The structure of these terms is discussed in the following subsections.

Inclusion of gain by rate equations
For the simulation of cw light amplification within the active medium, we extend a quasistationary approach which is known in literature [34,35] by a novel semi-analytical treatment of the gain saturation in Eq. ( 16).This semi-analytical approach allows an increased vBPM step size Δz and consequently a reduced computational time of the vBPM compared to the existing quasi-stationary technique.In this subsection we will first introduce the most important steps of the quasi-stationary approach in our nomenclature and then present our novel semi-analytical solution of Eq. ( 16).
For the interpretation of the complex value Δn gain αβ (r r r,V V V ) in Eq. ( 16), we use the Lambert-Beer law [15], which describes the light amplification of a cw plane wave inside an active medium by: dI dz = gI (19) with the gain coefficient g and with the light intensity defined by the time-averaged Poynting vector S S S , which can be approximated for paraxial fields by [36] Here c is the vacuum speed of light, Re( ñ0 ) the real part of the linear refractive index share of the active medium, V are the polarization components of the complex amplitude of the harmonic field and ε 0 is the vacuum permittivity.Formally Eq. ( 19) has the solution: which has a structure equal to that of Eq. ( 16).Consequently the gain coefficient g in Eq. ( 21) is related to Δn gain αβ (r r r,V V V ) in Eq. ( 16) by: for α = β .For α = β there is Δn gain αβ (r r r,V V V ) = 0, meaning that polarization cross-talk by the light amplification effect itself is neglected.However of course we separately include polarization cross-talk due to stress-induced birefringence by the term Δn birefring αβ .Please note that the additional factor 2 on the left hand side of Eq. ( 22) is due to the fact that Eq. ( 16) is formulated for fields V and Eq. ( 21) is formulated for intensities I. Fields and intensities are related with each other through Eq. (20).The gain coefficient g in Eq. ( 22) depends on the energy level diagram of the active medium and its corresponding system of rate equations.The most important rate equation systems (e.g.4-level-energy-diagram, 3-level-energy-diagram) can be solved analytically for cw laser operation, using a quasi-stationary approach.A good overview of stationary solutions of various rate equation systems can be found in [37].For example in the case of a 4-level-energy-diagram Eq. ( 22) leads to: for α = β with the small-signal gain g 0 and the saturation intensity I s .In principle now Eq. ( 23) can be substituted into the operator given by Eq. ( 16) which can be solved approximately in the same way as Eq. ( 18).However this approximated numerical solution requires a decreasing vBPM step size Δz for an increasing small-signal gain g 0 and for a decreasing saturation intensity I s .Therefore in the following we will suggest a semi-analytical solution of the integral given in Eq. ( 16).This semi-analytical solution in principle works for any step size Δz as long as the small signal gain g 0 can be assumed to be z invariant within this step.Consequently, this semi-analytical vBPM operator is not limited by the nonlinear gain saturation factor anymore.Only thermal lensing, birefringence and diffraction effects are finally limiting the maximum step size of the vBPM.The operator given by Eq. ( 16) can be expressed by the two first-order differential equations with where Δn gain αβ (r r r,V V V ) was substituted by the expression given in Eq. ( 23).As shown in Appendix A, the exact solution of Eq. ( 24) is with Please note that there is no explicit analytical expression for the LambertW function.Nevertheless there are several math libraries available which can calculate W (a) numerically (see e.g.[38]).However these math libraries are typically limited in their applicability for larger g 0 and/or small V 2 sat .Especially the exponential dependency of Eq. ( 26) on g 0 is crucial.The argument of W in Eq. ( 26) might exceed the maximum possible value which can be expressed by a number in the floating-point format.In the examples we will show in section 4, this was the case and we were not able to apply any math library directly on Eq. (26).However, especially for large g 0 and/or small V 2 sat , a semi-analytical inclusion of the non-linear gain is essential to increase the vBPM step-size Δz.Fortunately there is an approximate explicit expression for W (a) which is particular suitable for large arguments a [39,40]: Please note that for the calculation of the intensity I of the resonator mode the intensities of the forward and backward propagating waves I + and I − must be taken into account for stable or unstable standing wave cavities using the equation: Here we have used the assumption that several slightly different axial modes oscillate with the same transversal mode in the resonator [15,41].Since the intensity minima of different axial modes are located at different axial positions, spatial hole burning can be neglected.This assumption is reasonable for several realistic solid-state lasers.Furthermore the axial modes are incoherent to each other, due to their slightly different wavelengths.Consequently it is suitable to model only one axial mode to represent all oscillating axial modes in the cavity and neglect the interference between the counterpropagating waves in Eq. ( 29).If the resonator round trip operator is built up according to Eq. ( 2), two vBPM propagation operators given by Eq. ( 10) are required.One vBPM operator is used to propagate the light in +z direction ending up with the intensity distribution I + (x, y, z) within the active medium volume.This distribution is stored for the evaluation of Eq. ( 29) in the reverse propagation by the second vBPM operator calculating and storing I − (x, y, z).Now again the forward propagation can be performed using I − (x, y, z) of the previous step, ending up with a new I + (x, y, z).This procedure can be repeated iteratively till a convergence of the round trip operator result is reached.Expressions of g and the resulting semi-analytical formulation of the vBPM for other energy-level-diagrams in the stationary case can be found in a similar way.

Inclusion of thermal lensing
An inhomogeneous temperature distribution T (x, y, z) within the volume of a solid-state active medium induces thermal dispersion and birefringence [5,6,42].The temperature distribution is brought about by the generation of heat Q(x, y, z) due to non-radiative relaxation processes in an inhomogeneous pumped active medium in combination with heat flow to the outer periphery due to cooling [41].Due to the cw operation of the laser, we will assume in the following that the inhomogeneous temperature distribution is constant over time.In literature there is a huge variety of possibilities how the 3D temperature distribution T (x, y, z) in a solid-state active medium can be obtained.In principle these techniques can be classified into simulation and experimental measurement techniques and their useful application depends on the concrete situation of the resonator setup, pump distribution, active medium material and so on.We would like to refer interested readers to [43][44][45][46] for different simulation approaches and to [47,48] for experimental measurement techniques.
Once the stationary temperature distribution is obtained from one of the above given approaches, thermal dispersion and consequently thermal lensing, which is also sometimes called thermal beam distortion, can be calculated using vBPM.If we assume an isotropic thermal dispersion, the term Δn therm αβ (r r r) in Eq. ( 9) can be calculated from the temperature distribution by a Taylor series: where T 0 should be chosen close to the average operation temperature of the active medium.The thermal coefficients dn dT , d 2 n dT 2 , ... can be fitted from experimental measurements of the refractive index for different temperatures.For Nd:YAG the coefficient dn dT is around 7.3 × 10 −6 1/ • K for T 0 = 40 • C [49].Please note that in Eq. ( 30) the refractive index share due to thermal lensing is not a function of the resonator mode V V V , because a linear material response is assumed.In principle it is also possible to get Δn therm αβ (r r r) directly by measurements.Direct measurement techniques for the determination of the thermal lenses are described e.g. in [5,50].

Inclusion of stress-induced birefringence
Besides isotropic thermal dispersion, an inhomogeneous temperature distribution results in anisotropic mechanical stress in solid-state active media.Due to the photoelastic effect this stress induces birefringence which can be described according to [51] by an indicatrix under stress with the principal coordinate system axes x, y, z and the relative dielectric impermeabilities B 0 xx , B 0 yy , B 0 zz of the unperturbed active medium's indicatrix.The changes of the relative dielectric impermeability tensor ΔB i j = ∑ kl p i jkl εkl (32) with the indices i, j, k, l = 1, 2, 3, are caused by the strain tensor εkl and the photoelastic fourthrank tensor p i jkl .Dependent on the temperature profile, the strain tensor can be obtained from analytical formulas [46,51] or from a finite-element analysis [52,53].Please note that for the calculation of Eq. ( 32) the strain tensor and the photoelastic tensor must be in the same crystal coordinate system given by the crystal cut direction [51].If necessary, the tensors have to be transformed into the desired coordinate system using coordinate rotations [51].In this work the changes of the relative dielectric impermeability tensor are calculated by the photoelastic tensor.Furthermore it is possible to calculate the photoelastic tensor by p i jkl = ∑ mn q i jmn C mnkl (33) from the piezo-optic tensor q i jmn and the elasticity tensor C klmn .Tabulated tensor data for many crystals can be found in [54,55].Please note that in the above given literature the tensor data were measured for temperature regions which might be different from the operation temperature region of the actual active medium.In this case the photoelastic tensor or the piezo-optic tensor could be measured e.g. by interferometry [56].Before the contribution of n birefring αβ (r r r,V V V ) on the nonlinear vBPM step given by the operator C NL s z 0 + s − 1 2 Δz in Eq. ( 13) and Eq. ( 14) is derived, we would like to point out again that the operator C NL s z 0 + s − 1  2 Δz can be interpreted locally as a multiplication of the incident field with a polarization and positiondependent optical path length.This optical path length is obtained by adding up the optical paths of rays propagating parallel to the optical axis within one vBPM slab with thickness Δz.As illustrated in Fig. 3 different optical path lengths for different positions ρ ρ ρ = (x, y) T in the initial transversal plane z 0 + (s − 1)Δz are obtained due to the inhomogeneous ñαβ (z ,V V V ) distribution.To include the polarization dependency of the optical path length due to birefringence, each ray carries the local electric field vector (V 1 (r r r),V 2 (r r r)) T .For small slab thicknesses Δz it is allowed to assume that the direction of the rays is not changing within a single nonlinear vBPM step.Consequently all rays, which can be also interpreted locally as plane waves, propagate in the normalized direction û u u = (0, 0, 1) T and we can find for each ray principal axes of the indicatrix given by Eq. (31).Therefore the eigenvalue problem [51] has to be solved, ending up with the two eigenvalues 1/n 2 x and 1/n 2 y .The directions of the two orthogonal axes in the x-y-plane are given by the eigenvectors V V V x and V V V y .Before we evaluate the optical path length of each ray in the nonlinear vBPM step, it is convenient to project the polarization vector in the x-y-coordinate system onto the principal x -y -coordinate system.So we can rewrite the nonlinear vBPM operator given by Eq. ( 13) and Eq. ( 14) as: with Fig. 3. Physical interpretation of a single nonlinear vBPM operation C NL s within a single slab with thickness Δz: the result of the operation is obtained by mapping the complex optical path length (OPL) introduced by the complex refractive index modulation within the slab onto the initial field distribution (V 1 ,V 2 ) T , which is sampled on a transversal grid ρ ρ ρ n .Due to the small slab thickness and a paraxial initial field, the OPL can be evaluated along rays propagating parallel to the optical axis.Each transversal grid point is the origin of a ray/local plane wave with a certain polarization vector V V V in the transversal x-y-plane.Due to the stress-induced birefringence, the OPL for each ray/local plane wave has to be evaluated for the two orthogonal polarization components parallel to the principal axes α and β .Due to the position-dependent nature of the electric permittivity tensor of the active medium, for each ray another orientation and shape of the indicatrix must be considered in the OPL analysis.In the left part of the figure two example indicatrices for different rays with different polarization vectors V V V are given.
for α = 1, 2, and the polarization vector projection operator per position ρ ρ ρ = (x, y) T in the transversal plane in front of slab s.After the nonlinear vBPM step, the inverse projection operator projects the polarization vector back into the original x-y-coordinate system.Finally we can write Due to the assumption of an isotropic behavior of the thermal lensing and nonlinear gain effects, the terms n therm α β and Δn gain α β are equal to the corresponding terms in the x-y-coordinate system given by Eq. ( 30) and Eq. ( 23).If we assume a dielectric active solid-state crystal which has an isotropic behavior in the absence of stress, we can write B 0 xx = B 0 yy = B 0 and ΔB xy = ΔB yx .In this case the eigenvalue problem given by Eq. ( 34) has the two solutions [51]: Consequently we can express Δn birefring α β (r r r,V V V ) by [51]:

Numerical example
In the previous sections novel concepts for the calculation of the fully-vectorial transversal eigenmode of a cw solid-state laser were introduced including nonlinear gain, birefringence and thermal lens effects.In the following we will apply these techniques exemplarily to improve the beam quality and the output power of the cw flash-lamp-pumped Nd:YAG laser given in Fig. 4. The laser resonator is a linear cavity consisting of two spherical mirrors M 1 and M 2 with Fig. 4. Initial setup 1: solid state laser, consisting of two spherical mirrors M 1 and M 2 , a linear polarizer and a flash-pumped Nd:YAG crystal.To suppress higher order modes, an aperture is placed in the plane of the outcoupling mirror M 1 .The parameters of the Nd:YAG crystal are given in Table 1.
focal lengths f M1 = f M2 = 100 mm.Mirror M 1 is the outcoupling mirror and has a reflectance R M1 = 0.9.It is desired that the beam quality M 2 emitted by the laser be close to 1 and that the laser operates only in a single transversal mode.Therefore the circular aperture diameter of mirror M 1 has to be chosen small enough to suppress higher order resonator modes.In this example the circular aperture of mirror M 1 should be 250 microns.All light outside of this aperture is perfectly blocked.Resonator mirror M 2 has a sufficiently large aperture size and a reflectance of R M2 = 1 to ensure that it is not introducing any additional resonator losses.It is desired that the emitted light is perfectly linearly polarized under an angle θ 0 = arcsin(y/x) = 30 • .Therefore a linear wire-grid polarizer is placed directly in front of the outcoupling mirror.
The polarizer is orientated under an angle of 30 • generating the requested linear polarization direction.The thickness of the polarizer can be neglected, so that the length of the resonator is 76 mm.The active medium is a Nd:YAG crystal, with the same parameters and flash-pump as discussed in [5,57].The parameters of the active medium are given in Table 1.Please note that the pump efficiency is not equal to the pump quantum efficiency.The pump efficiency is the fraction of the electric power of the flash lamps which is converted into stimulated emission.
It is calculated by the multiplication of the first three efficiencies given in Table 1 of [57].
The active medium and two flash-lamps are placed in a double-elliptical cylinder to ensure a high pump-transfer efficiency and a homogeneous pump distribution in the active medium [5,57].For simplicity, in the following we will assume that the heat which is dissipated in the active medium is homogeneously distributed over the hole active medium and an annular cooling configuration is used.Also the pump power is homogeneously distributed over the whole Nd:YAG rod.Please note that these assumptions in principle are not necessary to apply the concepts described in the previous sections.However the assumptions are pretty realistic for the flash-lamp double-elliptical pump configuration, ending up with the following cylindrical temperature distribution within the active medium [57] Table 1.Parameters of the Nd:YAG active medium.If no extra citation is given for the parameter, its value was taken from Koechner [5,57].
with the temperature in the rod center T (0) and at the rod surface T (0.5d) given in Table 1.Plugging Eq. ( 42) into Eq.( 30) and using the values for dn/dT given in Table 1 results in a thermal lens of Furthermore Eq. ( 42) can be used to solve the eigenvalue problem given by Eq. ( 34), ending up with azimuthally stress-induced birefringence [5].Then the nonlinear operator in Eq. ( 35) has the form [41]: with the azimuthal angle θ and the contributions of the stress-induced birefringence on the refractive index: = 0. Due to the broad spectrum of pump light emitted by the flash lamps, several pump levels are used to fill the upper laser level.Therefore we can express the small-signal gain g 0 in Eq. ( 23) by [2] To calculate the dominant eigenmode of the above given example resonator setup, we have implemented the vBPM component operator in the optics design software VirtualLab [58] and solved the eigenvalue problem given by Eq. (1) using the minimal polynomial extrapolation method (MPE) [59] and a random phase distribution as initial condition.For the vBPM the active medium was discretized in z-direction in 30 slabs.The mirrors were simulated by a thin element approximation (TEA) [28] and the linear polarizer by a Jones matrix [41].The E x and E y polarization components of the dominant transversal laser resonator mode in two different planes in the resonator are given in Fig. 5.The output power of the beam is 0.19 W. Furthermore we have applied the second-order momenta method [41] on the transversal mode to calculate the beam quality in the plane of the outcoupling mirror.It is M 2 x = M 2 y = 2.8 in x-and y-direction, respectively.If we analyze the transversal mode structure and power at different positions in the laser cavity, we can see that there are significant depolarization losses at the linear polarizer and diffraction losses at the aperture of the outcoupling mirror.The corresponding values are given in Table 2.In a next step we improve the diffraction losses by compensating the thermal lens of the active medium.Therefore we have to rearrange our resonator setup.We replace the spherical outcoupling mirror by a plane mirror with reflectance R = 0.9 and use a Fourier lens in combination with an aspherical mirror as shown in Fig. 6.The Fourier lens in combination with a smaller mirror aperture suppresses the higher order resonator modes.The aspherical mirror compensates the wavefront deformations caused by thermal lensing.For the design of the aspherical surface of the mirror the vBPM operator is used by the following steps: 2. The vBPM operator and the other component and free-space operators mentioned below are applied to propagate the Gaussian beam to a plane in the desired aspherical mirror position z 0 .Then the phase Φ(x, y, z 0 ) of the propagated field is calculated in this plane.
3. After that the phase of the field is conjugated: Φ(x, y, z 0 ) → Φ * (x, y, z 0 ) 4. Due to the fact that the beam in the resonator has to pass the active medium twice in a full round trip, the transmission function of the aspherical mirror must be t(x, y, z 0 ) = exp(i2Φ * (x, y, z 0 )).
Once the aspherical surface is designed, we can apply again the round trip operator concept in combination with the vBPM to calculate the dominant resonator eigenmode.For the free-space propagation we have used the angular spectrum of plane waves operator [32].The thin lens, the  plane mirror and the aspherical mirror are simulated by TEA.The polarizer is simulated again by a Jones matrix.Figure 7 shows the dominant eigenmode at two different positions in the resonator.The output power of the beam could be improved to 10.9 W and the beam quality to M 2 x = M 2 y = 1.1.As shown in Table 2 the polarization losses are still pretty high.Consequently we have to modify the resonator setup further to reduce the polarization losses.To that end we use a 90 • quartz polarization rotator as described in [60] ending up with the resonator setup given in Fig. 8.The quartz rotator and the linear poarlizer are simulated by Jones matrices and the lenses and mirrors by TEA.Again MPE is used to solve the eigenvalue problem.The corresponding transversal eigenmode in two different planes is given in Fig 9 .Clearly the depolarization loss could be decreased by this resonator setup.The improved outcoupling power and beam quality of the laser are given in Table 2.

Summary
In this work the fully-vectorial formulation of the Fox and Li method was combined with the vectorial beam propagation method (vBPM), enabling the calculation of the dominant transversal eigenmode of realistic cw solid-state lasers.The vBPM was used to simulate light propagation through the active medium including the effects of thermal lensing, stress-induced birefringence and nonlinear gain saturation.To achieve a numerically efficient simulation, a semi-analytical expression of the vBPM was introduced for modeling the nonlinear gain effects.In the end the techniques were applied to optimize the output power and beam quality of a cw Nd:YAG laser.To that end an aspherical lens and a 90 • quartz rotator were added to compensate the thermal lens and the stress-induced birefringence.By these modifications of the resonator geometry the beam quality M 2 could be improved from 2.2 to around 1.3 and the output power from 0.19 W to 71.4 W.

A. Derivation of semi-analytical inclusion of nonlinear gain
In the following Eq.( 26) will be derived.Starting from the two coupled differential equations given in Eq. ( 24), we can write Adding of both equations leads to (48) which has the exact solution for g 0 being Δz invariant.Here W (a) is the LambertW function (also called product logarithm) which can be defined by a = W (a) exp [W (a)] [39].c 1 is a constant, depending on the initial value of the differential equation at Δz = 0.It can be written as: with ln being the natural logarithm in the base e.In the next step we separate the two terms on the left hand side of Eq. ( 49).We remind that the polarization of the electric field is maintained within the operator given by Eq. ( 16).Consequently we get: In addition to the real-valued gain g, an additional phase shift ig is introduced in reality by the active medium.Typically this additional phase shift only depends on the angular frequency ω 0 of the field oscillating in the resonator cavity.In a first approximation there is no dependence on the position, ending up with a constant phase shift for each wavelength [15].Consequently in a single nonlinear C gain αβ ,s step the wavefront of the field does not change it's shape.In this case we can rewrite Eq. (51) to In this work we will neglect the constant phase shift in our calculations, meaning that ig = 0, because we are just interested in the transversal mode profile of the laser.In this case we end up with Eq. ( 26).As discussed in our previous publication [27], in principle we can also analyze the axial modes of the laser resonator, by the calculation of the eigenvalue's phase term.Then the constant phase shift would be important for all wavelengths except of the center wavelength of the spectral gain profile, where the phase shift is zero anyway [15].

Fig. 1 .
Fig. 1.Exemplary resonator setup consisting of several different active and passive optical components.The transversal resonator mode can be calculated in a arbitrary plane using a sequence of diverse component operators C and free space propagation operators P describing a full resonator round trip.Each operator solves Maxwell's equations in a rigorous or approximated manner for a single subdomain of the resonator.

Fig. 2 .
Fig. 2. Schematic illustration of the symmetrized vBPM for light propagation through an active intracavity component: the component is subdivided into segments of width Δz.Within a segment only linear effects, like diffraction or linear absorption are included in the angular spectrum of plane waves operators P SPW s .The effect of nonlinearity is introduced by the nonlinear operators C NL s in the center of the segments (represented by dashed lines).

Fig. 5 .
Fig. 5. Dominant transversal eigenmode of the laser resonator setup given in Fig. 4. The upper row shows the V 1 (a) and V 2 (b) field components of the mode in the plane of the outcoupling mirror M 1 .The lower row shows the V 1 (c) and V 2 (d) field components of the mode in the plane between the active medium and the linear polarizer.In both rows the propagation direction of the mode is towards the outcoupling mirror and the fundamental transversal eigenmode has a non-Gaussian shape due to birefringence, diffraction and thermal lensing effects.

Fig. 6 .
Fig.6.Resonator setup 2 with thermal lens compensation: a stable Fourier transform resonator geometry is used to improve the beam quality and output power.To ensure that the thermal lens in the flash-pumped Nd:YAG laser is compensated an aspherical mirror was designed.The parameters of the Nd:YAG crystal are the same as in the initial setup.

Fig. 7 .
Fig. 7. Dominant transversal eigenmode of the laser resonator setup 2 given in Fig. 6.The upper row shows the V 1 (a) and V 2 (b) field components of the mode in the plane of the outcoupling mirror M 1 .The lower row shows the V 1 (c) and V 2 (d) field components of the mode in the plane behind the active medium, next to the linear polarizer.In both rows the propagation direction of the mode is towards the outcoupling mirror.

Fig. 8 .
Fig.8.Resonator setup 3 with thermal lens and birefringence compensation: a stable Fourier transform resonator geometry is used to improve the beam quality and output power.To ensure that the thermal lens in the flash-pumped Nd:YAG laser is compensated an aspherical mirror was designed.The initial Nd:YAG crystal is split into two rods with equal length and optical pump.Stress-induced birefringence is compensated by placing a 90 • quartz polarization rotator and a thin lens between the two Nd:YAG crystals.

Fig. 9 .
Fig. 9. Dominant transversal eigenmode of the laser resonator setup 3 given in Fig. 8.The upper row shows the V 1 (a) and V 2 (b) field components of the mode in the plane of the outcoupling mirror M 1 .The lower row shows the V 1 (c) and V 2 (d) field components of the mode in the plane behind the active medium 1 next to the linear polarizer.In both rows the propagation direction of the mode is towards the outcoupling mirror.

Table 2 .
Output power and beam quality of laser resonator setups 1-3.The depolarization loss and the diffraction loss are defined by the ratio of the beam power before and after the polarizer and mirror aperture, respectively.