Thermo-Optical Modeling of Double-Pass Type-II Second Harmonic Generation in KTP

1. Introduction

1.1. Background

Second harmonic generation (SHG) converts two photons at a fundamental frequency $\omega$ into a single photon at $2\omega$ inside a nonlinear crystal. The process is widely used to produce visible coherent light from infrared sources [1,2], with applications in material processing, biomedical imaging, spectroscopy, underwater communications [3], and display technology. Potassium titanyl phosphate (${\mathrm{KTiOPO}}_4$, KTP) is one of the most common crystals for converting $1064\mathrm{nm}$ Nd:YAG or Nd:YVO₄ output to $532\mathrm{nm}$[4,5], owing to its effective nonlinear coefficient ($d_{\mathrm{eff}}=7.3\mathrm{pm}/\mathrm{V}$) [4], high damage threshold [4], broad transparency range (350–$4400\mathrm{nm}$) [4], and favorable thermal conductivity [6].

In continuous-wave (CW) operation, tight Gaussian focusing raises conversion efficiency but also deposits heat through residual absorption at both fundamental and harmonic wavelengths. Two consequences limit performance. First, temperature-dependent refractive indices perturb the phase-matching condition

$$\Delta k=k_1+k_2-k_3=0,$$

(1)

producing thermally induced phase mismatching (TIPM) [7]. Second, the resulting index gradients create a thermal lens that distorts beam profiles [8,9]. Under CW pumping, steady-state temperature distributions develop over $0.1$–$1\mathrm{s}$ [10], making thermal management integral to the optical design.

Double-pass configurations, in which a mirror returns residual fundamental and generated second-harmonic waves for a second transit, effectively double the interaction length and improve efficiency over single-pass designs [11,12]. Type-II phase matching—where ordinary and extraordinary fundamental polarizations combine to generate an extraordinary second-harmonic field—is well suited to KTP and offers favorable walk-off and acceptance bandwidth characteristics [7,13]. Reported CW efficiencies range from $11\%$ to $28.5\%$ [1,2], with thermal effects identified as the principal limitation [5].

1.2. Modeling Challenges

Three classes of partial differential equations govern heat-coupled SHG: nonlinear wave equations for field propagation and frequency conversion, transient three-dimensional heat diffusion with absorption-dependent source terms, and auxiliary relations mapping temperature to refractive index and accumulated phase mismatch. Despite work by several groups on subsets of this physics, four gaps persist.

First, most published models are implemented as one-off scripts for specific configurations, without publicly available code or standardized inputs. Re-implementation from incomplete descriptions introduces unquantified discrepancies between nominally equivalent models.

Second, boundary-condition treatments are often oversimplified. Perfectly isothermal or fully adiabatic assumptions do not represent practical cooling. Mixed Robin conditions coupling conduction, convection, and radiative exchange are more realistic but computationally demanding and therefore frequently omitted [14], even though they can shift predicted temperatures by $30\%$–$50\%$.

Third, steady-state solvers cannot capture the efficiency transients observed during the first $0.5$–$1\mathrm{s}$ after pump turn-on, which are governed by KTP’s thermal diffusivity $\alpha =K/\left(\rho C\right)\approx {10}^{-6}{\mathrm{m}}^2/\mathrm{s}$ [15].

Fourth, comprehensive parametric sweeps over crystal length, pump power, beam waist, and coolant temperature [7,16] remain rare because of computational cost. This is particularly limiting for double-pass type-II configurations [11,17], which involve six propagating waves requiring careful treatment of mode overlap, polarization-dependent absorption, and mirror boundary conditions.

1.3. Prior Work

The present model builds on a series of studies of thermal and nonlinear-optical phenomena in KTP. Sabaeian et al. [18] obtained analytical steady-state temperature distributions for cubic laser crystals. Subsequent work [15,19] introduced temperature-dependent thermal conductivity and radiation boundary conditions, showing that $K\left(T\right)$ variations are significant while radiation terms are negligible. Rezaee et al. [20,21] solved the time-dependent heat equation for repetitively pulsed Gaussian beams using FDTD methods, and later [22,23] addressed spatiotemporal TIPM through experimentally fitted $n_i\left(T\right)$ dispersion relations. The undepleted-pump assumption was removed in [24,25] by solving coupled amplitude equations for double-pass SHG, and the framework was extended to pulsed Bessel–Gauss beams in [26,27].

The eight-equation coupled model combining TIPM and thermal lensing was formulated and experimentally validated in [28,29]. The present work provides a documented implementation of this model with reduced memory consumption (approximately 60% lower than the original research code), adaptive time-stepping, and spatially resolved diagnostics.

1.4. Contributions

This work makes three contributions. (i) A version-controlled Fortran implementation with standardized input files enables replication of the published results. (ii) Optimized memory management allows execution on hardware with less than $8\mathrm{GB}$ RAM, with runtimes below 30 minutes per parameter set on consumer-grade CPUs. (iii) Validation against experimental data [30] shows deviations below $4\%$ in conversion efficiency over pump powers from 22 to $53\mathrm{W}$ and beam waist sizes from 50 to $75\mu \mathrm{m}$.

2. Numerical Model

Table 1 lists the principal inputs, internally solved quantities, and reported outputs. The solver is organized into four modules.

The Core Solver integrates the eight coupled partial differential equations describing double-pass type-II SHG with thermal effects. It handles both transient and steady-state regimes through adaptive time-stepping. Wave propagation uses a shooting method with iterative boundary enforcement at the crystal entrance ($z=0$) and the mirror ($z=l$). Forward fundamental waves—ordinary (${\psi }_{1+}$) and extraordinary (${\psi }_{2+}$)—together with the forward second-harmonic wave (${\psi }_{3+}$) are advanced from $z=0$ to $z=l$ by a fourth-order Runge-Kutta integrator with adaptive refinement in regions of rapid phase variation. At the mirror, reflectivity coefficients $R_{f,\mathrm{FW}}$ and $R_{f,\mathrm{SHW}}$ seed backward waves (${\psi }_{1-}$, ${\psi }_{2-}$, ${\psi }_{3-}$) for the return transit. Transverse Laplacian terms ${\nabla }_T^2{\psi }_i$ are discretized on a cylindrical grid using second-order centered finite differences over $N_r\times N_z$ points. An optional paraxial mode omits transverse diffraction for beam waists below $100\mu \mathrm{m}$.

The Geometry and Materials module stores crystal-specific data: length l, radius a, phase-matching angles $\theta$ and $\varphi$, temperature-dependent thermal conductivity $K\left(T\right)=K_0/(T_0+\Delta T)$, specific heat C, density $\rho$, and polynomial fits for principal refractive indices $n_X(T,\lambda )$, $n_Y(T,\lambda )$, $n_Z(T,\lambda )$ from Eqs. (14)–(16) of [28]. Ordinary and extraordinary indices $n_o(\theta ,\varphi ,T)$ and $n_e(\theta ,\varphi ,T)$ are computed for arbitrary cut angles via Eqs. (17)–(19) of the same reference; their temperature derivatives feed the TIPM calculation. Default absorption coefficients are ${\gamma }_{1,2}=0.005{\mathrm{cm}}^{-1}$ and ${\gamma }_3=0.02{\mathrm{cm}}^{-1}$; the default nonlinear coefficient is $d_{\mathrm{eff}}=7.3\mathrm{pm}/\mathrm{V}$ for type-II matching at $\theta ={90}^{\circ }$, $\varphi =23.4^{\circ }$.

The Boundary and Cooling module implements thermal boundary conditions. End faces ($z=0$, $z=l$) satisfy mixed Robin conditions $-K\left(T\right)\nabla T·\mathbf{n}=h(T-T_{\infty })+\sigma \epsilon (T^4-T_s^4)$, where h is the convective coefficient, $T_{\infty }$ the coolant temperature, $\sigma$ the Stefan-Boltzmann constant, $\epsilon$ the emissivity, and $T_s$ the ambient radiative temperature. Presets include water cooling ($h=10\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$, $T_{\infty }=298\mathrm{K}$), air convection, and adiabatic limits. The lateral surface ($r=a$) accepts a Dirichlet condition $T(r=a,z,t)=T_0$ or a Robin condition. The $T^4$ nonlinearity is linearized via Newton-Raphson at each time step.

The Parametric Sweep module automates scans over pump power $P_F$, beam waist ${\omega }_f$, crystal length l, or coolant temperature $T_0$. Per-case outputs include $T(r,z)$ distributions, phase mismatch $\Delta \varphi \left(z\right)$, field amplitudes $|{\psi }_i^{\pm }{\left(z\right)|}^2$, conversion efficiency ${\eta }_{\mathrm{SHG}}=P_{\mathrm{SH},\mathrm{out}}/P_{F,\mathrm{in}}$, and peak temperature $T_{\max }$. Results are written in Tecplot ASCII format (.plt) with supplementary CSV summaries.

Execution proceeds in three stages. Initialization reads crystal geometry, material properties, beam parameters ($P_F$, ${\lambda }_{\mathrm{pump}}=1064\mathrm{nm}$, ${\omega }_f$), mirror reflectivities, boundary conditions, and numerical controls ($\Delta t$, mesh resolution, convergence tolerance). Computation iterates between wave and heat solvers: at each time step the wave equations are solved for ${\psi }_i^n(r,z)$ using the current temperature $T^n(r,z)$, the absorption source $S^n(r,z)$ is evaluated, the heat equation advances to $T^{n+1}(r,z)$, TIPM is updated to $\Delta {\varphi }^{n+1}\left(z\right)$, and convergence is tested via $|{\eta }_{\mathrm{SHG}}^{n+1}-{\eta }_{\mathrm{SHG}}^n|<{ϵ}_{\mathrm{tol}}$ (${ϵ}_{\mathrm{tol}}=0.01$). Output writes post-processed results to disk.

A typical build-and-run sequence is:

Optional flags include –transient, –no-diffraction, and –threads for OpenMP parallelization.

At the default resolution of $500\times 200$ points and $231,087$ time steps ($1\mathrm{s}$ of physical time), memory usage remains below $6\mathrm{GB}$. Loop ordering maximizes cache reuse, achieving approximately $75\%$ of peak throughput on x86-64 CPUs with AVX2. OpenMP threading over 4–8 cores gives near-linear speedup for $N_r>200$. Interprocedural optimization (-ipo) and fast floating-point arithmetic (-fp-model fast=2) reduce runtime by 30–$40\%$ with negligible accuracy loss, as confirmed by double-precision reference runs. For the baseline case ($P_F=53\mathrm{W}$, ${\omega }_f=50\mu \mathrm{m}$, $l=9\mathrm{mm}$), single-core execution takes approximately 28 minutes on a $3.2\mathrm{GHz}$ Intel Core i5; eight-core execution takes $4.5$ minutes. Reducing $N_r$ to 250 lowers peak memory below $2\mathrm{GB}$ at the cost of less than $2\%$ in accuracy. The code compiles with Intel Fortran (ifort 19.0+) or GNU Fortran (gfortran 9.0+) and requires no external libraries. Supported platforms include Linux (Ubuntu 20.04+) and Windows via WSL2.

Table 2 shows the coupling paths between subsystems: what each subsystem receives and what it passes onward.

3. Validation

The primary benchmark is the double-pass CW green laser of Bai and Chen [30]: a $9\mathrm{mm}$ KTP crystal ($\theta ={90}^{\circ }$, $\varphi =23.4^{\circ }$) in a water-cooled copper block ($T_0=298\mathrm{K}$) with a convex lens (M₃, $R=100\mathrm{mm}$) and a concave mirror (M₄, $R=50\mathrm{mm}$). For pump powers $P_F=22$–$53\mathrm{W}$ and ${\omega }_f=50\mu \mathrm{m}$, the reported efficiency was $25.5\%$ at $22\mathrm{W}$. Under identical conditions the model yields ${\eta }_{\mathrm{SHG}}=29.0\%$ at $22\mathrm{W}$ and $39.8\%$ at $53\mathrm{W}$, deviations of $3.5\%$ and below $4\%$ after accounting for unmodeled cavity losses. Predicted efficiency plateaus near $39\%$ over $P_F=37$–$65\mathrm{W}$, consistent with the saturation behavior shown in the left panel of Figure 4 of [28].

Transient behavior is tested at $P_F=53\mathrm{W}$ and ${\omega }_f=50\mu \mathrm{m}$. The model reproduces a monotonic efficiency decay from ${\eta }_{\mathrm{SHG}}(t=0)=50.0\%$ to $39.8\%$ over approximately $1\mathrm{s}$, matching Figure 6 of [28]. The $10.2$ percentage-point drop corresponds to on-axis heating from $298\mathrm{K}$ to $T_{\max }=303.5\mathrm{K}$ and an accumulated TIPM of $\Delta \varphi (z=l)\approx -1.1\mathrm{rad}$. The extracted thermal time constant ${\tau }_{\mathrm{th}}\approx 0.2\mathrm{s}$ is consistent with ${\tau }_{\mathrm{th}}\sim a^2/\left({\pi }^2\alpha \right)$, confirming correct implementation of transient heat diffusion. Runtime for $231,087$ time steps on the $500\times 200$ mesh is consistent with the figures reported in Section 2.

A crystal-length sweep at $P_F=53\mathrm{W}$ and ${\omega }_f=50\mu \mathrm{m}$ yields ${\eta }_{\mathrm{SHG}}=40.2\%$ at $l=20\mathrm{mm}$, only $0.4$ percentage points above the $39.8\%$ baseline at $l=9\mathrm{mm}$. This saturation is expected because the double-pass effective path ($18\mathrm{mm}$) exceeds the interaction length

$$l_0=\sqrt{{\epsilon }_0{c}^2\pi {\omega }_f^2{n}_o^{\omega }{n}_e^{\omega }{n}_e^{2\omega }/\left(4d_{\mathrm{eff}}^2{P}_F\right)}\approx 17\mathrm{mm},$$

(2)

as shown in the right panel of Figure 4 of [28]. Each sweep point runs in 25–35 minutes; a full $l\in [6,20]\mathrm{mm}$ scan completes in 4–6 hours.

Setting $R_{f,\mathrm{FW}}=R_{f,\mathrm{SHW}}=0$ reduces the model to the single-pass limit, recovering published single-pass curves to within truncation error ($<0.1\%$).

A constant-versus-variable conductivity comparison provides a further consistency check. Using

$$K\left(T\right)=K_0/(T_0+\Delta T),$$

(3)

with $K_0=13\mathrm{W}/(\mathrm{m}·\mathrm{K})$ and $T_0=298\mathrm{K}$, the model predicts ${\eta }_{\mathrm{SHG}}=39.8\%$ and $T_{\max }=303.5\mathrm{K}$. Fixing $K=K\left(T_0\right)$ gives $41.3\%$ and $302.8\mathrm{K}$, overestimating efficiency by $1.5$ percentage points and underestimating temperature by $0.7\mathrm{K}$. This confirms that constant-conductivity assumptions introduce systematic errors even for moderate temperature rises ($\Delta T\sim 5\mathrm{K}$). This test runs in under 45 seconds with coarse meshing ($N_r=100$, $N_z=50$, $50,000$ time steps).

4. Limitations

The wave solver assumes CW Gaussian beams with constant waist, ignoring diffraction-induced spreading. The approximation holds for $l\ll z_R=\pi {\omega }_f^{2}n/\lambda$; at ${\omega }_f=50\mu \mathrm{m}$ and $1064\mathrm{nm}$, $z_R\approx 7.3\mathrm{mm}$, so the $9\mathrm{mm}$ validation crystal lies near the validity boundary. Non-Gaussian or vectorial beam profiles must be computed externally and coupled through modified source terms. The materials database is limited to KTP in the type-II configuration at $1064\mathrm{nm}$; other crystals (LiNbO₃, BBO, LBO, periodically poled variants) require user-supplied property fits.

No GPU acceleration is provided. The finite-difference heat solver uses optional OpenMP spatial parallelization but remains CPU-bound. Memory bandwidth limits performance for meshes with $N_r\times N_z>{10}^6$. Only x86-64 processors have been validated.

CFL-violating time steps ($\Delta t>\Delta r^2/\left(2\alpha \right)$) produce divergent temperature oscillations; halving $\Delta t$ restores convergence. At extreme parameters ($P_F>200\mathrm{W}$ or ${\omega }_f<20\mu \mathrm{m}$), TIPM can exceed $-5\mathrm{rad}$ and the iterative coupling may fail to converge within the default 100 iterations.

5. Code Availability

Source code, documentation, and example configurations are available at GitHub. The release corresponding to this paper is tagged v1.0.2 and archived on Zenodo (DOI 10.5281/zenodo.17362470) under the MIT License.

6. Conclusions

An eight-equation coupled model for CW double-pass type-II SHG in KTP has been implemented and validated against experimental data, with conversion-efficiency deviations below $4\%$ over the tested range of pump powers and beam waists. Restructured memory management reduces RAM requirements by approximately $60\%$ relative to the original research code, allowing execution on standard desktop hardware. The implementation is openly available with version control and standardized inputs to support reproducibility.