1. Introduction

Modern optical nonlinear crystals when used with high peak power pulsed lasers allow efficient frequency conversion in a single pass through the crystal [1,2]. The efficient nonlinear interaction causes distortion of the optical wave fronts of the participating beams due to depletion, and the light is consequently diffracted away from the incident modes. The result is a change in the mode composition of the interacting beams and reduced efficiency of the conversion process. The problem has been treated in a plane wave approximation, where the conversion was calculated for each part of the wave front, but diffraction was not included [3]. Other calculations have employed a numerical wave propagation method, where the crystal was divided into thin slabs and the nonlinear differential equations describing the interaction was approximated by difference equations. After passage of each slab the beams were expanded in plane waves, and both nonlinear interaction and diffraction were included in the calculations [4]. This method allows calculation of the generated beam shape and of the efficiency of the conversion. In [5] general shaped input beams were synthesized from orthogonal Gauss-Laguerre modes to find the most efficient beam shape for sum- and difference frequency generation. However, nonlinear effects such as depletion and phase shift were not included, and the method only applies to very low power beams.

In this paper an alternative formulation is used, where the beams are expanded in higher order Gauss-Hermite (G-H) or Gauss-Laguerre (G-L) modes after passage of each thin slab. The choice of mode expansion depends on the symmetry of the problem. Since both G-H and G-L modes form complete orthogonal sets of functions, an expansion of the beams in these modes fully describes the result of the interaction, and diffraction due to the distortion of the wave front is included by the mix of the higher order modes. Besides forming an alternative way of performing a wave propagation calculation, this method has the advantage that it allows to follow the generation of higher order modes both spatially and temporally through the crystal. Furthermore it is possible to calculate the power generated in each of these modes. It is also possible to input higher order pump modes or superposition of modes and calculate the power conversion efficiency into the generated modes, a method that has recently been used in generation of squeezed light in higher order G-H modes [6]. In the calculations the phase matching may be adjusted at will; it is for example possible to adjust the phase matching to compensate at least partially for the Gouy phase shift in a particular mode and thus selectively generate light in one mode by pumping in another mode [7].

In the following, single pass second harmonic generation (SHG) in an efficient nonlinear crystal pumped by a pump beam of known mode composition is considered. As an example we shall look at a single mode Q-switched Nd:YAG laser pulse frequency doubled through a PPKTP crystal and investigate the formation of higher order G-L modes as the fields propagate through the crystal, as well as illustrate the temporal development of the pulse shape as a function of time.

2. Wave propagation with G-L expansion

The calculations are performed in the following way:

1. The crystal is divided into thin slabs normal to the optical axis.

2. The three interacting beams that make up the total field incident on a slab interact on their way through the slab. The slab is assumed so thin that the nonlinear differential equations governing the interaction (the coupled wave equations) can be approximated by difference equations, and diffraction effects on the way through the slab are unimportant (plane wave approximation within a slab). Losses can be included in each slab. The beams may be focused anywhere inside or outside the crystal. Collinear beams are assumed, although the simulations can easily be expanded to include non-collinear beams and walk-off.

3. After passage of a slab, the field incident on the slab has been distorted due to the nonlinear interaction. The fields leaving the slab are expanded in G-L modes.

4. The beam waist radius and phase of the individual modes of the expansion are recalculated for the position of the next slab. The different phase shifts of the different order G-L modes implicitly take care of diffraction.

5. The G-L modes with the new phases and beam waists are added for each beam, and these new beams are used as input for the next slab; the calculation is repeated from 2.

6. The output beams are found by summing the G-L modes leaving the last slab.

2.1. Nonlinear interaction

The coupled wave equations (CWE) for the 2^nd order nonlinear interaction between three beams, described by the normalized electric fields a_j(x,y,z) where j identifies the field, are in the difference equation form given by [8]

where n enumerates the slab and Δz is the thickness of the slab. g is the nonlinear gain coefficient, and Δk is the phase mismatch.

After passage of the slab, the fields are updated with the change,

2.2. Gauss-Laguerre expansion

Assuming cylindrical symmetry of the interaction and a cylindrically symmetric pump wave the fields are most efficiently described as G-L modes, u ^G-L _j,p,l(r,θ,z,t):

where w _0,j is the beam waist radius of the j’th beam and $w{\left(z\right)}_j=w_{0,j}\sqrt{1+{\left(\frac{z}{z_{0,j}}\right)}^2}$ is the beam radius at coordinate z with z _0,j=πw ² _0,j/λj being the Raleigh range and λ_j=_0,j λ/n(ω_j) and λ _0,j are the wavelengths of the beams in the medium and in vacuum, respectively. n(ω_j) is the index of refraction of the medium at angular frequency ω_j and R(z)_j=z(1+(z _0,j/z)²) is the local radius of curvature of the wave front. Further ($r=\sqrt{x^2+y^2},\theta =\arctan \left(\frac{y}{x}\right),z$) are cylindrical coordinates, while k_j=2π/λ_j is the wave number.

L^l_p(r) is the Laguerre polynomial of order p and index l. The normalization constants for the G-L modes are given by

where δ(l)=1 for l=0 and 0 otherwise. In the following it is assumed that there is no azimuthal dependence of the pump field, which means that l=0. By symmetry the generated modes will also have index l=0. The normalization constants then become equal for all p modes:

Diffraction is included in the expansion through the different phase dependence of the G-L components. The individual phase dependence of the different G-L modes is given by the so-called Gouy phase term: ${\phi }_{G,j}\equiv \left(2p+l+1\right)\arctan \left(\frac{z}{z_{0,j}}\right)$ and by the azimuthal phase term lθ. The quadratic phase term k_jr ²/2R(z)_j is the same for all G-L modes as are the propagation phase term (k_jz-ω_jt) and the phase mismatch term Δkz. The latter two are both independent of r and θ.

2.3. Computational procedure

The three fields at the position of the n’th slab are first allowed to interact during the passage of the slab. The field increment is calculated using difference equations and plane wave propagation. After passage of the slab the three fields are expanded in the G-L modes, and the expansion coefficient for each p - mode is calculated:

where the factor 2π stems from the θ-integration. z_n is the z-coordinate of the n’th slab.

The G-L modes making up the beam after leaving the slab are given by

After calculation of the mode components, the beam waist and phase for each mode function is updated to the values corresponding to the position at the end of the slab. The resulting beam after passage of the slab is found by summing over all modes up the highest order, pmax:

The procedure is then repeated for the next slab.

The output beams leaving the crystal are simply found as the beams exiting from the last slab:

where N is the number of slabs.

2.4. Second harmonic generation

In the following we consider SHG and the three interacting beams are designated the fundamental (F) beam (two identical fields, j=1, 2) and a generated second harmonic beam (SH): (j=3). Since the beam waist radius and position of the generated SH so far have not been defined (an arbitrary basis for the expansion of the SH beam can be chosen), a basis has to be assigned. The judicious choice is to choose beam waist parameters so that the fundamental pump mode couples to the fundamental mode of the generated SHG field; this is achieved by using the same location for the beam waists and choosing w _0,3=w ₀₁/√2. This choice results in the same centre and radius of curvature of the wave fronts for beams of all orders and will result in the simplest expansion of the generated SH beam in terms of G-L modes.

2.5 Power in the G-L-modes and efficiency of the SHG process

It is now possible to find the power in each of the modes as they propagate through the crystal

as well as the total SH output power by summing over all modes included in Eq. (8).

where P ⁽ⁿ⁾ _j,p,0 designates the power in mode p,0, and P _{(out) j} is the total power of the j’th beam at the end of the crystal.

The numerical simulations in section 3 are compared to the results obtained by Boyd-Kleinman [9]. As the B-K approach does not include depletion, a substantial difference is seen between our calculations and the B-K-theory for high incident power.

According to the B-K formula, the efficiency of SHG is given by Eq. (12).

where l_c is the length of the crystal, h(σ,B,ξ) is the B-K factor, σ=bΔk/2 is the phase mismatch parameter, $B=\rho \sqrt{\frac{l_c{k}_1}{2}}$ is the walk-off parameter, and ξ=b/z _0.1 is the focusing parameter. Here Δk is the phase mismatch, ρ is the walk-off angle, and b=2πw ² _0.1/λ ₁ is the confocal parameter.

2.6. Temporal pulse shape

The distortion of the temporal pulse shape may be found simply by stepping the incident pump field through incremental time steps. The method assumes that the dynamic response of the nonlinear process to the incident field is so fast that the internal nonlinear process follows the pulse shape, or in other words, an adiabatic assumption has been made regarding the response of the crystal to changes in the pump beam. This is undoubtedly true for an incident Q-switched laser pulse and probably also for ps pulses, but the assumption is not warranted for fs pulses.

The calculation as outlined above provides both the transverse spatial development of the pulse through the crystal and the temporal response to the pulse shape by plotting the output field as a function of time. Furthermore, the generated pulse energy in each mode can be found integrating over time.

3. Results of the numerical simulations

In this section some calculated results assuming two identical incident Gaussian beams are presented. In the nomenclature of Eq. (13) the Gaussian beams are given by

where z_i is the coordinate of the first slab. The following parameters were used for the nonlinear crystal: Crystal length l_c=1 mm or l_c=10 mm, passive losses α=0.01, beam waists between w _0,1=0.100 mm and w _0,1=0.019 mm, crystal refractive index n=1.7, and nonlinear coefficient d_eff=10 pm/V. These values corresponds approximately to single-pass SHG of a Nd:YAG laser using PPKTP as the nonlinear material.

The gain in the CWE is found to be: $g=2d_{\mathrm{eff}}\sqrt{\hslash {\omega }_1^3{\eta }^3}$, where $\eta =\frac{{\eta }_0}{n}=\sqrt{\frac{{\mu }_0}{{\epsilon }_0}}\frac{1}{n}$.

3.1. Weakly focused SHG

The weakly focused or quasi-plane wave regime corresponds to a situation, where the focusing parameter is small, and there is negligible phase mismatch throughout the crystal. Figure 1 shows the incident pump beam (a) and generated SH modes (b) at the end of the crystal (online movie shows the development of the fields as they propagate though the crystal). Assuming that the pump power is not high enough to cause significant depletion, the Gaussian shape of the pulse is not significantly distorted, but the power of the SH beam is low.

 

Fig. 1. Quasi-plane wave case, low power and short interaction length (1 mm); a-b: Fundamental and SH power (1 mm nonlinear crystal) and c-d: Fundamental and SH power (10 mm nonlinear crystal). Movie shows the field developing through the crystal.

Download Full Size | PDF

Still considering the weakly focused regime, but increasing the interaction length from 1 to 10 mm, allows some depletion of the incident pump beams. A significant generation of higher order G-L-modes is seen, and the pump fields at the output display the characteristic doughnut feature also reported in [3]. The transmitted pump beam and the generated SH beam at the output of the nonlinear crystal are seen in Fig. 1(c) and 1(d), respectively. The on-line movie shows the development of the fundamental and the SH as they propagate through the crystal. The incident fundamental power in Fig 1(c) and 1(d) was chosen high enough so that the emergence of back conversion of SH to F could be seen.

 

Fig. 2. The most powerful higher order G-L-modes at the output of the 10 mm long nonlinear crystal, using the same parameters as in Fig.1© and 1(d).

Download Full Size | PDF

The individual modes containing most power at the output of the nonlinear media are shown in Fig. 2. Parameters are the same as in Fig. 1. It is interesting that rather few modes (p _max<6) are needed to describe the distorted field even in the case of significant depletion.

 

Fig. 3. Development of power in F and SH beams through the crystal with optimal B-K focusing parameter ξ=2.84. Pump power is kept low to avoid back conversion.

Download Full Size | PDF

a) Zero phase mismatch, Δk=0, for the fundamental G-L modes.

b) Phase mismatch Δk=218 m^-1 to increase overall SH conversion efficiency.

3.2. Strong focusing, high nonlinear interaction and optimum phase mismatch

In case of tight focusing, many higher order modes are generated. Since the higher order modes all have different Gouy phase shifts, the conversion into different orders can be influenced by a careful control of the phase matching. The effect of changing the phase mismatch is clearly seen in the figure.

 

Fig. 4. Development of power in F and SH beams through the crystal with optimal B-K focusing parameter ξ=2.84. High pump power is considered to allow for back conversion.

Download Full Size | PDF

a) Zero phase mismatch, Δk=0, for the fundamental G-L modes.

b) Phase mismatch Δk=218 m^-1 increases overall SH conversion efficiency.

The phase match between different beams of different frequencies is typically controlled by slight angle tuning of a birefringent crystal or a change in temperature of the nonlinear media. In Fig. 3(a) and 3(b) two cases are shown, one in which the phase mismatch is zero for the fundamental orders and one, where the phase matching has deliberately been tuned to Δk=218 m^-1. This results in higher overall conversion. Note however, that this does not necessarily mean maximum conversion into the fundamental SH mode.

The phase mismatch of Δk=218 m^-1 corresponds to the phase mismatch denoted “the optimum phase mismatch” in the original paper by Boyd and Kleinman [9]. When depletion occurs, we find that greater overall conversion can be achieved with even higher phase mismatch. Figure 4 shows the same calculation at higher pump power, resulting in significant back conversion. It is seen from Fig. 4, that a crystal of about 5 or 6 mm length would be optimum depending on the phase mismatch.

3.3. Pulse distortion

The depletion of the spatial distribution also leads to a distortion of the temporal shape of the pulse. Assuming an initial Gaussian pulse shape, the distortion of the pulse can be calculated stepping the incident pump pulse through a series of time steps. As a first example, a weakly focused pump pulse, ξ=0.10, is considered in Fig. 5(a), (b) and (c). Figure 5(a) and (b) show the fundamental and second harmonic transverse field distribution at the peak of the pulse. The movie shows the fields as a function of time. Phase mismatch was set to zero. Figure 5(c) shows the time evolution of the power in the fundamental pump pulse, the transmitted fundamental pulse, and the generated SH pulse. Figure 6 shows the spatial profile of a) the fundamental and b) the SH at the peak of the pulse in the case of a strongly focused beam, ξ=2.84. The pump is strongly depleted and the conversion efficiency is higher than in the case of the weakly focused beam. ξ=2.84 corresponds to the optimum B-K focusing parameter for SHG when there is no walk-off. Figure 6(c) shows the time development of power.

 

Fig. 5. The transverse field distribution of a) the fundamental input pulse and b) the generated SH in the middle of the pulse, ξ=0.10. Movie shows time development of these fields. c) shows the development of power in pump input pulse, transmitted pump pulse and generated SH output as a function of time. The peak pump power is 400 W.

Download Full Size | PDF

 

Fig. 6. The transverse field distribution of a) the fundamental input pulse and b) the generated SH in the middle of the pulse, ξ=2.84. Movie shows time development of these fields. c) shows the development of power in pump input pulse, transmitted pump pulse, and generated SH output as a function of time. The peak pump power is 200 W.

Download Full Size | PDF

3.4. Comparison to Boyd-Kleinman theory

It is often interesting to know for what crystal length or pump power the pump beam gets significantly depleted, and how this affects the conversion efficiency, both the total conversion efficiency and the efficiency of conversion into the fundamental mode of the SH. It would be interesting to do this for the optimum focusing parameter according to the Boyd-Kleinman theory. The Boyd-Kleinman theory does not take into account the depletion of the pump wave and the diffraction due to mode distortion of the waves caused by the nonlinear interaction. The present model takes these effects into account, and the difference between the present model and the Boyd-Kleinman theory at higher powers will be considered in the following. Figure 7 shows a comparison between the present results for the conversion efficiency from fundamental pump wave to the SH and results obtained using the Boyd-Kleinman factor. The optimum focusing factor is ξ=2.84 according to the B-K approximation, and the corresponding optimum phase mismatch, Δk=218 m^-1, were used in these simulations. Furthermore, zero passive losses and focusing at the centre of the crystal were used. The results show that the B-K theory significantly overestimates the single-pass conversion efficiency already at conversion efficiencies around 20 %. Since the B-K theory does not include depletion, the results of course become absurd as the power or crystal length increases, resulting in conversion efficiencies exceeding 100 %. At low power levels, however, there is perfect agreement between the two methods.

 

Fig. 7. Comparison of SHG conversion efficiency using Boyd-Kleinman theory and the G-L mode expansion. Same parameters as in previous figure.

Download Full Size | PDF

4. Conclusion

A wave propagation calculation in which the three optical fields interacting in a nonlinear crystal are expanded in Gauss-Laguerre modes as they propagate through a nonlinear crystal has been developed and applied to both CW beams and a pulsed pump beam. The expansion in Gauss-Laguerre modes requires only few higher order modes to account for the distortion of high power beams resulting from the nonlinear interaction and the resulting diffraction into higher order modes. The calculation accounts for both nonlinear depletion, diffraction into higher order modes and deliberate or accidental phase mismatch due to e.g. temperature deviations or angular misalignment. The method allows one to go beyond the traditional Boyd-Kleinman theory to see the effect on conversion efficiency due to distortion of the beams and to calculate how much power is converted into the desired fundamental mode. The method also allows one to follow the development of the higher order modes as they propagate through the crystal. Since the Gouy phase shift due to the focusing of the beams depends on the mode order, it is possible by deliberately adjusting the phase mismatch to compensate at least partially for the Gouy phase shift in a selected mode and thereby optimizing the conversion efficiency into that mode.

The method also allows one to follow the temporal development of a pulse. This is illustrated with the example of a Q-switched pulse with an initial Gaussian temporal pulse shape. At high powers a typical “temporal hole burning effect” is seen in the shape of the pulse near the centre of the pulse.

The method would also allow treatment of crystals in resonators by feeding some of the output from the crystal back to the input. Walk-off or off-axis propagation can be accounted for by including the azimuthal modes or expand in Gauss-Hermite modes.