Exploring the dynamics of finite-energy Airy beams: A trajectory analysis perspective

: In practice, Airy beams can only be reproduced in an approximate manner, with a limited spatial extension and hence a finite energy content. To this end, different procedures have been reported in the literature, based on a convenient tuning of the transmission properties of aperture functions. In order to investigate the effects generated by the truncation and hence the propagation properties displayed by the designed beams, here we resort to a new perspective based on a trajectory methodology, complementary to the density plots more commonly used to study the intensity distribution propagation. We consider three different aperture functions, which are convoluted with an ideal Airy beam. As it is shown, the corresponding trajectories reveals a deeper physical insight about the propagation dynamics exhibited by the beams analyzed due to their direct connection with the local phase variations undergone by the beams, which is in contrast with the global information provided by the usual standard tools. Furthermore, we introduce a new parameter, namely, the escape rate, which allow us to perform piecewise analyses of the intensity distribution without producing any change on it, e.g., determining unambiguously how much energy flux contributes to the leading maximum at each stage of the propagation, or for how long self-accelerating transverse propagation survives. The analysis presented in this work thus provides an insight into the behavior of finite-energy Airy beams, and therefore is expected to contribute to the design and applications exploiting this singular type of beams.


Introduction
In 1979, Berry and Balazs found [1] that the free-propagation Schrödinger equation, a rather simple and even trivial equation, admitted a rather counter-intuitive nondispersive solution.Provided the initial ansatz is described by an Airy function, these authors proved that the time-evolution of the corresponding probability density remains shape-invariant during the whole propagation.At the same time, without the presence of any external potential, these wave packets undergo a self-accelerating motion pushing them forwards.When this behavior is analyzed within the framework of the equivalence principle, as it was done by Greenberger [2], it can be shown that Schrödinger's equation can be recast as the time-independent equation describing a free falling mass, of which the eigensolutions correspond precisely to Airy functions.Later on, Unnikrishnan and Rau [3] proved the uniqueness of these solutions to Schrödinger's equation as the only free-space self-accelerating nondispersive one.
The beams conjectured by Berry and Balazs have not become a reality until their experimental realization, formerly produced in 2007 by Siviloglou et al. [4], taking advantage of the isomorphism between the free-space Schrödinger equation for a non-relativist particle with a mass , and the paraxial Helmholtz equation for monochromatic light propagating in a linear medium with a refractive index .Ever since, Airy beams have become one of the relevant components in the puzzle of structured light [5,6], a very active and promising research field strongly developed over the last decade due to the interesting fundamental and technological applications brought in by conveniently designed light beams.In this regard, different control mechanisms for their spatial manipulation [7] and applications [8] have been discussed in the literature.
From the point of view of their experimental implementation, the production of Airy beams introduces an important technical challenge: how is it possible to reproduce such beams, with considered in the analysis.In Sec. 3, we present and discuss the mains results obtained for the three aperture functions considered.The simulations presented include both propagation properties and also information extracted from trajectory statistics.Finally, Sec. 4 summarizes the main conclusions from this work.

General aspects
Consider the paraxial monochromatic scalar field where r ⊥ = (, ) and  = 2/ 0 (with  0 = /).The complex-valued amplitude (r ⊥ , ) satisfies the paraxial Helmholtz equation, where ∇ 2 ⊥ =  2 / 2 +  2 / 2 is the transverse Laplacian.If this equation is multiplied on both sides by  * (r ⊥ , ), and then we subtract to the resulting equation its complex conjugate version, we obtain Rearranging terms on the right-hand side of this equation, it can be further simplified in the form of the following transport equation which describes the transversal changes undergone by the intensity distribution, || 2 , as it propagates along the -direction.Such changes are determined by the term on the right-hand side, where the flux density [14] is associated here with a transverse electromagnetic flux vector, Actually, according to the usual definition of a flux density, the second line of Eq. ( 5) is specified in terms of the intensity distribution and the transport field This field quantifies the local rate of change of the flux, analogous to a local velocity field describing the flow of energy on a given point r ⊥ at a given value of the  distance.Hence, from now on, we shall refer to v ⊥ as the effective transverse velocity field.
It is worth noting that, although here we have proceeded directly from the paraxial Helmholtz equation, Eq. ( 2), in order to reach the expression for the flux (5), in analogy to how we could have proceeded in quantum mechanics with Schrödinger's equation to obtain the quantum flux [15], it could have also be done considering the full vector nature of the electromagnetic field and then assuming paraxiality conditions, as formerly done by Allen et al. [16].In this latter case, one explicitly observes that the time-averaged Poynting vector consists of two contributions, one simply related to the forward propagation in the paraxial case, and another one dealing with the space variations undergone by the field amplitude, which, in turn, are connected to the field local phase variations.This convenient decoupling between the longitudinal and the transverse contributions in the Poynting vector has been considered more recently, for instance, by Broky et al. [17] to precisely investigate the self-healing properties of Airy beams in two dimensions.
As it can be noticed, if the complex field amplitude is recast in the following polar form where both (r ⊥ , ) and (r ⊥ , ) are real-valued fields, the above defined velocity field reads as This expression unambiguously shows that any change in the transverse flux is directly associated with the local (transverse) spatial variations undergone by the phase of the field phase during its propagation, which eventually leads to measurable changes in the intensity distribution.Accordingly, to better understand these physical consequences, one can compute ensembles of trajectories obtained from the integration of the equation of motion which will provide us with first-hand local information directly associated with the phase of the field (or, more strictly speaking, its spatial phase variations) rather than from the profile displayed by the intensity distribution (from which any phase effect is indirectly inferred a posteriori).
In this regard, these trajectories have already been used to describe light propagation through optical fibers [18], to account for the differences between Gaussian beams and Maxwell-Gauss beams [19], or, more recently, to describe the dynamics of ideal and finite-energy Airy beams [13] or to investigate light focusing [20].Furthermore, it has also been used to establish a connection between quantum and optical analogs given the isomorphism between the paraxial Helmholtz equation and the time-dependent Schrödinger equation [15].In analogy to their quantum-mechanical counterparts [15], the trajectories obtained after integrating Eq. ( 9) have three distinctive properties with practical interest concerning the problem that we are dealing with here: (i) The trajectories cannot get across the same point r ⊥ at the same value of .This is a direct consequence that arises from the single valuedness of the solutions of the paraxial Helmholtz equation (2), except for a constant phase factor.From the solution in polar form (7), we note that this phase factor corresponds to an integer multiple of 2, which does not affect the trajectory dynamics, as it can be inferred from the Eq. ( 8).
(ii) Consider a closed curve C(r ⊥ , 0) at  = 0 defined by a set of initial conditions r ⊥ (0) for Eq. ( 9).At any other value of , C(r ⊥ , 0) will uniquely evolve into another closed curve C(r ⊥ , ) consisting of the positions r ⊥ () of the trajectories started at r ⊥ (0).
(iii) Following properties (i) and (ii), the partial energy content within the surface Σ(r ⊥ , ) determined by the curve C(r ⊥ , ), defined as remains constant regardless of the value of  and the shape acquired by C(r ⊥ , ) along its propagation [21].Note, on the right-hand side of the second equality, that  Σ (r ⊥ ,) is directly related to the set of all trajectories laying within C(r ⊥ , ) and with initially distributed within C(r ⊥ , 0) according to Specifically, as it is shown below, here we use these properties to quantitatively analyze some aspects associated with the leading peak of the input finite-energy Airy beams.

Finite-energy Airy beams
Let us consider the particular case of the one-dimensional Helmholtz equation, although this will be done, for simplicity, in terms of the reduced coordinates x = / 0 ( 0 denotes an arbitrary typical value for the transverse coordinate) and z = / 2 0 =  0 /2 2 0 [15], (to simplify notiation, from now on tildes will be removed).As Berry and Balazs showed [1], for wave packets with an initial amplitude described by an Airy function, propagating according to Eq. ( 11), their amplitude remains shape invariant and undergoes a transverse displacement that increases rapidly with  2 : Alternatively, the propagation of these beams can be described at a local level in terms of trajectories associated with the transverse energy flux [19], Specifically, the trajectories arise from the integration along  of the equation of "motion" where (, ) is the local field describing the energy flow along the transverse direction () at a given value of the longitudinal coordinate (), and Because the field ( 14) is directly related to the beam local phase variations along the transverse direction, it becomes a suitable tool to explore any property associated with the phase and, more specifically, with effects due to possible variations induced on it (see, for instance, Refs.[13,20]).Note that for ideal Airy beams, (, ) depends linearly on the -coordinate and hence the associated trajectories are fully analytical: These flux trajectories are parabolas, in agreement with the renowned parabolic transverse acceleration undergone by the beam during its propagation along .
In the case of non-ideal Airy beams, that is, realistic beams with a finite energy content, because the infinite tail of the ideal Airy beam is not physically realizable in practice, the corresponding truncation introduces extra phase factors that eventually lead to the gradual dispersion of the beam.In other words, the beam will gradually lose both its shape invariance and its transverse self-accelerated propagation beyond a certain spatial range.In this regard, the above trajectory-based approach constitutes a suitable tool to detect and explore any change directly related to the phase at a local level.
The truncation can be introduced by directly acting on the tail of the beam, as it was formerly proposed by Siviloglou and Christodoulides [9], who considered an exponentially decaying factor acting on an input (ideal) Airy beam, i.e., (, 0) = ()   . (16) The propagated beam in this case can be analytically determined and reads as where the complex argument  =  −  2 /4 +  in the Airy function makes the latter to become also complex, thus leading to the appearance of another phase factor, apart from those arising from the multiplying exponential prefactor.More specifically, the expression for the velocity field becomes [13] (, ) which leads to a non-analytical equation of motion, where we clearly observe that the extra term (second term on the right-hand side) is responsible for the deviation of the trajectories (flux) from the typical parabolic shape associated with an ideal Airy beam.
In the above case, the field can be represented as a product of the ideal Airy beam and an exponential.Hence, the corresponding transverse spectrum, is the direct product of the transverse spectrum corresponding to the ideal Airy beam, and a complex Gaussian function.Taking into account this result and that, according to the convolution theorem, the Fourier transform of the convolution of two functions is the direct product of the Fourier transforms of such functions, we proceed now as follows to investigate the effects induced by known aperture functions on the flux dynamics.Consider that the finite-energy Airy beam results from the convolution function where () represents a given aperture function.From a formal point of view, here we are going to assume that () satisfies the following two properties: () it is  2 in order to warrant the finiteness of the energy content of (), and () under some control parameter it approaches a Dirac -function, which allows us to recover the original (ideal) infinite energy Airy beam.Thus, according to the convolution theorem, the Fourier transform of (21) reads as from which we already infer that any additional phase effect is going to be associated with the presence of the new term w(  ), as the Airy prefactor only governs the curvature of the flux.More importantly, in the current case, note that the transverse power spectrum of the ideal Airy beam is totally flat, so the parabolic curvature of the flux only comes from the cubic phase factor of the Fourier transform.The appearance of w(  ) in ( 22) is not only going to affect the phase, as seen in ( 18), but it may also induce a modulation of the spectral components, thus removing the flatness of the power spectrum.Both factor are therefore responsible for the eventual diversion of the propagation of the finite-energy beam from the one displayed by an ideal Airy beam.
Taking into account the functional form of a propagated Airy beam, Eq. ( 11), the amplitude (21) can be recast as where The amplitude (23) leads to a generalized equation of motion that now reads as containing the action of an extra term on the right-hand side.This extra term is directly connected with the spectrum of the aperture functions, which induces a new local phase configuration apart from the truncation of the beam, thus leading to the distortion of the beam with respect to its ideal counterpart.So far the trajectories obtained after integrating Eq. ( 25) provide us with an idea of the newly induced phase effects on the flux dynamics.However, we can also extract valuable quantitative information regarding the finite-energy beams if we consider particular ensembles of such trajectories and the above properties satisfied by the trajectories arising from Eq. ( 25).Here, in particular, we consider two quantities to investigate in a more quantitative manner the dynamics associated with the leading maximum of the input finite-energy Airy beans generated: (i) The escape rate,  (), defined as the partial energy content that remains confined within two parabolic trajectories delimiting the main maximum of an ideal Airy beam, given that the main maximum for the three input beams here considered basically overlaps with that of the ideal Airy beam.This quantity is computed according to Eq. ( 10) and the C(, 0) is determined by the above mentioned limiting trajectories.In order to produce a fair sampling in this single-event dynamics, the initial conditions are randomly chosen according to the weight |(, 0)| 2 (see Appendix A for additional technical details).Depending on the functional form of the aperture function, the escape rate will fall faster or slower, thus indicating either a faster or a slower deviation with respect to the ideal Airy beam.
(ii) The average position,  av (), determined at each -value as the direct average taking into account all trajectories in the ensemble.Again, like the escape rate, this quantity supplies local information about how fast the main peak disappears, and hence how fast the finite-energy beam looses the properties that characterize an ideal Airy beam.

Aperture functions
Specifically, to explore the feasibility of aperture functions with different functional forms and then to determine which one shows a better performance, we have considered the following aperture functions: In particular, for a comparison on the same footing conditions, we have considered   = 0 in the three cases, while  has been chosen in a way that the three functions have the same FWHM.In this latter regard, we have fixed the width of the Gaussian aperture function to be   = 26.5 m (0.5 in reduced units, in terms of  0 ), which renders a FWHM = 44.13m and the following width values for the two other aperture functions: where is the (positive) solution for sinc  arising from the Bhāskara I's sine approximation formula [22], for the equation sinc  = 0.5, with  = ( −  , )/  .It is worth mentioning that we have considered the same FWHM for the above three aperture functions instead of in their Fourier spectrum, because the truncation, in the end, is going to manifest along the spatial transverse section of the beam.Thus, regardless of the spectrum associated with the above aperture functions, if we set a given width for all of them, we will control the spatial limitation of the beam.However, if we choose instead to set the same FWHM for the corresponding Fourier spectra, we may end up with aperture functions with very different (non comparable) widths.
The three aperture functions are shown in Fig. 1a, all with the same FWHM, 44.13 m, which is nearly half the FWHM of the leading maximum of the ideal (infinite energy) Airy beam considered, 96.13 m.For a better comparison, the three functions are plotted without the constant prefactor, so that their maxima coincide with 1, while they asymptotically approach 0 at large .We notice that, although the three cases present nearly the same trend around   , the spatial extent of the sinc-type function is larger than in the other two cases.Nonetheless, if we inspect the transverse power spectrum, | ψ(  )| 2 , corresponding to the initial ansätze constructed with these aperture functions (see Fig. 1b), we notice that the three cases essentially span the same range of momenta, even though the profile of the respective spectra looks very different (again, for a better comparison, the three spectra have been normalized to the unity at their center).Given that the power spectrum of an ideal Airy beam is flat, with an infinite extension, we note that, while the Gaussian and Lorentzian aperture functions are going to introduce a gradually decreasing of the amplitude components involved in the Airy beam, the sinc-type The curves have been normalized to unity taking into account the maximum value of the input beam (black solid line).In the inset, the same curves but each normalized to unity with respect to their corresponding maxima.aperture function, on the contrary, is going to essentially preserve the same amplitude for the spectral range spanned.In the current case, as it can be seen in Fig. 1b, this spectral range is Δ  = 2/  ≈ 171.76 mm −1 , with all components, low and large spatial frequencies, essentially contributing the same.This is in sharp contrast with the other two aperture functions, where the low spatial frequencies will play a more prominent role than the large ones (more remarkable in the Lorentzian case than in the Gaussian one).

Numerical dynamical simulations
To check the feasibility of the aperture functions here proposed, we have conducted a series of numerical simulations due to the non-analiticity of the propagation process in the three cases.As it is mentioned above, the parameter values considered are analogous to those used in the original experiments [4], which enables a more direct comparison with such experiments.Accordingly, in Figs. 2, 3, and 4 we show both an overall view of the propagation process in terms of trajectories superimposed on the corresponding density plots of the intensity distribution (see part a in these figures), and a series of snapshots of the intensity distribution at given values of , which allow us to elucidate in a qualitative manner how fast the properties of ideal Airy beams are gradually lost (see part b in the figures).
In Fig. 2a we observe the propagation for a beam obtained with the Gaussian aperture function.As it is seen from the density plot, the beam behaves similarly to an ideal Airy beam for a relatively long distance, with a well-defined main leading maximum surviving up to nearly the largest -value considered ( = 30 cm).The snapshots displayed in Fig. 2b show clear evidence of such a behavior, particularly in the inset in this figure, where we have renormalized the three distribution for a better comparison.Nonetheless, unlike an ideal Airy beam, there is a progressive decrease in the intensity as well as a loss of the characteristics oscillations that form the tail of the beam.To understand this behavior, we resort to the trajectory description (see The curves have been normalized to unity taking into account the maximum value of the input beam (black solid line).In the inset, the same curves but each normalized to unity with respect to their corresponding maxima.Fig. 2a).Accordingly, we note that, as it was mentioned in [13], the limited extension of the tail releases the beam from the infinite energy influx coming from the back.Consequently, not only the beam propagates forward, but it can also get diffracted backwards, which is what the trajectories (white solid lines) clearly show.Actually, we can observe that, if we choose an initial condition at the position of the maximum of the main peak of the initial beam, the corresponding trajectory (black solid line) start deviating from a parabolic-type curve for distances beyond  = 15 cm, approximately.In other words, typical Airy beam properties are warranted for about 15 cm; afterwards, the beam starts getting distorted.
The results for the Lorentzian aperture function are shown in Fig. 3.In this case, we observe a faster degradation, at shorter distances (around 10 cm).In Fig. 3a the density plot reveals a large transfer of flux from the main peak to the rear secondary maxima, although with the seemingly net effect of getting stacked around  = 0 mm.If we inspect Fig. 3b, effectively we clearly notice an accumulation of intensity around  = 0 mm (see also inset), decreasing towards larger values of .From the trajectory description, we find an overall picture that resembles that of the Gaussian aperture, but with the subtle difference that now the trajectory starting at the position of the main maximum remains nearly around the same position, without following the foremost ensemble of trajectories that still show a parabolic-type shape.
The two previous cases can be understood in the light of the strong filtering produced by the aperture functions considered in each case, which remove dramatically the high frequencies from the spectrum (see Fig. 1a).In sharp contrast, the case of the sinc-type aperture function, represented in Fig. 4, shows a very different picture.First, except from some intensity bursts propagating backwards, the density plot shows that all maxima, both the main peak and the secondary maxima, exhibit the typical parabolic shape expected for an Airy beam.The results in Fig. 4b clearly indicate that, effectively, the beam behaves like an ideal Airy beam for longer distances.On the other hand, the behavior of the trajectories also reveals that there is a larger survival of the parabolic shape for longer distances.Unlike the two previous cases, note that here the trajectories show a back and forth displacement, which keeps confined along parabolic-like , and  = 30 cm (green solid line).The curves have been normalized to unity taking into account the maximum value of the beam at  = 10 cm (red solid line).In the inset, the same curves but each normalized to unity with respect to their corresponding maxima.channels for longer distances, until the backward energy influx reduces and they can escape (backwards).This behavior makes that a large portion of the foremost trajectories will keep showing their parabolic shape for rather long -distances.Indeed, note that the trajectory started at the position of maximum of the main peak is nearly parabolic for the full propagation.

Dynamics of the main maximum
Apart from a local description of the propagation process, the trajectories also render interesting information that cannot be directly obtained from the density plots.In particular, to further understand physically the performance of the three aperture functions, we can investigate the local dynamics associated with a piece of the beam without perturbing its propagation.As it is inferred from [23], the intensity content within two any initial conditions will remain constant along  provided we propagate those conditions along the corresponding trajectories.That is, if we consider two initial conditions delimiting a given intensity maximum (for instance, setting such initial conditions near the adjacent nodes, but not on top of them), we can determine at any subsequent -value how the corresponding intensity distributes and within which boundaries.
Accordingly, we have considered pairs of initial conditions between which we collect most of the intensity associated with the main peaks in the three cases (the ranges are given in Table 1 in Appendix A).Then, in order to understand how fast the corresponding beams deviate from an ideal Airy beam, we have computed two quantities.On the one hand, we have computed the escape rate,  (), which measures, along , the intensity confined within two boundaries given by two parabolic trajectories.It is clear that the better the preservation of the ideal Airy beam properties the larger the value of  () for larger and larger -values.This quantity is plotted in Fig. 5a.As it can be seen, the decay of  () is relatively fast for the Gaussian and the Lorentzian aperture functions, although the falloff for the latter is still slightly larger.While the energy (proportional to the number of trajectories) decays to about a half at  ≈ 18.8 cm for the Gaussian aperture function, in the Lorentizan case the same takes place at  ≈ 13.5 cm.This corresponds to the degradation of the intensity distribution profile observed in Figs.2a and 3a, respectively.On the contrary, for the sinc-type aperture function,  () remains nearly constant until  ≈ 19.8 cm and, extrapolating data, a decrease to about 50% is only expected for larger distances ( ≈ 33.5 cm), which is in correspondence with the behavior observed in Fig. 4a.
On the other hand, taking into account that those initial conditions have been randomly distributed according to the intensity distribution, their average,   , can provide an objective measure of the self-accelerated propagation of only this beam.Again, the better the preservation of the ideal behavior, the closer this average will be to a parabola.This quantity is represented in Fig. 5b (see solid lines), where we have also included the trajectories associated with the initial condition taken at the position of the maximum of the respective main peaks (dashed lines).To compare with, we have also included the trajectory of an ideal Airy beam with its initial position at the -value where the main peak of such a beam reaches its maximum (see green dotted line).As we can see, the trajectories corresponding to the maximum in the Gaussian and Lorentzian cases start deviating from the average value at short and medium distances, respectively, in agreement with what we already saw in Figs. 2 and 3.More specifically, such a deviation starts being remarkable at about 20 cm for the Gaussian aperture function, and about 10 cm for the Lorentzian one.These distances can be used as distinctive objective criterion to determine when the performance of the corresponding aperture function becomes poor.In therms of the average,   , this loss of efficiency is indicated by the transition from a nearly parabolic behavior, at the initial stages of the propagation, to a seemingly linear behavior at larger distances, which is produced by the gradual drain of trajectories (flux) backwards.
In the case of the sinc-type aperture function, again we observe a very different behavior, since both the average and the trajectory of the maximum basically propagate parallel one another.More interestingly, these two quantities overlap in a good approximation the trajectory of the maximum for an ideal Airy beam; only at large -values we start observing a discrepancy, which is a signature of the loss of the ideal behavior.

Discussion
We have introduced a trajectory-based methodology to investigate the effect of limiting the energy content of Airy beams by means of aperture functions with different shapes, and hence to elucidate which functional form better preserves the properties of an ideal Airy beam at the input.Thus, we have combined conventional propagation simulations, to determine the intensity distribution at difference distances from the input plane, with the above mentioned non-conventional approach based on an on-the-fly computation of the transverse energy flux, which provides us with information at a more local level.The transverse energy flux, described in terms of sets of trajectories, enables a description and analysis of the propagation process on a single-event basis, making more intuitive the properties of Airy beams (shape invariance and self-accelerated propagation) at the same time that it is also possible to perform some additional quantitative analyses of specific features.
From our analysis, we conclude that, among the three types of aperture functions investigated, sinc-type ones show a better performance, because they basically introduce a cut in the beam spectrum, without modifying the amplitude of its components, unlike other types of aperture functions, such as Gaussian or Lorentzian profiles.The trajectories show a very interesting picture in this regard.On the one hand, we note that the trajectories associated with the Gaussian or the Lorentzian aperture functions make evident a monotonous energy dissipation backwards.This continuous flow produces a rather fast lost of fringe visibility of the intensity distribution, as it is manifested by the corresponding escape rates, and hence the also fast loss of the typical properties of an Airy beam at relatively close distances from the input plane.In particular, it has been noticed that the self-accelerated displacement that should undergo the average position of the energy accumulated in the main maximum is gradually lost.
On the other, in the case of the sinc-type aperture function, the fact that all spatial spectral components play have nearly the same amplitude provokes that the presence of a wiggling displacement of the trajectories and hence a back-and-forth energy flow between adjacent intensity maxima.This oscillatory behavior enables the preservation of the typical Airy beam features for much longer distances, even though the gradual degradation observed in the pattern along them.Furthermore, unlike the two previous cases, here it has been shown that the average position of the energy accumulated in the main maximum follows very nicely the behavior expected for an ideal Airy beam.
In sum, the present analysis is aimed at providing new insights into the behavior of finite-energy Airy beams, and therefore also expected to contribute to the design and development of novel applications exploiting finite-energy Airy beams.Moreover, it is also worth noting that, by virtue of the tight analogy between the paraxial Helmholtz equation and the Schrödinger equation, the results shown and discussed here are also directly transferable to the generation and exploitation of Airy beams with nonzero mass particles, as it is the case of electron Airy beams, first produced in 2013 by Voloch-Bloch et al. [24].

A. Numerical and data analysis methods
The simulations conducted in this work to obtain the beam propagation and the corresponding trajectories, have been carried out following the same procedure formerly used in [13], based on the split-step for beam propagation and a spectral recast of the equation of motion for the trajectories [25].For the trajectories shown in Figs. 2, 3, and 4, we have considered an evenly-spaced distribution of 50 initial conditions, taken from a range that covers the first 10 intensity maxima of the input beam (the leading maximum and the nine subsequent secondary maxima).
Concerning the random distribution required to compute the quantities shown in Fig. 5, first we have fitted the main maximum of the input intensity distributions (see black line curves in Figs. 2, 3, and 4) to the biGaussian distribution, which captures fairly well the asymmetry of such a maximum.The values for the parameters  0 ,   , and   obtained from the corresponding fitting procedures for the three cases are given in Table 1.
Regarding the specific random choice, first, we select randomly from a [0, 1) uniform distribution a number.If this number is smaller than the ratio   /(  +   ), then we associate the initial condition to the left side of the biGaussian distribution (i.e.,  <  0 ); otherwise, the initial condition is associated with the right side.A number is then picked up randomly following a Gaussian distribution centered at  0 and with a width   , if the previous choice laid on the left side of the distribution, or   , if the choice led to the right side.It might happen that the number obtained from this procedure will lay on the opposite side where it is supposed to be; in such a case, the algorithm re-start the choice of the side and repeats the process.

Fig. 1 .
Fig. 1.(a) Spatial dependence of the three aperture functions considered here, all with the same FWHM, 44.13 m.The Gaussian aperture function is represented with black dotted line, the Lorentzian aperture function with red dashed line, and the sinc-type aperture function with blue solid line.For a better comparison, the three functions are normalized to the unity.(b) Transverse power spectrum, | ψ(  )| 2 , of the input convoluted wave functions arising after applying the convolution with the aperture functions displayed in part (a) [type and color lines are in correspondence with the respective figures in (a)].Again, for a better comparison, the three spectra have been normalized to the unity at their center.

Fig. 2 .
Fig. 2. (a) Flux trajectories (white solid lines) associated with an Airy beam convoluted with a Gaussian aperture function.The initial conditions cover homogeneous the range spanned by the first foremost 10 maxima (including the leading one).The thick black line denotes a trajectory with its initial condition at the same position where the main maximum lies.(b) Profiles of the intensity distribution at:  = 0 (black solid line),  = 10 cm (red solid line),  = 20 cm (blue solid line), and  = 30 cm (green solid line).The curves have been normalized to unity taking into account the maximum value of the input beam (black solid line).In the inset, the same curves but each normalized to unity with respect to their corresponding maxima.

Fig. 3 .
Fig. 3. (a) Flux trajectories (white solid lines) associated with an Airy beam convoluted with a Lorentzian aperture function.The initial conditions cover homogeneous the range spanned by the first foremost 10 maxima (including the leading one).The thick black line denotes a trajectory with its initial condition at the same position where the main maximum lies.(b) Profiles of the intensity distribution at:  = 0 (black solid line),  = 10 cm (red solid line),  = 20 cm (blue solid line), and  = 30 cm (green solid line).The curves have been normalized to unity taking into account the maximum value of the input beam (black solid line).In the inset, the same curves but each normalized to unity with respect to their corresponding maxima.

Fig. 4 .
Fig. 4. (a) Flux trajectories (white solid lines) associated with an Airy beam convoluted with a sinc-type aperture function.The initial conditions cover homogeneous the range spanned by the first foremost 10 maxima (including the leading one).The thick black line denotes a trajectory with its initial condition at the same position where the main maximum lies.(b) Profiles of the intensity distribution at:  = 0 (black solid line),  = 10 cm (red solid line),  = 20 cm (blue solid line), and  = 30 cm (green solid line).The curves have been normalized to unity taking into account the maximum value of the beam at  = 10 cm (red solid line).In the inset, the same curves but each normalized to unity with respect to their corresponding maxima.

Fig. 5 .
Fig. 5. (a) Escape rate,  (), as a function of  for the three aperture functions here considered.(b) Average position,  av (), over a sampling of 1,000 flux trajectories randomly distributed, at  = 0, within the area covered by the leading maximum for the three aperture functions.The distribution of initial conditions have been taken according to a biGaussian distribution (see Appendix A).To compare with, the -dependence of the leading maximum of an ideal (infinite energy) Airy beam is also shown, represented with green dashed line.

Table 1 .
Parameters used in the fitting of the main maximum to a biGaussian distribution.aperturefunction   (mm)  1 (mm)  2 (mm) i.c.range (mm)