{'fraction_non_alphanumeric': 0.05009137021160316, 'fraction_numerical': 0.022526599832630268, 'mean_word_length': 4.112100576217915, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-Paraxial Wave Analysis of 3D Airy Beams

The 3D Airy beam (AiB) is thoroughly explored from a wave-theory point of view. We utilize the exact spectral integral for the AiB to derive local ray-based solutions that do not suffer from the limitations of the conventional parabolic equation (PE) solution, and are valid far beyond the paraxial zone and for longer ranges. The ray topology near the main lobe of the AiB delineates a hyperbolic umilic diffraction catastrophe, consisting of a cusped double-layered caustic, but this caustic is deformed in the far range where the field loses its beam shape. The field in the vicinity of this caustic is described uniformly by a hyperbolic umilic canonical integral which is structured explicitly on the local geometry of the caustic as obtained from the initial field distribution. In order to accommodate the finite-energy AiB we also modify the canonical integral by adding a complex loss parameter. The canonical integral is calculated using a series expansion and the results are used to identify the validity zone of the conventional PE solution. The analysis is performed within the framework of the non-dispersive AiB where the aperture field is scaled with frequency such that the ray skeleton is frequency-independent. This scaling enables an extension of the theory to the ultra wide band (UWB) regime and ensures that the pulsed field propagates along the curved beam trajectory without dispersion, as will be demonstrated in a subsequent publication.


Introduction
Airy beams (AiB) are exact solutions of the time-harmonic parabolic wave equation (PE) in free space, that are generated by an Airy function initial condition in an extended aperture. The AiB's have attracted a considerable attention mainly because of their intriguing features, the most distinctive one is the propagation along curved trajectories [1][2][3][4]. Other intriguing properties are the weak diffraction [2] and self-healing (i.e., regeneration of the solutions if the beam is partially obstructed) [5]. Many interesting applications that utilize these properties have been suggested, e.g. [6][7][8][9][10].
In a previous paper [11], we provided a cogent wave-based description to the 2D-AiB's that clearly explains all these features. Using uniform asymptotic evaluation of Kirchhoff's radiation integral [12] and uniform geometrical optics (UGO) [13,14], it has been shown that the AiB is not generated by the main lobe of the Airy function in the aperture, but it is in fact a caustic of rays that emerge sideways from the side lobes of the extended aperture distribution. This implies that the AiB cannot be regarded as a "beam field" in the sense that the evolution of the main lobe along the curved propagation trajectory is not described by a local wave dynamics. It has also been established that ray theory in conjunction with the UGO provides an accurate solution for the AiB, which fully agrees with the conventional PE solution of [1,2] when the latter is valid. The UGO solution is accurate even at large observation ranges where the PE solution fails due to the curving beam axis. It has been noted that these ray-based tools can be used to construct other curving beam solutions in uniform and non-uniform media by back-projecting the rays from the caustic to the aperture and synthesizing the aperture distribution that focuses onto that caustic [11]. This concept has been applied recently in [15].
Using the ray description of [11] we have also constructed in [16] a class of non-dispersive Airy pulsed beams (AiPB). These solutions are obtained by synthesizing the aperture field distribution such that the ray skeleton is the same for all frequencies. This scaling of the aperture distribution, unlike previous suggestions [17], ensures that all frequency components focus on the same curved caustic, thus giving rise to a pulse that propagates along the curved trajectory with no dispersion. A closed-form time-domain expression for the AiPB has also been derived in [16] using the spectral theory of transients (STT) [18][19][20].
The main goal of this paper is to clarify the wave mechanism of the 3D-AiB by deriving uniform ray-based solutions that reveal the local features of the field. Furthermore, these solutions rely on the exact ray and caustic topology of the AiB and therefore do not suffer from the limitations of the conventional PE solution. We note that due to the curving of the beam trajectory, the PE solution is valid only for short propagation distances and in particular it does not predict the termination of the AiB at a certain range, as will be discussed here. In order to meet this goal we first analyze the ray topology of the 3D-AiB (Sec. 3) and then use this analysis to derive uniform solutions that are structured explicitly upon the local geometry of the caustic (Secs. [4][5]. As will be explained below, we emphasize here local uniform solutions rather than globally uniform solutions since the former yields explicit solutions that explain all the parameters and features of the field. The analysis is performed within the framework of the 3D non-dispersive AiB, defined in Sec. 2 such that the ray skeleton is frequency-independent. The non-dispersive formulation not only simplifies the ray analysis by making it more transparent, but it also enables the construction of the non-dispersive AiPB as noted above. This subject will be explored in a subsequent publication [21], where the results of the present work will be used to synthesize 3D-AiPB and to derive closed-form time-domain solutions. In Sec. 3 we study the exact ray and caustic topology of the 3D-AiB, concentrating on the region near the curved beam axis and deferring the discussion on the global caustic topology to Sec. 6. The analytical details are summarized separately in Appendix A. The topology near the beam axis has many intricate features and it is considerably more complicated than the approximate ray topology of the 3D PE solution in [7] or the 2D case in [11]. Specifically, we show that the main lobe of the AiB is described by a cusped double-layered caustic which is of the hyperbolic umbilic type [22][23][24][25][26] (see Figs. [4][5]. We also identify the range where the isolated hyperbolic umbilic caustic that describes the AiB field near its propagation axis, evolves into the more complicated parabolic umbilic caustic considered in Sec. 6, so that the AiB loses its form. It should be noted that these caustic evolutions (unfoldings) are expected from catastrophe theory principles (see e.g., Fig. 4.4 and Fig. 7.3. (3)(4)(5) in [24], respectively), yet, in this paper we provide explicit expressions for this overall topology and also use it to construct the field solution, as described below.
The field calculations are considered next in Secs. 4-5. We start from the exact spectral diffraction integral for the AiB, and then discuss its asymptotic evaluation in terms of geometrical optics (GO). The field near the beam axis, where the ray topology is of the hyperbolic umbilic type, is considered separately in Sec. 5 within the framework of diffraction catastrophe theory [24,27,28], discussed below.
Catastrophe theory [24,27,28] classifies the caustic topologies in stable polynomial forms, referred to as elementary catastrophes. This also provides a systematic framework for calculating the field near caustics with arbitrary topologies. The general methodology is to map the original integral representation of the field, expressed either spatially or spectrally (e.g., Kirchhoff or plane-wave integrals, respectively), onto the diffraction catastrophe integral corresponding to the elementary catastrophe that describes the caustic. This integral, also referred to as canonical integral, constitutes a standard special function [29,Sec. 36.2].
The mapping of the diffraction integral into the canonical catastrophe integral may be applied globally via a systematic method first introduced in [30,Sec. 5]. This method ex-tends to multi-dimensions the method of Chester et al. [12] (see also [31,32]) for uniform asymptotic evaluation of 1D integrals with several nearby stationary points and singularities [31]. The global mapping approach, while very elegant and systematic, leads to implicit representations where the field is described in a generalized coordinate system that has to be calculated numerically. Alternatively, this mapping may be performed locally around any given point in configurational space, leading to explicit field representations where the generalized coordinates are the local Cartesian coordinates. In this paper we utilize the local mapping approach since our goal is to explore the wave mechanism of the 3D-AiB by deriving explicit field solutions where all the parameters are directly related to the aperture field distribution. We start in Sec. 5.A by mapping the exact AiB spectral integral in the vicinity of the curved beam axis, ending up with a hyperbolic umbilic diffraction integral for the field. Next in Sec. 5.B we extend this integral to the case of finite-energy AiB's by adding a complex loss parameter into the phase of this integral. In Sec. 5.C we calculate this integral by extending the closed-form series expansion of [29,Sec. 36.8] to finite energy AiB's. Finally, in Sec. 5.D we use this solution in order to check the validity of the conventional PE solution, and demonstrate that the latter suffers from a translational error due to an error in the beam axis under the PE approximation.

Non Dispersive Airy Beams
In a previous paper [11] we have dealt with the ray interpretation to the 2D-AiB and demonstrated that its curved propagation trajectory is a caustic of rays that emerge sideways from an extended aperture. Based on this concept we have extended the 2D-AiB to the ultra wide band (UWB) frequency domain and to the short-pulse time domain, by employing an initial field distribution that gives rise to a frequency-independent ray skeleton, thus ensuring that all frequency components propagate along the same trajectories and focus onto the same caustic [16].
Extending this concept to the 3D case, the non-dispersive AiB field in the half-space z > 0 is generated by the aperture field distribution in the z = 0 plane where (x, y) are the transversal coordinates, Ai is the Airy function, and a harmonic timedependence exp(−iωt) is assumed and suppressed throughout. k = ω/c is the wave-number and the wave-speed c is assumed uniform. The exponential window in Eq. (1) is added in order to render the energy of this distribution finite. The parameters β x , β y , α x , α y are frequency independent scaling factors, with β x,y controlling the widths of the field distribution along the x and y axes, respectively, and α x,y controlling the decay rates along these axes (see Fig. 1).
It should be noted that the conventional expression for the AiB initial condition of Eq. (1) is U 0 (x, y) = Ai(x/x 0 )Ai(y/y 0 ) exp(a x x + a y y) [5,Eq. (2)]. The scaling parameters in this expression are related to those in Eq. (1) via x 0 = β 1/3 x k −2/3 , y 0 = β 1/3 y k −2/3 and a x,y = kα x,y . We prefer the frequency independent β x,y parameterization since it renders the ray skeleton and the propagation trajectory frequency independent (see Eq. (10) and Sec. 3).
The radiating field at r = (x, y, z) due to the distribution in Eq. (1) has a closed-form solution within the parabolic wave equation (PE) which can be expressed as [5] U (x, y, z) = F (x, z; ω, β x , α x )F (y, z; ω, β y , α y ) exp(ikz), where F is given by The beam trajectory is defined by the vanishing of the arguments of the Airy functions in Eqs.
Note that the β x,y parameterization implies that the beam trajectory is frequency independent.
For later use we introduce the coordinate system (x,ỹ, z) such that the beam trajectory is contained in the (x, z) plane (see Fig. 1). This system is obtained by rotating the (x, y) system at the angle ϕ = tan −1 (β x /β y ) about the z-axis, viz where the superscript T defines the transpose of a matrix and R ϕ is the rotation matrix For simplicity, but without loss of generality, we shall consider the symmetrical case 1 Actually, the peak trajectory is described by replacing the 0 on the right hand side of Eqs. (4)-(5) by the asymptotically small parameter (kβ) −2/3 a 1 where a 1 −1.0188 is the first zero of Ai (x) [29,  such that ϕ = 45 • (see Fig. 1). The parameters in the (x,ỹ) system are defined conveniently asβ andα, and the beam trajectory (the "beam axis") is given bỹ The general case β x = β y can be reduced to the symmetrical case by the scaling x + β 2 y ) 1/2 /β y ] 1/3 and applying the exponential windows in the new coordinate system.
The solution in Eqs.
(2)-(3) is derived within the PE (paraxial) approximation and is therefore inaccurate at large distances where the AiB propagation trajectory curves. The ray solution of Secs. 4-5, on the other hand, remains accurate even at large distances since it is structured upon the ray skeleton of the exact wave equation. The ray solution also provides a clear physical interpretation for the AiB and a tool for constructing other curving beam solutions in uniform and non-uniform media.

3.A. Ray System
The properties of the AiB are fully described by the ray system that emerges from the aperture plane. We therefore start with the ray-system topology. Position r = x, y, z T along a ray can be expressed as where r = x , y , 0 T is the ray initiation point in the z = 0 plane, σ is the distance along the ray, andŝ defines the ray direction viẑ with θ x,y,z being the ray angles with respect to the x,y,z axes, respectively, and In order to findŝ for a given initiation point (x , y ) we substitute in Eq. (1) the asymptotic form of the Airy function, and recast the initial field in the form where A 0s (x , y ) = iβ 1/6 4πk 1/3 e kα(x +y ) with the ± signs in Eq. (16) taken as ++, +−, −+, −−, for s = 1, 2, 3, 4, respectively. In Eq. (15), M 0s = 0 for s = 1, M 0s = 1 for s = 2, 3 and M 0s = 2 for s = 4. We assumed here that α is small so that the decay term exp[kα(x + y )] is slowly varying on a wavelength scale and can be included in the amplitude A 0s rather than in the phase term ψ 0s . This assumption is justified in Sec. 4.B.
One finds that each point (x , y ) in the aperture plane, with x , y < 0, emits four rays, corresponding to s = 1, .., 4 in Eq. (15). The ray directionsŝ s are obtained by local matching to the phase gradient viz where the "±" signs are taken as explained in Eq. (16). Note that the ray directions are frequency independent if β is frequency independent. The directions of the four ray species projected onto the z = 0 plane are depicted in Fig. 2. Finally, we note that the result in Eq. (18) has been derived asymptotically via Eq. (14), but as shown in Appendix A.1, it is also valid exactly. The ray trajectories in Eq. (11) are parameterized in terms of the ray coordinates (ξ, η, σ) where the ray direction (ξ, η) tags the initiation point (x , y ) via Eq. (18). The topology of the ray system is described by the Jacobian of the mapping from the (x, y, z) coordinates to the ray coordinates (ξ, η, σ), defined by where r σ = ∂ σ x, ∂ σ y, ∂ σ z T and denotes the determinant. An explicit expression for J is obtained by noting from Eq. (18) that the exit point coordinates are defined by the ray coordinates via (x , y ) = −β(ξ 2 , η 2 ). Substituting in Eq. (11) and calculating the derivatives in Eq. (19) we obtain The mapping is singular when J = 0. This condition defines the envelope of the ray system, also termed caustic surface or catastrophe [24,28]. The structure of the caustic near the beam axis is discussed next in Sec. 3.B, with further details in Appendix A. The complete caustic topology including also regions away from the beam axis is presented in Sec. 6.

3.B.1. Caustic topology within the PE model
The PE approximation for the beam axis in Eq. (10) is valid only for moderate ranges before the axis curves too much. In this range one may approximate in Eq. (13) ζ 1 − 1 2 (ξ 2 + η 2 ). Substituting into Eq. (20) one obtains a closed form solution for the caustic equation J = 0 (see Appendix A.3). Figure 3 depicts the caustic structure within this approximation. It consists of two curved sheets joined at an edge that delineates the beam axis of Eq. (10). Fig. 4 shows a cross sectional cut of the caustic in two constant-z planes, where it is seen to form a right-angled corner. Figure 4 also depicts the exact caustics to be discussed next.

3.B.2. Exact caustic topology along the beam axis
Employing exact ray theory (see Appendix A.2) reveals further details of the caustic, which are also evident from the canonical catastrophe integral to be discussed in Sec. 5. Here we discuss the structure of the caustic near the beam axis; the complete global structure beyond this zone will be discussed separately in Sec. 6.
The exact caustic is shown in Fig. 5. It consists of two surfaces, referred to as "caustic 1" and "caustic 2," each consisting of two sheets joined at an edge. Caustics 1 and 2 are very close and hardly distinguishable on the figure scale. The caustic within the PE approximation is a degenerate case in which caustics 1 and 2 coalesce. Figure 4 zooms on cross-sectional cuts of the caustics at two constant-z planes, showing that caustic 1 is a cuspoid, whereas caustic 2 is smooth, and that the distance between these caustics increases with z. This topology is identified as a hyperbolic umbilic diffraction catastrophe (see e.g., Fig. 4.4 in [24]) whose focal plane is located at the z = 0 plane where the two surfaces coalesce along the semi-infinite lines {x = 0, y < 0} and {y = 0, x < 0} and form a right-angled corner at the origin.
The explosion of the PE caustic into a hyperbolic umbilic caustic at z > 0 can also be expected from catastrophe stability theorems [27,28] as follows from the fact that the rightangled corner at the cross-section of the PE caustic is not a stable form (the only stable forms in 2D are the fold and the cusp; see e.g., Table 3.1 in [24]). Therefore, the corner will not survive the perturbation due to exact propagation 2 . Furthermore, the simplest stable singularity in 3D which contains a corner at a particular cross-section (termed focal plane) is the hyperbolic umbilic catastrophe (see also [33] where a similar argument has been used). Figure 4 also shows that the deviation of the PE caustic from the exact one increases with z, thus limiting the range of validity for the PE solution in Eqs. (2)-(3). In Fig. 6 we concentrate on the difference between the edges of the exact caustics 1 and 2 (Eqs. (A.16)-(A.19)) and of the PE caustic (Eq. (10)) since this edge delineates the propagation trajectory of the AiB. Figures 7 and 8 describe the ray formation of caustics 1 and 2, respectively, by depicting the rays emitted from points in the aperture plane and their point of tangency at the caustics. Actually, the figures show traces of the exit points in the aperture (thick lines), and the traces on the caustics of the corresponding points of tangency (thin lines with the arrows indicating the direction of the traces). To emphasize the symmetry, traces withỹ ≷ 0 appear as full or dashed lines, respectively. Lines of constantỹ and of constantx are plotted in magenta and blue, respectively.
As discussed in Eq. (18) and Fig. 2, each point in the aperture emits four rays, belonging to species s = 1, ...4. As seen in these figures, caustic 1 is formed by ray species 1,2,3, whereas caustic 2 is formed by species 1 only. Specifically, referring to Fig. 7, the rays of species 2 and 3 touch only caustic 1 (red points), whereas the rays of species 1 penetrate through caustic 1 (through the white circle), then touch caustic 2 (green point; see also Fig. 8), penetrate again through caustic 1 (white circle) and finally touch caustic 1 (red point). In the limitỹ → 0, the rays of species 1 tend to the symmetry planeỹ = 0 thus giving rise to the somewhat counterintuitive situation in Fig. 10 where the ray in the symmetry plane penetrates through the edge of caustic 1 and touches caustic 2, but then as it penetrates back through the cusped-edge of caustic 1 it is also tangent to the sheets that generate the cusp (see the blue solid-line trajectory in Fig. 7). This situation will be discussed again in connection with Figs. 9 and 10.

4.A. Spectral Representation
Using the spectral representation of the Airy function [29, Eq. (9.5.1)], the initial field in Eq. (1) can be expressed as where exp[iω(ξx + ηy)/c] is the Fourier kernel, (ξ, η) are the spectral variables associated with (x, y) and The spectral integral in Eq. (21) has been expressed in a normalized non-dispersive form [18] where ω appears as a linear parameter in the phase and τ 0 is frequency independent. As a result (ξ, η) have a pure frequency-independent geometrical interpretation in terms of the plane wave angle (see discussion below), and τ 0 is the departure time of the (ξ, η) plane wave from the z = 0 plane. We use this normalization not only because its physical transparency but also since it may be extended to the time-domain [21]. The radiated field, obtained by adding the spectral propagator, is where ζ = 1 − ξ 2 − η 2 is the spectral wave-number in the z-direction, chosen with Im ζ ≥ 0 for ω > 0. Consequently the integration extends along the real ξ, η axes from −∞ to ∞ such that the ξ contour passes above and below the branch points ξ = ∓ 1 − η 2 corresponding to ζ. Equation (24) expresses the field as a spectrum of plane-waves propagating in the direction s(ξ, η) defined in terms of (ξ, η) via Eq. (12). Note that in Sec. 3, (ξ, η) were parameters that define the ray in Eqs. (11)- (12), whereas in this section, they are spectral variables that define the spectral-plane-wave direction. For simplicity we use the same notations and distinguish between them from their context. Actually, the two are related since the ray coordinates of Sec. 3 are stationary points in the spectral (ξ, η) domain, as explained after Eq. (26) and also in Appendix A.1.
As noted after Eq. (21), the spectral integral Eq. (24) has a non-dispersive form where ω appears as a linear parameter in the phase and τ is frequency independent, describing the arrival time of the plane wave that passes through r. This special form will be used in [21] to evaluate the field of the Airy pulsed beam (AiPB).
The direct numerical evaluation of the spectral integral in Eq. (24) is difficult since the integrand is very oscillatory at high frequencies. It may be evaluated asymptotically using the steepest descent method [34] which yields the geometrical optic (GO). Equation (24) will be also used in Sec. 5 in order to calculate the AiB within the framework of Catastrophe theory.

4.B. Geometrical Optics
The main contributions in the integral of Eq. (24) come from the vicinity of the saddle points, defined by with τ given in Eq. (25). We denote the solutions of Eq. (26) as (ξ r , η r ) where r is the rayindex. For observation points far away from the caustic, the saddle points are distinct and the field can be expressed in terms of the isolated saddle point contributions [34] U with τ r = τ at (ξ r , η r ). H r is the determinant of the Hessian matrix evaluated at (ξ r , η r ), and µ r is determined by the eigenvalues h 1r and h 2r of H, evaluated at (ξ r , η r ), via Expressions for (ξ r , η r ) and H are given Eqs. As discussed in Eq. (15), α is relatively small and hence we consider first the case α = 0. In this case, the saddle points (ξ r , η r ) corresponding to an observation points r on the lit side of the caustic are real and equal to the ray parameters in Eq. (12) of the rays reaching r (see proof in Appendix A.1). The ray-index r should not be confused with the species-index s introduced in Eq. (15): depending on r, the rays described by the rth saddle points can be of any species.
The result in Eqs. (27)- (29) can therefore be expressed in the GO form (see details in Appendix A.4) where σ r is the ray length in Eq. (11) and J r is the Jacobian in Eq. (19), both evaluated at (ξ r , η r ). All squared roots are taken as positive. M r is the Maslov index that counts the number of times ray r has touched the caustic. τ 0r and A 0r in Eqs. (30)-(31) are, respectively, the initial delay and amplitude of the rth ray, given by Eq. (16)- (17), with τ 0r = ψ 0s /c, both are taken at (x r , y r ), the initiation point of the rth ray. Note that the GO amplitude in Eq. (31) is invalid at the caustics where J r (σ r ) = 0. Fig. 9 depicts an illustrative example of the rays reaching a given observation point r. For simplicity we consider a point in the vicinity of the beam trajectory where all rays are of species 1 (see Fig. 2(a)). A point on the lit side of the caustic is reached by four rays, tagged by the index r = 1, 2, 3, 4. The Maslov index of Eq. (31) in this example is as follows: Ray r = 4 has not touched the caustic hence M 4 = 0; Rays r = 2, 3 have penetrated through "caustic 1," then touched "caustic 2", and then penetrated again through "caustic 1" (see also ray 1 in Fig. 7), hence M 2,3 = 1; Ray r = 1 is has penetrated through "caustic 1", touched "caustic 2" and then touched the cusp at the edge of "caustic 1" while penetrating it (cf. Fig. 10), so M 1 = 2. In addition to the four rays described above, the equation for the rays has three more solutions representing rays that emanate from the aperture at distant points where the field is very weak and the ray angles are very shallow with respect to the z = 0 plane. These rays will not be considered therefore in the analysis.
Finally, we note that the GO field in Eq. (27) has been derived via an asymptotic analysis of the exact spectral representation, but for σ = 0 it reduces exactly to the asymptotic initial conditions in Eq. (15) that were derived from the asymptotic approximation of the Ai function in Eq. (14).
If α = 0 then the saddle points become complex. The saddle point analysis above still applies, but the physical interpretation is less transparent. However, since typically α is very small, one may expend τ r and H r of Eq. (27) into Taylor series about the real saddle points (ξ r , η r ) for α = 0. Keeping terms to leading order in α one obtains the result in Eq. (27) with Eqs. (30)-(31) except that here A 0r → A 0r e kαx r (cf. Eq. (17); see details in Appendix A.4).

4.C. The field Near the Caustic
Near the caustic, the ray-system Jacobian vanishes and the GO amplitude in Eq. (31) explodes. From a spectral perspective, the relevant saddle points come close together and the second order asymptotic analysis in Eqs. (27)-(29) is no longer valid. One should resort, instead, to higher order asymptotics that may accommodate two or more nearby saddle points. The simplest situation occurs near the smooth parts of the caustics which are fold catastrophe where two saddle points coalesce. In these regions, the field may be described uniformly via the uniform geometrical optics (UGO) [13,14] where the field is expressed explicitly in terms of the ray-parameters (i.e., the phase and the Jacobian), but these parameters are combined together so that the solution remains regular near and on the caustic while reducing to the conventional GO expressions further away from it. The simple UGO is not valid near the beam trajectory where the caustic topology is a hyperbolic umbilic as discussed in Sec. 3.B.2 and Figs. 4-8. The field in that region is considered next in Sec. 5.

The Field Near the Propagation Trajectory
Near the beam trajectory, the caustic topology is that of a hyperbolic umbilic, hence, as explained in the introduction, we refer to the general framework of catastrophe theory in order to calculate the field. In this approach, the field is represented by the canonical catastrophe integral that can be written as where (x,ȳ,z) are generalized coordinates in the vicinity of the caustic such thatz is a coordinate along the caustic and (x,ȳ) are transversal coordinates, and (ξ,η) are the spectral variables associated with (x,ȳ). Φ has a canonical polynomial form for each elementary catastrophe (see e.g., Eq. (42) for the hyperbolic umbilic). Near the caustic, the integral in Eq. (32) can be evaluated by standard numerical techniques such as series expansion [29,Sec. 36.8] (see also discussion in Sec. 5.C), whereas further away from the caustic, it reduces asymptotically via saddle point analysis to the GO field. The spatial and spectral coordinates in Eq. (32) are typically normalized to be dimensionless and C 0 is a normalization constant such thatŪ = 1 at (x,ȳ,z) = 0.
In the theory of diffraction catastrophes [23,30], the exact integral representation of the field, expressed either spatially (i.e., via a Kirchhoff-type integral) or spectrally, is mapped asymptotically via a method outlined in [30,Sec. 5] into the canonical integral in Eq. (32), thus leading to a globally uniform solution of the field. 3 As discussed in the introduction, in the globally uniform approach, the generalized coordinates in Φ are related implicitly to the geometrical coordinates near the caustic. Our goal, however, is to obtain an explicit mapping such that these generalized coordinates are directly related to the local Cartesian coordinates along the beam trajectory as shown in Fig. 10, and also to the parameters of the aperture field distribution. Such explicit expressions are obtained by the locally uniform approach of Sec. 5.A where we concentrate on the region near the caustic. In this approach, we transform the exact spectral integral of Eq. (24) to the local Cartesian coordinate system in Fig. 10, and then approximate the spectrum paraxially around the caustic direction. This leads directly to the spectral integral in Eq. (32) where the phase is recognized as hyperbolic umbilic. This integral describes the field in the transition layer around the caustic, but its validity range is much larger and in particular it reduces to the GO field away from the caustic. As discussed in the introduction, this direct approach also provides a convenient tool for the synthesis of other types of beams.
As was done in Sec. 3, we consider first in Sec. 5.A the case where α = 0 where the rays are real. The modifications for the finite energy AiB where α = 0 are done in Sec. 5.B. In Secs. 5.C-5.D we discuss the numerical evaluation of the integral.

5.A. Mapping the Exact Spectral Integral to the Canonical Form
The canonical integral is constructed about a a typical reference point along the beam trajectory, which is taken to be on the edge of caustic 2 (see Fig. 10; recall from the discussion in Sec. 3 that the caustic topology near the beam axis is hyperbolic umbilic, consisting of two sheets denoted as caustic 1 and 2, each has an edge in the symmetry planeỹ = 0, see also Figs. [4][5]. This reference point is parameterized, conveniently, by θ z , the angle of the tangent to the edge at this point with respect to the z-axis. The coordinates at this point are given by (x 2 (θ z ), 0, z 2 (θ z )) of Eq. (A.15) (see also Eqs. (A.17) and (A.19)). The Cartesian coordinates around this point, denoted as (x,ŷ,ẑ), are chosen such thatẑ is the tangent to the edge and (x,ŷ) are the normal and the bi-normal, respectively (see Fig. 10), and are given, therefore, by where R θ is the rotation matrix in Eq. (7) with ϕ = θ z . The spectral integral in the (x,ŷ,ẑ) system is obtained from the one in Eq. (24). We start by expressing the latter in the (x,ỹ) system of Eq. (6) by rotating the spectral variables such that ξ ,η T = R 45 0 ξ, η T with R given in Eq. (7). Substituting into Eqs. (24)- (25) we obtain where withζ = (1 −ξ 2 −η 2 ) 1/2 . Recall that here we assumeα = 0. The spectral variables (ξ,η,ζ) corresponding to (x,ŷ,ẑ) are related to (ξ,η,ζ) via withζ = (1 −ξ 2 −η 2 ) 1/2 . Substituting Eqs. (33) and (37) into (35), expandingζ to a Taylor series about (ξ,η) = (0, 0), and keeping terms up to third order, we obtain the following expression for τ where Concentrating only on theẑ = 0 plane perpendicular to the edge, and normalizing the coordinates, the integral in Eq. (34) with the phase τ of Eq. (38) can be written in the form where τ a is given in Eq. (39) andŪ is the canonical integral in Eq. (32) with Φ(ξ,η;x,ȳ;δ) =ξ 3 +ξη 2 +ξx +ηȳ +η 2δ , The normalized spatial and spectral coordinates in Eq. (42) arē and the normalized parameters arē where in the second equality inδ we substitutedδ from Eq. (39), while the last term applies for θ z 1. The spectral phase Φ in Eq. (42), which is the simplest polynomial that describes the AiB in the vicinity of the caustic, is of the hyperbolic umbilic type [24, pp. 66], [28, pp. 41]. It follows that the field-scale near the caustic is governed by the factor (kβ) 2/3 that scales the local coordinates (x/β,ŷ/β) in Eq. (44). Note that the k 2/3 scaling is expected in view of the general theory of diffraction catastrophes [29,Sec. 36.6].
The caustic topology is controlled by the parameterδ of Eq. (46) which represents the normalized distance between caustics 1 and 2 such that it scales with frequency like (kβ) 1/3 and also grows with the observation range z. Forδ 1, caustics 1 and 2 are far apart and may be described by separated cusp and fold catastrophes, respectively. In this case, the hyperbolic umbilic integral may be reduced to the Pearcey and Airy integrals [29, Sec. 36.2(ii)] that represent, respectively, the field near these catastrophes.
Equation (41) describes the AiB field in a plane normal to the beam axis. The solution in the entire domain along the axis is obtained by tracking the local (x,ŷ,ẑ) coordinate system of Eq. (33) along the beam axis. The parameters in Eqs. (39)-(40) are therefore expressed in terms of θ z (z) which describes points on the beam trajectory (see Eq. (A.15)). Since this local spectrum representation is constructed entirely using the ray system topology near the caustic, it follows that the ray approach can be used to construct AiB-type solutions that follow any prescribed convex trajectory by back-projecting the rays to construct the relevant aperture distribution (e.g., the 2D example in [15]).
There are several alternative forms for the hyperbolic umbilic spectral phase Φ, all sharing the "germ", i.e., the highest order terms of Φ which depend only on (ξ,η) and define the catastrophe. The phase Φ in Eq. (42) contains the germξ 3 +ξη 2 , but there is another form used in [29,Sec. 36.3] which is discussed in Appendix B. The rest of the phase terms, called perturbation terms (or unfolding terms), can have different forms [35, chap. 4] without changing the catastrophe type, as long as they are small. For example, the phase used in [28, pp. 41] is the same as in Eq. (42) withη 2δ replaced byξ 2δ , but since these are perturbation terms, the two forms are equivalent. In fact, the latter can be obtained by setting the reference point of the local coordinates on the edge of caustic 1 and taking thê z axis along the direction of the ray that touches the cusped edge there (see pointẑ =δ in Fig. 10). We prefer the form in Eq. (42) since in the alternative one,ẑ is not tangent to the edge of caustic 1 (see dashed line in Fig. 10) and the tracking along this edge becomes more complicated. Finally, in Sec. 5.B below we extend Φ to the finite energy AiB case.

5.C. Numerical Evaluation of the Canonical Integral
The canonical integral in Eq. (32) with Φ in Eq. (42) or (54) can be evaluated by a straightforward numerical integration. To accelerate the convergence, the integration contours for |ξ|, |η| 1 can be deformed into the complex plane so as to reach infinity along the asymptotic valleys of exp(iΦ). Fastest convergence is obtained if the ξ-integration starts asymptotically along the line arg ξ = 5π/6 and ends along the line arg ξ = π/6. Likewise, the η integration may start asymptotically along the line arg η = −3π/4 and ends along the line arg η = π/4. Another approach is to use the 1D integral representation for the hyperbolic umbilic as explained in [29,Sec. 36.15(iii)].
Using (iξ) ↔ ∂x,(iη) ↔ ∂ȳ in Eq. (57) and taking the derivatives outside of the integral, we rewrite Eq. (57) asŪ I(x,ȳ) 1 (2π) 2 dξdη e iΦ(ξ,η;x,ȳ) with (ū,v) given by where R is given in Eq. (7). Equations (59)-(61) can readily be calculated using standard numerical routines. Note that in the limit¯ → 0, the series representation in Eq. (59) reduces to the standard series expansion in [29, Eq. 36.8.3]. In Fig. 11 the canonical integral solution and the PE solution are in an excellent agreement on both the lit and the shadow sides of the caustic. As expected, the GO solution breaks on the caustic, but otherwise it agrees well with the other solutions. Further away on the lit side of the caustic, the GO solution starts to deviate, but in this region it is, in fact, more accurate than the other solutions since it does not rely on the local spectrum approximation or on the PE approximation.

5.D. Numerical Example
In Fig. 12, on the other hand, on observes a significant translational error of the PE solution (notice, on the other hand, the perfect agreement between the canonical integral and the GO solutions, except on the caustic). The error of the PE solution is due to the deviation of the PE-based beam trajectory in Eq. (10) from the true ray-based trajectory, and can be parameterized byx p which is the normalized location of the PE caustic in thē x-system (recall from Fig. 10 that thex-system is centered on the edge of caustic 2, see also Fig. 4).x p can be found using Eq. (10) in conjunction with Eqs. (33), (A.14)-(A.15) and (44). Here we consider only the small θ z approximations ofx p , obtaininḡ In Fig. 11 the frequency is relatively low such thatx p = −2.7 × 10 −2 (i.e.,x p = −6 × 10 −2 λ), hence the error of the PE solution is hardly noticeable on the figure scale. In Fig. 12, on the other hand, the frequency is higher givingx p = −0.59 (i.e.,x p = −6λ). Sincex p is O(1), the translational error of the PE solution is readily discerned in the figure.
In addition to the translational error of the PE solution, the exact canonical integral has some intricate features which are readily discerned in the main-lobe of Fig. 12 (notice the peak 1% error of the PE-solution there). Recalling the discussion after Eq. (46), the field structure of the canonical integral is parameterized byδ which is the normalized distance between caustics 1 and 2. As noted there, this parameter grows with both frequency and observation range. Forδ 1 (near the aperture or small kβ), the distance between these caustics is small and the canonical integral reduced to the Ai function field structure in Eqs. (2)-(3), but for largerδ, the field structure of the canonical integral describes the separation between the two caustics and is more complicated than the PE solution. This is seen in Figs. 11 and 12 whereδ = 7.6 × 10 −3 and 3.5 × 10 −2 , respectively. In these figures, the cusped edge of caustic 1 is located atx 1 = −1.7 × 10 −4 andx 1 = −3.7 × 10 −3 , respectively (i.e.,x 1 = −4 × 10 −4 λ andx 1 = −4 × 10 −2 λ).
Finally, to help the reader we map the parameters used in these numerical examples to the conventional ones used for example in [2,5] (see also discussion after Eq. (1)). Noting that the "Fresnel distance" z 0 = kx 2 0 used there is related to ourβ parameter via kz 0 = 2 1/3 (kβ) 2/3 , it follows that the range at the caustics in Figs. 11 and 12 is z = 2.7z 0 and z = 12.5z 0 , respectively.

Discussion on the Complete Caustic Topology
We consider the global caustic topology which includes far away parts of the caustic which were not included in the analysis of Sec. 3.B. The field in these regions is formed by contributions of rays that emanate sideways from far away points in the aperture. This global picture is presented here in connection with the discussion in Sec. 6.C regarding the limitations of the AiB. Since the caustic structure is rather complex, we discuss separately the topologies of caustic 2 and 1 in Secs. 6.A and 6.B, respectively.

6.A. Global Structure of Caustic 2
The global structure of caustic 2 is depicted in Fig. 13, as a zoom-out view of Fig. 5. In addition to the sheet shown in Fig. 5, tagged here "caustics 2A", there is another caustic surface which is tagged "caustic 2B." These two surfaces are connected and form a pyramidlike structure which terminates at a sharp edge at the top. Caustic 2B is formed by rays exiting the aperture at far away regions, as illustrated in Fig. 16 which shows a cross section of the caustic surface in theỹ = 0 (the cross sectional line is also depicted in Fig. 13 by a red dashed-line).

6.B. Global Structure of Caustic 1
Figs. 14(a) and 14(b) depict the complete structure of caustic 1 as a zoom-out view of Fig. 5, with Fig. 14(a) showing the caustic parts formed by ray species s = 1 (see Fig. 2 for the species definition), and Fig. 14(b) showing the rest of the caustic. The surfaces in Figs. 14(a) and 14(b) are connected at the margins that are delineated in the figures. The connections of these surfaces are clarified in Fig. 15 which depicts cross-sectional views of the complete caustic at two constant-z planes. A cross sectional view in the symmetry planeỹ = 0 is depicted in Fig. 16.
As can be discerned from Fig. 14(a), the part of caustic 1 that is shown in Fig. 5, tagged here as "caustics 1A", is connected to a hat-like structure, tagged as "caustic 1B." The latter is formed by rays of species 1 exiting the aperture from far away points, as in the case of caustic 2B of Fig. 13. Caustic 1B is a cuspoid (as may also be seen in Fig. 15) with its edge in theỹ = 0 plane (red dashed line in Fig. 14(a)).

6.C. Discussion: Limitations of the AiB Solution
As discussed in Sec. 3.B, the AiB has been identified as a hyperbolic umbilic diffraction catastrophe. However, as can be seen from Figs. 13-14(a), near the termination point of caustic 2, the hyperbolic umbilic caustic is no longer isolated, and the caustic structure evolves into the more complex parabolic umbilic catastrophe [24,Fig. 7.3]. From this global analysis it then follows that the propagation range of the AiB is limited to the termination point of caustics 2A (see also Fig. 16) where the AiB loses its beam shape. The maximal propagation range is of orderβ and hence it can be controlled by choosingβ for a desired application. The range is also controlled byα which controls the decay rate along the beam axis (see Eqs. (56) and (50)).
The exact ray-based solutions of Secs. 4 and 5 may be used up to the maximal propagation range discussed above. In the near zone, one may also use the PE solution in Eqs. (2)-(3), but the accuracy of this solution reduces as the distance between the exact and the PE approximated caustics grows (see Fig. 6), giving rise to the translational error of the PE field in Fig. 12.

6.D. Comments on the Range Limitations of the 2D AiB
In [11] and [16] it has been shown that the 2D-AiB is supported by a smooth caustic which is a fold catastrophe [28]. Those references, however, have concentrated on the field near the beam axis. A global analysis reveals that this 2D caustic is actually a cusp whose form is very similar to the cross section of caustic 2 in Fig. 16 (caustic 1 does not exist in the 2D case). This global structure was not discussed in [11] and [16] since the rays forming the missing parts have weak contributions and late arrival times beyond the pertinent time window [16]. The analysis in Sec. 6.C regarding the limitations of the AiB can be similarly applied to the 2D case.

Conclusions
We presented a wave-theoretic analysis of the 3D-AiB and derived uniform solutions for the transition layer near the propagation axis. The formulation is based on exact spectral and ray representations and therefore does not suffer from the limitations of the conventional PE approximation. As in the 2D-AiB case, the 3D-AiB is described by sideways radiating rays that focus onto a caustic that delineates the beam trajectory. However, the 3D ray topology is considerably more complicated than the 2D topology, and it has many intricate features. Specifically, near the curved beam axis, the ray topology of the 3D-AiB consists of a cusped double-layered caustic identified as the hyperbolic umbilic catastrophe (Sec. 3). At a certain range away from the aperture, the hyperbolic umbilic caustic is no longer isolated and evolves into the more complex parabolic umbilic caustic, where the AiB loses its beam shape.
To clarify these issues we have investigated the ray topology of the 3D-AiB using the exact spectral representation for the AiB. We concentrated first on the region near the curved beam axis (Sec. 3), but we have also explored the global caustic topology (Sec. 6). Specifically, we have thoroughly explored all the ray paths that form the double-layered hyperbolic umbilic caustic noted above (Figs. 7-10). These contributions are essential in order to understand the time-dependent Airy pulsed beam (AiPB) that will be presented in a subsequent publication [21].
The field calculations in Secs. 4-5 were based on an asymptotic evaluation of the exact spectral integral of the AiB. In order to calibrate the final results, we derived first the lowest order geometrical optics (GO) solution (Sec. 4.B) which, however, fails near the beam axis. The field near the beam axis has been considered separately in Sec. 5 within the framework of catastrophe theory [24,27,28]. As discussed in the introduction (see also Sec. 5), a globally uniform field solution can be obtained via the procedure in [30,Sec. 5]. This approach, while very elegant and systematic, leads to implicit representations where the field is described in a generalized coordinate system that has to be calculated numerically. We therefore used in Sec. 5.A an alterative local approach where the spectral integral of the AiB has been approximated near the caustic using a local coordinate system, and the final result has been expressed as a canonical hyperbolic-umbilic diffraction integral [29,Sec. 36.2]. This more direct approach does not have the global uniformity of the general implicit procedure in [30], yet it has led to a solution that is expressed explicitly in the local coordinate system of the caustic, where all the field parameters are expressed in terms of the geometrical parameters of the caustic and of the aperture field distribution. This solution is valid in the transition layer near the caustic and reduces correctly to the GO field away from the caustic, but far away, one has to switch to the GO field of Sec. 4.B noted above. Using this approach we also modified the hyperbolic-umbilic integral by introducing a complex loss parameter that describes the finite energy AiB.
The final field solution in Eqs. (56) with (32) and (54)) provides an explicit solution in the transition layer near the beam axis. The width of the transition (or boundary) layer is of order k −2/3β1/3 (seex,ȳ coordinates in Eq. (44)) where k is the wave-number and β is a frequency-independent parameter describing the beam geometry (see Eq. (10)). The caustic topology is controlled by the dimensionless parameterδ of Eq. (46) appearing in Φ of Eq. (54). This parameter represents the normalized distance between the two layers of the hyperbolic umbilic caustic, scaled with the normalized frequency like (kβ) 1/3 . Asδ grows with the observation range z, the structure of the boundary layer field described by the canonical integral becomes more complex while forδ 1, the two layers are separable and the field can be described by separate cusp and fold catastrophes, as discussed after Eq. (46).
The canonical hyperbolic umbilic integral can be evaluated efficiently via contour deformation in the complex plane. Here, however, we calculated it using the series expansion of [29,Sec. 36.8] which has been modified to accommodate the finite energy AiB (Sec. 5.C). We used it to calculate the field near the beam axis under several parameter regimes (Sec. 5.D). The accuracy of the canonical integral has been demonstrated by noting in Fig. 11 that away from the beam axis it reduces to the GO field (the GO field near the hyperbolic umbilic caustic is a sum of four rays). Far away from the transition layer, however, the field is described more accurately by the GO solution of Sec. 4.B, since it is structured upon the exact ray skeleton and not on the local ray skeleton near the beam axis used here for the construction of the transition function (note that the solutions obtained via the globally uniform approach in [30,Sec. 5] reduce correctly to the GO field even very far from the transition layer, yet, as noted earlier, these solutions are implicit functions of the space coordinates). Comparing the canonical integral to the conventional PE solution one observes in Fig. 12 a significant translational error of the latter which is due to the error of the PE approximated beam trajectory. This error is parameterized by the dimensionless parameterx p of Eq. (62) which describes the distance between the exact and the PE caustics.x p grows with z (see Figs. 4 and 6) and scales with the normalized frequency like (kβ) 2/3 . As discerned from Figs. 11-12, the error is small forx p 1, but it becomes significant forx p ∼ 1. In summary, it has been demonstrated that the local canonical integral in conjunction with GO provides an accurate and convenient solution to the AiB, as an alternative to the conventional PE solution. This solution is structured upon the exact ray skeleton and therefore has a wider validity range than the PE solution. Furthermore, this solution provides a convenient design tool for synthesizing other types of beams propagating along any prescribed convex trajectory in uniform or non uniform medium. An application of this idea in 2D configurations has been presented in [15].
The analysis in this work has been performed within the framework of the non-dispersive AiB where the aperture field is scaled with frequency such that the ray skeleton is frequencyindependent. We use the non-dispersive formulation because it simplifies the ray analysis by making it more transparent geometrically. Moreover, this formulation also extends the AiB solution to the ultra wide band (UWB) frequency regime by ensuring that the pulsed field propagates along the curved beam trajectory without dispersion. This non-dispersive formulation will be used in a subsequent publication [21] to introduce a new class of 3D Airy pulsed beams (AiPB).

A. Appendix A: Derivation of the GO Field from the Spectral Representation
In Sec. 3 we introduced the ray representation of the AiB based on the approximate initial field distribution in Eq. (15). In this appendix we derive the ray representation based on the exact spectral representation of the initial field distribution in Eq. (24), and show that the results are the same as those of Sec. 3.
As discussed in Eq. (15), α is relatively small and hence we consider first the ray and caustic equations for α = 0 (Secs. A.1-A.3) and then we consider the GO field for the case α = 0 (Sec. A.4).

A.1. The Ray Equation
We start by deriving the ray equation in Eq. (11) from the saddle point condition of the exact spectral representation in Eq. (24).
Substituting the spectral delay function τ of Eq. (25) with α = 0 into the saddle points condition in Eq. (26) we obtain where ζ is defined in Eq. (13). which are, in fact the expressions for the ray trajectories in Eq. (11) with σ = z/ζ. It follows that the ray skeleton derived in Sec. 3 using the approximation in Eq. (14) is the same as the one derived here from the exact spectral representation.

A.2. Ray and Caustic Topology -Exact Analysis
There are two alternative routs to derive the expression for the caustic of Sec. 3. The first is to calculate J of Eq. (19) using the ray trajectories in Eq. (11), and then to solve for its zeros. The second relies on the spectral formulation of Sec.4.A. Here we start with the latter since it also complies with the formalism of catastrophe theory which is discussed in Sec. 5.
To calculate the caustic surface we note that the elements of the Hessian matrix H in Eq. (28) are given by and hence, the determinant H is given by At the caustic, two or more rays (saddle points) coalesce, implying that H = 0 there. This equation should be solved in conjunction with Eqs. (A.1)-(A.4) since the caustic is composed of rays. From Eq. (A.10) a ray with parameters (ξ, η) is tangent to the caustic at two points The (x, y) coordinates of these points are found by substituting Eq. where sin θ z = (ξ 2 + η 2 ) 1/2 , θ z being the ray-angle with respect to the z axis (see Eq. (12)). Alternatively, z 1,2 in Eqs. (A.14)-(A.15) can be expressed as function ofx where all square roots taken with a positive real part. At moderate ranges such thatx/β 1 these expressions may be approximated by (10), hence one finds for a givenx that z 2 < z 1 < z p (see also Fig. 6). Finally we note that the same results could have been obtained via the ray-system analysis. The Jacobian in Eq. (20) can be expressed explicitly as . σ c1,2 = z c1,2 /ζ define the ray lengths to the caustic.

A.3. The Caustic within the Parabolic Equation Approximation
In the paraxial regime one can approximate in Eq.

B. Appendix B: Additional Canonical Form for the Hyperbolic Umbilic Catastrophe
In order to prevent possible confusion we note that the Hyperbolic Umbilic catastrophe has another canonical form in addition to the one in Eq. (42), which has the following form [ where (x,ȳ) andδ are defined in Eqs. (44)-(46), (ξ,η) are defined in Eq. (45), and R is given in Eq. (7). Note that some terms which do not depend on (ξ,η) have been omitted since they do not affect the catastrophe type.   Fig. (a) along the x and y axes are given, respectively, by x = a n β(βk) −2/3 , y = a n β(βk) −2/3 , where a n with n = 1, 2, ..., are the zeros of Ai given in [29, 12) and (18) onto the z = 0 plane (white arrows). Also shown for reference the intensity |U | 2 from Fig. 1(a). Figures (a)-(d) show the four ray species defined in Eq. (18) such that (a) θ x , θ y ∈ [0, π/2]; (b) θ x ∈ [π/2, π] and θ y ∈ [0, π/2] (c) θ x ∈ [0, π/2] and θ y ∈ [π/2, π] (d) θ x , θ y ∈ [π/2, π].  Fig. 3. Gray and green lines: the exact caustic, which consists of two surfaces denoted as "caustic 1" and "caustic 2" respectively (see Sec. 3.B.2 and Fig. 5). Thex axis in this figure is centered aboutx 1 (z) the edge of caustic 1 in Eq. (A.16). The insets zoom on the area near the origin where caustic 1 has a "cusped edge" while caustic 2 has a smooth corner. One observes that the distance between caustics 1 and 2 and the PE approximated caustic increases with z (see also Fig. 6).    The notations are the same as in Fig. 7. Note that caustic 2 is formed only by ray species 1. The parts of the rays before and after r are depicted as solid and dashed lines, respectively, and the parts between caustic 1 and 2 are depicted as dotted lines. The ray initiation points (x r ,ỹ r )/β are: ray 1 (−0.14, 0)10 −2 ; rays 2,3 (−0.95, ∓.81)10 −2 ; ray 4 (−1.8, 0)10 −2 . The ray direction parameters (ξ r , η r ) are: ray 1 (2.7, 2.7)10 −2 ; ray 2 (2.7, 9.4)10 −2 ; ray 3 (9.4, 2.7)10 −2 ; ray 4 (9.5, 9.5)10 −2 . The paths of rays 2 and 3 are similar to that of ray s = 1 in Fig. 7: they penetrate caustic 1 (through the white circles), touch caustic 2 (green points), penetrate caustic 1 again (upper white circles) and then touch caustic 1 (red points) after passing through r (blue point). Rays 1 and 4 propagate in theỹ = 0 plane. Ray 4 touches the caustic after passing through r: It penetrates caustic 1 (through the white circle), touches the edge of caustic 2 (green point) and then penetrates back through the cusped-edge of caustic 1 and at the same point it is also tangent to the cusp (red point; see discussion in the last paragraph of Sec. 3.B.2 and also Fig. 10). The path of ray 1 is similar to that of ray 4, only it already touched the caustics before passing through r. When r moves toward the cusped-edge of caustic 1 from its lit side, rays 1-3 coalesce. When r is on the shadow side of edges 1 such that it is between the edges of caustics 1 and 2, it is reached by 2 rays which penetrate caustic 1, one which results from the coalesced rays 1-3 and has already touched caustic 2, and ray 4 which has not touched caustic 2 yet. These two rays coalesce as r moves toward the edge of caustic 2.  Fig. 4, caustic 1 is actually a cuspoid but only its edge is shown here. A ray (dashed black line) penetrates through caustic 1 before the point of tangency with caustic 2. After that point it penetrates again through caustic 1 (pointẑ =δ in the figure). At that point it is also tangent to the cusp of caustic 1 (see see discussion in the last paragraph of Sec. 3.B.2). Also shown, the local Cartesian coordinate system (x,ŷ,ẑ) centered around a point on caustic 2, such thatẑ is the tangent to the edge at that point, whilex andŷ =ỹ are the normal and the bi-normal there, respectively. This system is parameterized by the angle θ z of theẑ-axis with respect to z-axis, and it is used to construct the canonical integral of the field in Sec. 5.   The global structure of caustic 2 (green surface). "Caustic 2A" is the part that was shown in Figs. 5 and 8, whose edge follows the beam axis. This part of the caustic is formed by rays leaving the aperture at points close to the caustic. "Caustic 2B" is an additional part of caustic 2 not shown in Fig. 5 and it is formed by rays leaving the aperture sideways from faraway points (see Fig. 16). At some range, caustic 2 terminates at a pyramid-like edge. Red dashed-line: cross section of the caustic atỹ = 0 which is also shown as red dashed-line in Fig. 16.  , respectively). In (a), caustic 1A is the caspoid shown in Figs. 5 and 7 which is formed by rays leaving the aperture at points near the beam axis, while caustic 1B is formed by rays leaving the aperture sideways at faraway points. In (b), caustic 1C (blue) and caustic 1D (magenta) are two overlapping sheets connect to caustic 1B at margins 3 and 4 and to caustic 1A at margins 1 and 2 (see also cross-sectional cuts of the caustic in Fig. 15). . Caustics 1C and 1D of Fig. 14(b) are not formed by rays contained in theỹ = 0 plane and are therefore omitted here.