Iterative 3D modeling of thermal effects in end-pumped continuous-wave HO3+:YAG lasers

We present a highly accurate model for end-pumped continuous-wave Ho3+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Ho}^{3+}$$\end{document}:YAG laser resonators, based on a vectorial beam propagation method (BPM) algorithm. Thermal effects in the laser crystal like thermal lensing and stress birefringence are taken into consideration. The implementation of these effects is based on an iterative method, which updates the temperature and displacement in the crystal at specific roundtrip intervals. An experimental Ho3+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Ho}^{3+}$$\end{document}:YAG resonator cavity in continuous wave operation is used to validate the model. At higher pump powers, the resonator reaches its stability limit and shows beam distortions, which is well replicated with the presented model.


Introduction
High-power lasers operating in the 2μm spectral region can be used for a variety of applications, ranging from medical treatments and countermeasures to the pumping of midinfrared optical parametric oscillators (OPOs). For the latter application, solid-state lasers have shown great potential given the high pulse energies that these laser sources allow to reach. Thulium-based lasers are often used in this spectral region, because they can be pumped with highly efficient laser diodes at 800nm and have the advantage of high conversion efficiencies. However, Ho 3+ doped media like Ho 3+ :YAG are efficient for emitting at wavelengths greater than 2μm and, therefore, may be preferred for some applications, such as ZGP OPO pumping, by taking advantage of lower background absorption for wavelengths around 2.1μm and beyond [1]. The increasing availability of high-power thulium-based fiber lasers has made Ho 3+ :YAG more attractive for these applications. Since the absorption cross section of holmium is at a maximum at around 1.9μm depending on the host material, thulium fiber lasers are ideal high power pump sources at this spectral region. One of the main limitations for power-scaling solid-state laser resonators is the loss of resonator stability due to thermal lensing in the laser crystal. This can get even more problematic in holmium crystals, due to high energy transfer upconversion (ETU) depending on the doping concentration of the crystal. Consequently, the optimization of such a system can be much more difficult. A simulation model can greatly assist in the design and optimization of Ho 3+ :YAG resonators for high power, efficiency, and good beam quality by allowing more insight into the internal processes of the system.
In this work, we present a powerful model developed for accurately describing end-pumped Ho 3+ :YAG lasers. The general approach is based on a sequential field propagation method. A modular simulation environment allows for the inclusion of a variety of different element operators, like mirrors and laser crystals. By combining these operators, a linear or ring resonator can be modeled. To calculate the laser gain of the system, rate equations are used to describe the radiative and non-radiative processes in a spatially resolved crystal. The model for holmium includes all the relevant manifolds and their respective energy transfer processes. By applying the rate equation model, not only the laser gain but also the heat load of the crystal is accurately modeled. Subsequently, the temperature distribution as well as the displacement distribution and resulting mechanical stresses are calculated. The respective differential equations are solved with finite difference methods (FDMs), which are chosen for their simplicity. Resulting thermal effects like thermal lensing and birefringence can then be calculated for the crystal. In the last section, a continuous-wave Ho 3+ :YAG resonator setup is used to validate the model by comparing the laser output for the model and experiment. All the procedures described in this manuscript where implemented using the MATLAB language, in some cases incorporating the underlying theory and algorithms described in the references.

CW simulation approach
The fundamental cw approach is based on the sequential field propagation of a two-dimensional complex field. Using a modular environment, different element operators are placed in between the propagation steps to accurately model a laser resonator. Each element operator uses unique algorithms to transform the field for their appropriate application. They include mirrors, lenses, glasses, polarizers, and laser crystals. Using these operators, a linear or ring resonator can be modeled by following the approach of Fox and Li [2].
The resonator gain is recalculated iteratively over a certain amount of roundtrip steps until an end condition is met. The flowchart depicted in Fig. 1 shows the modular structure and basic operation of the model. When the laser gain in the resonator has stabilized, the thermal steady state of the laser crystal is calculated in between the roundtrip steps. This process is repeated until the thermal recalculation reaches a convergence criteria. Resulting thermal lensing and birefringence effects can have a large influence on the generated power output and beam quality of the laser.

Propagation method
Two-dimensional field propagation is usually achieved by solving Maxwell's equations for transverse electromagnetic waves. Here, a variety of different solvers can be used. Direct Maxwell solvers [3,4] like finite element methods (FEM) or finite difference time domain methods (FDTD) provide the highest accuracy at the cost of high simulation effort. More simplified approaches like Gaussian beam propagation [5] use specific field distributions and typically have less accuracy, but may provide lower simulation efforts. However, the beam propagation method (BPM) provides a good compromise between accuracy and numerical effort. For this model, a split step BPM approach developed by Fleck et al. [6] is used. The method is based on the slowly varying envelope approximation, which simplifies Maxwell's equations by assuming the complex amplitude only varies slowly in time and space. This may cause problems for some applications where a broad-band spectrum is to be propagated, but has proven to be very compatible with this model. Maxwell's wave equation is solved in the parabolic form where n is the refractive index, E is the slowly varying field amplitude with the Intensity I given by I ∝ |E| 2 and k is the wave number k = 2 ∕ . For the numerical solution, equation (1) is solved in a discretized form as a finite twodimensional Fourier series wherein a square of the size L = N x = N y is discretized in m and n pixels in x-y space. With this, the initial field E mn (z) is propagated along the distance z . In between the propagation step, at half the step length the phase delay through a refractive medium is added to the complex field as shown in Fig. 2. To include polarization effects in the later simulation, this propagation model can be extended vectorially by propagating one field for each polarization component and coupling the fields at specific element operators. In this model, the split step BPM is used to propagate the laser field as well as the pump field, so that all relevant refraction and diffraction effects are considered.

Ho 3+ rate equations
To calculate the laser gain and internal heat load of the active laser crystal, a rate equation model is used to determine the energy state population of the relevant manifolds. For quasi-three-level laser ions like Ho 3+ , the simplified laser process by stimulated emission can be described by just the upper laser level and the ground level, since the lower laser level gets depopulated quickly due to fast thermal transitions. The temporal population dynamic of the upper laser level can be expressed by the following rate equation in differential form [7] where M ij stands for the rate constant of phonon relaxation from the i to the j manifold, 1 i,sp describes the rate constant of spontaneous emission and W L and W P are the rates of stimulated laser emission and pump absorption, respectively. In this simplified model the upper laser level population N I7 and the ground state population N I8 can be calculated analytically in the steady state. The heat load in this model only depends on the quantum defect between the pump wavelength and the laser wavelength, as well as heat generation due to phonon transitions of the upper laser level. A more accurate rate equation model also has to take into account energy transfer upconversion (ETU) and cross-relaxation processes, since especially ETU is an important contributing factor to heat loss in holmium, as mentioned earlier. By including these transfer processes, more manifolds have to be considered. For the laser transition at 2.09μm from the I7 → I8 manifold, the I8, I7, I6, and I5 manifolds now have to be taken into account. The expanded set of rate equations is based on a 2D model for Er 3+ :YAG described by Eichhorn [8] and is changed to fit the specific requirements for holmium at this laser transition. In this case, the rate equations have to be solved for all four relevant energy levels resulting in the following set of nonlinear equations:   Figure 3 shows the schematic laser level diagram of the relevant manifolds for holmium, including all energy transfer processes of interest. The constants describing ETU k ij and cross-relaxation C ij quantify the ion-pair transitions taking place at different energy levels and are typically dependent on the doping concentration of the material [9]. Another level of complexity stems from the branching ratios for spontaneous emission ij , which determine the probability of the spontaneous decay from i to a certain manifold j [10]. Many of the constants describing the population dynamics and spectral properties of holmium can be found in literature [11][12][13]. However, as the relevant literature shows variation in values, it may be useful to validate some of the properties in experiments to achieve more precise simulation results. The spectroscopic material parameters used for an example simulation with this model are presented in Appendix A ( Table 1).
Using the population dynamics of the spatially resolved crystal, the laser gain and absorption of the electromagnetic field E can be calculated by applying the Lambert-Beer law to consider the attenuation in the gain medium. Subsequently, the heat load in the crystal can be determined by the following equation: [8] The first term in equation (8) accounts for the heat load by the absorption or emission of photons. The change of the intensity distributions in z direction I ± P,L is accounted for the pump and laser fields respectively in negative and positive propagation direction. To include the radiative power loss, the spontaneous emission from the different manifolds is subtracted in the second term. Here, E i→j = hc i is the average photon energy for a specific i → j manifold transition. The rate equations (4-7) are unique to each active laser material at a specific laser transition and may be modified accordingly, if a different material has to be simulated.

Temperature distribution
The heating of the crystal can lead to a variety of detrimental effects, for example, stability loss due to thermal lensing, power losses by depolarization due to birefringence, and (8) the temperature dependence of different material parameters like the emission and absorption cross sections. Solid-state laser crystals for high power applications typically have a heat sink in place, leading to a specific thermal distribution. Therefore, the crystal temperature has an important influence on the overall resonator performance and should not be neglected in a precise simulation model. With the heat load in the crystal accurately calculated by the rate equation model in a three-dimensional grid, the temperature distribution can be solved numerically with high accuracy. The following second-order partial differential equation (PDE) describes the thermal diffusion where k is the heat transfer coefficient, C is the heat capacity, and is the density of the material. The heat load calculated by the rate equation model is included through the internal energy deposition H. There are several numerical solving methods suitable to this problem. For this model, the finite difference method (FDM) was chosen for its relative simplicity. More specifically, the three-dimensional alternating-direction implicit method (ADI) is used, as it has high stability and accuracy [14]. The method is based on an iterative approach, splitting the individual time steps in the three spatial directions as depicted in Fig. 4. To achieve an accurate representation of the real temperature in a crystal, the boundary conditions play an important role and have to be chosen correctly. In the case of the temperature diffusion, the boundary can either be defined by a Dirichlet or a Neumann condition. For the Dirichlet condition, the temperate values at the boundary are set to constant values T lat and T end for the lateral surface and the endfaces respectively. In cases where the temperature at the surface is not kept constant, a Neumann boundary condition can be chosen for higher accuracy. Here, the boundary is defined by the effective heat transfer coefficients k e,lat and k e,end , calculated for the thermal resistances at the surface.

Displacement distribution
Besides the temperature, the displacement due to thermal expansion also plays an important role in laser performance, since some contributions to thermal lensing and birefringence depend on the stress distribution in the crystal. With the temperature distribution calculated in a three-dimensional mesh, the displacement distribution can be calculated numerically by solving the respective differential equations. However, while the temperature can be calculated by a single parabolic PDE, the displacement distribution in this model is described by a set of coupled second-order hyperbolic PDEs, which adds more complexity to the respective model. The solution is based on Hooke's generalized law in the three spatial dimensions given by the following set of equations [15] wherein xx , yy , zz are the normal stresses and xy , xz , yz are the shear stresses between the discretized mesh grid points as shown in Fig. 5. The displacements in the three spatial directions are given by u, v, w. E describes the elastic modulus, the Poisson's ratio of the material. The initial thermal expansion is included by the thermal expansion factor and the temperature T . By equating all the volume forces in the system and solving the resulting set of differential equations, the displacement is calculated. In this model, the set of equations is yet again solved by the ADI FDM, which is also suitable for this problem. However, to model more complex crystal shapes, a FEM algorithm should be considered. A control volume-based FDM approach using a staggered grid is detailed in [16]. As a simple solution for the boundary conditions, the degrees of freedom at a given surface can be set for the displacements in x,y, and z. For a cylindrical crystal rod, the boundary at the lateral surface DoF lat and the endfaces DoF end are defined to fit the realistic mounting condition of the crystal. With the resulting displacement distribution, the relevant stresses in the crystal can be calculated analytically with equations (10-15).

Thermal lensing
Due to the heating of the laser crystal and resulting thermal effects, the crystal in a solid-state resonator typically displays lens-like behavior as depicted in Fig. 6. This is usually unwanted, limiting the power scalability of the system. Minimizing thermal lensing is, therefore, an important part of resonator design. In this model, thermal lensing is implemented by including the thermooptic effect, the photoelastic effect, and endface bulging [7].
Typically, the greatest impact to the thermal lens can be expected from the thermooptic effect, which describes the temperature-dependent change of the refractive index quantified with the material constant n T . By taking the highly resolved temperature distribution calculated prior, the change in refractive index is implemented with [17]  where n 0 is the initial wavelength-dependent refractive index, T o is the initial crystal temperature, and T(x, y, z) describes the temperature distribution due to the heat load in the crystal. Depending on the sign of the constant n T , this effect can contribute to a focusing or defocusing in the crystal, for Ho 3+ :YAG a focusing effect can be expected. Another contribution stems from the photoelastic effect, which describes the change in refractive index with the mechanical stresses due to thermal expansion. This effect is tied to the stress induced birefringence discussed in the next section. The last contributing factor in this model comes from the endface bulging of the crystal, which causes a lensing effect due to the deformation caused by thermal expansion at the endfaces. This can be implemented by adjusting the crystal length with the highly resolved displacement distribution calculated prior and using the intrinsic refraction functionality of the BPM.

Stress-induced birefringence
Birefringence describes the optical property of a material, where the refractive index depends on the polarization direction incident on the material. Depending on the crystal structure, anisotropic materials can be inherently birefringent. Ho 3+ :YAG is an isotropic material, which typically does not exhibit this behavior. However, mechanical stresses in the crystal due to thermal expansion can deform the local optical indicatrix as shown in Fig. 7. This changes the refractive properties of the crystal to act birefringent, depending on the specific temperature profile. For end-pumped laser crystals, this will usually result in a radial stress profile, which acts birefringent towards an optical field incident along the optical crystal axis.
The implementation of stress-induced birefringence in this model is based on an approximate approach by Koechner and Rice [18]. The polarization dependent refractive index n ± = 1 √ B � ± that describes the resulting optical indicatrix in this approach, is expressed by the dielectric impermeability B.
The specific impermeability values used in equation (17) are calculated with the stress data derived from the displacement distribution and the photoelastic tensor p ij , which is a material tensor dependent on the specific crystal structure.
The birefringent behavior in YAG is investigated in greater detail by Ostermeyer et al. in [19]. Unlike the slab geometry described in their work, cylindrical Ho 3+ :YAG rods are typically cut in the < 111 > crystal axis at a cut angle of cut = relative to the [101] direction. This has to be considered when choosing the photoelastic tensor values for equation (18).

Experimental setup
The combination of the aforementioned algorithms results in a powerful model for simulating cw Ho 3+ :YAG resonators. To validate the model, a resonator is set up in the laboratory and the experimental results are compared to the model results. Figure 8 shows a schematic depiction of the resonator used. The active medium is pumped with a Tm 3+ -fiber laser operating at around 1908nm, with a maximum cw output power of 100W and a M 2 < 1.1 . A two-lens telescope is used to focus the pump beam in the crystal at a D4 waist diameter of around 1000μm . The diameter of was measured at the crystal position without thermal effects and is expected to change in operation due to thermal lensing. The resonator itself is Z-shaped, where the pump beam passes through two dichroic mirrors at a 45 • angle with high reflection at 2090nm and high transmission at 1908nm. This Z-shaped cavity was chosen to keep the resonator simple by utilizing a single pass pump design. The high reflector (HR, Mirror 1) has a concave surface with a curvature of 100cm. The plane outcoupling mirror (OC, Mirror 2) has a reflectivity of 50% at 2090nm. A 1.1 at.%-doped Ho 3+ :YAG crystal rod with a diameter of 4mm and a length of 18mm acts as the active medium of the system. Both end faces are anti-reflection coated for the pump and laser wavelength. For cooling, the crystal was mounted into a copper heat sink to dissipate the heat load. The crystal is placed at half the resonator length, while the whole resonator from HR to OC has a length of 29cm. To produce a linear polarized output, a Brewster's angle thin film polarizer is placed in the resonator in front of the OC. This basic setup allows for a good fundamental comparison between model and experiment.

Comparison between model and experiment
For many thermal and spectroscopic material properties, different values can be found in literature. The parameters used for the validation were therefore chosen in such a manner as to assure agreement with the literature and different experiments, like the measured absorption through the crystal at varying pump powers. The specific spectroscopic and thermal material parameters of the Ho 3+ :YAG crystal  Fig. 9a. The comparison shows a similar slope of 44.2% for the experiment and 42.9% for the model. It is also notable, that at around 60W of pump power, the output power drops significantly. This is most likely due to the thermal lens in the resonator getting too strong and higher order modes getting more dominant, and it is also very well replicated in the model. The residual pump power over the incident pump power is plotted in Fig. 9b.
The measured data as well as the simulation show an irregular pump absorption, a rise in the residual pump power is notable at around 60W of pump power. We believe that this is due to the same reason as the drop in output power. Since the mode overlap between pump and laser mode likely changes due to thermal effects and the stronger higher order mode content, this will have a direct impact on the pump absorption efficiency.
Another important validation can be made by comparing the output field distributions at certain pump powers, as shown in Fig. 10. Here, the upper row shows the simulated output fields, while the lower row shows the output fields generated in the laboratory. The field distributions and D4 diameters at the centroid were measured at a distance of 35cm behind the OC mirror. The evolution of the fields over the output power shows a clear degradation of the fundamental mode into higher order modes, with an especially strong degradation between 50W and 65W, where the drop in output power takes place. The simulated fields show a very similar evolution, further validating the accuracy of the model. One thing to note is that the Z-shaped Ho 3+ :YAG resonator used in this experiment was not optimized for a gaussian beam output. This approach was deliberately chosen to test the model at suboptimal beam qualities, where a resonator reaches its stability limit.
A further comparison between model and experiment can be made by taking a look at the depolarization field ejected at the polarizer in the resonator, where the power losses due to depolarization take place. Figure 11 shows the depolarization field measured in the experiment (a) and in the model (b). They both show a four-lobe intensity distribution, as expected from an end-pumped cylindrical laser crystal with isotropic properties. While the general fields look similar, there are some differences in the specific intensity distribution on the lobes. Some accuracy issues are to be expected here, since the birefringence implementation is an approximation and some material parameters may not be chosen correctly. Nonetheless, this comparison shows that in general, the birefringence algorithm does work as intended.

Conclusion
We developed a powerful model for simulating end-pumped Ho 3+ :YAG resonators in continuous-wave operation. The fundamental model is based on a modular approach, where the propagation between different element operators is achieved with the split-step BPM field propagation algorithm. Thermal effects in the active laser crystal are included by a combination of different algorithms. The heat load in the crystal is calculated with the rate equations considering the population dynamics of holmium. The resulting temperature and displacement distributions are included by solving the respective differential equations with an FDM approach. These distributions allow for the inclusion of the most important thermal effects like thermal lensing and birefringence. The model is then compared to an experimental Ho 3+ :YAG linear cavity, where the model shows a good agreement in the evolution of output power and the resulting field distribution.
To achieve a high simulation accuracy, finding the right set of material parameters for the active laser medium can be a challenge, since literature values can vary between publications and may be incomplete. However, if the parameters are chosen in a suitable range, the model can be an powerful tool for simulating resonators with high precision.

Appendix A
See below Tables 1 and 2.