## Abstract

We derive a new model and simulation technique called “Dynamic Multimode Analysis (DMA)” to simulate the 3-dimensional dynamic behavior of a laser. A Gaussian mode analysis is used to obtain resonator eigenmodes taking into account thermal aberrations. These modes are coupled by a set of rate equations to describe the dynamic behavior of the individual modes for cw and Q-switched lasers. Our approach analyzes mode competition and provides a detailed description of the laser beam in terms of output power, beam quality factor M^{2}, and pulse shape. Comparison of experimental data with our simulation results provides new insight into the effect of mode competition and the operation of Q-switched lasers.

©2009 Optical Society of America

## 1. Introduction

The performance of a laser depends on the complicated interaction of a number of physical effects such as dynamic behavior of population inversion, thermal lensing, and mode competition. For successful resonator design, accurate simulations which provide a quantitative understanding of these effects are important. So far, mainly two methods are used to model solid state laser resonators: the Gaussian mode approach [1,2] in combination with scalar rate equations [2,3], and the beam propagation method (BPM) based on the Fox and Li principle [4].

The Gaussian mode approach is most common, but in its present use it is not able to predict mode superposition and competition. The latter means that different transverse modes simultaneously overlap with the gain distribution in the active medium and are amplified depending on their relative intensity and the resulting gain saturation (see Chapter 25.4 in [2]). Therefore, for approaches which do not include this effect, it is necessary to know the beam structure in advance to solve the rate equations. In less complicated cases, such as fundamental mode operation, such approaches have been used to predict pulse duration and peak power of a Q-switched laser [2,5,6] including passive Q-switching and intra-cavity frequency conversion [7]. However, a detailed time-dependent analysis of pulse shape and beam structure could not be carried through.

The BPM approach basically is a time-independent approach. Since just one wave front is propagated forth and back in the cavity, the full 3D interaction of the electrical field with the laser active medium cannot be analyzed dynamically. Attempts to propagate several wave fronts simultaneously are computationally elaborate.

In this paper we present a new approach using a “Dynamic Multimode Analysis” to compute mode superposition and competition, laser output power, beam quality, and pulse shape. For this purpose the common set of rate equations is extended to a set of multimode rate equations, which involves specific rate equations for the individual transverse modes being excited in the cavity. Using this extended set of rate equations, we describe the optical wave as a dynamic superposition of several eigenmodes. For the computation of these eigenmodes we considered two options, both of them taking into account thermal effects by a thermal finite element analysis (FEA) using absorbed pump power density and cooling configuration as input.

The most accurate approach to compute the eigenmodes is a FEA of the electromagnetic field in the cavity. In [8,9] it was explained how such a FEA can be accomplished. However, this approach is connected with a high numerical effort. Therefore, to reduce the computational time, we used a Gaussian mode approach in combination with a transverse parabolic fit applied to the 3D results of the thermal FEA. The results, presented in Section 4, show that this already leads to reliable results for end pumped laser systems. A set of multimode rate equations was first presented in [10] where the required assumptions were discussed on a more mathematical level. The current work is a further development of this idea applied to real laser cavities, including thermal lensing and Q-switching. All computations have been implemented and performed in LASCAD [11].

## 2. Model

#### 2.1 Definitions

Let us assume that $\Omega ={\Omega}_{2D}\times [0,{l}_{res}]$ is the domain of a laser resonator with length ${l}_{res}$(see Fig. 1 ). Furthermore, let ${\Omega}_{a}$ be the domain of the active medium, i.e. the laser crystal. The length and the refractive index of the active medium are denoted by ${l}_{cryst}$ and ${n}_{cryst}$.

In case of a Q-switched laser, let ${\Omega}_{q}$ be the domain of the Q-switch, ${l}_{q}$ its length, and ${n}_{q}$ its refraction index. To describe our model, we use the following abbreviations:

• speed of light in vacuum: c,

• optical path length of the resonator: ${l}_{opt}=({n}_{cryst}-1){l}_{cryst}+({n}_{q}-1){l}_{q}+{l}_{cryst}$,

• round trip time: ${\tau}_{r}=\genfrac{}{}{0.1ex}{}{2{l}_{opt}}{c}$,

• reflectivity of the output mirror: ${r}_{r}$,

• output reflection rate: $R=-\mathrm{log}({r}_{r})$,

• round trip loss: ${L}_{rtrip}$,

• (logarithmic) resonator loss: ${L}_{Res}={L}_{rtrip}+R$,

• stimulated emission cross section: *σ*,

• resonator life time: ${\tau}_{c}=\genfrac{}{}{0.1ex}{}{{\tau}_{r}}{{L}_{Res}}$,

• fluorescence decay time of the upper laser level: ${\tau}_{f}$,

• doping density of ions which are responsible for the laser process: ${N}_{tot}$,

• space-dependent pumping rate: ${R}_{p}(\overrightarrow{x})$.

Note that in case of Q-switching, ${\tau}_{c}$ is a time-dependent function. Furthermore, ${\tau}_{c}$ can be used to describe different kind of losses in the resonator, like losses caused by apertures or Gaussian output mirrors. This will be explained in a subsequent paper.

#### 2.2 Space-dependent rate equations

### 2.2.1 Rate equations for general field distributions

In case of a 4-level laser system, the dynamic behavior of the laser usually is described by the following rate equations for the population inversion density $N(t,\overrightarrow{x})$ and the total photon number in the cavity $\mathrm{\Phi}(t)$ [3]:

The interaction between the optical wave and the active laser material is determined by the stimulated emission $W(t,\overrightarrow{x})$ given by

The rate Eqs. (1) and (2) provide a comprehensive description of laser dynamics if the energy density $\rho (t,\overrightarrow{x})$ is known. However, it is not possible to compute $\rho (t,\overrightarrow{x})$ by these equations.

In this paper, we propose to solve this problem by a “Dynamic Multimode Analysis” assuming that the electric field consists of an incoherent superposition of transverse resonator eigenmodes. For this approach, we cannot use so-called “cold resonator modes” [12], but we must include aberrations due to temperature dependent modification of the refractive index and deformation of the active crystal. These effects are known as thermal lensing, and have a major influence on resonator stability and beam quality power [13,14].

In the next section we describe how the thermally induced phase distortions are taken into account by a Gaussian mode approximation.

### 2.2.2 Computation of resonator eigenmodes

The resonator eigenmodes needed for the “Dynamic Multimode Analysis” can be computed by different methods. As already mentioned in the introduction, the most accurate, but numerically elaborate, method is a finite element analysis (FEA) [8,9]. This FEA approach can take into account aberrations in a general way using the full 3D information about the thermally induced distortions of the refractive index and the deformation of the crystal. This information can be computed by a finite element analysis of the thermal load and the resulting deformation of the crystal.

However, for this paper we used a much more straightforward method as described in the LASCAD program manual [11]. This method uses a FEA for the computation of the thermal distortion in combination with a parabolic fit. Such a fit is necessary because an ABCD analysis cannot be applied to general refractive index distributions and interfaces. Therefore, we approximate the crystal by a set of short Gaussian ducts together with curved interfaces at both ends of the crystal. To this end, the crystal is subdivided into short sections along the axis, as shown in Fig. 2 , and every section is considered to be a Gaussian duct (see Chapter 20.3 in [2]).

The elements of the ABCD matrices describing each short duct are computed by the use of a parabolic fit for the thermally induced refractive index distribution. Similarly, ABCD matrices are computed for the deformed end faces of the crystal. In this way it is possible to compute an ABCD matrix for the full round trip which takes into account thermal aberrations in a good approximation [8].

### 2.2.3 Accounting for anti-reflection coatings in the computation of the energy density

The ABCD matrix algorithm delivers the transverse shape of the modes. In order to use them in the rate equations, it is important to consider how the propagating intensity changes at interfaces between different media. Since in laser resonators those interfaces usually are anti-reflective coated, the Poynting vector passes unchanged. Accounting for the different speed of light in the media, this delivers at the interface position ${z}_{I}$

where ${\rho}_{cryst}$ is the energy density in the crystal, and ${\rho}_{out}$ denotes the energy density outside the crystal. Thus, by Eq. (5), the energy density is proportional to the refractive index $n(z)$ along the resonator axis. To take this fact into account, let $u(\overrightarrow{x},t)$ be the transversely normalized field shape as it is obtained by a standard Gaussian mode analysis. That means that the integral of $|u(\overrightarrow{x},t){|}^{2}$ over the transverse directions equals 1. Then, the energy $\rho (t,x)$ can be written asHere, Eq. (4) is used to compute the effective mode volume${V}_{eff}$ by

Note that the transverse normalization of the field shape $u(\overrightarrow{x},t)$leads to a value ${V}_{eff}$ independent of time.

Due to the anti-reflection coating, both Eq. (6) and Eq. (7) are proportional to the refractive index and not to its square value as in the case of an uncoated interface. Therefore, at an anti-reflective coated interface, the transverse electric field is discontinuous which must be taken into account in the derivation of rate equations. In the literature, sometimes this fact is not considered explicitly [3].

#### 2.3 Dynamic Multimode Analysis

Let us now derive the rate equations of the “Dynamic Multimode Analysis”. For this, we rewrite Eq. (2) using Eq. (3) and Eq. (6):

To approximate the optical wave in the cavity, a set of orthogonal resonator eigenmodes is coupled by space-dependent rate equations. Then, the main model assumptions of the DMA are as follows:

## Assumption 1

Each eigenmode oscillates independently. This has the following consequences:

1. Rate Eq. (8) independently holds for each eigenmode.

2. The optical wave is considered to be an incoherent superposition of the excited eigenmodes. Therefore, the overall energy distribution $\rho (\mathrm{t},\mathrm{x},\mathrm{y},\mathrm{z})$ is obtained as the sum of the energy distributions of the individual eigenmodes.

3. Integrating the above equation and using Eq. (4), we obtain that the overall number of photons is given as the sum of the photons oscillating in each eigenmode:

It was explained in [10] that Assumption 1 holds, at least in the sense of a spatial average or a temporal mean value. This guarantees correct prediction of the average output power in the cw case. Deterioration of laser beam quality due to degenerated mode frequencies and coherent mode superposition was discussed in [12] and [15]. According to [12], coherent resonances can be avoided in well designed lasers.

In addition to mode competition, in this paper, we are mainly interested in actively Q-switched lasers with pulse interval shorter than the thermal relaxation time or in cw lasers. For these lasers an analysis of the transient thermal lens shows that one can assume a steady-state thermal lens [13,16]. Therefore, we assume the following:

## Assumption 2

The Gaussian mode shapes ${u}_{i}(t,\overrightarrow{x})={u}_{i}(\overrightarrow{x})$ do not change in time. Therefore the energy density ${\rho}_{i}$ of the individual modes can be computed by

An immediate consequence of Assumption 2 is that any changes of the laser beam shape occur due to mode competition.

Using Eq. (9) and Eq. (11), the stimulated emission $W(t,\overrightarrow{x})$ defined by Eq. (3) can be rewritten as

The above assumptions yield the main equations of the DMA using *M*modes

*Φ*can be computed by Eq. (10).

#### 2.4 Numerical discretization

A time-dependent and 3-dimensional simulation of several Gaussian modes requires the numerical discretization of Eq. (13) and Eq. (14). To discretize these equations, we apply a finite volume discretization on a finite volume grid of cells ${c}_{I}$, where

Here ${m}_{x,y}$ and ${m}_{z}$ are the number of grid points in $xy-$direction and $z-$direction, respectively [17]. The finite volume discretization leads to the following set of ordinary differential equations:

#### 2.5 Application to active Q-switching

Active Q-switching is a common technique to generate short pulses of high energy [2,3,13]. The time scheme comprises two different periods which are repeated with a certain frequency. In period I, the active material is pumped to a very high level of population inversion. For this purpose, the lasing process is prevented by a high intra-cavity loss.

After the load period, the Q-switch becomes transparent in period II, i.e., the intra-cavity loss is reduced significantly, and the high population of the upper laser level is consumed by a short and highly energetic laser pulse. For numerical reasons it is useful to subdivide period II as shown in Fig. 3 . Period IIa contains the laser pulse, and time period IIb is a relaxation period without any noticeable output power. Whereas period IIa requires a high numerical time resolution, period IIb can be simulated using much larger time steps.

According to the explanations of the principle of Q-switching given above, we introduce a time-dependent Q-switch loss ${L}_{QS}(t)$ which results in a time-dependent resonator life-time${\tau}_{c}(t)$:

Equation (18) and Eq. (19) allow us to adapt our model for Q-switched lasers without further modifications of the rate equations. Now, let us explain the output quantities of our simulation, which are used to characterize the laser beam.#### 2.6 Computation of output power and beam quality

Besides the characterization of the pulse shape, output power and beam quality are the main information obtained by the “Dynamic Multimode Analysis”. According to [3], the output power ${P}_{out,i}$ of the *i*-th mode can be calculated using the photon number ${\mathrm{\Phi}}_{i}$ as follows:

The beam quality of a laser can be characterized by the beam quality factors ${M}_{x/y}^{2}$. Let ${p}_{x/y,i}$ be the order of th i-th Gaussian mode in x/y-direction. Then, according to [19] and using Eq. (22), the beam quality factors ${M}_{x/y}^{2}$ can be calculated by

Equation (24) can be used to obtain the average beam quality factors [2,19] over the simulation time T:

Note that in Eq. (25) we have to take a weighted mean value of the time-dependent beam quality factors to compute the time-averaged beam quality factors ${M}_{x/y}^{2}$ in agreement with ISO 11146 measurements [20]. Let us now describe our experimental setup.

## 3 Experimental setup

Experimental results were obtained using an actively Q-switched Nd:YAG laser at $1064nm$. The resonator of length ${l}_{res}=1074mm$contains a $0.3\%$ doped slab crystal of length ${l}_{cryst}=26mm$and an output mirror with a reflectivity of $95\%$. The laser crystal is end-pumped by two laser diodes at $805nm$ with a pump power between $7W$ and $14W$ each. The absorbed pump power distribution is shown in Fig. 4 . The laser is actively Q-switched with a frequency of $1kHz$ and is designed for pulses of several hundred ns and high pulse stability. With a measured beam quality of ${M}^{2}<1.1$, the Q-switched output beam is almost diffraction limited.

According to the pulse repetition frequency of $1kHz$, in our simulations a load period of $990\mu s$and a combined pulse and relaxation period of $1000ns$are used (see Fig. 3). For the material parameters, well-known values are used [21].

## 4. Results and discussion

In this section we compare the results of our simulations with experimental data. Figure 5 shows the output power of a laser pulse as a function of time. The maximum output power predicted in this way is helpful to avoid damaging the crystal and other optical elements in the cavity.

In Fig. 6 , the average output power, computed by the DMA according to Eq. (23), is compared with experimental data.

All output power measurements are time-averaged values of the laser in Q-switched mode. As shown in Fig. 6, for low pumping power our simulation overestimates the output power. The reason is that we used the same absorption coefficients at $805nm$ and 0.3% doping density [21] for all our computations. However, in the experiment the laser is optimized for high output power and the emitted wavelength of the pump diode decreases for low pump power. Due to the sharp absorption spectrum at $805nm$the effective absorption efficiency also decreases. In our simulation this behavior of the pump diodes is neglected. As a consequence, accurate simulation results for high output power and deviations for low output power are obtained as expected. Furthermore, the decreasing slope efficiency due to beginning instability of the resonator can be noticed in the simulation for $13.9W$ pump power per diode. For higher pump power, the resonator becomes unstable and does not contain any stable Gaussian resonator eigenmodes.

The dynamic evolution of the laser pulse in Fig. 5 provides information about pulse length and symmetry. The pulse length, depending on the pump power, is compared with experimental data in Fig. 7 .

In the range of high pump power for which the laser is designed, Fig. 7 shows that the simulation slightly overestimates the pulse length, but predicts the correct dependency on the pump power. However, for low pump power, we notice a strong increase in the measured pulse length. This can be again explained by the behavior of the pump diodes.

Table 1 shows the simulated output power of each mode in cw and Q-switched operation. Furthermore, the time-averaged beam quality factors ${M}_{x,y}^{2}$ computed by Eq. (25) are given. The dynamic behavior of output power and beam quality which is needed for Eq. (25) is shown in Fig. 8 for cw and in Fig. 9 for Q-switched operation. Note that for the computation of time-averaged output power and beam quality in cw mode, we only integrate over the time period when the laser has almost reached its steady-state behavior (see Fig. 8). Furthermore, note that we chose the same initial value for all considered modes. Therefore, Fig. 8, and Fig. 9 show a high ${M}^{2}$ value at the beginning. This, however, does not influence the simulation, because as soon as the output power increases, the beam shape is determined by mode competition.

According to Table 1, Fig. 8, and Fig. 9, the simulated beam quality factor ${M}^{2}$ is about 1.7 for cw and 1.0 for Q-switched operation. Thus, the predicted single-mode operation of the laser in Q-switched operation agrees with the laser specification of ${M}^{2}<1.1$ (see Section 3).

Similar results for output power, pulse energy, and pulse shape can be obtained by an analytic model [2,22] or the solution of rate equations with a fixed mode shape [5–7,23]. However, DMA additionally computes the beam profile and the beam quality. This is not possible by common rate equation models. An interesting result is that DMA correctly predicts differences in mode competition for cw and Q-switched operation, which cause different beam qualities.

Let us analyze mode competition for cw and Q-switched operation in more detail. For this we use the dynamic 3-dimensional simulation results for beam profile and population inversion as shown in Fig. 10 and Fig. 11 .

In case of cw operation, Fig. 10(a)(Media 1) shows the dynamic behavior of the laser beam and the occurring mode competition. The evolution of the population inversion density is shown in Fig. 10(b)(Media 2). Observe that cw operation leads to a smooth transverse distribution of the population inversion. Since the fundamental mode does not overlap perfectly with the population inversion density, from time to time, higher-order modes are needed to consume the population inversion at the peripheral regions.

In case of Q-switched operation, Fig. 11(a)(Media 3) shows a short TEM00 pulse. Since the laser was pumped far beyond threshold, the strong laser pulse burns a spatial hole into the center of the population inversion, which cannot be filled up during the current Q-switch period (see Fig. 11(b)(Media 4)). High-order modes would be needed to consume the population inversion at the outer regions. However, oscillation of these modes is impeded by their overlap with the hole in the center of the population inversion. The first spikes of the transient oscillation in Fig. 8 as well as the laser pulse in Fig. 9 show this effect. In cw operation, Fig. 8 shows that these spikes are followed by several high-order mode spikes before the laser finally reaches its steady state. In Q-switched operation, Fig. 9 shows an increasing beam quality factor ${M}^{2}$ at the end of the pulse period. However, before any high-order modes can gain significant power, the current pulse period ends and the cavity loss is increased by the Q-switch again. This leads to the different values for the beam quality given in Table 1.

Figure 12 shows the dynamic behavior of the population inversion. Here, the DMA provides useful information about the generated population inversion during the load period. In our case, it can be seen that the upper laser level already begins to saturate before the Q-switch is opened.

On the one hand, this leads to a decrease of the output power due to increased losses by fluorescent decay. On the other hand, pumping up to saturation guarantees an almost identical population inversion distribution after each load period. This promises excellent pulse stability in agreement with our experimental data. Thus, the time-dependent analysis of the DMA can help to find a good trade-off between output power and pulse stability.

## 5. Conclusion

We presented a “Dynamic Multimode Analysis” which provides a comprehensive 3-dimensional simulation tool for solid state laser cavities. Using an eigenmode analysis that takes into account thermal aberrations, we could derive an appropriate space-dependent set of rate equations to compute shape and dynamic behavior of the laser beam. Comparison with experimental data shows that the DMA correctly predicts resonator stability, power output, pulse shape, and beam quality for cw lasers as well as for actively Q-switched lasers. It could be shown that, due to mode competition, the beam quality is improved by Q-switching. Thus, the “Dynamic Multimode Analysis” turns out to be a flexible and powerful new simulation tool to improve the design of solid-state laser cavities.

The combined analysis of different spatial and dynamic effects applied in this paper to model Q-switch operation and beam quality allows various generalizations. Modelling of apertures, Gaussian output couplers, passive Q-switches, and intra-cavity frequency conversion by a “Dynamic Multimode Analysis” will be discussed in subsequent work.

## Acknowledgments

The first author gratefully acknowledges funding of the Erlangen Graduate School in Advanced Optical Technologies (SAOT) by the German National Science Foundation (DFG) in the framework of the excellence initiative. The work presented in this paper was supported by the InnoNet project SOL of the Federal Ministry of Economics and Technology in Germany.

## References and links

**1. **P. A. Bélanger, “Beam propagation and the ABCD ray matrices,” Opt. Lett. **16**(4), 196–198 (1991). [CrossRef] [PubMed]

**2. **A. E. Siegman, *Lasers*, (University Science Books, Mill Valley, 1986).

**3. **O. Svelto, *Principles of lasers* (Springer, New York, 1998).

**4. **A. G. Fox and T. Li, “Resonant modes in a maser interferometer,” Bell Syst. Tech. J. **40**, 453–458 (1961).

**5. **X. Zhang, S. Zhao, Q. Wang, B. Ozygus, and H. Weber, “Modeling of diode-pumped actively Q-switched lasers,” IEEE J. Quantum Electron. **35**(12), 1912–1918 (1999). [CrossRef]

**6. **O. A. Louchev, Y. Urata, and S. Wada, “Numerical simulation and optimization of Q-switched 2 mum Tm,Ho:YLF laser,” Opt. Express **15**(7), 3940–3947 (2007). [CrossRef] [PubMed]

**7. **G. Li, S. Zhao, H. Zhao, K. Yang, and S. Ding, “Rate equations and solutions of a laser-diode end-pumped passively Q-switched intracavity doubling laser by taking into account intracavity laser spatial distribution,” Opt. Commun. **234**(1-6), 321–328 (2004). [CrossRef]

**8. **K. Altmann, C. Pflaum, and D. Seider, “Three-dimensional finite element computation of laser cavity eigenmodes [corrected],” Appl. Opt. **43**(9), 1892–1901 (2004). [CrossRef] [PubMed]

**9. **B. Heubeck, C. Pflaum, and G. Steinle, “New finite elements for large-scale simulation of optical Waves,” SIAM J. Sci. Comput. **31**(2), 1063–1081 (2009). [CrossRef]

**10. **B. Heubeck and C. Pflaum, “Numerical simulation of multiple modes in solid state lasers,” Proc. SPIE **6190**, 61900X (2006). [CrossRef]

**11. **LAS-CAD GmbH Munich, Germany, http://www.las-cad.com/.

**12. **R. Paschotta, “Beam quality deterioration of lasers caused by intracavity beam distortions,” Opt. Express **14**(13), 6069–6074 (2006). [CrossRef] [PubMed]

**13. **W. Koechner, *Solid-state laser engineering*, (Springer, Berlin, 2006).

**14. **T. Y. Fan, “Heat generation in Nd:YAG and Yb:YAG,” IEEE J. Quantum Electron. **29**(6), 1457–1459 (1993). [CrossRef]

**15. **Q. Zhang, B. Ozygus, and H. Weber, “Degeneration effects in laser cavities,” Eur. Phys. J. Appl. Phys. **6**(3), 293–298 (1999). [CrossRef]

**16. **S. Wang, X. Wang, T. Riesbeck, and H. J. Eichler, “Thermal lensing effects in pulsed end pumped Nd lasers at 940 nm,” Proc. SPIE **7194**, 71940J (2009). [CrossRef]

**17. **R. J. LeVeque, *Finite volume methods for hyperbolic problems* (Cambridge Texts in Applied Mathematics, 2007).

**18. **E. Hairer, and G. Wanner, *Solving ordinary differential Eqs. (2): Stiff and differential-algebraic problems*, (Springer, Berlin, 1996).

**19. **A. E. Siegman and S. W. Townsend, “Output beam propagation and beam quality from a multimode stable-cavity laser,” IEEE J. Quantum Electron. **29**(4), 1212–1217 (1993). [CrossRef]

**20. **International Organization for Standardization EN ISO 11146:2000 and EN ISO 11146:1999,” Laser and laser-related equipment. Test methods for laser beam parameters. Beam width, divergence angle and beam propagation factor,” (ISO, Geneva, 2000).

**21. **F. Träger, *Springer Handbook of Lasers and Optics* (Springer, Berlin, 2007), Chap. 11.2.

**22. **J. J. Degnan, “Theory of the Optimally Coupled Q-Switched Laser,” IEEE J. Quantum Electron. **25**(2), 214–220 (1989). [CrossRef]

**23. **R. B. Kay and D. Poulios, “Q-Switched Rate Equations for Diode Side-Pumped Slab and Zigzag Slab Lasers Including Gaussian Beam Shapes,” IEEE J. Quantum Electron. **41**(10), 1278–1284 (2005). [CrossRef]