2.1 Spectral properties of the SrF₂ crystal

The 0.5%Nd,5%Gd:SrF₂ and 0.5%Nd,5%Y:SrF₂ disordered crystals used in this study can be grown using the vertical Bridgman method^[ Reference Lu, Zhang, Jiang, Kou, Zhang, Xu, Zhao, Wu and Su^26]. Their absorption and emission spectra are shown in Figure 1. The strongest absorption peak for both 0.5%Nd,5%Gd:SrF₂ and 0.5%Nd,5%Y:SrF₂ occurs around 797 nm, which is well-suited for commercial laser diode (LD) pumping. The two emission peaks for 0.5%Nd,5%Gd:SrF₂ are at 1052 and 1059 nm, while the two emission peaks for 0.5%Nd,5%Y:SrF₂ are at 1050 and 1057 nm. When the two crystals are used together, the emission spectrum is approximately flat at the top, which helps to mitigate the gain narrowing effect during amplification^[ Reference Gaul, Martinez, Blakeney, Jochmann, Ringuette, Hammond, Borger, Escamilla, Douglas, Henderson, Dyer, Erlandson, Cross, Caird, Ebbers and Ditmire^27– Reference Cao, Lu and Fan^29]. The key parameters of the crystals are summarized in Table 1.

Figure 1 Emission and absorption spectra for (a) 0.5%Nd,5%Gd:SrF₂ crystal and (b) 0.5%Nd,5%Y:SrF₂ crystal.

Table 1 The parameters of 0.5%Nd,5%Gd:SrF₂ and 0.5%Nd,5%Y:SrF₂.

2.2 Configuration of the amplifier

There are two methods to achieve high gain in an amplifier: one is to increase the small signal gain coefficient, and the other is to increase the optical path length. The former is influenced by the crystal’s intrinsic parameters and the pump energy density, while the latter is determined by the crystal geometry and the beam trajectory. In the pursuit of miniaturization, the design of the crystal shape becomes particularly important.

The proposed compact amplifier configuration is shown in Figure 2(a), where two identical isosceles right-angled triangular prisms are cut on a cuboid crystal, and the two cross-sections of the prisms serve as the incident and exit planes, respectively. This shape can also be considered a combination of a cuboid and a quadrangular prism. The bottom surface of the cuboid is square, while the bottom surface of the quadrangular prism is an isosceles trapezoid. The lower base of the trapezoid and the side of the square have the same length, with the height of the trapezoid being one-fourth of the side length of the square. Three LDs are placed on three sides of the crystal for pumping. The refractive index of the SrF₂ crystal is about 1.4328@1.05 μm^[ Reference Li^30], so the angle for total-internal reflection (TIR) is 44.2615°. The beam undergoes normal incidence and is reflected seven times within the designed amplifier due to TIR.

Figure 2 (a) Configuration of the compact amplifier. (b) Coordinates of the crystal.

Multiple parallel surfaces can easily lead to parasitic oscillations, which significantly affect the gain performance^[ Reference Chen and Xu^31]. To mitigate this issue, the bottom and top surfaces of each crystal are roughened, and in the x–y plane, most sides are tilted by 0.5° or 1° from the original position. Figure 2(b) shows the coordinates of the crystal in the x–y plane. The side length of the square is about 12 mm, which is compatible with the bar size of commercial LDs.

2.3 Setup of the two-stage two-pass amplifier

To increase the gain of the laser amplification system, we designed a two-stage, two-pass amplifier, as shown in Figure 3. Firstly, the seed pulse is p-polarized by passing through a thin-film polarizer (TFP) and injected into amplifier 1 and amplifier 2. Secondly, the beam is directed to a quarter-wave plate (QWP) and then to the end-mirror, which rotates the polarization by 90° after two passes. Thirdly, the returned s-polarized pulse passes through amplifier 2 and amplifier 1 again. Finally, the amplified beam is reflected by the TFP and exits the system. It should be noted that thermally induced depolarization, caused by temperature gradients, could disrupt the operation of the amplifier. Therefore, it is recommended to use the ‘zero depolarization’ orientation of the crystals, which helps minimize the energy loss^[ Reference Snetkov, Yakovlev and Palashov^32].

Figure 3 Schematic of the two-stage two-pass amplifier.

3.1 Calculation model of ASE

In the simulation, the crystal is divided into N parts, where each part is treated as an energy point with the same volume, except for the two angles. The energy supplied to the crystal by the LDs can be expressed as follows^[ Reference Koechner^34]:

(1) $$\begin{align}E={Pt}_{\mathrm{p}}{\eta}_{\mathrm{T}}{\eta}_{\mathrm{A}},\\[-24pt]\nonumber\end{align}$$

where $P$ is the pump power, ${t}_{\mathrm{p}}$ is the pump time and ${\eta}_{\mathrm{T}}$ is the pump transmission efficiency. ${\eta}_{\mathrm{A}}$ is the absorption efficiency, with its calculation model varying based on the object of calculation. When the small signal gain coefficient is calculated, ${\eta}_{\mathrm{A}}$ refers to the total absorption efficiency of the crystal. When the energy distribution is calculated, ${\eta}_{\mathrm{A}}$ is the absorption efficiency of a part of crystal. The expression of ${\eta}_{\mathrm{A}}$ is as follows:

(2) $$\begin{align}{\eta}_{\mathrm{A}} = \left\{\begin{array}{@{}cc@{}} 1-\mathrm{exp}(-\alpha_{\mathrm{D}}l) & (\mathrm{calculate}\ g_{0})\\\mathrm{exp}(-\alpha_{\mathrm{D}}l_2) - \mathrm{exp}(-\alpha_{\mathrm{D}}l_1)& (\mathrm{calculate\ the \ energy\ distribution})\end{array},\right.\\[-24pt]\nonumber\end{align}$$

where ${\alpha}_{\mathrm{D}}$ is the absorption coefficient of the crystal, which is $3.45\ \mathrm{cm}^{-1}$ , l is the total length of the LD beam passing through the crystal, and ${l}_1$ and ${l}_2$ represent the farthest and nearest propagation distances, respectively, that the LD beam travels to reach the boundary of the calculated part of the crystal. The pump time is divided into $t\hbox{'}$ parts, and the increased energy in time interval ( $\Delta t$ ) is ${E}/t\hbox{'}$ . In two adjacent time intervals, the earlier one is defined as the ‘previous time’, and the later one as the ‘current time’. The stored energy at the previous time ( $E\hbox{'}$ ) consists of two parts: the energy supplied to the crystal by the LDs and the energy consumed by the ASE effect (the calculation model for ASE is described later). The stored energy at the current time is as follows:

(3) $$\begin{align}{E}_{\mathrm{S}}=\left(E\hbox{'}+\frac{E}{t\hbox{'}}\right){\eta}_{\mathrm{S}}\exp \left(-\frac{\Delta t}{\tau_{\mathrm{f}}}\right),\end{align}$$

where ${\eta}_{\mathrm{S}}$ is the Stokes efficiency and ${\tau}_{\mathrm{f}}$ is the fluorescence lifetime. The small signal gain coefficient in $\Delta t$ is as follows:

(4) $$\begin{align}{g}_0=\frac{E_{\mathrm{S}}\sigma }{Vh\nu},\end{align}$$

where $\sigma$ is the peak emission cross-section, $V$ is the pump volume, $h$ is Planck’s constant and $\nu$ is the light frequency. The energy produced by spontaneous emission (SE) is as follows:

(5) $$\begin{align}{E}_{\mathrm{SE}}=\left(E\hbox{'}+\frac{E}{t\hbox{'}}\right){\eta}_{\mathrm{S}}\left[1-\exp \left(-\frac{\Delta t}{\tau_{\mathrm{f}}}\right)\right].\end{align}$$

Since the generation of SE is random, it is assumed that the probability of generating SE rays at different angles in $4\pi$ space is uniform. The number of rays per part is denoted as ${N}_{\mathrm{ray}}$ . The total SE energy of each part is denoted as ${E}_{\mathrm{SE}}^{\hbox{'}}$ , so the initial energy of a ray ( ${E}_{\mathrm{ray}}$ ) is equal to ${E}_{\mathrm{SE}}^{\hbox{'}}/{N}_{\mathrm{ray}}$ . The parameters used in the ASE simulation are listed in Tables 1 and 2.

Table 2 Values of the parameters in the ASE simulation.

Figure 4 shows the flow chart of the ASE calculation model. The first step is to select the time slice. In the second step, energy distribution in the crystal is calculated using Equations (1)–(5). The third step involves selecting the starting point and the direction of the trace. The coordinates (x, y, z) in millimeters represent the position of the point. To save computational time, only one z-value is considered, which is half the height of the crystal, since the crystal has the same shape in the z-direction. The direction in $4\pi$ space is given by two angles: the angle between the ray and the z-axis ( $\theta$ ), and the angle between the projection of the ray in the x–y plane and the x-axis ( $\varphi$ ). The fourth step involves comparing the ray energy with the energy threshold. The energy threshold serves as the cutoff condition for ray tracing. If the threshold is set too high, the simulation accuracy decreases. Conversely, if the threshold is set too low, the simulation time increases. After comprehensive consideration, the threshold is set to 5% of the initial ray energy. When the ray energy is lower than the threshold, the simulation returns to the third step to select another point. Otherwise it proceeds to the next step. The fifth step involves ray tracing. The amplified ray energy can be expressed as follows:

(6) $$\begin{align}{E}_{\mathrm{A}}={E}_{\mathrm{ray}}\exp \left[\left({g}_0-\alpha \right)L\right]R,\end{align}$$

where $\alpha$ is the loss coefficient of the crystal (0.0015 cm^–1) and $L$ is the optical path length. The reflections ( $R$ ) on the two rough surfaces of the crystal are neglected, while the reflection coefficients of the remaining surfaces are calculated using the Fresnel equations. The sixth step is to record the position and energy of the ray, and then return to the fourth step. After the energy of all rays has been calculated, the simulation returns to the first step to select the next time slice. The model continues until the energy for all time slices has been calculated.

Figure 4 Flow chart of the ASE calculation model.

3.2 Calculation results of ASE

To simplify the calculation, the divergence angle of the LD is set to zero. This assumption can be fulfilled using collimation lens arrays, which can control the fast and slow axis divergence angles of LD bars to within 0.35°^[ Reference Guo, Gao, Wu, Lv, Li and Li^35– Reference Liu, Zhao, Xiong and Liu^37]. The energy distribution is shown in Figure 5(a) with normalized values. In the x–y plane, more energy is concentrated at the two upper corners of the crystal, while the energy is almost zero at the bottom of the crystal. The relation between the small signal gain coefficient ${g}_0$ and pump time is plotted in Figure 5(b). The dashed lines indicate the results when ASE is not considered, and the solid lines indicate the results when ASE is included. The value of ${g}_0$ gradually increases, with the growth rate slowing down as the simulation time progresses. In addition, the influence of ASE on ${g}_0$ increases with higher pump power. For example, when the pump time is set to 320 μs, the presence of ASE reduces ${g}_0$ by 0.0764 cm^–1 when the pump power is 3.6 kW, and by 0.1602 cm^–1 when the pump power is 7.2 kW. Figures 5(c) and 5(d) show the ASE energy distribution of the output surface x’O’z’ when the pump power is 6 kW. The ASE energy increases slightly along the x’-axis due to the asymmetrical shape in the x–y plane. In the z’-direction, the ASE energy is the highest in the middle, with a symmetrical distribution on both sides. This is because only the beam from the starting points in the middle of the crystal (z = 3 mm) is tracked. The ray is reflected most of the time within the crystal, and the gain is highest when $\theta =0^{\circ}$ . The ray path to both ends of the z’-axis is longer, so the gain is also relatively higher when $\theta \ne 0^{\circ}$ . In summary, the higher the pump energy, the higher the gain, and the stronger the ASE effect. Therefore, to achieve the desired ${g}_0$ , the pump energy must be increased beyond the value required to reach the target ${g}_0$ without considering ASE.

Figure 5 (a) The normalized energy distribution. (b) The relation between the small signal gain coefficient and the pump time. (c) The energy distribution of ASE along the x’-axis. (d) The energy distribution of ASE along the z’-axis.

4.1 Calculation model of broadband pulse propagation

The broadband characteristics are simulated using MATLAB software. The broadband pulse propagation in the amplifier is derived from the resonant-dipole equation, the population inversion equation and Maxwell’s equation. During the derivation process, several methodologies, including Fourier transformation and inverse transformation, the slowly varying envelope approximation and the fourth-order Runge–Kutta method, are employed. The resulting expression is given as follows^[ Reference Chuang, Zheng and Meyerhofer^38]:

(7) $$\begin{align}&\frac{\partial {E}_0\left(z,t\right)}{\partial z}=-i\frac{\omega_0}{2\varepsilon c}{P}_0\left(z,t\right)+i\frac{\beta^{\hbox{'}\hbox{'}}}{2}\frac{\partial^2{E}_0\left(z,t\right)}{\partial {t}^2}-i\frac{\beta_2}{2}{\left|{E}_0\left(z,t\right)\right|}^2{E}_0\left(z,t\right),\nonumber\\[4pt] &\frac{\partial {P}_0\left(z,t\right)}{\partial t}+\frac{\Delta {\omega}_{\mathrm{a}}}{2}\left[1+i\frac{2\left({\omega}_0-{\omega}_{\mathrm{a}}\right)}{\Delta {\omega}_{\mathrm{a}}}\right]{P}_0\left(z,t\right)=i\frac{K}{2{\omega}_0}N\left(z,t\right){E}_0\left(z,t\right),\nonumber\\[4pt] &\frac{\partial N\left(z,t\right)}{\partial t}=i\frac{2^{\ast }}{4\mathrm{\hslash}}\left[{E}_0^{\ast}(z,t){P}_0(z,t)-{E}_0(z,t){P}_0^{\ast}(z,t)\right],\end{align}$$

where ${E}_0\left(z,t\right)$ is the electric field intensity, ${\omega}_0$ is the carrier frequency, $\varepsilon$ is the dielectric constant of the amplifier, $c$ is the velocity of light in the amplifier, ${P}_0\left(z,t\right)$ is the resonant polarization, ${\beta}^{\hbox{'}\hbox{'}}$ is the second-order derivative of the propagation constant at ${\omega}_0$ , the constant ${\beta}_2$ is equal to $2\pi {n}_{2\mathrm{E}}/{\lambda}_0$ , ${n}_{2\mathrm{E}}$ is the nonlinear refractive index, which is 1.4628×10^–20 m²/W^[ Reference Adair, Chase and Payne^39], ${\lambda}_0$ is the laser wavelength in vacuum, $\Delta {\omega}_{\mathrm{a}}$ is the full atomic linewidth, ${\omega}_{\mathrm{a}}$ is the center frequency of the line, the constant $K$ is equal to $\varepsilon c\sigma \Delta \omega_{\mathrm{a}}$ , $N\left(z,t\right)$ is the population inversion, ${2}^{\ast }$ is a dimensionless population saturation factor with values between 1 and 2 and $\mathrm{\hslash}$ is equal to $h/\left(2\pi \right)$ . The seed pulse waveform is defined as a fourth-order ultra-Gaussian distribution. It has a beam diameter of 2 mm, a wavelength of 1053 nm, a pulse duration of 3.5 ns, an energy of 3 nJ and a bandwidth of 18 nm. The beam passes through the two-stage two-pass amplifier, as depicted in Figure 3, where the pump power is 6 kW. The reflectance of the TFP for the s-polarized pulse and the reflectance of the mirror are both 99%. Meanwhile, the transmittance of the TFP for the p-polarized pulse and the transmittance of other optical elements reach 99%. The thicknesses of the TFP and QWP are 3 and 1.6 mm, respectively.

4.2 Calculation results of broadband pulse propagation

The two-stage two-pass amplifier features three configurations based on the type of crystal utilized: (1) both crystals are 0.5%Nd,5%Gd:SrF₂; (2) both crystals are 0.5%Nd,5%Y:SrF₂; and (3) the first crystal is 0.5%Nd,5%Gd:SrF₂ while the second one is 0.5%Nd,5%Y:SrF₂. When both crystals are of the same type, the pump time is set to be equal and is incremented by 5 μs at a time until the output energy exceeds 4.8 mJ. In cases where the crystal types differ, the increment step for the pump time of Nd,Gd:SrF₂ is 20 μs, and the pump time for Nd,Y:SrF₂ is increased by 5 μs at a time until the output energy exceeds 4.8 mJ.

Figure 6 Normalized distributions of the output pulse under three configurations: (a), (b) both crystals are 0.5%Nd,5%Gd:SrF₂; (c), (d) both crystals are 0.5%Nd,5%Y:SrF₂; and (e), (f) the first crystal is 0.5%Nd,5%Gd:SrF₂ while the second crystal is 0.5%Nd,5%Y:SrF₂. (a), (c) and (e) are time waveforms, while (b), (d) and (f) are spectra.

Table 3 The broadband characteristic under various conditions.

The normalized distributions of the output pulse are shown in Figures 6(a)–6(d) when the crystal type is the same. The black, red and blue lines represent the distribution of the input pulse, the pulse after the first amplification and the output pulse, respectively. ‘ $\Delta$ Time’ and ‘ $\Delta$ Wavelength’ are the offset time and wavelength relative to the reference time and wavelength, respectively. During the amplification process, the width of both the time waveform and the spectrum narrows, and the normalized intensity difference between the main peak and the secondary peak increases. This narrowing of the output spectrum bandwidth leads to a large Fourier transform limited pulse duration, making it more challenging to achieve high peak power. To increase the output bandwidth, Nd,Gd:SrF₂ and Nd,Y:SrF₂ are used in combination. Figures 6(e) and 6(f) represent the time waveform distribution and the spectrum distribution of the output pulse, respectively. The broadband characteristics under diverse conditions are presented in Table 3. As the pump time of the first crystal increases and the pump time of the second crystal decreases, the normalized intensity of the secondary peak gradually approaches or even exceeds that of the main peak. Moreover, the bandwidth is more than doubled compared to using the same crystal type. This phenomenon highlights the advantage of simultaneously using crystals with different emission wavelengths. Furthermore, different emission wavelengths of crystals can be achieved by tuning the doping concentrations or co-dopant species, a strategy that can be applied not only to SrF₂ crystals but also to other substrate materials, such as calcium fluoride^[ Reference Ma, Zhang, Jiang, Zhang, Kou, Strzep, Tang, Zhou, Zhang, Zhang, Zhu, Yin, Lv, Li, Chen and Su^40].

The additional phase acquired by the nonlinear effect during laser transmission is defined as the B-integral, and the expression is given by the following^[ Reference Perry, Ditmire and Stuart^41]:

(8) $$\begin{align}B=\frac{2\pi {n}_{2\mathrm{E}}}{\lambda_0}\int {E}_0{\left({z},{t}\right)}^2\mathrm{d}z.\end{align}$$

The B-integral of the two-stage two-pass amplifier is significantly lower than 1, which ensures the output pulse quality.

When the ASE is not considered, the pump power is 6 kW, the pump time for Nd,Gd:SrF₂ is 320 μs and the small signal gain coefficient of Nd,Gd:SrF₂ is 0.615 cm^–1. As shown in Figure 5(b), when the ASE is considered, the expression for the green solid line, which represents a pump power of 7.2 kW, is approximated by the following:

(9) $$\begin{align}{g}_{0\;\left(7.2\mathrm{kW}\right)}=0.66504-0.66616\exp \left(-\frac{t_{\mathrm{P}}}{188.49189}\right).\end{align}$$

When ${g}_{0\;\left(7.2\mathrm{kW}\right)}$ is 0.615, the corresponding pump time is 488 μs. The pump power of a commercial LD bar operating at a pump time of less than 500 μs and a repetition rate of less than 1 Hz can reach 500 W at 797 nm, and five bars can be placed within 2 mm. This design configuration is therefore feasible and should be realizable based on the given parameters.

5.1 Calculation model of temperature distribution

The temperature distribution of Nd,Gd:SrF₂ is simulated using COMSOL Multiphysics software. The transient-state equation to describe heat transfer in solids is given by the following^[ Reference Bavil, Safari and Mech^42]:

(10) $$\begin{align}\rho {{C}}_{\mathrm{p}}\frac{\partial T}{\partial t}+\nabla \cdot \left(-k\nabla T\right)=Q,\end{align}$$

where $\rho$ is the solid density, ${{C}}_{\mathrm{p}}$ is the solid heat capacity at constant pressure, $T$ is the temperature, $t$ is the simulation time and $k$ is the thermal conductivity. Here, $Q$ is the heat source, which can be written as follows:

(11) $$\begin{align}&{Q}_{\left(y\approx 12\ \mathrm{mm}\right)}\nonumber\\&\quad=\left\{\begin{array}{@{}cc}\frac{\alpha_{\mathrm{D}}\left(1-{\eta}_{\mathrm{S}}\right){P}_{\mathrm{s}}}{S_{\mathrm{LD}}}\exp \left[-{\alpha}_{\mathrm{D}}\left(l-y\right)\right]&\left(2\;\mathrm{mm}\le z\le 4\;\mathrm{mm},1\;\mathrm{mm}\le x\le 11\;\mathrm{mm}\right)\\ {}0& \left(\mathrm{else}\right)\end{array}\right.\!\!\!\!\!\!\!\!,\end{align}$$

(12) $$\begin{align}&{Q}_{\left(x\approx 0\ \mathrm{mm}\right)}\nonumber\\&\quad=\left\{\begin{array}{@{}cc}\frac{\alpha_{\mathrm{D}}\left(1-{\eta}_{\mathrm{S}}\right){P}_{\mathrm{s}}}{S_{\mathrm{LD}}}\exp \left(-{\alpha}_\mathrm{D}x\right)& \left(2\;\mathrm{mm}\le z\le 4\; \mathrm{mm},1\;\mathrm{mm}\le y\le 11\; \mathrm{mm}\right)\\ {}0& (\mathrm{else})\end{array}\right.\!\!\!\!\!\!\!\!,\end{align}$$

(13) $$\begin{align}&{Q}_{\left(x\approx 12\ \mathrm{mm}\right)}\nonumber\\&\quad=\left\{\begin{array}{@{}cc}\frac{\alpha_{\mathrm{D}}\left(1-{\eta}_{\mathrm{S}}\right){P}_{\mathrm{s}}}{S_{\mathrm{LD}}}\exp \left[-{\alpha}_{\mathrm{D}}\left(l-x\right)\right]& \left(2\;\mathrm{mm}\le z\le 4\; \mathrm{mm},1\;\mathrm{mm}\le y\le 11\; \mathrm{mm}\right)\\ {}0& (\mathrm{else})\end{array}\right.\!\!\!\!\!\!\!\!.\end{align}$$

Here, ${S}_{\mathrm{LD}}$ is the area of the pump light and ${P}_{\mathrm{s}}$ is the pump power of a single LD, which is equal to $P/3$ . The other parameters have been described in Section 3.1. Equations (11)–(13) represent the distribution of heat sources when three LDs pump the crystal at $y\approx 12$ mm, $x\approx 0$ mm and $x\approx 12$ mm, respectively.

The boundary conditions for establishing the model of the Nd,Gd:SrF₂ crystal include two aspects, the ambient temperature ${T}_0$ and the air convection mechanism as follows:

(14) $$\begin{align}\left\{\begin{array}{@{}l}{k}_x\frac{\partial T}{\partial x}={h}_{\mathrm{a}}\left(T-{T}_0\right)\\\\[-6pt] {}{k}_y\frac{\partial T}{\partial y}={h}_{\mathrm{a}}\left(T-{T}_0\right)\\\\[-6pt] {}{k}_z\frac{\partial T}{\partial z}={h}_{\mathrm{a}}\left(T-{T}_0\right)\end{array}\right..\end{align}$$

Nd,Gd:SrF₂ is an isotropic material with uniform thermal conductivity of $k$ in all directions ( ${k}_x={k}_y={k}_z=k$ ). Since the amplifier has no cooling device, it is assumed that all sides of the crystal are in contact with air, and the heat transfer coefficient is ${h}_{\mathrm{a}}$ . The values of the parameters in the thermal simulation are presented in Table 4.

Table 4 Values of the parameters in the thermal simulation.

5.2 Calculation results of temperature distribution

Figure 7 depicts the temperature distribution of Nd,Gd:SrF₂ over time. The output time step in Figures 7(a)–7(c) is set to 50 μs with a total simulation time of 5 s. Due to the hardware limitations, an ‘adaptive timestepping scheme’ is employed in COMSOL Multiphysics software to optimize memory usage and computational time. This scheme allows the software to automatically adjust the time step to maintain the tolerance, which is governed by physics^{[ 43]}. In the y–z plane (x = 0 mm, Figure 7(a)) and x–y plane (z = 3 mm, Figure 7(b)), the highest temperatures are concentrated in the area where the Nd,Gd:SrF₂ crystal is pumped. As time progresses, the temperature in the pumped area gradually decreases, while the temperature in the surrounding area increases slightly due to heat conduction.

Figure 7 Temperature distribution of Nd,Gd:SrF₂ crystal. The perspectives are as follows: (a) y–z plane, x = 0 mm; (b) x–y plane, z = 3 mm; (c), (d) temperature distribution at three coordinate points.

Figure 8 Schematic diagram of water cooling.

The temperature curves for the coordinate points (–0.055, 6, 3), (6, 6, 3) and (0, 0, 0) are shown in Figures 7(c) and 7(d) with all coordinates provided in millimeters. Figure 7(c) focuses on the temperature change within 5 s. At 0.005 s after the start of the pumping laser pulse, the highest temperature reached 293.254 K at (–0.055, 6, 3), and then decreased. A similar temperature trend is observed at (6, 6, 3). At (0, 0, 0), the temperature slowly increases from the ambient temperature. After 3 s from the start of the pumping laser pulse, the temperature at all three coordinate points stabilizes. To more accurately capture the highest temperatures at the three coordinate points, the output time step in Figure 7(d) is set to 10 μs with a simulation time of 0.03 s. The highest temperature is 296.109 K at (–0.055, 6, 3). Further reduction in the time step shows that the highest temperature remains virtually unchanged.

It can be postulated that the implementation of cooling devices, such as water cooling and heat sinks, may facilitate the operation of the laser at a repetition rate exceeding 1 Hz. Figure 8 shows a schematic diagram of a water cooling structure that can be referenced. Two cooling plates are in contact with the bottom and top surfaces of the crystal, respectively. Inside the cooling plates are many microfluidic channels, and the blue arrows indicate the direction of coolant flow. The heat deposited in the crystal can be dissipated by setting appropriate cooling temperatures and flow velocities^[ Reference Nagymihaly, Cao, Papp, Hajas, Kalashnikov, Osvay and Chvykov^44]. The output performance of the laser remains almost unaffected as long as the temperature remains within tolerable limits over a period of 1 h.

Figure 9 (a) Temperature as a function of time for three coordinate points. (b) Von Mises stress distribution in the x–y plane. (c) Displacement at different positions along the line segment from coordinate point (0, 0, 3) to (0, 12, 3). (d) Displacement at different positions along the line segment from coordinate point (0, 12, 3) to (12, 12, 3).

5.3 Calculation model of stress distribution

The heat load creates a thermal gradient that can lead directly to stress in the crystal, and stress can cause adverse effects such as birefringence and thermal lensing. The stress distribution of SrF₂ is also simulated using COMSOL Multiphysics software. In three dimensions, the stress ( ${\sigma}_{\mathrm{s}}$ ) and strain ( ${\varepsilon}_{\mathrm{s}}$ ) can be written as follows^{[ 45]}:

(15) $$\begin{align}{\sigma}_{\mathrm{s}}=\left[\begin{array}{ccc}{\sigma}_{xx}& {\sigma}_{xy}& {\sigma}_{xz}\\ {}{\sigma}_{yx}& {\sigma}_{yy}& {\sigma}_{yz}\\ {}{\sigma}_{zx}& {\sigma}_{zy}& {\sigma}_{zz}\end{array}\right],\end{align}$$

(16) $$\begin{align}{\varepsilon}_{\mathrm{s}}=\left[\begin{array}{ccc}{\varepsilon}_{xx}& {\varepsilon}_{xy}& {\varepsilon}_{xz}\\ {}{\varepsilon}_{yx}& {\varepsilon}_{yy}& {\varepsilon}_{yz}\\ {}{\varepsilon}_{zx}& {\varepsilon}_{zy}& {\varepsilon}_{zz}\end{array}\right],\end{align}$$

where each element in ${\varepsilon}_{\mathrm{s}}$ is defined as a derivative of the displacements:

(17) $$\begin{align}\left[\begin{array}{c}{\varepsilon}_{xx}\\ {}{\varepsilon}_{yy}\\ {}{\varepsilon}_{zz}\\ {}{\varepsilon}_{xy}\\ {}{\varepsilon}_{yz}\\ {}{\varepsilon}_{xz}\end{array}\right]=\left[\begin{array}{c}\frac{\partial u}{\partial x}\\ {}\frac{\partial v}{\partial y}\\ {}\frac{\partial w}{\partial z}\\ {}\frac{1}{2}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\\ {}\frac{1}{2}\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\\ {}\frac{1}{2}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)\end{array}\right],\end{align}$$

where u, v and w are the displacement vectors. The thermal strain ( $\Delta {\varepsilon}_{\mathrm{s}}$ ) is given by the following^[ Reference Harish, Bharatish, Murthy, Anand and Subramanya^46]:

(18) $$\begin{align}\Delta {\varepsilon}_{\mathrm{s}}={\alpha}_{\mathrm{e}}\Delta T,\end{align}$$

where  ${\alpha}_{\mathrm{e}}$ is the thermal expansion coefficient, which is 18.4×10^–6 K^–1 for SrF₂, and $\Delta T$ is the temperature difference. The three dimensions generalization of Hooke’s law can be written as follows:

(19) $$\begin{align}\left[\begin{array}{c}{\sigma}_{xx}\\ {}{\sigma}_{yy}\\ {}{\sigma}_{zz}\\ {}{\sigma}_{xy}\\ {}{\sigma}_{yz}\\ {}{\sigma}_{xz}\end{array}\right]={D}\left[\begin{array}{c}{\varepsilon}_{xx}\\ {}{\varepsilon}_{yy}\\ {}{\varepsilon}_{zz}\\ {}{\varepsilon}_{xy}\\ {}{\varepsilon}_{yz}\\ {}{\varepsilon}_{xz}\end{array}\right].\end{align}$$

For the isotropic SrF₂ crystal, ${D}$ is a function of Young’s modulus ( ${E}_{\mathrm{Y}}=99.91\times {10}^9$ Pa) and Poisson’s ratio ( ${\nu}_{\mathrm{P}}=0.25$ ):

(20) $$\begin{align}{D}=\frac{E_{\mathrm{Y}}}{\left(1+{\nu}_{\mathrm{P}}\right)\left(1-2{\nu}_{\mathrm{P}}\right)}\left[\begin{array}{@{}c@{\ }c@{\ }c@{\ }c@{\ }c@{\ }c@{}}1-{\nu}_{\mathrm{P}}& {\nu}_{\mathrm{P}}& {\nu}_{\mathrm{P}}& 0& 0& 0\\ {}{\nu}_{\mathrm{P}}& 1-{\nu}_{\mathrm{P}}& {\nu}_{\mathrm{P}}& 0& 0& 0\\ {}{\nu}_{\mathrm{P}}& {\nu}_{\mathrm{P}}& 1-{\nu}_{\mathrm{P}}& 0& 0& 0\\ {}0& 0& 0& \frac{1-2{\nu}_{\mathrm{P}}}{2}& 0& 0\\ {}0& 0& 0& 0& \frac{1-2{\nu}_{\mathrm{P}}}{2}& 0\\ {}0& 0& 0& 0& 0& \frac{1-2{\nu}_{\mathrm{P}}}{2}\end{array}\right].\end{align}$$

Since the bottom and top surfaces of the SrF₂ are expected to be clamped, their displacement changes are constrained to be zero. The values of the parameters related to temperature simulation are shown in Table 4.

5.4 Calculation results of stress distribution

The seed pulse is injected into the crystals at the end of the pump time, so the simulation time of stress distribution is set to 500 μs. The output time step is set to 50 μs. To verify the effectiveness of this time step, we simulate the temperature changes over time at three coordinate points, (–0.055, 6, 3), (6, 6, 3) and (0, 0, 0), as shown in Figure 9(a). The highest temperature of 295.458 K is obtained at (–0.055, 6, 3), which is consistent with the result at the same moment in Figure 7(d). The top view of the stress distribution in the x–y plane is shown in Figure 9(b), and the maximum stress is 6.6×10⁶ Pa. To make the stress-induced deformation of the crystal more visible, the value is magnified by a factor of 6000 in Figure 9(b). It can be observed that the deformation occurs mainly in the three pumped areas, and the degree of deformation on the left- and right-hand sides is essentially the same.

Figure 9(c) shows the displacement at different positions along the line segment from coordinate point (0, 0, 3) to (0, 12, 3). The displacement size increases with the pump time, and the maximum displacement size is 1.198×10^–7 m at 500 μs at coordinate point (0, 10.6, 3). Figure 9(d) shows the displacement at different positions along the line segment from coordinate point (0, 12, 3) to (12, 12, 3). There are two peaks near the positions x = 1.3 mm and x = 10.7 mm at different times. The maximum displacement size is 1.246×10^–7 m at 500 μs at coordinate point (10.7, 12, 3). This work focuses on the simulation of the transient state with the backward differentiation formula^{[ 47]}. In future steady-state analyses, the effects of thermally induced birefringence and thermal lensing can be explored in more detail.