1. Introduction

The Gaussian beam is a basic concept of wave optics and it has many, diverse applications in the areas such as laser resonators, waveguides, nonlinear optics, optical information processing, etc. From the theoretical point of view the waist region of a (weakly-focused) Gaussian beam constitutes perhaps the simplest realization of the plane wave concept. A lateral confinement required for a beam to be realizable (means not to possess an infinite energy) causes a minimal uncertainty-principle-determined angular spread in the case of the Gaussian beam. Also, the phase surfaces are planes in its waist region.

However, the Gaussian beam is not an exact solution to the free-space wave equation, instead – as a result of paraxial diffraction theory – it obeys the parabolic truncated version of the equation. Therefore the paraxial approximation is insufficient in a high-aperture case – it starts breaking down for numerical apertures of about 2^-½ (i.e. f/0.5), which is far from being a limit for various devices nowadays. That is why there has been a considerable and growing interest in the extending of the model of the Gaussian beam beyond the paraxial regime and a variety of infinite-series and more or less approximate solutions have been proposed (see [1–3] and references therein). Unfortunately, as the beam convergence/divergence angle grows, the mathematical expressions for the models become either computationally onerous or increasingly inaccurate. Moreover, the elegant mathematical way to get the common Gaussian beam as the field of a point source placed at an imaginary displacement along the optical axis, which was introduced by Deschamps [4], tends to give nonphysical singularities if applied beyond the paraxial case.

The advances in ultrafast optical technology have made possible the generation of electromagnetic pulses of femtosecond duration, which comprise few cycles (in the optical range) or even less than one cycle (in the terahertz range). Such ultrashort pulses exhibit – among other interesting features – directly the Gouy phase shift through a temporal reshaping and polarity reversal in the course of propagation, as demonstrated by thorough study of various models in the paraxial limit (see [5–10] and references therein). Detailed treatments of the single-cycle pulses in Refs. [5, 8] start from an exact solution of the wave equation, still a majority of useful results and relations have been obtained in the paraxial approximation.

Ultrawideband pulses contain the Fourier components with frequencies far below the mean one, the lowest one being the dc-component. For the low-frequency components the product of the wavenumber and the Rayleigh range (a half of the confocal parameter) or, equivalently, the ratio of the beam waist diameter to the wavelength, may easily turn out to be of the order of unity. This means that for wide-band and pulsed Gaussian beams the paraxial model is unsatisfactory even in a modest-aperture case, especially if the low-frequency tail of the spectrum has not been suppressed or cut off.

A surprisingly simple closed-form expression for Gaussian-like monochromatic focused beams, which is an exact solution to the scalar Helmholtz equation and in the paraxial limit reduces to the common Gaussian beam mode, has been found to be possible in oblate spheroidal coordinates [11]. This can be understood as a result of matching between the beam shape and one of the oblate spheroidal coordinate (henceforth the acronym OSC will be used) surfaces, which both are hyperboloids of one sheet. However, the solution suffers from the above-mentioned singularities and discontinuities. These nonphysical peculiarities originate from the description of the field as due to a source, which – despite its imaginary-coordinate location – is inherently contradictory, as physically the field is source-free, i.e. it obeys the homogeneous wave equation (HWE) in the spatial region under consideration.

In order to overcome these difficulties, recently a modification of the Deschamps complex point-source theory has been developed [2], according to which the sources are accompanied by sinks “absorbing” the external “incoming” field. This can be recognized as a version of Huygens’s principle. This approach – the addition of an incoming wave associated with a sink, which results in a standing-wave nature of the total field (in the monochromatic case), has also been fruitfully applied to the description of the nonparaxial beams in OSC [3]. A rigorous time-domain deduction of the sink-source method from the basics of electrodynamics can be found in [12, 13], particularly for the description of some recently-discovered localized pulse-shaped exact solutions of the free-space wave equation (see [13–16] and references therein).

In the present paper, we give a simple but an exact treatment of the behavior (illustrated by animated 3D plots) of ultrashort pulses under arbitrarily tight focusing. The ingredients in building the appropriate model in the next section are: (i) the source-and-sink representation possible for any free field, i.e. for any solution of the HWE, (ii) the Deschamps complex source-point technique, and (iii) OSC – the oblate spheroidal coordinates. No matter how wide the pulse bandwidth or the solid angle subtended by the focusing optics are, the model remains exact in the sense of obeying the scalar HWE, yet it reduces to the common Gaussian beam in the paraxial and narrow-band limits. Naturally, the vectorial nature of electromagnetic waves becomes essential at high apertures. However, for simplicity in this paper we restrict ourselves to the scalar treatment in order to better convey the overall physical picture. Possible extensions of the model are discussed in the 3^rd section.

2. The model

For beginning, let us consider an impulsive point source of unit strength, which is placed at the origin of coordinates. Thus, let the source density distribution be expressible via the Dirac delta functions ρ(t,r)=δ(t) δ(r). As it is well known from electrodynamics (or – mutatis mutandis – from acoustics), such source emits a δ-pulse-shaped wave, spherically expanding from the origin at the positive times t. This elementary wave – the retarded Green function of the wave equation, which in the Gauss units reads as G₊=(c/4πR) δ(R-ct), R being the distance between the field and the source points – allows one to express any emitted field as the 4D-convolution with the corresponding source distribution ρ(t,r). Due to the temporal symmetry of the HWE, one can in the same way deal with an elementary wave, spherically collapsing onto the origin – with the advanced Green function G_-=(c/4πR) δ(R+ct), in which case the point looks like a sink rather than a source. The difference between these Green functions, which in relativistic field theories is usually denoted by D:

is a specific elementary solution of the source-free wave equation (i.e., the HWE) and represents a spherical delta-pulse-shaped wave, first (at negative times t) converging to the origin (the right term) and then (at positive times t) diverging from it.

 

Fig. 1. The elementary free field D(t,R). The first 0.7MB video clip shows its temporal evolution. The field is nonzero only on the transparent spherical surface, which has been made blue to mark the negative values and red – the positive values of the wave function.

Download Full Size | PDF

What is important for our purpose, is that the free-field Green function D(t,R) allows one to express any solution of the HWE as the 4D convolution with an appropriately chosen “charge” distribution ρ(t,r), whereas the quantity ρ(t,r) expresses strength of both the source and the sink at the same point. The minus sign between the two terms in Eq. (1), which follows from the requirement that a source-free field cannot have a singular point at R=+0, is crucial and it guarantees the vanishing of the function at t=0. It is just the resulting change of the sign of D(t,R), which takes place when the elementary wave goes through the collapsed stage at the focus, that is responsible for the 90-degrees phase factor known from the Huygens-Fresnel-Kirchhoff principle and for the Gouy phase shift peculiar to all focused waves. Generally a given field does not determine uniquely the distribution of the charge generating the field. Consequently, the sinks and sources need not to cover a distant surrounding surface – the picture commonly associated with Huygens principle, but may be spatio-temporally localized in an appropriate way provided they generate the same given field.

The above-given was an abridgement of a rigorous deduction of the sink-and-source method (see also Refs. [12, 13]). Now, let a point-like “charge” at the origin function in time according to the Lorentzian curve with the width parameter Δ (HWHM), i.e. let the distribution be

where q is the total sink-and-source strength, the parameter t ₀ will be specified later on, and where the Lorentzian expressing the temporal dependence has been introduced in the form of the complex analytic signal for mathematical convenience. When expressing the source-free field through the 4D-convolution of the functions D(t,R) and ρ(t,r) like it is commonly done for obtaining source-induced retarded potentials, the spatial integration is trivial and one obtains:

Eq. (3) describes two different smoothly impulsive waves, still spherically symmetric and collapsing/expanding through the origin like in Fig. 1: the real part represents a unipolar (at large times |t|>Δ) pulse which is odd with respect to time and changes its sign at t=0, while the imaginary part represents a bipolar pulse which is temporally even and therefore unipolar at the origin (such tranformation of pulses will be illustrated by animated plots later on).

Now we make use of the Deschamps complex source-point technique, viz., we shift the “charge” from the origin to an imaginary location z ₀=id on the axis z. Naturally, this makes the expression of the distance R between a field point and the “charge” location rather complicated. Fortunately, we can avoid the complications by introducing the system of oblate spheroidal coordinates (OSC) instead of the Cartesian x,y,z by the following relations [3,11]:

The system of OSC might be imagined as system of spherical coordinates where the origin point has been expanded into a disk of radius d placed in the equatorial plane (see Fig. 2). The coordinate surfaces φ=const are – just as in the case of spherical coordinates – the meridian planes containing the z axis. The surfaces η≡cosθ=const≠±1 are hyperboloids of revolution of one sheet – just matching the shape of the Gaussian beam waist -whose asymptotes pass through the origin inclined at an angle θ with the z axis, while two degenerate surfaces η=±1 constitute the z axis. The surfaces ξ=const≠0 are oblate ellipsoids having an interfocal spacing of 2d, while the surface ξ=0 is the circular disk.

 

Fig. 2. Intersection of a surface ξ=const (red) and a surface θ=const (blue) with the coordinate plane φ=±π/2 in the case of the oblate spheroidal coordinate system with the parameter d=0.5, see also Fig.3. (0.4 MB video clip shows the tranformation from the spherical coordinate system).

Download Full Size | PDF

In OSC the square of the distance that we need for Eq. (3) is simply R ²=d ²(ξ-iη) ², while the length d turns out to be the Rayleigh range of the field in the paraxial limit [3]. For convenience, we set t ₀=d/c and Δ=t ₀+a/c. The latter definition introduces an arbitrary contribution a>0 to the pulselength but assures the fulfilling of the condition Δ>t ₀ required for finiteness of the field at r=t=0. Thus, from Eq. (3) the model scalar field reads

The real and the imaginary part of the relatively simple expression of Eq. (5) give us two different exact solutions of the HWE, which are studied in the next section.

3. Results and discussion

To visualize the pulsed solutions in nonparaxial conditions, we set the pulselength parameter a=0.05, where the confocal parameter 2d has been taken for the unit of length. The same relative unit is also used on the length scales of the forthcoming figures. In an absolute scale this unit might possess values from a few centimeters or less – to relate our results to typical terahertz pulses – down to, say, ten microns for a frontier femtosecond optical experiment. Consequently, the ratio a/d=0.05/0.5=0.1 is rather large stressing the nonparaxiality. Eq. (5) does not depend on the angle φ, i.e. the field is axisymmetric around the z axis and it suffices to study the behavior of the pulse in the plane y=0, which is depicted in Fig. 3 with coordinate lines of ξ,θ shown in it. Thus, by making use of the third axis of a 3-D plot for depicting the strength of the field and relating the time t to the frame number of animation, we can visualize all spatio-temporal dependences of Eq. (5) by two video clips (Fig. 4 and 5). The field values have been computed on a 60×60 grid over the range (ξ=0…4.8)×(θ=0…π). Thus, both plots consist of 7200 data points with their density increasing towards the center of the ellipse-bounded region in the z-x plane.

 

Fig. 3. The z-x plane (or two meridian planes φ=0 and φ=π) and lines of intersection with the coordinate surfaces ξ=0.5, 1, 1.5, 2,…. and θ=0, 0.1π, 0.2π,… Note that the pair of blue θ-hyperbolas, nearest to the axis z, together with the arcs of red ξ-ellipses between them, are reminiscent of a scheme of the (paraxial) Gaussian beam.

Download Full Size | PDF

 

Fig. 4. The frame t=0 of the animated plot (2.2 MB) of the pulse given by the real part of Eq. (5). The propagation axis z has been directed from left to right. For a better visualization a color “lighting” of the surface plot has been used and a colored contour plot of the same data is shown at the bottom. The Rayleigh range d=0.5, the pulsewidth a=0.05 (values given relative to the confocal parameter, the same unit of length has been used in the scales on the basal plane). The vertical scale of the amplitude has been normalized to the unit “charge” q=1. Fast start and end frame jumps have been caused by cutting the clip shorter in order to meet the file length limit. The “shakes” of the scales in the animation indicate the time moments of the penetration of the basal plane by the pulse peak.

Download Full Size | PDF

 

Fig. 5. The frame t=0 of the animated plot (2.5 MB) of the pulse given by the imaginary part of Eq. (5). Other parameters and units are the same as in Fig. 4. The pulse peak has been partially cut off, as its highest value (~19) falls outside the vertical scale.

Download Full Size | PDF

 

Fig. 6. A frame of the animated 2-D plot (1 MB) showing the on-axis behavior of the real part (blue line), the imaginary part (magenta), and the modulus (dotted) of Eq. (5). Other parameters and units are the same as in Fig. 4. The frame depicts the distribution of the field along the z axis at the converging stage (at instant 2d/c prior to the focus).

Download Full Size | PDF

We can see that both pulses undergo a considerable temporal reshaping in the course of propagation – a spectacular manifestation of the Gouy phase shift, which is an almost imperceptible finesse in the case of monochromatic or narrow-band focused beams. In passing through the focus the real pulse is both time reversed and polarity inverted while the imaginary pulse is only time reversed. This is explained by the following space-time symmetries of the field: one pulse of the Hilbert-conjugated pair (the real one in our case) is antisymmetric under the reflection of both the z axis and the time, while the other pulse is symmetric – the rules that had been earlier established for the on-axis behavior [6, 8]. Our results as well as an inspection of Eq. (5) show that these symmetry rules remain also valid beyond the paraxial approximation.

For a better comprehension of Eq. (5) and a comparison with the results of other authors, we introduce a local dimensionless time τ=ξ-ct/d=ξ-t/t ₀ for the pulse as well as two dimensionless width parameters δ-=(1-η)+a/d and δ₊=(1+η)+a/d, and extract distinctly the amplitude and phase in Eq. (5):

Here it is more convenient to have in mind an alternative choice of the regions of the variation of the spacial variables: -∞≤ξ≤∞ and 0≤η≤1. In the paraxial limit (a/d≪1, |η|≈1), the quantity 2ξ-τ can be replaced by 2ξ and we obtain

The first multiplier in the denominator shows that the pulse intensity (Lorentzian) profile remains invariant along any given hyperbola η=cosθ=const as far as the second one expresses nothing but the change of the pulses’s overall intensity depending on the distance from the focus. The constant 900 phase factor is unimportant and it results from our particular choice of the phase in Eq. (2). If the last term in the exponent were absent, the phase – and hence the profiles of both the real and the imaginary pulse – would also remain invariant along the hyperbolas. The pulse reshaping is introduced by the last term tan^-1(ξ), which is the Gouy phase shift that any finite wave field obeys upon passing through a focus.

Our results can straightforwardly be extended for the few- and many-cycle pulse cases by taking a temporal derivative of a corresponding order from Eq. (5–6A). In this regard it is appropriate to point out that one reaches the same final expression Eq (5) when working in the Fourier (wavenumber k or frequency ω=kc) domain. Namely, the calculation of a k-domain superposition of the monochromatic nonparaxial Gaussian beams expressed in OSC [3] reduces to taking an integral which is nothing but the Laplace transform of S(k)sin(Rk), where the complex number R=d(ξ-iη) and S(k) is the spectrum without the exponential factor exp(-ak). There is a general formula in the transform tables for choice S(k)=k^p, which means that the time-derivative pulses of any order p can be readily obtained. The setting p=0 reproduces Eq. (5), yet many other interesting exact solutions to the HWE in OSC can be found via the Laplace transform tables.

Eq. (6A) – rather compact thanks to OSC – can be related to Eqs. (2.10–2.15) of Ref. [10], taking there p=0 in the Fourier spectrum ω^p exp(-ωτ₀) of the pulse. By choosing the simplest temporal dependence in Eq. (2) for our model we have intentionally dealt with the spectrum of the simplest form exp(-ak), that is setting S(k)=const or p=0, in which case the low-frequency tail of the spectrum is not suppressed.

For the generalization of the results to electromagnetic vector fields one can use different ways [5, 8, 17]. In order to get linear field polarization which corresponds to the usual TEM₀₀ mode in the paraxial limit, we should replace the “charge” with an electric dipole and a magnetic dipole, oriented along the x and the y axis, respectively, and combined with the sinks possessing the same dipole properties [17]. This procedure, aimed at nonparaxial expressions for the pulses of the fields E and H (or B), is perhaps the most convenient to be carried out by making use of the Hertz vector technique (like it is done in [8]), taking the scalar solution Eq. (5) as the magnitude of the Hertz vector. As a matter of fact, when rewriting Eq. (5) in cylindrical coordinates, it turns out that it coincides with Eq. (2.5) of Ref. [8] – an exact two-parameter (q ₁ and q ₂) solution to the scalar HWE – if one takes the the pulsewidth parameter q ₁ equal to a of our model and sets q ₂=2d+a, which in the paraxial limit a≪d means that the confocal parameters of both model scalar fields concide. Thus, by Eq. (5) we have expressed the so-called modified power spectrum pulse – a particular version of “electromagnetic directed energy pulse trains” (see [14] and references therein) in the oblate spheroidal coordinate system. Also, the procedure of deducing the fields E and H from the scalar generating field has been accomplished in [8] already and a number of interesting results concerning single- and few-cycle pulses have been obtained, yet in the paraxial limit.

However, since E and H are second-order spatio-temporal derivatives of the Hertz vector, the spectra of their components acquire the k ²-dependence, thus suppressing the nonparaxiality of the initial scalar field. Fortunately, for the role of the generating scalar field we could take the indefinite integral of Eq. (5), in other words – choose the negative power p in the spectrum k ^p of the generating field, which means an arctan-type temporal dependence instead of the Lorentzian one introduced into our model by Eq. (2). Moreover, the above-mentioned general formula in the Laplace transform tables proves the existence of the transform of k ^p sin(Rk) even for the case p→-2. Hence, it should be possible, by appropriately changing the initial time dependence in Eq. (2), to reach finally the fields E and H possessing spectra, where the low-frequency tail is not suppressed and which therefore represent essentially nonparaxial vector fields.

If one wants to describe the pulses whose temporal profile cannot be deduced from the Loretzian via differentiation/integration and/or the Hilbert transform, one has to insert an appropriate temporal dependence into Eq. (2) or, equivalently, to use an expression for the spectrum S(k), other than the simple powers of the wavenumber or the frequency. Yet in a number of such cases closed-form final expressions can be found with the help of the Laplace transform tables.

One might question the applicability of the term “the Gaussian beam” for the pulses whose transverse profile is clearly non-Gaussian. A justification is that the energy of the pulsed beam diffracts in the same manner as a monochromatic Gaussian beam [8]. One could also make an issue of the backward-propagating components obviously introduced into the field by the sink-and-source pair. Fortunately, as it is known for focus wave modes [14, 16] – the Gaussian constituents of the modified power spectrum pulse – the forward-propagating plane wave components can be made dominant and the backward-propagating ones neglected, which is just the case for the pulsed fields E and H, as shown in [8]. Naturally, in the case of in-cavity propagation, there is no such issue at all.

Finally, the wave given by Eq. (5), as it does not depend on the asimuthal angle φ, is merely the lowest-order solution. Our model can be extended to non-axisymmetric higher-order pulsed solutions like it is done in the case of monochromatic Hermite-Gaussian or Laguerre-Gaussian modes [3].

4. Conclusion

By making use of the oblate spheroidal coordinate system to represent free field as if generated by a combined point-like source and a sink placed at a complex-number coordinate, a simple but exact treatment of the spatiotemporal behavior of sub-cycle pulses under arbitrarily tight focusing has been developed and illustrated by animated 3D plots. The model and the results obtained are not restricted by the common paraxial approximation and are applicable to ultrawideband isodiffracting ultrashort pulses in high-aperture resonators, such as a microcavity, and they are usable in femto- and attosecond optics in general.