Stretching and compressing of short laser pulses by chirped volume Bragg gratings: analytic and numerical modeling

Abstract. The theory of stretching and compression of short light pulses by chirped volume Bragg gratings (CBGs) is reviewed based on spectral decomposition of short pulses and on wavelength-dependent coupled wave equations. Analytic theory of diffraction efficiency of CBG with constant chirp and approximate theory time-delay dispersion is presented. Comparison of approximate analytic results with exact numeric coupled-wave modeling shows excellent agreement for smooth heterogeneities of CBGs.


Introduction
Chirped pulse amplification 1 is the main approach for highpower/high-energy ultrashort laser pulse' generation. It is achieved by stretching ultrashort laser pulses. Those stretched pulses have a much lower value of instantaneous power and, therefore may be amplified by broadband laser amplifiers without damaging the medium. Amplified chirped pulses may then be directed to the compressing element, which allows collection of all spectral components back into ultrashort pulse, but with higher energy. The conventional way for doing this is the use of a pair of diffracting gratings for both stretching and compression (Treacy stretchers and compressors 1 ). An especially elegant idea was suggested and implemented by Galvanauskas et al. 2,3 with the use of chirped fiber Bragg gratings, which are gratings with a gradually varied period along the fiber. It is possible to use the same grating for compression, and stretching stage, but to illuminate it from the opposite end. In this case, the influence of smooth heterogeneities of time-delay dispersion (TDD) at the stretching stage, T str ðλÞ, is compensated for by those of the compression stage, so that T str ðλÞ þ T compr ðλÞ ≈ const. A similar performance was demonstrated with the use of chirped volume Bragg gratings (chirped VBGs or CBGs), which are produced by holographic recording in the bulk of photo-thermo-refractive (PTR) glass and have dramatically higher apertures when compared with those of fiber gratings. 4 Recent development of VBGs based on PTR glass 5 allowed for operation of stretching-compression schemes at much higher values of power, see for example, Ref. 6.
The present paper is devoted to the development of numerical and analytic tools for the study of stretching and compression by volume CBGs. While this analysis is applicable for any volume CBG, the examples shown in this paper are based on CBGs recorded in PTR glass.

Basic Scheme and System of Equations,
Definition of Chirp, and Time-Delay Dispersion The basic scheme of stretching, amplification of stretched pulse, and subsequent recompression back into short pulse by means of a volume CBG is presented in Fig. 1. Due to gradual variation of CBG period in the z direction, different spectral components of an incident pulse are reflected from different parts of the CBG and, therefore, have different delays. After amplification, the stretched pulse is launched to the same CBG from the opposite side and compressed back to its original width. In this work, we consider the amplifier as a linear device that does not affect any parameters of a laser pulse but power. Therefore, we assume that a CBG-reflected stretched pulse is transmitted from the cross-section at the front surface of the CBG (z ¼ 0) to the opposite end of the same CBG (z ¼ L), without any additional distortions.
We assume that the dielectric permittivity and magnetic permeability of CBG depend on coordinate z as where δnðzÞ is the variation of refractive index. with k ¼ ωn∕c; n 0 ≡ nðωÞ is where the terms of the order of OðδnÞ 2 were ignored. We also assume that the spatial variation of the refractive index δnðzÞ has the form δnðzÞ ¼ n 2 ðzÞ þ n 1 ðzÞ cos½Qz þ ψðzÞ: Here, smooth heterogeneity n 2 ðzÞ and spatial modulation n 1 ðzÞ are small in comparison with n 0 and are relatively slow varying functions of z; ψðzÞ is a differentiable function, and Q (rad∕m) has the value close to 2k, see Eq. (13). (It should be emphasized that our notation n 2 ðzÞ for the smooth correction to the background refractive index has nothing to do with the nonlinearity of refractive index.) In the standard approach of slow varying envelope approximation (SVEA), one should ignore second derivatives of slowly varying amplitudes AðzÞ and BðzÞ, keep the resonant terms ∝ expðikzÞ and ∝ expð−ikzÞ only, and equalize each of them to zero. As a result, one gets the following system of coupled equations: It is convenient now to introduce new amplitudes aðzÞ and bðzÞ, defined by Then for these small amplitudes aðzÞ and bðzÞ, one gets our system in the form Here νðzÞ is the local value of detuning, measured in (rad∕m), and κðzÞ is the local strength of coupling, of dimensions (1∕m) κðzÞ ¼ ω c n 1 ðzÞ 2 : Possible background absorption in the medium with coefficient αð1∕mÞ, for intensity IðzÞ ∼ expð−αzÞ, may be described by purely imaginary n 2 ðzÞ, namely n 2 ðzÞ ¼ iλα∕4π.
Solution of the linear system of ordinary differential equations [Eq. (7)] may be represented as the solution of the Cauchy problem, where both boundary values are known at the same point, e.g., at z ¼ 0. Linearity of the system allows writing. Mð0Þ Actual boundary conditions for the reflection problem are usually prescribed at the opposite ends of CBG; for example where we assumed unit amplitude of incident wave a at the beginning z ¼ 0 of CBG and no incident wave b at the end z ¼ L of CBG. Then the reflection and transmission amplitudes are found by imposing the conditions in Eq. (10).
The determinant of the matrixMðzÞ from Eq. (9), which is present in the numerator of Eq. (11), is unity as the consequence of Eq. (7), even in the presence of absorption. This allows additional checking of the accuracy of numeric integration of those equations. The amplitude of reflection and transmission for the opposite boundary conditions may be expressed through the elements of the same matrixMðLÞ.
which is especially convenient for the case when one uses the same CBG both for stretching and for recompression. Subsequently, we neglect the n 2 ðzÞ term. Then, equating the local value of detuning νðzÞ to zero, one can find the position z B ðλÞ of the point of exact Bragg resonance for the given vacuum wavelength λ.
Here λ 0 is the central wavelength defined by basic parameter of CBG: Q ¼ 4πnðλ 0 Þ∕λ 0 , with Q defined by Eq. (4). This condition allows connecting spectral chirp parameter dλ∕dz B of CBG with the second derivative of the phase correction ψðzÞ. Indeed, one has to assume λ ¼ λðzÞ in Eq. (13) and take z-derivative of the left-and right-hand sides of that equation. As a result, one gets Here the expression for n gr is an alternative form of denoting group velocity, v gr ¼ dω∕dk ¼ c∕n gr and n gr ¼ nðωÞþ In the approximation of unchanged group velocity v 0 ¼ c∕n gr , the value of time-delay in the stretching process as a function of wavelength λ is We will see below ( Fig. 7 in Sec. 5.3) that this formula yields a reasonable result for CBG with low efficiency only.
3 Analytic Expression of Diffraction Efficiency of CBG Consider first a weak CBG for which the first order of perturbation theory is valid. In the zeroth approximation and at n 2 ≡ 0, one can take AðzÞ ≈ 1 for all z. Then for the reflection amplitude coefficient, one gets from Eq. (5) which, for ψðzÞ ¼ β · ðz − 0.5LÞ 2 , yields explicit value of An analytic expression has been derived in Ref. 8, which is valid for any strenth of coupling. It may be simplified for constant chirp profile of the grating, with the use of Eq. (15) to: Derivation of this expression, see Refs. 7 and 8, is similar to the calculation of quantum-mechanical transmissivity of parabolic potential barrier, see Ref. 9, and will not be discussed here. To go over the rather heavy derivation of Eq. (19) from Refs. 8 and 7, we suggest using the known structure of Eq. (19) and checking all the coefficients in it via the first-order perturbation result [Eq. (18)]. Actual numerical modeling (see below) confirmed the validity of Eq. (19) with great accuracy, especially for apodized CBG, where spatial refractive index modulation comes smoothly to zero at the ends of CBG in this particular modeling, see Eq. (40).

Approximate Expression for Time-Delay
Dispersion Approximate evaluation of TDD may be done on the basis of assumption of slowly varying behavior of the coefficients κ ¼ kðzÞ and ν ¼ νðzÞ in Eq. (7). We take the SVEA equations in the symmetric form.
In the assumption of ν ¼ const; κ ¼ const, the solution of the form ða; bÞ ¼ expðpzÞða 0 ; b 0 Þ exists. From that we get the result for two eigenvalues of matrixV.
We define the positive root as the one for which, in the region jκj < jνj, the real square root μ ¼ þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ν 2 − jκj 2 p satisfies the condition μ∕jμj ¼ ν∕jνj; in other words, the sign of μ is the same as the sign of ν in that region where propagation is not forbidden. Eigenvectors for p 1 ¼ þiμ and for Group velocity of the wave a in the absence of coupling, jκj ¼ 0, is taken as v a ¼ v 0 , while group velocity of the wave b is taken as v b ¼ −v 0 . The general expression for the group velocity in the presence of both constant values ν and κ for a given mode may be intuitively written as v gr ¼ P½Watt∕m 2 ∕U½Joule∕m 3 ; (23) where P is the Poynting vector and U is the energy density.
In the approximation of constant (or reasonably slow varying) ν and κ, one gets for mode #1 (a predominantly a-wave) and for mode #2 (a predominantly b-wave) We take those expressions as approximations for group velocities of incident and reflected waves in the use of CBG for the stretcher-compressor scheme and verify the degree of their validity by comparison with results of direct wave modeling of CBG, Sec. 5.
In this approximate approach (Fig. 2), the incident wave a (the one to be stretched by CBG) propagates from the point z ¼ z 1 of the left boundary of CBG to the left point z ¼ z 2 inside CBG, where decisive reflective transformation of a into b takes place. The points z ¼ z 2 and z ¼ z 3 are defined by the condition νðz 2;3 ; λÞ ¼ AEjκj. In CBG approximated by constant chirp, νðz; λÞ ¼ β · ½z − z B ðλÞ, those points are symmetrically positioned around z B ðλÞ, i.e., around perfect Bragg-resonant point.
After the decisive reflection from the point z ¼ z 2 , incident wave of type a gets transformed into reflected wave of type b. Total time-delay due to propagation from the entrance z ¼ z 1 to reflection point z ¼ z 2 and back is where vðzÞ ¼ v 0 · jμðzÞ∕νðzÞj as a consequence of Eqs. (22) and (24). At the compression stage, the wave enters CBG through the cross-section z ¼ z 4 and propagates in the bform down to z ¼ z 3 , and after the reflection, it propagates back to z ¼ z 4 . The corresponding round-trip delay time of this stage is given by where again vðzÞ ¼ v 0 · jμðzÞ∕νðzÞj. For a general z-dependence of νðzÞ and κðzÞ, one should numerically find z B ðλÞ; z 2 ðλÞ; z 3 ðλÞ and numerically calculate the corresponding integrals. Calculation can be done analytically for constant chirp and constant κ.
In the case when both z 0 ≪ z B ðλÞ and z 0 ≪ L − z B ðλÞ, one can expand T cycle ðλÞ in terms of small ratio z 0 ∕L. (32) One should take the value of T cycle from Eqs. (31) or (32) above and then average the expression over the spectral content of the pulse in question. As an even more crude estimate, one may try to take the value of T cycle at jz B − L∕2j ≪ z, and then the formula for time-delay of the peak of recompressed pulse becomes It is worth discussing separately two effects that influence the delay time of the cycle stretching-recompression.
The first one is that the round-trip length of the stretching process is shorter than 2½z B ðλÞ − z 1 by the thickness of forbidden zone: 2z 0 ¼ 2jκ∕βj, see Fig. 2. This leads to a shorter delay time of the stretching-recompression cycle: shorter by approximately ð4z 0 ∕v 0 Þ, i.e., of the first order in z 0 .
The second effect is due to considerable (∼50% and more) decrease of group velocities in the vicinity of reflection points z ¼ z 2 and z ¼ z 3 . The thickness of this vicinity is of the order of z 0 . Thus, this second effect results in a longer delay time, also of the first order in z 0 .
What is truly remarkable is that these two effects compensate each other in the first order in z 0 . The resulting delay time does decrease (in comparison with 2L∕v 0 ), but in the second order in the coupling constant, i.e., proportionally to z 2 0 . A monochromatic wave stretches in time from t ¼ −∞ to t ¼ þ∞. One can define time delay dispersion TDDðλÞ for the quasi-monochromatic packet with the wavelength λ. Many sources advise to one calculate the phase of reflection coefficient φðλÞ ¼ argðrðλÞÞ and postulate that Here L is the total thickness of CBG; z ¼ z 2 is the point where the input wave a of the given frequency ω ¼ 2πc∕λ. hits the region forbidden for propagation (marked by curved dashed line).
dφðλÞ dλ : (34) Our numeric modeling shows that this is valid for large modulation of φðλÞ only, jφj ≫ 1.

Parameters of Numerical Modeling
We have developed the program for numerical modeling of the scheme presented in Fig. 1 with the use of Mathematica software package. In our approach, the fields of all pulses were decomposed into time Fourier series via discrete Fourier transform (DFT) subroutine. The discrete index j in DFT programs takes only non-negative values in the range 0 < j < N − 1, where N is the total number of points (either in time domain or in frequency domain). For that reason the frequency of an individual component was related to that index j by In that manner, we account for both positive and negative values of frequency detuning. Meanwhile discrete time points were numbered just as Here T total is the total time interval under consideration. Average refractive index of PTR glass was taken for each wavelength λ j ¼ 2πc∕ω j as either constant value n ¼ 1.4891 (for central wavelength λ central ¼ 1.06 μm) or from Sellmeier formula for that glass. 10 The system of ordinary differential equations was integrated numerically for each frequency ω j . Sometimes it was advantageous to divide all the integration length L (thickness of CBG) into four separate regions, with accuracy goal achieved for each region independently. Examples below are demonstrated mostly for L ¼ 0.1 m (≡10 cm) and constant chirp parameter β ¼ 8.10 4 rad∕m 2 from Eq. (15). Value of coupling constant jκj was chosen via dimensionless parameter S ¼ jκjL, with S varying from 15 to 150.

Modeling of Stretching-Compressing by CBG with Imperfections
It should be noted that CBGs with large apertures and long stretching times 4,6,11 show some spatial variations of phase resulting from optical aberrations in a holographic recording system and optical homogeneity of a recording material (PTR glass). This is why, among other results of our numerical modeling of stretching-compression process, we would like to discuss here the influence of small and very inhomogeneous variations of phase of the grating. Equation (13) allows finding the position z B ðλÞ, where formal Bragg condition is satisfied for a given wavelength λ. We took a small oscillatory addition to the phase, δψ ¼ a cos γz, where γ½1∕meter ¼ 2π∕Λ 1 and Λ 1 is the period of perturbations.
Since the term dψ∕dz is present in the right-hand side of Eq. (13), the equation for detuning z B ðλÞ becomes If aγ 2 is larger than 2jβj, then there is a possibility that Eq. (37) has several solutions ðz B Þ 1 ,ðz B Þ 2 , . . . for the given wavelength λ. Then application of the simplest formula yields multivalued TDDðλÞ and gives a suspicion of lowquality recompression, to say nothing about peculiar oscillations, which are even more severe in the stretched pulse. However, the results of our numerical modeling show for phase modulation that if the amplitude 2a (top to bottom) is moderately small, e.g., 2a ≤ 0.4 rad, the influence of formally multivalued feature of TDDðλÞ from Eq. (38) is not important. In particular, Fig. 3 shows z-dependence of Bragg resonant wavelength detuning (in nanometers) for particular perturbation of the phase given by Eq. (39). Actual wave modeling of stretching by perturbed CBG was done for L ¼ 0.1m and β ¼ 8 · 10 4 rad∕m 2 , coupling coefficient jκj ¼ S∕L, and S ¼ 30.
Additional modulation δψ from Eq. (39) was taken with the parameters 2a ¼ 0.4 rad, γ ¼ 1405 rad∕m. so that aγ 2 ∕ 2β ¼ 2.47, and the condition of multivalued function TDDðλÞ from Eq. (38) was satisfied. Stretching and compression of incident Gaussian pulse with τ 0 ðHWe −2 IMÞ ¼ 3.0 ps was modeled and depicted in Fig. 4. It shows spectra of the incident pulse (upper curve) and of the recompressed pulse (lower curve). One can see that for a CBG with spectral width ΔλðFWÞ ≈ 0.1 nm equal to that of a laser pulse at the level of 1∕e 2 , the spectral width of a recompressed pulse is practically the same. However, there is some modulation of the spectrum in the vicinity of the maximum, which is caused by ripples in the dispersion curve depicted in Fig. 3; cutting off short-and long-wavelength wings of the spectrum and the spectral width of the CBG is due to finite spectral width of CBG. Figure 5 shows the intensity profiles of the input pulse, stretched pulse (multiplied by factor 80 for illustrative There is a small precursor in the recompressed pulse (arrow in Fig. 5) containing~6% of recompressed energy. Nevertheless one observes rather good quality of recompression: diffraction efficiency of full cycle of stretching-compression for energy was 0.92 and ratio of peak intensities was 0.61. It means that steep perturbation of phase δψ with amplitude 2a≲0.4 rad is not harmful for recompression. Meanwhile, our modeling for perturbations with amplitude 2a≳0.4 rad had shown considerable hindering recompression.

Study of Time-Delay Dispersion
For the evaluation of our analytic (but approximate) results for time-delay dispersion, TDD(λ), we had to choose a special approach in numerical modeling. Namely, illumination of CBG by a very short incident pulse (< 1 ps and therefore with a very broad spectrum) resulted in a stretched pulse with very flat top, for which it was difficult to find TDD at the stretching stage. If any distortions in the temporal profile of stretched pulse were present, the arrival time of the peak of reflected pulse was difficult to determine. On the other hand, if the incident pulse was relatively long (and hence was relatively narrow-band), the reflected pulse was not stretched to full possible duration 2L∕v gr . Here v gr is the group velocity at λ 0 in unexposed glass. Still, determining the delay time for stretching stage was difficult due to large duration of the input pulse.
For that reason we calculated from numerical wave modeling the delay time of a recompressed pulse in comparison with the arrival time of an incident pulse. The delivery of the stretched pulse to the back of CBG for recompression was considered to be instantaneous in the modeling. Two types of CBG were considered. One type was with jκj ¼ S∕L taken as a constant through the whole thickness L of CBG; hereafter, we call this type of CBG as uniform or nonapodized. The other type had the profile This apodization function suppressed not-quite Bragg contributions of the ends of CBG; the latter contribution leads to Fresnel function-like oscillations in the reflection spectrum of nonapodized CBG. For both cases we calculated TDD (at λ 0 ¼ 1.06 μm) for different values of dimensionless parameter S (and for constant chirp β ¼ 8 · 10 4 rad∕m 2 ). The results for TDD slightly varied as a function of duration of input Gaussian pulse τ 0 (HWe −2 IM). The shortest pulse duration τ 0 , for which this CBG still reflected almost all the spectrum, was τ 0 ¼ 3.0 ps. The stretched pulse in both cases had a duration of ∼1 ns, i.e., 300 times longer than the input one. Figure 6 shows the comparison of numerical modeling of time-delay T cycle with our analytic approximation, Eqs. (27), (28), and (33). The square points in Fig. 6 show full duration of stretching-recompression process versus dimensionless coupling coefficient S ¼ jκjL for unapodized CBG.    as the dashed curve for unapodized CBG. We see quadratic deviation of those curves from 2L∕v 0 (form horizontal line) versus coupling jκj ¼ S∕L. All data are for the input pulses with τ 0 ¼ 12.0 ps at the central wavelength λ 0 ¼ 1.06 μm of our CBG. We see that the approach used in the derivation of the integrals [Eqs. (27), (28), and (33)] yields very reasonable correspondence with the results of numerical modeling.
Here is an important (and originally unexpected) observation. Too large a value of coupling (S ≥ 70) for that particular CBG and τ 0 ¼ 3.0 ps, while yielding 99.4% diffraction efficiency for energy of recompressed pulse, resulted in relatively poor quality of recompression: peak intensity of recompressed pulse constituted 52% of the incident one. Our interpretation of that decrease of recompression quality is the following. Different spectral components of incident pulse exhibit somewhat different values of the sum T cycle ðλÞ ¼ T 1 ðλÞ þ T 2 ðλÞ, from integrals [Eqs. (27) and (28)] with different values of z B ðλÞ, see Fig. 7 for apodized CBG with S ¼ 90. The presence of shorter values of T cycle ðλÞ points to the formation of precursor in the recompressed pulse. This effect leads to deteriorated recompression. Meanwhile, the fact that two considerably distant wavelengths have those shorter values of T cycle ðλÞ explains interference-type oscillations in the precursor.

Conclusions
We have developed a detailed Mathematica-based numerical tool for modeling the process of stretching-recompression by CBG with arbitrary profiles of the grating's phase and coupling coefficient. To better understand the results of numerous variants of that modeling, we developed approximate analytical model of time-delay dispersion TDDðλÞ. An unexpected result of that analytic model is that dependence of TDD on the coupling constant jκj ¼ S∕L starts with terms proportional to jκj 2 . We have shown that excessively large coupling in CBG leads to deterioration of recompression quality. We show that perturbations of grating phase with small (< ∼0.4 rad top-to-bottom) amplitude do not hinder the recompression quality much, enabling >90% of energy in the recompressed pulse.