1. Introduction

Bulk solid-state Master Oscillator Power Amplifier (MOPA) systems have, for several decades, been of considerable interest to power scale low-power laser beams from a master oscillator to higher powers [1]. In such a MOPA configuration, the amplifier laser gain medium is excited using a pump source that is typically a narrow-wavelength diode-laser (single or multimode) in either a side- or end-pumping configuration. In a side-pumped configuration, pump light from horizontally stacked diode laser bars is directed perpendicular to the optical axis while the seed beam (output from the oscillator), which is to be amplified, is propagated along the optical axis [2,3]. Conversely, in end-pumped systems, the input pump source, which is typically a multimode fiber-coupled (MMFC) diode-laser, impinges a laser gain medium in parallel to the optical axis to ensure the seed beam is co-axially propagated with the pump beam for optimal spatial overlap. In applications that employ high-brightness laser beams, which require both the output laser beam power and its beam quality to be maximized, both design configurations may be implemented, however, a prominent advantage of end-pumped systems is that it may be applied to a variety of laser gain media geometries such as thin-disks [4], cylindrical rods [5], fiber [6], and slab architectures [7]. To ensure that both the output power and beam quality are maximized, a stringent modeling approach is required to optimize the design strategy. For this purpose, we consider cylindrical rod geometries as they prove to be cost-effective, controllable, robust, easy to simulate, and are still widely used as pre-amplifier stages for the low to intermediate power scaling of low-power seed beams [8].

Modeling of MOPA systems using contemporary finite element methods (FEM) has been well demonstrated in literature [9–11]. However, there are still advantages to using an analytical modeling approach, such as short computational times, easy-to-interpret data, and detailed control over the physical input variables. There have been several seminal papers that advanced the analytical modeling of MOPA systems, particularly end-pumped systems. These include the simulation of Laguerre-Gaussian beam propagation in Raman, dye, and free-electron lasers [12], consideration of the transverse intensity distribution of the pump and seed beam during amplification [13], and derivation of an analytical expression for the temperature distribution of the pump beam [14]. The temperature distribution inside the gain medium induces both mechanical thermal expansion and thermo-optical variations, which can aberrate the seed beam resulting in beam distortions [15]. In recent years, highly accurate analytical expressions have been derived for end-pumped isotropic materials, specifically YAG. These analytical models incorporate steady-state rate equations [16], temperature distributions of various pump beam profiles [17], and the effects of stress-strain on seed beam distortions [18], employing the beam propagation method [19,20]. However, in these modern approaches, the treatment of the multimode fiber-coupled diode pump beam in 3D still poses a non-trivial modeling problem, often resulting in a reduced 2D static approximation. End-pumped systems typically utilize a fiber-coupled diode laser as the pump source, which is focused onto the crystal face or at a specific position along the crystal length. In the case of multimode sources, the 2D intensity profile of the light output undergoes a rapid transformation from a circular flat-top (FT) shape at the fiber exit face to a Gaussian or bell-shaped intensity profile along the propagation axis. Since the pump beam plays a fundamental role in a MOPA system and greatly influences the overall amplifier performance, it is essential to characterize not only the spatial profile of the pump beam at its highest fluence but also its dynamic propagation behavior throughout the volume of the gain medium. Current 2D or 3D end-pumped models consider a pump beam output that is static in its spatial intensity profile and is approximated as being either a uniform beam or top-hat [20–23], Gaussian beam [9,14,24,25], super-Gaussian beam [19,26–29] or an arbitrary polynomial distribution [30]. Another crucial aspect that requires attention in modeling is the interplay between the spectral emission characteristics of the pump light source and the absorption properties of the gain medium. It is well-known that the emission wavelength of a diode pump laser varies with the applied driving current due to temperature changes. However, most current amplifier models oversimplify this behavior by considering a single value for the pump absorption cross-section, typically corresponding to the maximum absorption wavelength of the gain material. This oversimplification fails to account for the non-static nature and non-linear variation of the pump absorption cross-section across different current levels. Incorrectly selecting this parameter can result in significant errors in determining the temperature distribution and overall accuracy of the amplifier performance.

In this paper, we present a novel approach for analytically modeling end-pumped Nd:YAG rod amplifiers in three dimensions, utilizing a multimode fiber-coupled diode laser as the pump source. We begin by describing the experimental setup used to calibrate and validate our 3D model. In the subsequent theory section, we discuss and mathematically formulate the phase-only transmission function for achieving Gaussian-to-FT beam shaping. Additionally, we highlight the importance of employing a weighted calculation for the absorption cross-section. Moving forward, we provide analytical expressions for the gain and absorption coefficients, the thermal gradient, and the resulting refractive index change specific to the 4-level [1 1 1]-cut Nd:YAG crystal. To integrate these expressions effectively, we develop an iterative Fourier beam propagation method that applies amplitude and phase modulations to the seed beam. In the results and discussion section, we showcase the correlation between the free-space propagation of the multimode fiber-coupled diode laser and the Gaussian-to-FT pump beams. Moreover, we compare the 3D temperature and refractive index gradient profiles obtained through the static super-Gaussian approximation and the dynamic Gaussian-to-FT methods, respectively. Finally, we evaluate the predicted output power of our unified 3D model against traditional modeling approximations and experimentally amplified Gaussian beams, demonstrating the accuracy of each approach.

2. Experimental setup

To validate and calibrate the analytical 3D model, a single-pass amplification experiment was conducted using a compact end-pumped Nd:YAG MOPA system in a contra-propagating pump and seed configuration, as depicted in Fig. 1. The physical properties of the Nd:YAG crystal used are summarized in Table 1. The end-pumping was achieved using a multimode fiber-coupled (MMFC) diode laser with a fiber core diameter of $\phi _{\mathrm {core}} = 400$ $\mathrm {\mu m}$, with an emission wavelength range of $\lambda _{0} = 803-808$ $\mathrm {nm}$, a spectral bandwidth of $\Delta \lambda = \pm 1.85$ $\mathrm {nm}$ and maximum continuous-wave output power of $P^{\mathrm {max}}_{\mathrm {p}} = 38$ $\mathrm {W}$. The spatial profile of the MMFC beam was analyzed using a CCD at positions $z$ = L, 3L/4, L/2, L/4, 0 along the crystal length, as indicated in Fig. 1. To achieve a uniform intensity flat-top (FT) profile with a full-width at half-maximum (FWHM) radius of $\omega ^{\mathrm {p}}_{f} = 200$ $\mathrm {\mu m}$ at the face of the Nd:YAG crystal ($z$ = 0), the output fiber face of the MMFC diode laser was imaged in a 1:1 fashion using a pair of $f_{2} = 100$ $\mathrm {mm}$ lenses in a 4f configuration. The pump beam radius (with a Gaussian-like intensity) upon exiting the crystal ($z$ = L) was measured to be $\omega ^{\mathrm {p}}_{i} = 1.7$ $\mathrm {mm}$. A $45^{\circ }$ dichroic mirror (D1) with a high-reflective and anti-reflective coatings at wavelengths of $1064$ $\mathrm {nm}$ and $808$ $\mathrm {nm}$, respectively. This allowed the Gaussian seed beam $\Psi _{\mathrm {seed}}$, with a wavelength of $\lambda {s} = 1064$ $\mathrm {nm}$, power of $P_{\mathrm {s}} = 100$ $\mathrm {mW}$, and beam quality of $\mathrm {M^{2}} = 1.04$, to be easily coupled out of the MOPA system after amplification. The seed beam was focused using an $f_1 = 200$ $\mathrm {mm}$ lens, resulting in a beam radius of $\omega _{s} = 181$ $\mathrm {\mu m}$ at the entrance face of the crystal ($z$ = L) and beam waist radius of $\omega ^{s}_{0} = 175$ $\mathrm {\mu m}$ at the pumped-face of the crystal ($z$ = 0). This produced an optimal pump and seed beam size overlap ratio of 0.8 at position $z$ = 0. The power of the MMFC pump and Gaussian beam was measured using a power meter (Gentec model: UP19K-30H-H5-D0).

 

Fig. 1. Schematic diagram of the single-pass end-pumped experimental configuration.

Download Full Size | PDF

Table 1. Summary of Nd:YAG properties: The parameters marked with superscripts ’a’ and ’b’ were obtained from Refs. [31] and [18], respectively, while parameter ’c’ was calculated using the method outlined in Ref. [17].

View Table

3. Theory

3.1 Pump beam modeling

The emission of multimode fiber-coupled (MMFC) diode pump lasers lacks a well-defined phase structure, resulting in an output beam that is incoherent and highly divergent. As such, simulation of its free-space propagation or propagation in complex media is a non-trivial task. The spatial profile of the MMFC pump beam, as it exits the fiber cable, can be approximated as flat-top (FT) shaped and is typically relay-imaged using lenses to the face of the gain medium, as shown in Fig. 1. However, as illustrated in Fig. 2, the FT beam diverges and transforms gradually into a Gaussian-like intensity profile as it propagates over the length of the gain medium. Through the stationary phase approximation, this behavior can be modeled by utilizing a lossless phase-only beam shaping technique which transforms a Gaussian beam to a FT beam, at the Fourier plane of a lens [32–34]. This approach enables the MMFC pump beam to be a complex field that mimics the divergence and spatial intensity transformation in free-space propagation. The phase-only transmission function $\Phi _{\mathrm {FT}}\left (r\right )$ and Fourier lens function T$_{\mathrm {lens}}(r)$ to perform a Gaussian to FT beam transformation is expressed as:

 

Fig. 2. Illustration of the spatial evolution of a typical MMFC pump beam (red) propagating through a crystal material (magenta cylinder), with co- and contra-propagating seed beams (green arrows). The two surfaces of the crystal are labeled A and B, respectively.

Download Full Size | PDF

3.1.1 Contra-propagation

The contra-propagation configuration consists of a convergent seed beam that is directed toward the surface B of the crystal adjacent to surface A where the pump beam enters the crystal. At surface B, the pump beam has a large Gaussian-like profile and as the seed decreases in size as it traverses the length of the crystal, it will arrive at its desired width on surface A where the pump beam is a tightly focused FT profile. To simulate this scenario, we create the phase transformation $\Phi _{\mathrm {FT}}\left (r\right )$, by specifying the desired FT beam radius $\omega _{f}$, Gaussian beam radius $\omega _{i}$, and the Fourier lens focal length $f$ to produce the transmittance of a single phase-element $\mathrm {DOE_{1}}$

For a Gaussian beam of radius $w_{i}$ and power $P_{\mathrm {p}}$, the corresponding beam will be modulated with the phase-transformation of $\mathrm {DOE^1}$

The modulated Gaussian beam $\Psi _{\mathrm {G}}(r)$ will transform gradually to a FT beam $\Psi _{\mathrm {FT}}$ of radius $w_{f}$ during free-space propagation over the distance of the Fourier lens $f$, as depicted by Fig. 3. It is important to note that in the simulation of the contra-propagating case, the propagation direction of the pump beam is reversed. As a result, the pump power $P_{\mathrm {p}}$ of the Gaussian beam $\Psi _{\mathrm {G}}(r)$ is determined by experimentally measuring the power after the beam has propagated through the crystal. Therefore, it is necessary to change the sign of the absorption coefficient in the simulation from an exponentially decreasing function to an exponentially increasing function, this will be discussed in more detail later.

 

Fig. 3. Initial Gaussian beam modulated by the phase of $\mathrm {DOE_{1}}$ that transforms into a FT beam $\Psi _{\mathrm {FT}}(r)$ at the Fourier plane of lens $f$, used for a contra-propagation simulation.

Download Full Size | PDF

3.1.2 Co-propagation

In the co-propagating configuration, the seed and pump beams propagate in the same direction, at surface A of Fig. 2. Here, the seed beam, at its desired beam width, is injected into the crystal where the pump beam is a tightly focused FT profile and exits the crystal at surface B where the pump beam has a large Gaussian-like shape. To simulate this configuration, the pump shaping process must be reversed from a Gaussian-FT to a FT-Gaussian transformation. To achieve this, the complex-conjugate phase of the FT beam $\Psi _{\mathrm {FT}}(r)$ at the plane of the Fourier lens $f$, shown on the right of Fig. 3, is extracted, and labeled as $\varphi _{\mathrm {FT}}(r)$. The complex-conjugate field $\varphi _{\mathrm {FT}}(r)$ is then used to create the reverse transformation transmission function, defined as $\mathrm {DOE_{2}}$.

The phase profile $\mathrm {DOE_{2}}$ is then applied to the FT beam $\Psi _{\mathrm {FT}}(r)$ twice, once to cancel out the existing phase and a second time to modulate the beam with the reverse phase transformation

Now, $\Psi _{\mathrm {FT^{*}}}(r)$ will transform from a FT beam back to the initial Gaussian beam $T_{\mathrm {G}}(r)$ during propagation over the length of the Fourier lens $f$, as depicted in Fig. 4.

 

Fig. 4. FT beam $\Psi _{\mathrm {FT}}(r)$ modulated by the phase transformation of $\mathrm {DOE_{2}}$ that reverses the transformation process back to the original Gaussian beam $\Psi _{\mathrm {G}}(r)$ at the Fourier plane of lens $f$, used for a co-propagation simulation.

Download Full Size | PDF

In the co-propagating scenario, where pump and seed beams enter the crystal at the same crystal surface, the pump power $P_{\mathrm {p}}$ of Eq. (4) should be equal to the total input pump beam power before it enters the crystal and the sign of the absorption coefficient should be exponentially decaying. Now, we have mathematically defined two transmission functions, namely $\mathrm {DOE_{1}}$ and $\mathrm {DOE_{2}}$, allowing for same direction propagation to simulate either co- or contra-propagating pump and seed fields during amplification. To achieve the highest level of simulation accuracy, it is necessary to obtain experimental measurements of the input parameters $\omega _{i}$, $\omega _{f}$, and $P_{\mathrm {p}}$.

3.2 Absorption cross-section

The rate at which pump photons are absorbed by the neodymium (Nd$^{3+}$) ions inside the yttrium aluminum garnet (YAG) host lattice is determined by the absorption cross-section $\sigma _{\mathrm {abs}}$ (i.e. the probability of an absorption event to take place). The $\sigma _{\mathrm {abs}}$ spectrum for 0.5$\%$ atm. doped Nd:YAG is depicted by the solid black curve in Fig. 5. The original data for $\sigma _{\mathrm {abs}}$ spectra in Fig. 5 was measured for a 1$\%$ atm. doped Nd:YAG sample and had a peak value of $10\times 10^{-20}$cm$^{2}$ at wavelength $808$ $\mathrm {nm}$ [35]. At low doping concentrations, $\sigma _{\mathrm {abs}}$ varies linearly with the dopant $\%$ and therefore, the original data can be scaled to represent a 0.5$\%$ doped Nd:YAG crystal [36]. The spectral line shape of a typical MMFC diode laser emission can be approximated as a standard Lorentzian distribution function, indicated by the red curve in Fig. 5, having a well-defined central wavelength $\lambda _{0}$ (vertical red dotted line) and bandwidth $\Delta \lambda$.

 

Fig. 5. Absorption cross-section of Nd:YAG crystal (black) and emission spectrum of the pump diode (red) with resulting weighted overlap integral (blue shaded area).

Download Full Size | PDF

Experimentally, the central wavelength $\lambda _{0}$ and spectra bandwidth $\Delta \lambda$ of the MMFC emission spectrum varies with the driving current due to a temperature change in the diodes. This is the typical behavior of diode lasers and manufacturers generally provide a value for the expected $\mathrm {nm/C^{\circ }}$ or $\mathrm {nm/A}$ central wavelength shift as a function of temperature or current. The spectral emission characteristics of most $808$ $\mathrm {nm}$ wavelength MMFC diodes vary between central wavelengths of $802$ $\mathrm {nm}$ at threshold current to $808$ $\mathrm {nm}$ at the maximum driving current with a spectral bandwidth of $<2$ $\mathrm {nm}$. Nevertheless, most amplification models use the peak absorption cross-section $\mathrm {\sigma ^{\mathrm {p^*}}_{\mathrm {abs}}} = 5\times 10^{-20}$ $\mathrm {cm^{2}}$ (vertical magenta line in Fig. 5) across all driving currents and do not consider the spectral overlap of the emission and absorption curves. This is a significant oversight that leads to an overestimation of the pump absorption and seed amplification behavior. Ideally, $\lambda _{0}$ and $\Delta \lambda$ of the diode emission spectrum should be measured experimentally, as a function of the pump power (i.e. driving current). With this information, an accurate $\sigma ^{\mathrm {p}}_{\mathrm {abs}}$ value can be calculated by performing an overlap integral between a Lorentzian function $f(\lambda )$, using the experimentally determined $\lambda _{0}$ and $\Delta \lambda$ values, and the absorption cross-section data $g(\lambda )$ of the Nd:YAG crystal

The Lorentzian function $f(\lambda )$ is normalized to the peak absorption value. As an example, the overlap between the two curves, depicted as the magenta shaded region in Fig. 5, is the weighted value of $\sigma ^{\mathrm {p}}_{\mathrm {abs}}=9.1\times 10^{-21}$ $\mathrm {cm^{2}}$ which is $5.5\times$ smaller than the peak value of $\sigma ^{\mathrm {p*}}_{\mathrm {abs}}$.

3.3 Gain and absorption coefficients

In a lossless 4-level Nd:YAG gain medium under steady-state continuous pumping conditions, the time-independent gain coefficient $G(r,z)$ that describes the exponential growth rate of a seed beam can be written in cylindrical coordinates as [16,29]

Analogous to the gain coefficient, which defines the increase of the seed intensity $I_{\mathrm {s}}$, the time-independent absorption coefficient $\alpha (r,z)$ describes the strength of the longitudinal exponential absorption of the pump intensity $I_{\mathrm {p}}$ by the crystal. Similarly, for a lossless 4-level Nd:YAG gain medium under steady-state continuous pumping conditions, the absorption coefficient may be written as

3.4 Temperature-gradient

The difference in energy between the pump absorption and laser emission energy levels contributes to the heat generated inside the gain medium. This is caused by an effect known as quantum defect heating $\eta _{\mathrm {h}}$, as only a fraction of pump photons that are absorbed reach the laser emission band and the rest are converted to heat through non-radiative processes ($\eta _{\mathrm {h}} = 0.24$-$0.3$ for $0.5\%$ atm. Nd:YAG). The analytical expressions of gain $G(r,z)$ and absorption $\alpha (r,z)$ coefficients, in steady-state, are calculable at any point during the Gaussian-FT transformation since they require only the intensity profile of the pump and seed beams. However, an analytical expression for the temperature distribution induced by the Gaussian to FT pump inside the gain medium is non-trivial. Recently, an improved analytical solution was derived for the temperature $T(r,z)$ produced when end-pumping a cubic crystal (such as Nd:YAG) with a generalized n$_{\mathrm {th}}$ order super-Gaussian pump beam $\mathrm {SG^{n}}(r,z)$ [17]. The boundary conditions in the derivation accounted for active cooling by a surrounding heat sink, schematically depicted in Fig. 6, and solved using the steady-state heat conductance equation in cylindrical coordinates with temperature-dependant thermal conductivity $k(T)$

Equation (12) can then be solved to obtain the analytical expression for the temperature distribution $T(r,z)$ induced by the generalized super-Gaussian pump beam $\mathrm {SG^{n}}(r,z)$ of order $\mathrm {n}$, using the temperature dependant thermal conductivity $k(T)$ of Eq. (11). Rather than use the absolute temperature $T(r,z)$, it is convenient to define the temperature gradient relative to the boundary interface of the crystal rod and heat-sink material $T(R,z)$

 

Fig. 6. Schematic diagram of an end-pumped (blue arrows) laser rod (magenta) of radius $R$ and length $L$ mounted inside a copper heatsink (brown) with red arrows indicating the radial direction of heat flow.

Download Full Size | PDF

Analytically solving Eq. (12) at each $N_{z}$ interval to obtain a unique $\Delta T(r,z)$ for the Gaussian to the FT pump transformation is not possible. Therefore, we use the generalized solution of Eq. (13) as is, and select a $\mathrm {SG^{n}}(r,z)$ approximation that best fits the spatial profile and radius of Gaussian to FT transformation at each $N_{z}$ position. This approximation is acceptable since the parabolic shape of $\Delta T (r,z)$ is similar for super-Gaussian beams of order $n =2, 4, 6$ & $8$, what is more important is the change in the pump beam size due to divergence which is accounted for at each $N_{z}$ interval.

3.5 Change in refractive-index

The thermal gradient $\Delta T(r,z)$ that develops inside the gain medium causes the material to experience stress-strain effects leading to mechanical bulging of the crystal end-faces defined by the thermal expansion coefficient $\mathrm {C}$. Additionally, thermo-optical changes occur in the material which causes variations to the refractive index $n$ defined by the coefficient $\beta$(i.e. $\frac {\mathrm {d} n }{\mathrm {d} T}$).

In most amplifier models, the change in refractive index is simply defined as $\Delta n(r) = \frac {\mathrm {d}n }{\mathrm {d}T} \Delta T(r)$, where $\frac {\mathrm {d}n }{\mathrm {d}T}$ is the thermal dispersion coefficient with values ranging between $7-10\times 10^{-6}$ $\mathrm {K^{-1}}$ for Nd:YAG at room temperature $\sim$ 300 K [31,37]. The change in refractive index produced by the quadratic to logarithmic variation in the thermal gradient $\Delta T(r,z)$ of Eq. (13), will simulate both thermal lensing and spherical aberration effects [29]. This simple method for calculating $\Delta n(r)$ is a reasonable approximation, depending on how well $\Delta T(r)$ is defined but neglects thermal expansion caused by stress and strain of the material.

A holistic $\Delta n(r)$ calculation for end-pumped amplifiers, especially when high-pump powers are involved, must include the effects of stress and strain and use temperature-dependant coefficients so that $\mathrm {C} \equiv \mathrm {C} (T)$ and $\beta \equiv \beta (T)$. Extensive research has been conducted on the topic of thermally-induced stress-strain effects in cubic crystals, particularly Nd:YAG. Detailed explanations on the use of temperature-dependent coefficients $\mathrm {C} (T)$ and $\beta (T)$ have been provided in several publications [18,31,38]. The mechanical properties of cubic crystals, such as [1 1 1] cut Nd:YAG, can be estimated as an isotropic solid media where only two elastic coefficients (Poisson’s ratio $\nu$ and Young’s modulus $E$) are needed for a description of induced stress-strain and the resulting change in refractive index $\Delta n(r)$

The terms $\beta _{\sum }(r)$ and $\mathrm {C}_{\sum }(r)$ are the temperature-dependant thermo-optic and linear-expansion coefficients, respectively, which have been integrated with respect to temperature. The coefficients $\mathrm {C}^{R}_{\sum }$ & $\mathrm {C}^{r}_{\sum }$ are obtained from an additional integration over the crystal radius $R$ and radial coordinate $r$, respectively. The above-mentioned integration steps to arrive at Eq. (14) are non-trivial, therefore, the reader is referred to the detailed discussion and derivation by D. Bričkus et. al. (2017) [18]. The coefficients $a_{1}$ and $a_{2}$ contain the Poisson’s ratio $\nu$, the initial reference refractive index $n_{0}$ and the $p_{ij}$ coefficients of the elasto-optic matrix, which are provided in Table 1. The (+) and (-) signs in front of the last term in Eq. (14) describe the radial and tangential components of the change in refractive index. For low pumping fluence ($\mathrm {W/cm^{2}}$), the thermal dependencies of $\beta$ and $\mathrm {C}$ can be neglected, however, for high-power pumped crystals the temperature-independent equations may underestimate the actual stresses by 1.5-1.8 times, which could result in crystal fracture if the thermal-fracture limit is exceeded ($235 \pm 16$ MPa for YAG).

3.6 Beam propagation method

Using the scalar angle-spectrum theory of diffraction, based on the Fast Fourier Transform (FFT) method, the complex fields of the pump $\Psi _{\mathrm {p}}(r,z)$ and seed $\Psi _{\mathrm {s}}(r,z)$ beams are propagated using a split-step procedure, whereby the crystal volume of length $L$ is divided into $N_{z}$ intervals each of length $\delta z = L/N_{z}$. Let $\Psi _{i}(r,z)$ be an arbitrary complex field that is an eigenmode of Maxwell’s equation, the propagation of $\Psi _{i}(r,z)$ over a $\delta z$ increment is as follows:

3.6.1 Initialisation of pump and seed fields

As discussed in section 2.1, the seed field $\Psi _{\mathrm {s}}(r,z)$ can be either co- or contra-propagating with respect to the pump field $\Psi _{\mathrm {p}}(r,z)$. For the sake of completeness, we describe the starting pump fields $\Psi _{\mathrm {p}}(r,z)$ for both cases, however, in our work, the experiment could only be performed in the contra-propagating configuration. To simulate the co-propagating case, the starting pump field is set as the modified FT beam of Eq. (6)

Once $\Psi _{\mathrm {p}}(r,z)$ has been correctly defined, the spatial transformation of the pump (i.e Gaussian to FT or FT to Gaussian) will occur during the iterative propagation procedure. In this work, the seed field $\Psi _{\mathrm {s}}(r,z)$ was a Gaussian beam of radius $\omega _{\mathrm {g}}$ with an initial power $P_{\mathrm {s}}$ which is defined as

3.6.2 Iterative procedure

The beam propagation method (BPM) of Eq. (15) is a tool that enables the seed and pump fields to interact in an iterative manner with the analytical expressions for the gain coefficient $G(r,z)$ of Eq. (8), absorption coefficient $\alpha (r,z)$ of Eq. (10), change in temperature gradient $\Delta T(r,z)$ of Eq. (13) and the change in refractive index $\Delta n(r)$ of Eq. (14), providing a complete solution for 3D modeling of Nd:YAG end-pumped amplifiers.

The BPM iterative procedure for the contra-propagation configuration can be explained with the aid of Fig. 7. The simulation was performed in MATLAB using a mesh grid of dimensions $N_{x,y} = 3000$ and $N_{z}=36$ so that each slice of the crystal was of width $\delta _{z} = 0.69$ $\mathrm {mm}$. Increasing $N_{z}>36$ did not improve on the model accuracy but only increased the computational time. The position $z=0$ marks the entrance face of the crystal, where the starting seed field $\Psi _{\mathrm {s}}(r,0)$ of Eq. (18) and pump field $\Psi _{\mathrm {p}}(r,0)$ of Eq. (17) are propagated midway into the first slice $N_{z=1}$ of the crystal (i.e. distance $\delta z/2$) using the propagation Eq. (15). At $z = \delta z/2$ (blue plane in Fig. 7), the pump $I_{\mathrm {p}}(r,\delta z/2)=|\Psi _{\mathrm {p}}(r,\delta z/2)|^{2}$ and seed $I_{\mathrm {s}}(r,\delta z/2)=|\Psi _{\mathrm {s}}(r,\delta z/2)|^{2}$ intensities are calculated and used to determine the absorption coefficient $\alpha (r,\delta z)$. The pump field $\Psi _{\mathrm {p}}(r,\delta z/2)$ is then modulated by the calculated absorption coefficient $\alpha (r,\delta z)$ as follows

The sign in front of the absorption coefficient is $+$ for contra- and $-$ for co-propagating beams, respectively. The resulting pump intensity after absorption $I_{\mathrm {p}}(r, \delta _{z}) = |\Psi _{\mathrm {p^*}}(r,\delta z/2)|^{2}$ is used in the calculation of the gain $G(r,\delta z)$ coefficient and temperature gradient $\Delta T(r,\delta z)$. When calculating $\Delta T(r,\delta z)$, it is important to substitute the value of $P_{\mathrm {p}}$ corresponding to the absorbed pump $I_{\mathrm {p}}(r,\delta z)$ and not the starting power $P_{\mathrm {p}}$. The seed field $\Psi _{\mathrm {s}}(r,\delta z/2)$ is then modulated in amplitude by the gain distribution $G(r,\delta z)$ and in phase by the change in refractive index $\Delta n(r)$ as follows

 

Fig. 7. The schematic diagram illustrates the split-step iterative procedure for the pump beam (red) and seed beam (green) in a contra-propagating configuration, where the pump beam travels from right to left and the seed beam travels from left to right.

Download Full Size | PDF

4. Results and discussion

4.1 MMFC diode vs FT-Gaussian spatial transformation comparison

The accuracy of modeling an in-line end-pumped amplifier relies on experimental verification steps, as the spatial behavior of MMFC diode lasers varies based on the numerical aperture and core diameter of the delivery fiber. The data points in Fig. 8 represent the experimentally measured MMFC pump beam intensity profiles (normalized to unity) at positions $z = 0$, $L/4$, $L/2$, $3L/4$, and $L$. The inset of the figure displays the 2D intensity profiles corresponding to $z = 0$, $L/4$, $3L/4$, and $L$ as $a$, $b$, $c$, and $d$, respectively. It is evident that there is a significant deviation in the intensity profile and radius of the pump beam between the entrance ($z=0$) and exit ($z=L$) faces of the crystal. To ensure the physical representation and accuracy of the FT to Gaussian beam shaping process, the experimentally measured MMFC pump beam radius at the entrance ($z=0$) and exit ($z=L$) faces of the crystal are used as inputs $\omega _{i}$ and $\omega _{f}$, respectively, in Eq. (4). The resulting FT to Gaussian transformation is plotted on the same axis (dotted lines) and shows an accurate correlation with the experimental MMFC data with spatial overlap percentages of $98\%$, $95\%$, $99\%$, $94\%$ and $96\%$, at positions $z = 0$, $L/4$, $L/2$, $3L/4$ and $L$, respectively. It is important to reiterate, that only the input ($z=0$) and output ($z=L$) MMFC beam sizes were used as inputs to Eq. (4), while the intermediate profiles ($0<z<L$) are produced naturally during Fourier propagation of the FT to Gaussian transformation. Upon closer inspection of the experimental data, we find that at $z=0$ the MMFC pump beam is slightly more rounded than the simulated FT to Gaussian curve. The fit of the simulated curve can be slightly optimized by varying the position of the crystal relative to the length of the Fourier lens $f$ in Eq. (1), to adjust the roundedness of the FT at $z=0$. We found that for a crystal of length $25$ $\mathrm {mm}$ a Fourier lens of length $f=25.5$ $\mathrm {mm}$ produces an FT beam that is slightly rounded and maintains the correct steepness.

 

Fig. 8. The figure illustrates a comparison between the FT-Gaussian beam shaping profiles (represented by dash-dot lines) and the measured MMFC diode output (indicated by markers) at various points along the free-space propagation corresponding to the crystal length, specifically at $z=0,$ $L/4$, $L/2$, $3L/4$, and $L$. The inset of the figure displays the experimental 2D intensity profiles of the MMFC pump beam measured at $z = 0$, $L/4$, $3L/4$, and $L$, labeled as $a$, $b$, $c$, and $d$, respectively.

Download Full Size | PDF

4.2 Absorption cross-section calibration

Model calibration involved measuring input and transmitted pump powers, central wavelength ($\lambda _{0}$), and spectral bandwidth ($\Delta \lambda$) at seven calibration points. The input power was simply measured experimentally (black asterisk) and set into the model (magenta), as illustrated in Fig. 9 c). The transmitted pump power after absorption through the crystal (black asterisk) is shown in Fig. 9 b) and illustrates non-linear behavior due to the shift in the pump emission central wavelength $\lambda _{0}$, shown in Fig. 5, which causes $\mathrm {\sigma ^{p}_{abs}}$ and $\alpha (r,z)$ in Eq. (10) to vary. To obtain the correct transmitted powers, adjustments to $f(\lambda )$ Lorentzian function used in Eq. (7), at each calibration point were made. This is done by experimentally measuring the central wavelength $\lambda _{0}$ (left axis) and spectra bandwidth $\Delta \lambda$ (right axis) used to generate $f(\lambda )$ as shown Fig. 9 a) at each of the 7 calibration points. The magenta points on a) represent the optimal values of $\Delta \lambda$ that produce the best fit (magenta circles) in b) for each of the calibration points. With the calibrated absorption cross-section and associated FT to Gaussian pump transformation, the intensity of the pump beam absorption along the length of the crystal is shown in Fig. 10.

 

Fig. 9. Plot a) shows the experimentally measured ($\pm$ 0.1 nm error bars) central wavelength $\lambda _{0}$ (left axis, black) and spectra bandwidth $\Delta \lambda$ (right axis, magenta), with the optimal values used in the simulation (asterisk). Plots b) and c) are experimentally measured (black asterisk) and simulated (magenta circles) values of the transmitted and input pump beam powers, respectively.

Download Full Size | PDF

 

Fig. 10. 2D simulation, using the calibrated absorption cross-section calculation method, of the absorption into the crystal for the Gaussian-FT evolving pump beam at maximum pump power $P_{\mathrm {p}}$= 38 W.

Download Full Size | PDF

4.3 SG of order ’n’ fitted to FT-Gaussian for $\Delta T(r,z)$ calculation

The analytical expression for the temperature gradient $\Delta T(r,z)$ is derived for a generalized super-Gaussian heat-load of order $\mathrm {SG^{n}}$ for a specified radius $\omega _{\mathrm {p}}$, absorption coefficient $\alpha (r,z)$ and pump power $P_{\mathrm {p}}$. Since it is not feasible to solve the heat-conduction differential equation analytically at each $N_{z}$ interval for the FT to Gaussian pump transformation, a suitable $\mathrm {SG^{n}}$ distribution was fitted to the FT to Gaussian transformation at each $N_{z}$ interval. Importantly, only the order "n" of the $\mathrm {SG^{n}}$ was varied to establish the best fit, while the radius $\omega _{\mathrm {p}}$, absorption coefficient $\alpha (r,z)$ and pump power $P_{\mathrm {p}}$ were extracted at each $N_{z}$ interval from the BPM procedure and substituted into Eq. (13) to obtain $\Delta T(r,z)$. Fig. 11 shows the normalized $\mathrm {SG^{n}}$ (dotted lines) distributions for the corresponding FT to Gaussian spatial transformation (solid lines) at various positions over the length of the crystal. At position $z=0$, a super-Gaussian $\mathrm {SG^{20}}$ matches the steepness of the pump beam distribution accurately, but slightly overestimates the plateau flatness, as seen by the purple dotted line in the enlarged inset in Fig. 11. In the literature, we find that $\mathrm {SG^{n}}$ pump beams having equal sizes but slightly different peak intensities and orders of ’n’, produce similar spatial distributions and peak values of $\Delta T(r,z)$ [17]. What is more important is the correlation of the steepness between the super-Gaussian estimation and the simulated pump beam, especially since this defines the area where the highest thermal gradient is experienced. Therefore, a slight overestimation of the pump beam flatness is traded to achieve an accurate estimation of the beam size at $z=0$. For positions $z>0$, the "best-fit" values of ’n’ vary rapidly over the first few $N_{z}$ intervals from $\mathrm {SG^{20}}$ at $N_{1}$, $\mathrm {SG^{6}}$ at $N_{2}$, $\mathrm {SG^{4}}$ from $N_{3-6}$ and $\mathrm {SG^{2}}$ for the remaining $N_{7-36}$ slices.

 

Fig. 11. Gaussian-FT pump beam (solid) compared to the "best fit" super-Gaussian SG$^{\mathrm {n}}$ (dotted) used in the calculation of $\Delta T(r,z)$.

Download Full Size | PDF

4.4 FT-Gaussian vs SG thermal gradient

The thermal gradient over the intervals $N_{z}$ is obtained by substitution of the corresponding best fit $\mathrm {SG^{n}}$ parameters: "n", $\omega _{\mathrm {p}}$, $\alpha (r,z)$ and $P_{\mathrm {p}}$ into Eq. (13). Fig. 12, shows the dynamic temperature gradient evolution inside the crystal as a function of $z$ under maximum input pump power of 38 W, using the calibrated absorption cross-section value. To illustrate the differences between the modeling approach presented here and that of a typical static and non-diverging $\mathrm {SG^{n}}$ approximation, we consider a $\mathrm {SG^{20}}$ pump beam at $z=0$. The 3D plots in Fig. 12 a) $\&$ d), have a similar form with the pumped end of the crystal at $z=0$ having the highest peak temperature gradient that dissipates as the pump beam is absorbed into the crystal. As illustrated, the thermal gradients of the two approaches decrease at different rates along the length of the crystal which is evident in the 2D and 1D plots of the respective images. This is primarily due to the non-diverging nature of the static $\mathrm {SG^{20}}$ pump beam approximation whereas, in the Gaussian-FT dynamic pump beam approximation, the temperature gradient rather exhibits a more natural dispersion over the crystal length. We can further quantify these differences by analyzing the 1D traces through the crystal, as shown in Fig. 12 c) $\&$ f). At $z=0$, the static and dynamic pump beams produced thermal profiles of the same shape and peak value of $\Delta T \sim 17.5$ $\mathrm {K}$, as expected, since they are both the same $\mathrm {SG^{n}}$ approximations at $z=0$. Additionally, we see can clearly see the quadratic radial dependence of the temperature gradient in the region of the pump beam and the transition to a logarithmic dependence so that $\Delta T(R,z) = 0$ at the crystal boundary. For $z>0$ the quadratic thermal gradient of the non-diverging super-Gaussian approximation is retained and remains centralized throughout the crystal compared to the dynamic approach which smooths out gradually. The peak $\Delta T$ values for the static pump beam exceed the dynamic pump beam by $1.4\times$, $1.8\times$, $2.2\times$ and $2.6\times$ at $z=L/4$, $L/2$, $3L/4$ and $L$, respectively.

 

Fig. 12. 3D (top row), 2D (middle row) and 1D (bottom row) visualization of the thermal gradient $\Delta T(r,z)$ inside the Nd:YAG crystal as a result of end-pumping by a non-diverging super-Gaussian beam (left column) versus the FT-Gaussian transformation pump beam (right column).

Download Full Size | PDF

4.5 Change in refractive index

The thermal gradient gives rise to a change in the refractive index $\Delta n (r)$ as shown in Fig. 13. Here, we compare the corresponding $\Delta n (r)$ produced by the static $\mathrm {SG^{20}}$ (dotted lines) and FT-Gaussian pump beams (solid lines), at various z-positions along the length of the crystal and it is clear that the $\Delta n(r)$ distributions follow the same spatial shape as that of the thermal gradient. $\Delta n(r)$ is calculated using Eq. (14) which includes the temperature-dependant linear expansion $\mathrm {C}(T)$ and thermo-optic $\beta (T)$ coefficients. The inclusion of the higher order expansion terms introduces non-parabolic variations to $\Delta n(r)$, however, these effects only become apparent when $\Delta T(r,z)$ is significantly large. In this case, due to the relatively low maximum pump power $P_{\mathrm {p}} = 38$ $\mathrm {W}$, the peak $\Delta T(r,z) = 17.5$ $\mathrm {K}$ is too small to induce higher-order thermal variations and therefore $\Delta n(r)$ is reduced to the simplified $\Delta n = \frac {dn}{dT} \Delta T$ expression, where the thermal gradient and change in refractive index are proportional. Under conditions of high-intensity pumping conditions or poor spatial overlap between the pump and seed beam sizes in the crystal, the effects of aberrations induced by the refractive index variation, $\Delta n(r)$, become more pronounced in the seed beam. In this work, the pump-to-seed beam size ratio was optimized to minimize aberrations and maximize power scaling. Figure 13 illustrates the parabolic shape of the $\Delta n(r)$ curve for $z\geq 0$, encompasses the region occupied by the seed beam radius ($\omega ^{s}_{0} = 0.175$ mm). As a result, the seed beam experiences only uniform parabolic lensing effects. If the seed beam were larger and extended beyond the parabolic region of $\Delta n(r)$, non-uniform spherical aberrations would be observed. These aberrations would be further magnified when using a static super-Gaussian (SG) approximation, as it significantly overestimates $\Delta n(r)$. Since the analytical expressions in this model are derived for a circularly symmetric system, non-circular phase aberrations such as astigmatism and coma will not be observed.

 

Fig. 13. Comparative plot showing the change in refractive index $\Delta n (r,z)$ caused by the non-diverging super-Gaussian beam of order n=20 (dotted lines) and the FT-Gaussian transformation pump beam (solid lines), respectively, at $z=0$, $L/4$, $L/2$, $3L/4$ & $L$ positions of the crystals.

Download Full Size | PDF

4.6 Amplification performance

The individual components presented here are all combined in the determination of the power scaling performance of the end-pumped amplifier. As a fair point of comparison, we consider 6 cases to model where 3 are based on a fixed value for the absorption cross-section $\mathrm {\sigma ^{p}_{abs}}$ and 3 consider the dynamic weighted absorption cross-section approach presented above to obtain $\mathrm {\sigma ^{p}_{abs}}$ across the pump power levels. For each subset, we consider a 2D modeling approach where only a single slice in the crystal, $N_{1}$, is evaluated. Here the pump beam is a super-Gaussian of order 20, $\mathrm {SG^{20}}$, which holds for the FT-Gaussian transformation at the first slice. The remaining cases consider modeling in 3D for the static $\mathrm {SG^{20}}$ and dynamic FT-Gaussian pump beams, respectively as shown in Fig. 14. To supplement the modeling results, the experimentally obtained results (blue asterisks) for the amplification of a Gaussian beam are plotted on the same curve. Here the Gaussian beam had an input power of $P_{\mathrm {s}} = 100$ $\mathrm {mW}$ and was amplified at pump powers $P_{\mathrm {p}} = 5.8$ $\mathrm {W}$, $12.4$ $\mathrm {W}$, $18.5$ $\mathrm {W}$, $25$ $\mathrm {W}$, $31.5$ $\mathrm {W}$ $\&$ $38$ $\mathrm {W}$. The error bars correspond to a $\pm 4\%$ error deviation to account for the power meter and master oscillator fluctuations. In Fig. 14, the insets on the right of the graph show the simulated and experimentally measured 2D intensity profiles of the Gaussian seed beam after the MOPA stage during maximum pumping power $P_{\mathrm {p}}=38$ W, and show that the spatial profile was well preserved - in good agreement with the numerical simulation.

 

Fig. 14. The red, magenta and black plot colors represent the 3D SG, 2D SG and FT-Gaussian pump beam models, respectively. The inset on the right graph shows the 2D spatial intensity of the simulated (top) and experimentally measured (bottom) Gaussian seed beam after amplification at $P_{\mathrm {p}}=38$ W.

Download Full Size | PDF

There are several interesting observations that can be made for the data presented in Fig. 14. Firstly, the 3D models considering the static $\mathrm {SG^{20}}$ pump beam (red curves) drastically overestimate the power scaling potential of the amplifier as compared to the experimental data. For the case where the dynamic absorption cross-section was used, the slope of the output power follows an almost asymptotic curve which deviates strongly from the anticipated exponential behavior in the small signal gain region. While this behavior was partly mitigated by considering a fixed absorption cross-section, it is clear that 3D modeling of an end-pumped amplifier using a static pump beam gives rise to significant errors. The 2D models (magenta curves) produce more reliable results which are to be anticipated as the majority of the pump absorption occurs within the first few $\mathrm {mm}$ in the crystal. However, as the pump power increases, the absorption of the pump beam in the crystal extends the region of gain along the crystal length due to pump saturation and the 2D models start to deviate from the experimental data more rapidly. By considering the dynamic FT-Gaussian pump beam and a fixed value of the absorption cross-section (black dashed curve) in 3D, the output power at the maximum pump power is underestimated by $2.4\times$ in comparison to the experimental data. These results are circumvented by considering the full model specifications with a dynamic FT-Gaussian pump beam and a weighted absorption cross-section over the pump power levels (solid black curve). This gives rise to the most accurate determination of the power scaling potential as is evident with the high fidelity of the overlap with the experimental data. This model and associated calibration steps offer a holistic approach to determining the power scaling potential of in-line end-pumped Nd:YAG amplifiers.

5. Conclusion

In this manuscript, we presented a novel analytical method for modeling in-line end-pumped Nd:YAG amplifiers using an MMFC diode laser as the pump beam. Our approach incorporated contemporary analytical expressions and introduced a dynamic treatment of the MMFC pump beam’s transverse characteristics, allowing it to undergo a metamorphosis from a flat-top to a Gaussian beam, inspired by lossless laser beam shaping theory. This innovative approach enabled us to treat the pump beam as a complex field and seamlessly integrate it into a Fourier propagation split-step procedure - previously not possible for analytical modeling. This approach demonstrated a remarkably high correlation with the actual shape transformation of the MMFC pump beam throughout the crystal length. To further enhance the robustness of our model, we considered the absorption cross-section as a weighted variable, determined by the output spectrum of a physical diode laser under various pump powers. Accurately representing the behavior of both the pump beam evolution and the absorption cross-section required the use of experimentally determined physical inputs and calibration methods. By combining these calibrated steps with analytical expressions for the gain and absorption coefficients, thermal gradient, and change in refractive index, we developed a comprehensive 3D model, unified by a Fourier beam propagation method. To validate the model, we compared its predictions with the experimental amplification of a Gaussian beam at a wavelength of 1064 nm. The results demonstrated excellent agreement, surpassing previous approximations. Therefore, our analytical model offers an accurate and computationally efficient alternative to complex FEM modeling, making it highly valuable to the laser development community.