## Abstract

An analytical expression for the time evolution of the diffraction pattern of an ultrashort laser pulse passing through a circular aperture is obtained in the Fresnel regime. The diffraction is not constant in time as the pulse travels through the aperture. This may have implications in experiments involving fast dynamics. Examples of the evolution of the diffraction pattern are given.

©2003 Optical Society of America

## 1. Introduction

Pulsed laser beams have many applications in areas such as communications, spectroscopy, chemistry, plasma fusion and microscopy. The pulse length has continuously been reduced to the level of only a few cycles, or even of less than one cycle in the case of terahertz beams. The usual theory for short pulses, which assumes that all spectral components of the pulse can be considered in a manner that is independent of the wavelength, is valid only for pulses of many cycles. For ultrashort pulses, however, this theory breaks down, and we must consider how the individual spectral components propagate. Thus, one cannot treat short pulses as continuous beams; in order to fully describe the behavior of ultrashort pulses in a system, a deeper understanding of the propagation, imaging properties and diffraction effects of these pulses is needed. In the literature, many aspects related to propagation and diffraction of ultrashort pulses have been studied [1–9]. However, although diffraction studies of the integrated intensity of the whole pulse have been done, there are as far as we know no detailed studies on the evolution in time of the diffraction pattern. This is mainly due to the fact that because the light field oscillates so fast there is no detector in the optical regime which could resolve the temporal evolution of the diffracted pattern. But this evolution may play an important rule in understanding the fast dynamics of the interaction of light and matter. Also, a direct detection of the evolution of the diffraction pattern may be achieved with terahertz pulses with one or less cycles, since the time domain is considerably larger than in the optical regime and thus the time dependent detection should be realizable.

In this paper, we analyze the evolution in time of the diffraction pattern of an ultrashort pulse through a circular aperture in the Fresnel regime. In the case of a very short pulse, the diffraction pattern is not a simple Airy disc, nor is it constant in time; it has a temporal evolution as the pulse hits the aperture. This effect is mainly due to the difference in the propagation time of the wavelets that originate from different points of the aperture to a point at the observation plane. In a recent paper, the time dependent behavior of the on-axis diffraction has been discussed [2]. Here, we derive an analytical expression for the whole diffraction pattern as a function of time, so that the build-up of the diffraction pattern at a screen behind the aperture can be seen starting at the arrival of the pulse at the aperture and following it until the pulse ends.

In the first part of this article we show the calculation based on the Fresnel diffraction formula [10,11]. In the second part we show the profile of the diffraction patterns as a function of time using as an example pulses in the terahertz regime. A film showing the evolution of the 2D pattern is also available.

## 2. Theoretical aspects

In order to calculate the diffraction pattern of an ultrashort pulse we follow the lines of Ref. [1]. Under paraxial approximation and for incident illumination of wavelength λ, the Fresnel diffraction formula of the spectrum of the diffracted field from a circular aperture of radius a in radial co-ordinates in the observation plane is given by:

where ρ_{1}=*r*_{1}*/a* and ρ=*r/a* are the normalised radial co-ordinates with *r*_{1}${\mathit{=}\mathit{(}x}_{\mathit{1}}^{\mathit{2}}$${\mathit{+}y}_{\mathit{1}}^{\mathit{2}}$*)*^{1/2}
and *r=(x*_{2}*+y*_{2}*)*^{1/2}
being the real co-ordinates in the aperture and the observation planes, respectively. Here, the Fresnel number is defined as *N*=π*a*^{2}*/(*λ*z)* and *z* is the distance between these two planes. *Ẽ*_{1}
is the spectrum of the pulse in the aperture plane given by [12–14]:

where ω_{0} is the frequency shift, which gives an estimate of the average wavelength of the spectrum, ω_{1} is the spectrum width, and s is a positive parameter (here we are interested in the small values of *s*, such as 0, 1 or 2, which corresponds to very short pulses in the time domain). Γ is the Gamma function, and particularly Γ(*s*+1)=*s*! if *s* is a positive integer. This model for the spectrum of the pulse satisfies the physical requirement that it includes only positive-frequency components and the pulse has a simple analytical form. Note also that as s increases, the spectrum tends to a Gaussian [12].

Since we consider the spectrum of a plane wave in the diffraction plane, according to Eq. (1) the field *Ẽ*
_{1}(ρ_{1}, ω) does not depend on ρ_{1}. As a consequence it is straightforward to calculate the integral of Eq. (1) analytically, using the Lommel functions [11]:

*U*_{1}
and *U*_{2}
in Eq. (3) are the *U* Lommel functions defined by:

By these considerations we can rewrite the spectrum by combining Eqs. (2) and (3) into Eq. (1):

$${e}^{-\left(s+1\right)\frac{\omega -{\omega}_{0}}{{\omega}_{1}}}{e}^{-i\frac{{a}^{2}\omega}{2cz}}\frac{2cz}{{a}^{2}\omega}\left({U}_{1}(2N,2N\rho )+i{U}_{2}(2N,2N\rho )\right)$$

$$=i{\alpha}^{s+1}\frac{{\left(\omega -{\omega}_{0}\right)}^{s}}{\Gamma \left(s+1\right)}{e}^{-\alpha \left(\omega -{\omega}_{0}\right)}{e}^{-i\omega \beta}\left[{U}_{1}(\frac{{a}^{2}\omega}{\mathit{cz}},\frac{{a}^{2}\mathit{\omega \rho}}{\mathit{cz}})+{\mathit{iU}}_{2}(\frac{{a}^{2}\omega}{\mathit{cz}},\frac{{a}^{2}\omega \rho}{\mathit{cz}})\right],$$

where $\alpha =\frac{s+1}{{\omega}_{1}}$ and $\beta =\frac{1}{c}\left(z+\frac{{a}^{2}}{2z}\left(1+{\rho}^{2}\right)\right).$

We then separate *Ẽ* in *Ẽ*
_{1}+*Ẽ*
_{2}, each field *Ẽ*_{j}
being associated with the Lommel function *U*_{j}
. The temporal evolution of the electric field *E*_{j}
is obtained by taking the Fourier transform of Eq. (5) (note that the integral starts at the point ω_{0} because of the function *f*_{s}
):

This integral can be evaluated analytically if ω_{0}=0. Thus, by the change of variables ω'←ω-ω_{0} and ω←ω' we obtain

$$\mathrm{where}\phantom{\rule{.2em}{0ex}}K=\frac{i}{\Gamma \left(s+1\right)}{\alpha}^{s+1}{e}^{-i{\omega}_{0}\beta}.$$

This integral was evaluated by following the steps:

- First, we write the Lommel function as a sum of Bessel functions. Then, by interchanging the sum sign and the integral sign (which is allowed when the Lommel function converges, i.e., for the shadow region ρ>1)[11], we obtain:

In the final formula we have to avoid the negative indices for the Bessel functions (because of the Gamma functions). In order to do so we rewrite Eq. (9) using *J*
_{-ν}(*x*)=(-1)^{ν}
*J*
_{ν}(*x*).

Thus, *J*
_{ν}(*a+b*) can be rewritten as

where the Bessel functions on the right hand side have only positive indices (if ν is positive).

Finally we have three double-sums and the evaluation of the integral can be done (for instance it is written here only according to the first expression of the addition theorem):

where ${N}_{0}=\frac{{a}^{2}{\omega}_{0}}{2cz}.$

The last integral is calculated as follows [15, p. 384]:

This integral converges provided *s* is positive, (which is always verified in our case), and provided |*b*|<|*a*| (either the two parameters are real or complex).

This leads to the final expression for the integral:

$$\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}{\left(\frac{\rho {a}^{2}}{2cz}\right)}^{m+2k}\frac{\Gamma \left(s+1+m+2k\right)}{\Gamma \left(k+m+1\right)}{\left(\alpha +i\left(\beta -t\right)\right)}^{-\left(s+1+m+2k\right)}$$

We rewrite now Eq. (13) using Eq. (10) in order to avoid the Gamma function Γ(*n*) with n<0. Then the second form of the addition theorem using only positive indexes for the Bessel and Gamma functions becomes:

$$K{e}^{i{\omega}_{0}t}\sum _{p=0}^{\infty}\sum _{k=0}^{\infty}\sum _{m=0}^{\infty}\frac{{\left(-1\right)}^{p+m+k}}{k!}{\rho}^{-j-2p}\times {J}_{m+j+2p}\left(\rho {N}_{0}\right)\times {\left(\frac{\rho {a}^{2}}{2cz}\right)}^{m+2k}\frac{\Gamma \left(s+1+m+2k\right)}{\Gamma \left(k+m+1\right)}{\left(\alpha +i\left(\beta -t\right)\right)}^{-\left(s+1+m+2k\right)}+$$

$${Ke}^{i{\omega}_{0}t}\sum _{p=0}^{\infty}\sum _{k=0}^{\infty}\sum _{m=1}^{j+2p}\frac{{\left(-1\right)}^{p+k}}{k!}{\rho}^{-j-2p}\times {J}_{j+2p-m}\left(\rho {N}_{0}\right)\times {\left(\frac{\rho {a}^{2}}{2cz}\right)}^{m+2k}\frac{\Gamma \left(s+1+m+2k\right)}{\Gamma \left(k+m+1\right)}{\left(\alpha +i\left(\beta -t\right)\right)}^{-\left(s+1+m+2k\right)}+$$

$${Ke}^{i{\omega}_{0}t}\sum _{p=0}^{\infty}\sum _{k=0}^{\infty}\sum _{m=j+2p+1}^{\infty}\frac{{\left(-1\right)}^{p+m+k-j}}{k!}{\rho}^{-j-2p}\times {J}_{m-j-2p}\left(\rho {N}_{0}\right)\times {\left(\frac{\rho {a}^{2}}{2cz}\right)}^{m+2k}\frac{\Gamma \left(s+1+m+2k\right)}{\Gamma \left(k+m+1\right)}{\left(\alpha +i\left(\beta -t\right)\right)}^{-\left(s+1+m+2k\right)}$$

Finally we have *E*(ρ, *z, t*)=*E*
_{1}(ρ, *z, t*)+*iE*
_{2}(ρ, *z, t*). The instantaneous intensity (considering that the detector has a delta function response in time) is given by *I* (ρ, *z,t*)=|*E*(ρ, *z,t*)|^{2} [1,2].

We make a couple of remarks about the computation of Eq. (14). First, it does not converge for ρ<1 because of the *U* Lommel functions in Eq. (7), which diverge as ρ tends to 0. We verified numerically however that this divergence is very slow. If we use Eq. (14) to compute the diffraction patterns, this problem does not seem to be disturbing since the diffraction pattern extends to an area which is much bigger than ρ=1. In order to avoid the centre (ρ=0) we simply substitute ρ=0 by a small number. The problem in using the *U* Lommel function will appear if one tries to calculate the temporal evolution of the intensity of the diffracted light exactly on the axis. Then an algorithm involving for example the *V* Lommel function will be necessary [11]. Here we concentrate mostly in looking at the diffraction pattern on a large area, so that the calculation of the field exactly at the axis was avoided. Actually, the calculation on the axis only can be done directly and has been already analysed for example in Ref. [2].

Another problem already evoked is the integral of Eq. (12) which may diverge in function of the modulus of its parameters: formally the integral converges only if the condition $\frac{\rho {a}^{2}}{cz}\le \sqrt{{\alpha}^{2}+{\left(\beta -t\right)}^{2}}$ is satisfied. In the computed diffraction patterns shown here the condition of convergence is always satisfied.

Finally, this formula takes the time propagation between the two planes into account; as the pulsed beam passes the aperture plane at *t*=0, it arrives around *t=z/c* in the observation plane. We can observe the diffraction pattern around *t*=0 by simply substituting *t* by *t+z/c* in the formula, i.e., by just rewriting
$\beta =\frac{{a}^{2}}{2cz}\left(1+{\rho}^{2}\right).$
The formulas computed here use this expression.

We used the formula given by Eq. (14) in a Matlab program to compute the diffraction pattern. We have made movies which show the temporal evolution of the 2D diffraction pattern using absolute values in order to simulate the actual diffraction pattern seen at the screen. Also, cross-sections of one of the axis using normalised values are shown in the figures below. Although the formulas given by Eqs. (13) and (14) are very general, we have limited ourselves to calculate the diffracted field in the low Fresnel number regime. For practical applications this is more interesting and numerical convergence problems are avoided.

## 3. Simulations

In this section we present simulations of the temporal evolution of the diffracted pattern from a circular aperture in the Fresnel regime. We have chosen to calculate the patterns in the terahertz domain. The investigation of diffraction properties in the terahertz regime is interesting since pulses with one or two cycles are routinely being generated in laboratories and there are interesting applications for imaging. Also, the time scale of the terahertz pulses is not so short as for visible light so that detection of the evolution is feasible. Here we show a few examples of some normalised cross-sections for different physical parameters calculated at consecutive times. Using the same parameters as for each series of cross-sections we include a film of the 2D pattern and the cross-sections calculated without normalisation.

The graphics presented below show a number of examples on how the diffraction pattern of an ultrashort pulsed laser beam evolves as a function of time. In the first series, the radius of the diffracting pupil *a* is 0.5 mm and the distance between the two planes *z* is 1 m. The spectrum width is 6.10^{12} rad/s (about 1 THz), which corresponds to a pulse whose duration is about 1 ps, and the frequency shift is 10^{14} rad/s (about 16 THz or a wavelength of 20 µm). With these parameters we have a Fresnel number which is equal to 0.04. The parameter s is set to 1. The non-diffracted pulse is at its maximum when *t*=0.

In Fig. 1(a) the diffraction pattern of the pulse at *t*=- 0.5 ps is shown. This pattern is approximately the same up to *t* almost equal to 0 (time in the observation plane minus the propagation time between aperture and observation planes). When *t*=0, the peak widens at the top and narrows a little bit at the bottom, as shown in Fig. 1(b). At t=0.1 ps, two peaks (in fact a ring in the whole surface of the plane) start to appear from both sides of the centre whose relative value is decreasing (Fig. 1(c)). When the pulse is vanishing, the value in the centre decreases faster than the two lateral peaks as shown in Fig. 1(d).

In the second series, the radius of the aperture a, the frequency shift and the spectrum width remains the same as in Fig. 1 while the distance z is decreased to 250 mm. With these parameters we have a Fresnel number of 0.17. Figure 2(a) shows the pattern for time t=0. Before this time, the form is about the same but with lower intensity. As t=0.2 ps (Fig. 2(b)) the peak starts to divide as in the previous case, and it follows it up to about t=0.5 to 0.6 ps. For later times, the peaks to centre ratio diminishes quickly, and two peaks appear on both sides, which correspond to the first secondary peaks of the Airy disc.

At *t*=0.8 ps, the first two peaks vanishes and there is a big flat centre. The secondary peaks also grow compared to the centre. At *t*=1 ps, the peaks which appeared at *t*=0.2 ps have totally disappeared; it remains only the secondary peaks whose relative value is still increasing. Finally, at *t*=1.5 ps the centre is vanishing, but not the side peaks. The latter will still spread away from the centre while their value decreases to zero.

From these two examples one can see that the obtained results agree with the argument given at the beginning of the paper that light will take longer time to arrive at the outer region of the observation plane than at the region around the centre. If one watches the evolution in time of the pattern one can see that at earlier times there is light only at the centre of the screen and only at later times the pulse will spread out far from the centre. At the latter it happens also that the pulse has already passed the central region, resulting in a ring configuration. In the second example the ring is less evident than in the first example (as can be seen by watching the films of the 2D patterns) since the observation plane is nearer than in the first example and light stays more concentrated in the central region.

## 4. Conclusions

By analyzing the diffraction of an ultrashort pulse through a circular aperture we have shown that the diffracted pattern is not constant in time but evolves as the pulse hits the aperture. We derived a general analytical expression of the time evolution of the diffracted field which under certain conditions (related to the stability of the integrals and stability of the computer algorithm) allows us to calculate the diffracted pattern for a given ultrashort pulse and a given Fresnel number. In two examples in the terahertz domain, we have shown that this evolution depends on the different input parameters, and particularly of the Fresnel number *N*_{0}
.

In order to make the solution more complete some general problems in the algorithm of the program, as suggested in the first part should be solved. First, there are problems at the center of the diffracted pattern (divergence of the *U* Lommel function), and second, is the problem of the integral given by Eq. (12) which does not converge for a certain choice of parameters. We also would like to note that although the analytical result given by Eq. (14) seems to be very complicated, the convergence of the sums is very rapid (only a couple of terms are necessary) which makes the evaluation quite efficient. Alternatively one may also try to compute Eq. (6) directly numerically, but the main reason that we have chosen for the analytical solution is that this gives more insight on the time dependence and the shape of the diffracted pattern. With the analytical solution it is also possible to explore certain limiting parameters which may be intersting.

Finally, we emphasize the importance of looking at the evolution of the diffraction pattern of ultrashort pulses. As seen in Figs. 1 and 2, the diffracted pattern takes different shapes as a function of time. This phenomenon may have consequences in experiments involving dynamics of systems due to interaction of light with ultrashort pulses. This evolution could be observed in actual experiments involving a few cycles terahertz pulses.

## Acknowledgments

We acknowledge the financial support from the Dutch Foundation for Fundamental Research on Matter (FOM). The authors appreciate the helpful discussions with Prof. Colin Sheppard and Prof. Joseph Braat, and the assistance of Dr. Erik Koelink of the Department of Applied Mathematics of the Delft University of Technology in evaluating the integrals analytically.

## References and links

**1. **M. Gu and X. S. Gan, “Fresnel diffraction by a circular and serrated apertures illuminated with an ultrashort pulsed-laser beam,” J. Opt. Soc. Am. A **13**, 771–778 (1996). [CrossRef]

**2. **Z. Ziang, R. Jacquemin, and W. Eberhardt, “Time dependence of Fresnel diffraction of ultrashort laser pulses by a circular aperture,” Appl. Opt. **36**, 4358–4361 (1997). [CrossRef]

**3. **H. Quan-Sheng and Zhu-Zhen-he, “Diffraction of ultrashort pulses,” Acta Phys. Sin. **37**, 1432–1437 (1988).

**4. **M. Kempe and W. Rudolph, “Analysis of confocal microscopy under ultrashort light pulse illumination,” J. Opt Soc. Am. A **10**, 240–245 (1993). [CrossRef]

**5. **J. Anderson and C. Roychoudhuri, “Diffraction of an extremely short optical pulse,” J. Opt. Soc. Am. A **15**, 456–463 (1998). [CrossRef]

**6. **M. Kempe and U. Stamm. “Spatial and temporal transformation of femtosecond laser pulses by lenses and lens systems,” J. Opt. Soc. Am. B 9, 1158–1165 (1992). [CrossRef]

**7. **H. Ichikawa, “Analysis of femtosecond-order optical pulses diffracted by periodic structure”, J. Opt. Soc. Am. A **16**, 299–304 (1999). [CrossRef]

**8. **J. Cooper and E. Marx, “Free propagation of ultrashort pulses” J. Opt. Soc. Am. A **2**, 1711–1720 (1985). [CrossRef]

**9. **R. W. Ziolkowski and J. B. Judkins, “Propagation characteristics of ultrawide-bandwidth pulsed Gaussian beams,” J. Opt. Soc. Am. A **9**, 2021–2030 (1992). [CrossRef]

**10. **W. Goodman, *Introduction to Fourier optics* (Mc Graw-Hill, NY, 1968).

**11. **M. Born and E. Wolf, *Principles of Optics* (Pergamon, NY, 1980)(p.437–439, chap. 8.8.1).

**12. **C. F. R. Caron and R. M. Potvliege, “Free-space propagation of ultrashort pulses: space-time couplings in Gaussian pulse beams,” J. Mod. Opt. **46**, 1881–1891 (1999). [CrossRef]

**13. **S. Feng, H. G. Winful, and R. Hellwarth, “Gouy shift and temporal reshaping of focussed single-cycled electromagnetic pulses,” Opt. Lett. **23**, 385–387 (1998). [CrossRef]

**14. **C. J. R. Sheppard, “Bessel pulse beams and focus wave modes”, J. Opt. Soc. Am. A **18**, 2594–2600 (2001). [CrossRef]

**15. **G. N. Watson, *A treatise on the theory of Bessel functions* (Cambridge University Press, Cambridge, 1966).