Post-acceleration of electron bunches from laser-irradiated nanoclusters

In this paper the energy gain of attosecond electron bunches emitted during the interaction of intense, few-cycle linearly polarized lasers with nanoscale spherical clusters is determined. In this case electron bunches are emitted from the rear side of the cluster and are then further accelerated while co-propagating with the laser. A previous study has shown how this two-stage process readily occurs for clusters whose radii lie between the relativistic skin depth, δ r = γ 1/2 c/ω p , and the laser spot size σ L (Di Lucchio & Gibbon, Phys. Rev. STAB 18, 2015). An analytical model for focused light waves interacting with compact, overdense electron bunches in vacuum is derived heuristically from world-line equations of motion of an electron. The functional integral approach is followed under the mathematical point of view of integration with respect to a stochastic variable. The resulting picture of the laser wave crossing the electron’s trajectory leads to a finite energy gain of the electron in light–matter interaction in vacuum. The analytical theory is compared with three-dimensional PIC simulations from which trajectories of the electron bunches can be extracted. The effective increase in bunch energy is determined under realistic conditions both for the peak (mode) and the cutoff energy of the emitted bunch, in order to make quantitative comparisons with theory and the experimental findings of Cardenas et al , Nature Sci. Reports 9 (2019).


Introduction: direct laser acceleration of electrons in free space
The physics of laser acceleration of single electrons in vacuum has been widely studied, and despite the Lawson-Woodward theorem seemingly excluding useful electron energy gain due to phase slippage between the electron and the longitudinal electric field in the linear, non-relativistic regime, it is now well established that significant acceleration is possible due to nonlinear processes at highly relativistic intensities [1][2][3][4][5], where the dimensionless transverse momentum of the electron is bigger than unity, namely p ⊥ /mc > 1. One fundamental process, nonlinear ponderomotive scattering, assumes a collinear propagation geometry and has been studied making use of the paraxial approximation for the wave propagation, with laser pulse durations of 20 fs and above. However, when laser pulses are shorter in time than 20 fs, the rapid carrier-wave amplitude variations must also be taken into account [2].
Acceleration of sub-micron size electron bunches to relativistic energies in vacuum by means of radially polarized electromagnetic pulses has been demonstrated experimentally and computationally [6,7]. These processes can be efficient in an unbounded free space acceleration process; however, they require the electron bunch to be confined along the laser propagation axis. This mechanism foresees no angular emission and the laser duration is not specified. The axial confinement of the laser pulse is also the basis for the simplest schemes so far proposed for achieving a linear energy gain for a relativistic electron bunch in vacuum. As an example, a current state-of-the-art facility can provide a linac as a mechanism for producing an electron bunch [8] holding a non-zero initial energy , typically of the order of few keV. The bunch can experience a very high electromagnetic field gradient while co-propagating phase-locked with an ultrashort laser pulse [9], thus acquiring relativistic energies through a linear energy gain process. Here, much effort was also dedicated to isolate the second-order quadratic energy gain typical of the ponderomotive scattering. Linear behaviour of the energy gain in an ultrashort, ultrarelativistic process can only be considered in a qualitative way due to the lack of a simple relationship between conservation of energy during photon-electron interactions and the general theory of ultrashort lasers acceleration in vacuum. Charged electron bunches ejected from an obliquely irradiated solid slab have been modelised by means of electromagnetic work done in terms of emitted and reflected field in [10]. In this study, short laser pulses are used, and the variation of momentum appears to be more complicated than it was described in [11], due to the decomposition of the laser wave. There holds still an analogy between [10] and [11] in the basic mechanism in which the electrons acquire energy first from the charge separation in the plasma following the irradiation, and then by traveling together with the laser in vacuum. In addition to that, the bunches traveling together with the reflected wave are phase-locked like the ones emitted from a cluster. In spite of the promising similarity between the two pictures, differences emerge after a second deep analysis. The mechanism proposed in [11] includes an ultraintense few-cycle laser pulse and a spherical geometry which allows to select exact emission directions according to scattering laws, a very important improvement in precision and also a simplification in the evolution of the bunches, which appear to be packed together into a nanometric size dense plasma region very soon. The novelty in [11] also relies in the fact that the particle distribution inside the bunches is narrowed and overdense. The latter result also sheds new light on gaussian nonlinear optics on the attosecond time scale. It is one of the recent advances in this field.
Few-cycle laser technology has now reached the point where relativistic focal intensities can be achieved, allowing for new accelerator target designs which take advantage of small time and space scales. Laser-plasma based acceleration schemes such as laser wakefield acceleration [12] are already well established for longer pulses and gas targets, but high-current injection is still a challenge. Energetic ion acceleration from laser-illuminated droplets [13] or other mass-limited targets [14] has also received interest because of the potential for generating beams with narrow energy bandwidth.
In our previous work on laser-irradiated nanodroplets [11] we demonstrated that there is the optimal choice of droplet size depending on the pulse intensity. Using 2-dimensional simulations it was possible to predict emission angles for electron bunches generated by solid-density, nanometer-sized droplets after irradiation by a few-cycle, ultraintense, linearly polarized, Gaussian laser pulse. Special attention was paid to deviations in the electron emission pattern in the relativistic regime with respect to the classical Mie theory's predictions in the attosecond time domain, when the electron acceleration is naturally sub-cycle. It was shown that emission angle depends strongly on the laser intensity, roughly following the trend expected for ponderomotive scattering of single electrons [15], but with much higher final energies. The reason for this was that after emission from the cluster surface, the bunches propagate in the laser polarization plane, phase-locked with the laser electromagnetic field.
In the present paper a fully three-dimensional analysis is presented in order to reconstruct this second 'postacceleration' phase of the bunches' interaction with the laser pulse. The paper is organized as follows: first the interaction dynamics at the surface is briefly examined, with a special focus on the electric field enhancement at the surface. The first step, namely the initial energy acquired by the electron cloud formed immediately after the passage of the laser pulse through the droplet, is a prerequisite for the onset of the subsequent ultra-relativistic scattering process. Then, the energy evolution of the bunch is numerically tracked using PIC simulations in order to determine the peak and cutoff energy values of a typical electron bunch with a finite energy spread. In section III a three-dimensional analytical model is formulated and compared with the particle-tracking analyses of corresponding PIC simulations performed in the context of a recent experiment [16].

Three-dimensional PIC simulation results
Following our initial work examining the electron bunch emission characteristics [11] and subsequent numerical studies accompanying the experimental campaign reported in [16], we have carried out additional three-dimensional simulations to investigate the long-time evolution of the post-acceleration of attosecond electron bunches emitted from sub-micron scale spherical targets due to irradiation by a few-cycle laser pulse.
It has already been shown in [11] that significant departures from Mie theory are found for electron bunch emission for droplets whose radii satisfy the condition δ r < R < 10δ r , where δ r = γ 1/2 c/ω p is the plasma relativistic skin depth. To examine the formation and propagation dynamics of electron bunches created during such interactions, we have made use of the electromagnetic, 3-dimensional EPOCH code [17] on the JURECA cluster [18] at the Jülich Supercomputing Centre. The droplet was modeled in simulations by a solid 200nm diameter hemisphere with density n 0 = 100n c placed on the top of a cone, in order to mimic the configuration used in the experiment [16]. This was irradiated by a 2-cycle, 4.5 fs laser pulse with 1.22 μm × 1.22 μm focal spot size, linearly polarized in the plane perpendicular to the cone axis, so that a quasi-2D bunch emission geometry is achieved. The laser intensity I was varied between 2 × 10 19 W cm −2 and 10 21 W cm −2 . The sphere size is chosen to match the empirical criterion above to avoid complete relativistic transparency and efficiently produce the overdense attosecond bunches with experimentally feasible target dimensions.
Despite the short interaction time, the target still undergoes some expansion, apart from likely being subject to some laser prepulse which will also serve to pre-ionize the target. Field and impact ionization may well alter the details of the density profile, but this is beyond the scope of the present study. An example of the effective plasma profile seen by the main pulse is shown in figure 1 for an intensity I = 5 × 10 19 . Here the density rapidly falls towards its critical value within the focal dimension. This is the actual cloud with which the laser pulse will interact for the post-acceleration mechanism, with some oscillations due to the electric field value-the stronger Figure 1. Plasma profile expansion soon after the interaction of the needle with the laser.The cluster size is 200 nm, the laser intensity is I = 5 × 10 19 W cm −2 and its duration is 4.5 fs. The plasma hydrodynamic expansion due to prepulse effects has been reproduced by simulating an initial temperature of about 5KeV. The blue curve is an exponential fit with scalelength L = 0.069 μm Figure 2. Enhanced vector electric field (red arrows) surrounding the three-dimensional sphere in the emission plane: (a) after at t = 12.5 fs and (b) at t = 18.75 fs after the start of the laser-droplet interaction. The color scale on the right represents the electron density, n n log e c ( ). The inside part of the droplet remains almost unperturbed.
the electric field, the faster the depletion. However, the electron mean cloud density lies between n = n crit inside the focus for I = 2 × 10 19 and n ≈ 0.5n crit for I = (2 − 8) × 10 19 .
The acceleration follows a two-stage process: attosecond electron bunches are first emitted from the target surface and are injected into the electric field null phase of the laser pulse with relativistic energies [16]. There then follows a second post-acceleration step in which the electron bunch gains energy while co-propagating with the transmitted part of the laser pulse for at least one Rayleigh length beyond the original target location. This is a very important aspect which distinguishes the present study from the conventional 'vacuum focus' in which phase slippage tends to cause any energy gain to be cancelled out due to a transition between different half-cycles of the laser. Later we further analyze the linear ponderomotive energy gain in more detail, under the special circumstances where the interaction of the bunches with the laser takes place immediately after the laser-droplet illumination.
To analyze both the bunches' emission and their subsequent post-acceleration in the electromagnetic field, a reduced set of particles from the 3D PIC simulation (around 3% of the total) were tracked. Emission occurs primarily in the laser polarization plane to the left and right of the propagation axis , as it is visible in figure 2. The asymmetry between left or right here does not compromise the analysis since the same statistical distribution can be applied to both sides. It turns out that forward trajectories with a very reduced set of particles, that is, below 10, are not representative of the true phase space evolution of the bunch since space charge effects must be properly taken into account.
Before drawing an analogy between a first-step accelerator and the near-field effects in the vicinity of the droplet-where soon after the laser is incident an enhanced electric field sustains the initial phase of the acceleration (see figure 2)-we split the process in two steps. The electrons are initially stationary on the surface, so it is assumed that the enhanced field around the droplet soon after the particles' ejection enables the attosecond bunch thus formed to acquire an initially relativistic energy within a laser subcycle. This assumption will be examined more closely shortly. A visual representation of the enhanced electric field structure in which the electron cloud is bound to the plasmonic initial energy is shown in figure 3. A dense electron cloud will be formed after the laser-plasma interaction when the longitudinal field island vanish. However, initial energy for cloud electrons and for bunches' electrons are the same at this point of space and time, where transverse field is about to provide confinement effect for bunches.
There are two consequences of this early interaction phase: first, the bunches will acquire a relativistic injection energy which is a consequence of high electric field O(TV/m) generated at the droplet surface. Second, the bunches are overdense, and very much concentrated in space-a few hundreds of nm-and with very little divergence with respect to the laser beam. These properties of the bunches have practical advantages, but makes it challenging to analyze the bunch dynamics in detail via PIC simulation.
In this case it is important to have sufficient statistics to resolve a range of relativistic trajectories and collect enough single particle interactions to reasonably represent the final energy spread of the bunch. The peak energy value is statistically calculated as the mean of the Gaussian distribution outside the filtered electron bunch. Any Gaussian distribution has a peak, an aperture and two tails. We are mainly interested in the peak value, where most of the electrons are, and in the extreme tail value, namely the maximum reachable energy gain which is allowed for the electron bunch. The relationship between optimal energy gain and phase shift of the bunches due to electron trajectory position crossing with laser wave, can also be examined analytically in the quasiparticle approximation.
Trajectories collected with a time-resolution of 250 attoseconds are summarized in figure 4. Filtering is performed dynamically with the help of a simple artificial intelligence selection technique [19], as a postsimulation analysis on the bunches selected by the PIC code. This happens because the dynamics is non-linear and therefore electron trajectories in the PIC selected bunch are allowed to cross each other. Previous nonlinear studies have so far dealt with longer pulse schemes, for example high-intensity VLA schemes [20] or plasma mirror injectors [21]. In these cases the mean direction and energy emerge progressively as part of a temporally limited process whose duration is strictly connected to the pulse duration itself. This is not applicable in our context.

Analytical model of the electron bunches' energy gain
One recent example of a linear energy gain in a few-cycle laser accelerated bunches when the laser focus is small or of the order of the laser wavelength is to be found in [9]. However, when considering both theoretical [11] and experimental results [16], bunch dynamics treatments at nonrelativistic laser intensities are insufficient to explain our findings. No direct evidence for longitudinal acceleration of the bunch by the laser in vacuum could be identified here. Instead the oblique acceleration of the bunches happens far way from the focus and at two symmetric orientations relative to ponderomotive scattering.
Moreover, most of the previous studies have devoted their efforts into distinguishing the paraxial approximation, which is valid for long pulses, from non-paraxial treatments which allow for a shorter pulse duration; but little attention has been devoted to the case of a few-cycle laser pulse to date.
Here we will first assume that our 4.5 fs laser pulse is in the short-pulse regime according to the subdivision made by Hartemann et al [2]. The multi-electron dynamics are studied by investigating the very high-energy gain observed in three-dimensional PIC simulations (section 2).

Basic equations
According to [2] the maximum energy in a finite duration laser wave can be simply inferred from the temporal evolution of the γ of the electrons, which in turn corresponds to the third member of the relativistic energymomentum transfer equations. We are interested in analyzing the process immediately after the bunch extraction from the cluster, that is, the laser-vacuum linear acceleration in the relativistic regime. Therefore we recall the already established analysis in [5] and take this model as a starting point. The electron 4-velocity u μ , 4-momentum p μ and the 4-vector potential A μ can be defined as a function of the proper time τ of the electron particle along its world line as follows:   1 In equations (1) the laser phase f is expressed as a function of the 4-vector wave number k μ and of the axis coordinates 4-vector x μ . The scalar ω 0 is the wave frequency. The vector (1,0,0,1) indicates that the laser travels along z direction. The light-cone invariant k k ) and the transverse momentum vector u ⊥ (τ) are defined respectively by the following two equations: It is straightforward to see how u ⊥ (τ) is a function of proper time τ through the wave phase f. Given that g = + + u u 1 z 2 2 2 the following general result for the evolution equation of γ in the context of focused plane waves holds-for details see [5] and notes in appendix A: In the Lorentz gauge and for a solenoidal vector potential, ∇ · A = 0, the vector potential and the electric field transformed as in equation (3b) and equation (3c). Equation (3a) can be regarded as the relativistic ponderomotive scattering version of the Lawson-Woodward theorem. In principle, the electron cannot acquire energy in the plane wave, which means that, when the electron is accelerated, it also radiates. In the final state, the laser wave attenuates, and a permanent destructive interference between the laser wave and the electron wave holds. The electron radiates waves that hold a destructive interference with the laser pulse. Even for two-cycle pulses the laser electric field can still be mathematically represented as the product of a complex envelope function a(t − z/c) and an optical carrier wave we j t z c 0 ( ) . For our specific case, a slowly varying envelope would fit the phase invariant picture when put in the form

Functional integral
The following treatment departs somewhat from Hartemann's equations of motion. In fact what we want to get here is an explicit expression for the vector potential A μ (f). The functional derivative rule can be used to state the vector potential variation in the most general way, as a function of f Here , f(τ) is the laser envelope phase, which can be expressed as a function of the proper time. Due to the non-invariance of the proper time itself, the functional representation will be proven to be useful because it allows simplifying the picture by considering only the phase variation in the vector potential. As long as the gauge conditions are maintained, the integral must be considered throughout the path along the effective phase displacement between the quasi-particle representing the electron bunch and the carrier envelope in the ultraintense laser, in order to express the vector potential variation. The high density of the bunch and the relativistic energy of its electrons in turn induce a strong correlation within the overdense bunch structure, whose nature is new and linked to the motion of the bunches in the laser. The second integral of the identity for the vector potential can be integrated by parts.
has been used, apart from the general definition of the functional derivative. The envelope function is now expressed as with g(f) as in equation (4).
The extra 2ω 0 factor in equation (7) can be absorbed in a dimensionless picture where c = 1 and κ μ = ω 0 The function g(u) inside the integral can be developed using the trigonometric identity, cos cos sin sin ( ) , and thus yields The shift of the phase between the electron trajectory and the laser field crossing it is already known to be a consequence of the switching on of the laser electric pulse also in the simplified ponderomotive picture [22] and may be mathematically expressed in three different ways: The expression in equation (9c) will be used in the following integration to represent the effective phase shift in energy gain final calculation. The kind of variations considered in 9a, 9b will be used in appendix A to figure out how the functional can be applied in this specific case. All of these variations do not take into account directly proper time, because it is not an invariant. Our aim is not to model time in all systems of reference, since the phase is enough to get a description of the process. We reconstruct the path along which the γ value actually increases up to the maximum and cutoff energies recognized as the main features of the electron bunch by solving the integral in equation (4) explicitly, instead of calculating a parametric solution as it is done in [2], in order to determine the effective variation of the γ during a finite time in which the electron remains in the laser wave. To do this we take our lead from the particle tracking results for the energy gain in figure 4. Equation (7) is integrated by parts thus: Integrating the 2nd term a further two times leads to a recurrence relation which can be solved to yield the final result: an analytical solution explicitly dependent on the wave phase shift with respect to the electron trajectory: ( ) The sin 2 envelope is perturbed by the oscillatory terms, thus generating in principle a different numerical result for the gain if the phase shift was not zero-namely, if the bunches were not phase-locked with the laser [23]. The function in equation (10) is displayed in figure 5 for the realistic values in the case of a two-cycle laser pulse (k = 0.27). In order to state the mathematical difference between the x integration variable in figure 5 and the phase variable (non necessarily a function of proper time) used in the treatment by Hartemann et al [5], we refer to the definition of functional as a number associated to a integral defined over a vector set of functions as it is presented in [24].
In the integral expression obtained in equation (10) a physical discontinuity must be avoided since when the general shift is zero, the integral is also zero. This means that when there is no phase difference between the electron trajectory and the laser field, then there is no gain (no assumption is made on the initial laser profile by setting the lower extreme of the integral equal to zero). However, this is not supposed to really happen for any nonzero shift that is within the allowed interval and the periodicity of the sine and cosine factors, namely whenever the parenthesis of the factors yield a null result for (Δf + π). In order to remove the discontinuity, we directly consider the square of the integral-i.e. the true acceleration yielded apart from the corrective terms, in the critical phases that are an integer multiple of π. A straightforward purely mathematical calculation, including a macLaurin development for the generic function coming from the expression for the square of the integral I defined as in equation (7) depending on the phase when f = 0, yields the result where E| f+π represents the electric field convolution expression whose cycle-limited Fourier transform gives the product of the two signal envelopes f p + The above calculation provides a linear fit, and therefore, a linear dependence in E 0 , when the actual value of the phase shift is given. This is especially relevant for the points in which exact gain is reached, namely, the allowed multiples of Δf = π in the 0, 3π phase shift range. In particular, the dimensionless measure for the product ωΔτ can be set to Δf = 2π for the integration on the cone. Keeping the argument of the sin function within the one laser period, one can easily observe how the choice of multiples of π/2 for Δf yields a null contribution in the cosine factor within the first term of equation (10), and therefore only integer multiples of π should be considered. If Δf = π the factor k -or the third term in equation (8)makes the first term in equation (10) on the cone surface reduce to which means that the square integral yields a dimensionless gain equal to E 4 0 . It is found that for Δf = 2π the first consistent gamma increase takes place and the calculation of energy leads to the statistical peak value of the Gaussian bunch corresponds to a dimensionless angle p 2 and a gain equal to

From initial energy to peak value
In order to confirm the results coming from the integral resolution in equation (10), the energetic pseudoparticles have been tracked in the code since when the sub-cycle dynamics appears. Since the electron cannot either be emitted or be injected with a non-zero initial energy , an approximation for the initial energy γ 0 in equation (3a) must be found.
Previous plasmonics experiments with lasers and nanotargets have been carried out for much lower laser intensities and in principle different electron and light emission models have been considered [25,26]. An analytical treatment based on exact solutions of Dirac's equation for monochromatic linear waves at relativistic intensity in plasmas has been developed in [27], where the main difference is that the Drude free electron model of a plasma medium only extends to the interface of the droplet with the vacuum.
In [11] the dielectric permittivity ( ) is assumed for the droplet, where ω p is the plasma frequency and ω is the laser frequency. However, in contrast to the assumptions used in [27], the considered plasma is not the one of the medium through which the laser is traveling initially. It must be noticed that, in any case, the considered electron density n e would be more than critical, even applying the empirical criterion in [11] for droplets, because, under this condition, plasma transparency reduces the initial electron density (100n c ) by a factor equal to the γ of the laser, which is still not enough to get an underdense medium (see density log scale in figure 2 as an example). According to the free electron model, which is widely used in nanoplasmonics and nanophotonics , if we can consider the outcoming electron cloud as a plasmonic surface with respect to the positively charged solid cluster left behind, we always get a density less than critical or around critical. This hypothesis has been verified by means of experimental measurements in a recent work showed in [28], where PIC calculations were carried on in the context of an intense laser irradiation for ion acceleration from submicron droplets. This scheme reproduces the conditions for the Debye length to be of the same order of the droplet diameter, as in [13]. In addition to that, it is shown how the polarization choice associated to this mathematical treatment, namely, the usage of 2ω 0 irradiation source, is responsible for the formation of a hot electron cloud with reduced temperature, which is propaedeutic to the development of a subsequent ion cloud generated by operating a charge separation mechanism. The 2ω 0 factor yielded by equation (7) constitutes the link between our model and the experimental 2ω polarization that can be chosen for an ion acceleration that happens in two steps. This is analogous to the process described in [13].
The plasmonic energy peak in the vicinity of a spherical nanoobject has been demonstrated to be tunable, according to a dimensionless generalized plasmonicity index [29]. The evolution of the resonance energy shows a square root trend as a function of electron density. This is symmetrical with respect to the change of plasma frequency in the droplet as assumed in [11], where w w g = p r reflected the relativistic transparency effect. The electron density shows a quadratic dependence on the plasma frequency, therefore in this context, if the plasma frequency is not to be changed, the effective plasma density will be obtained from a multiplication by the factor γ. The time dependence of the γ factor is not yet discussed in these previous studies, due to the lack of dedicated calculations and/or measurements with the existing detection methods.
In [29], the plasmonic behaviour just depends on sphere size and material. It is important to remark that plasmonicity is not field enhancement; it is a way of parametrizing resonance that constitutes a good theoretical approach both in quantum theory and in classical electrodynamics approach, and it constitutes a second-order correction to the field enhancement calculations as it is effectively done by means of electromagnetic equations and consistent numerical modeling. Therefore, our Mie calculations for the near-field enhancement that have been discussed in [11] do not need a specific correction by means of the quantum density of states, neither in the Mie code not in the initial setting of the PIC code. Due to the relativistic intensity of the laser, the initial configuration vanishes rapidly and a specific determination of the density of states as in [30] would be washed away by the time evolution of the γ factor in the density expression. In the end, the only possible approximation for the initial mean energy value from the artificial intelligence analysis over a dense sample in the cloud can only be determined with a general formula for the plasmonic energy. The so far accomplished PIC calculations in our specific work do not allow to identify the plasmonic energy of the transition or the transition damping, two factors which actually change the effective plasmonicity. The shift in the energy as a function of the electron density of the cloud is coming out as a collective process in the initial determination of the average energy. The prevalent energy in the statistical population of the cloud is found, and its fluctuations are only formally connected to the radiation damping changes. The frequency of the plasmonic oscillation is an artificial process here, since it comes from further interaction of the cloud with the laser pulse. The cloud itself represents a medium with its proper permittivity and its characteristic resonance, and undergoes functional treatment according to its strong cross-correlation effects (see appendix B). This mechanism is new outside metals and other well-known plasmonic materials, and it is for the first time identified here at relativistic light intensities. When the laser energy increases, the plasmonic cloud is more and more depleted by the extraction process and its density is lowered down to n = 0.25n cr .
Due to the complexity of this process, only the results for I = 5e19 W cm −2 will be shown, although a consistent table could be built for all intensities.
The initial value of the electron energy in the cloud as presented in figure 4(c) in [16] is reconstructed by means of filtered tracking of the bunches along the scattering directions. The results for different intensities are presented in table 1.
Higher intensities before the upper limit of I = 10 23 W cm −2 have been also investigated and included as upper graph curves in figure 4. Their behaviour, although similar to the lower energy curves, shows a saturation process. The upper laser intensity limit specified above is found when the droplet material becomes optically transparent.
The saturation process leads to a final peak value for the bunches' electron energy. This value can be compared with the results from equation (10) ,provided that the right phase difference values are chosen. In principle equation (10) can cover all transient energy values which are obtainable in the laser wave crossing the Table 2. table of electron bunches' peak energies, as calculated theoretically in equation (11) and as inferred from tracking (see figure 4).

Laser Irradiance
Predicted electrons' trajectory. In practice, when substituting Δf with physical values which can be related to the cycles of the laser, we obtain the favoured phase values at which optimal energy gain takes place. A straightforward calculation for these cases shows that all terms in equation . The second factor shift can be neglected if we make abstraction out of asymmetry effects. At the same time, the first term of the integral expression can be expressed in terms of the same initial electric field, only calculated at a point in which the phase is seen as shifted by π. Therefore the total maximum available shift will be the one of an entire laser pulse plus the shift, or 3π.

Interpretation of tracking results
Before presenting the full comparison of three-dimensional tracking with fits from the mathematical treatment, we shall mention how the presence of a finite region of acceleration highly simplifies this process both in the theory and in the numerical calculation. A close look at the analytical function inside the integral immediately shows such an advantage. It is important to notice how, on the cone, the argument of the laser envelope function g(u) is smaller than π, therefore the quasi-particle only sees the first laser cycle on the light cone surface. When Δf = 4π, corresponding to a dimensionless shift of π then the sin term goes to zero, that is, no gain after the maximum phase shift. This dumping is possible because of the mathematical form of the laser envelope function, which in turn implies a finite duration in the physical picture of the laser. Therefore, the two-cycle interaction of length 2π reduces to a one-cycle interaction in the argument of the g(u) function. This is essentially a natural consequence of the two-cycle duration of the laser pulse. Taking into account the value of n = 2k, with k defined as in equation (8), and its dependence on ω and Δτ = 4.5fs, the final result can be inferred for the different laser intensities provided that the initial value of the gamma is found by means of the initial ponderomotive condition g = + a 1 0 0 2 ( )and the subsequent application of equation (4) where the square of the integral is substituted by its MacLaurin development as in equation (11). The analytical results hold physical significance for a shift in phase on the cone equal to π as it is determined by the difference with the initial point for Δf and the approximated calculation in Δf + π in the development in equation (11). The selected results for different laser intensities are compared with tracking outcomes in table 2). This result has been verified by means of tracking of the bunch at longer times, namely up to 165 fs and beyond, and it is shown in figure 9. At this time the energy gain saturates and the Gaussian form of the energy spectrum of the single bunch is represented in figure 10 for I = 5 × 10 19 W cm −2 , while still keeping the global contribution of background particles coming from the cluster as a strict comparison. The selection of energies within a single bunch is evident. The contribution to the highest possible apparent energy comes from the onaxis RPA electrons, which are not bunches. It is important then to notice how the bunches'filtering allows  figure 1(d) in [16]. The peak bunch energy is the result of dedicated tracking as it is represented in figure 4. On the right side the maximum angle in equation (14) corresponds to the denser bunches' regions, whereas the small dense emission close to the x-axis follows the lowest direction in equation (14) and corresponds both to the blue triangle shown on the left picture of figure 2 in [16] and to the minimum angle in equation (14). perfect agreement with theory and with experimental findings. These results hold for different focus sizes, for example, from just below to double its original value w 2 = 1, 22 × 2 = 2, 44 μm. In both cases the initial cloud and the short-range scattering evolution exhibit the same behaviour. The long range tracking statistical and dynamical analysis shows a transition in the statistically determined overall mean energy. As long as time scale approaches the 3π gain phase, some electrons start moving toward the cutoff energy and this process shifts the apparent peak energy. However, intermediate values between peak and cutoff energy do not represent special optimal gain themselves, but are instead an effect of spreading of the acceleration process between different phases as seen by the electrons. A single monoenergetic electron bunch whose central direction is deviated according to optical transparency effects [11] is nevertheless formed and shows a continuum distribution (see figures 9, 10).
To check the intensity dependence of this process, we perform several more simulations at higher intensities. The electric field analysis for droplets irradiated by few-cycle lasers of intensity I = 2.4 × 10 20 W cm −2 and 1 × 10 21 W cm −2 respectively are presented in figure 7. In the same picture it is also visible how the number of extracted attosecond electron bunches also increases with increasing laser intensity, leading to more overdense slices in between the zones of the zero electric field. This process happens at the expense of the electron cloud surrounding the droplet.
The maximum energy tracking has been carried out in an analogous way with respect to the mean energy tracking, allowing for efficient comparison at long simulation times. Substituting Δf = 3π in equation (10) the cutoff energy value of the Gaussian distribution comes out from direct analytical calculation.The corresponding values for the absolute amplitude of the integral value can also be visualized in figure 5. The critical values of Δf = 2π, Δf = 3π correspond numerically to the foreseen critical phases for scattering in [2] and long range additional energy-gaining interactions are to be excluded in advance. Therefore only two laser cycles are the optimal choice for a positive gain. Negative cycles, and corresponding negative energy solutions on the light cone, are not allowed in this treatment (see methodological notes).

Understanding of the phase effects in terms of the spectrum of the bunches
The limitation in the final phase shift value and the specific periodicity of positive energy gain are consistent with the existence of a finite interaction region, although the bunches can travel for several Rayleigh lengths before being separated from the laser pulse. The phase is an invariant and therefore only the phase spectrum that the particles see before the slippage is important. On the other hand, the phase is related to the number of entire cycles which are in principle needed to ensure the energy gain , and that cannot exceed 3π. The explicit calculations for the corresponding space are addressed in the following. The space covered for Δψ = 2π, 3π respectively corresponds to the saturation length for the peak and the cutoff value of the energy for the bunch. This has been verified against the bunches'tracking in the simulation for I = 5 × 10 19 W cm −2 . In particular, it turns out from formula in equation (13) that at this laser intensity the distance at which the peak is formed is 6.27 μm, and the distance for the cutoff is 9.41 μm. Both of them have been observed in the corresponding tracking curves in figures 4 and 8. The oscillations in the mean and maximum energies that start to appear at this distance points, are a reflection of how the tracking has been built. In practice, the logical functions that 'cut' the bunch region in the post-simulation analysis only identify as a starting point the mass and charge density. However, it is possible to set up independent calculations for identifying the bunch charge statistical distribution on a virtual detector that catches the bunches at different time instants. According to the results of both three-dimensional simulations and experimental findings in [16], the power distribution of the spectrum can be taken as Gaussian from early times of the post-acceleration. By calculating separately the basic frequency and the inverse of the standard deviation of the spectrum, one can separate the coherent region of the signal-which takes place before the saturation lengths are reached; and the incoherent region-showing up after the saturation length. For a bunch of a sub-micron length these regions are given by w st 1  (coherent) and w s - [31,32]. The value of σ t is given by the standard deviation of the inverse of the frequency interval , or temporal width of the bunch.
The incoherent region of the spectrum is normally associated with no correct information of the bunch, unless the whole analysis is changed, and shifted to an imaginary spectrum. A detailed Fourier analysis of the bunches'spectrum goes beyond the scope of this study (see appendix C).

Discussion of methods
As long as of both longitudinal and transverse laser pulse components are present, the electrodynamic currents generated by the laser must play a role in the energy increase. This is also in good agreement with the geometrical momentum decomposition in the two-dimensional description in [11].
Due to the lack of a specific theoretical treatment which states an absolute tuning of the few cycle laser interaction with matter at nanoscale, three-dimensional particle-in-cell simulations must be used to determine threshold for bunch production, spacing and the effects of ultra-tight focusing or ultrashort pulse duration. In [9] it is well established how the shortness of the laser pulse effectively changes the energy scaling. However, in the original work by Wong et al [9] the initial conditions appear to be different and a complete matching with the laser theory is not carried on. Here, the results of three-dimensional simulations clearly indicate the limitation on any approximation neglecting the charged particle trajectory in the laser field. Therefore the treatment in [9] represents a statistical approach which shall not be used as a substitute for theoretical analysis. The realistic description based on the effective duration of the laser-bunch interaction can be only provided by threedimensional full analysis.
The determination of bunches' cutoff energy on a statistical point of view-with post-processing and basic artificial intelligence techniques is represented by graph lines in figure 8. When starting with a initial value of γ = 2.1, as it is set in subset tracking in the three-dimensional EPOCH code, with the aim of following the initial condition for the formation of an after-expansion cloud as it is shown in figure 1, the time evolution maximum attainable value of γ can be easily calculated in the limit of a short plane wave (see equation (16)), from hence the cutoff energy value can be estimated. The cutoff bunches' energies corresponding to different applied laser intensities are summarized in table 3.

Discussion of emission angles
The angles of emission in tables 2, 3, may be explained with the maximum and minimum emission angles as a function of the initial and final γ and β of the electron as foreseen in [2]. The corresponding formulas are here recalled and solved by substituting the same initial conditions as for the subset tracking in EPOCH and for determining the maximum energy: For a resume of the cutoff energy angle dependence on laser intensity, see table 3. The formula for maximum energy gain in short plane waves from [5] reads and is also in agreement with figure 8 and table 3. The limit of plane waves sets the maximum energy which is achievable far from the focus, and that is consistent with a many Rayleigh length extended acceleration. The current treatment sets no overcoming of this limit, therefore respecting the already known limitations for the scattering process. The only way to increase further energy of the emitted bunches is to explore all the possible highest laser intensities. Once again for the determination of the cutoff energies, we have applied the plasmonic initial condition, which reads γ 0 = 1.1MeV. The highest intensities, until the upper limit, are represented in figure 8. The maximum energy from equation (10) and equation (16) reads  figure 8. Those are very high energy values and therefore especially long simulation time and tracking are required. The non-invariance of the proper acceleration time -which cannot be stated before-forces the tracking up to the computational limits; the matching at 1 × 10 21 W cm −2 is nevertheless reached within the computational capabilities of the resources. Therefore, intensities above 1 × 10 21 W cm −2 have not been investigated. However, the energy curves in [11,16] are nevertheless fully covered.

Angles of emission
The absence of second-order effects in energy gain, and consequently the linear dependence of the γ increment on the electric field, is due to the laser pulse's very short duration , which is compatible with the absence of significant double Doppler up-shift and quivering down-shift effects. Both radiating effects can occur when relativistic electron beams interact with lasers via multiphoton processes, under the classical view of Compton scattering or in Thomson scattering at some interaction angle [33]. However, the initial electron γ must be much greater that the laser a 0 , whereas here γ 0 is smaller than a 0 . Therefore no quivering effect must be considered and Hartemann's three dimensional treatment can be followed, for determining the angle at which the electron bunch should be traveling in order to satisfy energy conservation. For plane waves only, a strict correlation is drawn between the electron trajectory angle q =ârctan u u z (|( ) ( )|) and its energy gamma. This leads to the following set of mathematical relationships: Substituting the final mean value of γ = 8.3 in equation (17), one obtains Θ = 0, 32rad = 18°for I = 5 × 10 19 W cm −2 . This should be compared to the deflected angle in figure 3 in [11]. Substituting the final cutoff value of γ = 22.56 the same formula yields Θ = 0.2rad = 11.44°, which is to be compared with the results for later time tracking in figures 9, 10. Substituting γ = 35.32 results in Θ = 0.16rad = 9.21°. It is interesting to notice that the final correspondence between equation (17) and equation (14) reinforces the three-dimensional validity of the previous treatment. For the lowest intensity I = 2 × 10 19 W cm −2 the slightly different initial condition in γ 0 = 2.0 yields a foreseen angle of Θ = 0.24rad = 13.7°for the final mean energy-see table 2.
We can analyze the relationship of the two expressions for the angles 14 with the three-dimensional general case by making use of the linear development already used in 11 for Δf = 0, keeping the square away from the calculated integral and considering only the plain result: The we get the expression (16) in [34] We are aware that the factor π introduces an extra rotation around the plane and we wish that further work will be dedicated in order to get more insight into the theme of the connection with three-dimensional geometrical picture of this configuration in laser-plasma interaction.

Conclusions
We have analyzed both analytically and with three-dimensional PIC simulations the dynamics of the postacceleration mechanism which allow attosecond electron bunches emitted by a few-cycle laser irradiated droplet to be injected into the laser pulse and therefore be accelerated up to MeV energies. We have found how the ultrarelativistic regime can be applied for intensities starting from the threshold of 5 × 10 19 W cm −2 , for a fewcycle laser pulse. A new integration method has been applied and verified for the determination of the gain in energy by the electron bunch in the laser pulse wave. Three-dimensional and focus effects have been taken into account. We have demonstrated that acceleration in vacuum is possible while the bunches are copropagating with the pulse, and that the predicted energy gain agrees well with PIC simulations for intensities for which the cluster remains opaque, or satisfying the empirical criterion in [11]. We have also calculated the saturation distance at which the bunches acquire their peak and cutoff energies respectively. All numerical results have been successfully compared with post-simulation tracking. The transient region provides a useful interpolation for determining the linear gain function in the case of the few-cycle laser pulse. This interpolation has been carried on analytically and then efficiently compared to the local decomposition of the spectrum in a new way, which settles an immediate connection with the more traditional approach presented in appendix B. It must be noticed that there is some elasticity and that initial conditions are chosen according to equation (10). However, the final results must be calculated for phase-locked bunches and under the choice of the corresponding parameters on the ultrarelativistic cone the MacLaurin development must be used for the fit of the peak energies (see equation (11)). This result actually yields the required linear fit, as it must be inferred from simulation tracking findings.
Both mathematical and numerical calculations are consistent with the previously published experimental results in [16], for the special case of I = 5 × 10 19 W cm −2 intensity. It is important to notice that for the first time a finite integral calculation is carried on for determining the mean-or peak-energy an electron bunch copropagating in vacuum, overcoming the difficulties presented by the fully parametric treatment in [2].
Former analytical results have been confirmed and extended to the three-dimensional case, with the significant advantage that the integral in equation (10) can be solved with proper initial conditions coming from the plasmonic energy of the produced electron cloud, and therefore elevated to a systematic calculation of the energy gain. All results are shown to be dependent on the electric field of the laser, as well as the number of emitted bunches. It is interesting to ask whether there are fundamental limits to scaling this acceleration scheme to higher intensities, which are currently restricted to the maximum reachable intensity with the present laser technology, namely 10 22 W cm −2 . To test this, it is clear that further experimental and simulation work is needed within the context of this novel scheme.
result for a rapidly varying phase; this allows us to choose a selected path for the stochastic approximation where the gain is stable. The reduction, with respect to variational principles already applied in stochastic approximations [36], of the N degrees of freedom down to one possible choice for the values of Δf giving the exact experimental gain was made possible by adapting the functional expressions actually exposed in appendix B. Furthermore, the lagrangian density must be left generic, so that the Euler-Lagrange principle only serves as a confirmation that the field f can be nevertheless extracted, provided that useful conditions are chosen for the Wiener function variation dW(Appendix D). This way, we neglect the number of situations that can be represented mathematically speaking by a stochastic process individuation of any proper kind. In particular, we try to neglect a specific expression for the composition of functions in the functional integral and we choose the general form in appendix B. The quasi-particle approximation actually supports the extraction of information for single electrons outside the bunch as it is stated in appendix C. The second mathematical problem would be constituted by the validity of the choice of the extremes of the integral as a unique path between a initial time and a final time provided that the time flow is relativistic and we are choosing a geometrical measure on the cone. Of course in this case the quantum picture is neglected, also because the material in appendix B is exhaustive enough. In other words, we should not consider correlation between states and therefore the measure in 10 is not a quantum measure. The laser phase can therefore be conceptually described by means of well-known solutions of the differential diffusion equations: the most complete transportation mechanism for the functional integral to the problem of a dynamical system was examined in [37]. Both the original mathematical and physical treatments as they are resumed in classic books about stochastic differentials and relevant stochastic differential equations are explained by means of a special choice in the appendix D, where the passage to the relativistic regime is very well found.

Appendix B. Methodological notes for the functional integral
In this section we aim at justifying the functional approach leading to 10 , by referring to an already known mathematical context. The original literature will be therefore cited outside the relevant physical context because we are interested in the general method. We proceed as following, ignoring the details and keeping everything as rigorous a possible, but without specifying the actual fields and variables involved, in order not to superpose with the text. While keeping the analogy with the formalism of the classical field theory, we recall the same gauge choice in [5], including the condition ∂ μ A μ which allows for the equivalence between the classical Lagrangian and the field theory picture. That allows for the light cone invariant to be chosen in a suitable way, keeping into the classical electrodynamics picture as it is done in [5]. However, this analogy is not enough to our case as we have to take into account a dependence of the integration variable f on proper time. A proper electromagnetic field development as a function of the invariant variables is found in [38]. The mathematical analogy between the definition of functional derivative and functional integral proposed in [38] and ours is based on the dependence of the vector potential on the phase, which in turn is a function of the proper time. It will be shown the usefulness of the definition of functional and functional derivative in this specific context, being the specific approach allowing, in principle, for any envelope variation function. This way, the proper envelope can be chosen at a later time, when the path to explicit calculation has been simplified. In addition to that, the slowly varying envelope can be stated independently, allowing for a specific choice of the laser waveform. The following treatment is thus in principle separated by the pulse approximation, which is flexible, but not from the polarization choice. The polarization choice can affect significantly the gain calculation and will not be analyzed here. The general Euler-Lagrange equation generalized to field theory expressed in relativistic covariant notation reads [38] ¶ ¶F - Here, the four-dimensional gradient was written as usual as ∂Φ r /∂x μ ≡ ∂ μ Φ = (∂ t , ∇)Φ r . The fields Φ r (not to be confused with the laser phase f) can have several independent degrees of freedom Φ r , r = 1,....,N. The variational principle in terms of the Lagrangian reads ò ò We choose as gauge fields the electromagnetic fields. The argument x represents all space coordinates Φ(x 0 , x 1 , The Lagrange density may depend on the field function Φ, on its time derivative F  , and also on the gradient ∇Φ. In the Lorentz gauge ∂ μ A μ = 0. Here we adopt the Coulomb gauge ∇ · A = 0, and according to the functional derivative rules δ∇Φ = ∇δΦ = 0. This gauge requires that the three-dimensional divergence vanishes ∇ · A = 0. The field Φ(x) is assumed to be localized at the space-time point x = {x 0 , x 1 , x 2 , x 3 }. The condition u μ u μ = −1 shows that the velocity vector is time-like, namely, the charged particle lies within the light cone and positive energy trajectories are separated from the negative energy trajectories. Therefore, if the particle starts gaining energy, it will remain on the light cone surface which ensures the positive gain is not perturbed by any deviation. This picture is symmetric under space reversal, but not under time reversal.
The general methods of field quantization in more advanced form can be found in various textbooks on relativistic quantum mechanics and/or field theory. For the field theoretical foundations of photon-electron interaction in vacuum the reader may want to explore the purely quantum electrodynamics approach in [39]. The mathematical treatment presented here has been directed towards applications-as it is widely addressed in [40]with a focus on canonical integration methods. With respect to the discrete systems of point masses , we turn to the study of fields by associating to every point in our finite region of space with a continuous field variable, which is Φ(x, t). The dynamical variables of the theory are now the values of the fields Φ at each point of space at a proper time instant. The Lagrange function now is a functional of the field, i.e. a mapping from a space of functions to the real numbers. Following Hamilton's principle we define the variation of a functional F[Φ(x)] as The variational principle is equivalent to the Euler-Lagrange equations. From the already known quantization procedure of the electromagnetic fieldjust as an example, one can see related treatments in [38,41], the canonical variables can be chosen to be the vector potentials A and P from the ordinary Maxwell equations. However, it must be reminded that the vector potential is not a true physical quantity, namely, it cannot be measured as well as the electric and magnetic field. In electromagnetic field quantization procedures it is convenient to set ∇ · A = ∇ · P = 0. Subsequently, also divE = 0 can be set, in order to eliminate mathematically all the inhomogeneities in the electric field geometrical setup. The condition divE = 0 is a supplementary condition, which allows for avoiding field canonical variables commutation issues. Both divE and divH -where H is the magnetic fieldare constants of motion and therefore can be set equal to zero in every point of space. The resulting integration procedure associated with the variational principle and the Maxwell equations yields the equality where only the time integration is necessary. The Gaussian focus existence can be allowed by creating a divergence-free vector potential, to be substituted explicitly to the vector potential A in the expression equation (B3). The representation chosen by Hartemann [5] reads As long as h(t) is the temporal form of the laser wave, we can consider the phase f = κ μ x μ (τ) for its analytical expression and by a simple substitution we can apply the instrument of the functional derivative to the divergence free vector potential. The Coulomb gauge allows to select the spatial focus terms out of the derivation process and the correspondence h(t) → h(f(τ)) leads to a straightforward application of the chain rule in the differentiation procedure. For the purpose of strict calculation, we use a not entirely rigorous approach in which the functional derivative is defined according to its properties as a mathematical operation.
where one of the two initial functions can be reversed accorded to time. The difference (t − τ) can be shifted according to the initial point of definition of the functions, and τ is the new integration variable; if τ > 0, the upper extreme of the integral would actually be shifted to the right. If α < t < β, with α, β real numbers , the integral must be calculated between t − β and t − α and solved by parts. The convolution of a periodic function with a rectangular function (merge of step functions) yields an ascending ramp with a saturation point, depending on the choice of extremes. In our treatment this choice is due to a bond coming from the transient form of the resolution of the integral. The cross-correlation of the convolution of two functions reads and y * (t) is the complex coniugate of y(t). In the case of the function = y t Arect t T ( ) ( ) the complex conjugate is also real, and coincides with the function itself. The expression in C2 is valid for τ > 0. The analogous expression for τ < 0 is not analyzed here, because of the physical constant on the phase shift imposed by the analytical resolution of integral 10. Generally speaking, integrating with respect to t would be of little physical meaning in every transient regime, since time dependence is not known. Therefore we will not examine any superposition in the phase space integration, even in the transient regime calculation. The present work makes use of the purely mathematical adaptation from the original Arect(), due to the advantageous statistical properties of this operator and and its form, which is tunable on a finite impulse duration.
The statistical filtering of the bunches was made on an assumption of uniformity inside the bunches' oscillations. The modulation in C2 is in principle not a realistic function for an experimentally-like constructed bunch or series of bunches, but allows us to interpret the negligible statistical deviation inside the bunch structure in a more precise way than just assuming a small transverse direction in a gaussian form. It also represents a good fit for the shape of the ascending ramp in figure 4. In the case of a Arect() function shape, filtering easily reduces to mean and max calculations in the first order approach. For a more complete time series approach, one can refer to [42,43]. More insight in the dedicated field is out of the scope of this work. In [44][45][46], first principles for any kind of ultrashort bunch detection have been set up. In the filtered bunches the frequency ω has been compared by means of standard deviation coefficients calculated in a straightforward way, so as to separate the coherent part from the incoherent, complex overall spectrum. The incoherent part represents the easiest subset to be fit, since a real analysis is sufficient. In the following treatment, there are the basic formulas for the dedicated approximation: where  k is the wave number, λ is the wavelength and F(λ) the longitudinal form factor. In turn, F(λ) depends on the three-dimensional charge distribution S(z). The wave-vector orientation with respect to the angle of measurement is not due in this treatment, as long as a correct spatial filtering of the bunches has been carried on.
The shape of the correct electron density distribution has been selected among the two tolerance angles aside the central direction of emission. This is possible by means of both a spatial and temporal filtering, which can be made within the mathematical approximation directly on the particle vectors in python language. The process is dynamical, namely, the size of the arrays can be squashed while maintaining the bunches' profile. This is done at least until the ascending energy ramp is finished. After that, the coherent part of the spectrum begins to apply and, more correctly, an imaginary numbers' spectrum analysis should be carried on. However, the statistics are much heavier here and slight oscillations with the present method should be taken into account. A specific treatment for separating two statistical regimes is out of the scope of the present study.
All of the computational data can be made available from the 3D PIC simulations, through a purely electromagnetic calculation, and most of the findings have been presented in [16]. In particular, the two terms adding up in the F(λ) expression represents, respectively, the incoherent and coherent part of the spectrum, proportional respectively to the total number of particles N and N(N − 1) ∝ N 2 . In particular, if index 1 refers to a reference particle, E i (t) = E 1 (t + Δt i ) is the electric field of the other particles in the bunch affecting the socalled quasi particle travelling and gaining energy in the electric field. This picture does not change in the case in which the electric field actually crosses the trajectory, because we are interested in the frequency spectrum.
A complete calculation of the spectrum after the saturation length has been reached should include also current evaluation outside of the transient treatment and therefore has been avoided here. Tracking itself would be not enough to ensure an overall description, since so far only a description based on the phase is in agreement with invariance on the cone surface where energy is acquired. Therefore only the coherent part of the spectrum was introduced in the post-processing calculations, isolating the frequency in the single bunches after the mean and standard deviation. The quantity W(t) represents the one-dimensional Wiener process whose average is zero and whose stopping time is related to the discretization of several instants Δt j = t j+1 − t j setting the limit of D max t j j close to the final time of the process τ. Analogously, the Wiener process gives the Riemann-Stieltjes measure of the integral ΔW j = W(t j+1 ) − W(t j ). Those instants named as t j can be put in succession as long as we consider this sequence to be tending to a limit of the process time τ constituting the stopping time of the dynamical system. In our case, the stopping time can be made coincident with the few-cycle laser duration as long as we consider acquiring energy from the laser to be a stochastic process whose random average can be simplified by reducing the formula D2 to its first term after considering the mean quantities in the bunch and stating < W(t) > = 0. Consequently the average value of the stochastic integral in dW(s) must be calculated for our case by assuming that the finite stopping time duration. This is compatible with the track analysis that has been performed. The white noise W(t) cannot be directly compared with small variations in mean energies and in cutoff energies in both tables 2 and 3 because such a detailed connection would require a separate research work which is out of the scope of the present study. However, as long as we make the stochastic variable X(t) be the necessary fit for the phase variable f(τ) seen as a function of the proper time by the many particles of the bunch, the parallelism with the variation of the results in the gain calculations in equation (10) appears to be perfectly feasible. Due to the special nature of our calculation for the energy gain, which is based on average quantities and on collective effects from dense correlation phenomena, it is to be expected to get a small error for the fit of the theoretical formula within 5%, compatibly with the assumed tolerance in a Wiener process. In addition to that, it is interesting to consider that the higher extreme of the integral calculation in equation (10) is the phase difference and therefore the initial condition for the equation (D2) can be neglected in our effective calculation. By assuming the phase difference seen by a phase-locked bunch is a random process within the ultrarelativistic cone associated with the very high energy transfer process, we set a limit to the structure of the bunch and an effective connection with the quasiparticle theory as presented in appendix A. The quantum field theory approach in [38] can be integrated by assuming that the integration variable extremes x can be substituted by an arbitrary kind of defined variable. Therefore also the variable of integration can become a function, under the proper assumptions. It is straightforward to notice that the functional integral -first presented in early seventies in [50] has not been invented in a relativistic context, but it is transported to relativistic theory here in this treatment. The passage to the proper time is assumed to be made continuous by the definition for the functional integral as an ultimate limit to sum over small definite time intervals, as it happens for the nonrelativistic treatment. In order to express the formula D2 as a one-dimensional solution, we must make use of a proper approximation. The following can be applied to a gaussian few-cycle laser pulse, being not in contrast with the assumptions made in the paper. The definition of the laser phase f(τ) can be physically related to the diffusion equation known as the solution of the differential stochastical approximation of the system of equations. Hereby we assume that the definition of the second term in equation (D5) is exact and not necessary, as long as we are going to neglect the Wiener process part , after stating that analogous conditions for the expectation values observable quantities with respect to the previous treatment are still valid. Once we have got rid of the Wiener process by means of averaging, it is easy to recognize the usual phase approximation as a function of the laser frequency. As in our treatment we are dealing with a monochromatic wave, no further decomposition for the phase is necessary and the differential variation of the phase can be easily considered as in equation (9a). Therefore we end up with the integral expression in 10. Regarding the point of view of analogous laser-plasma physical processes including stochastic representation for the gain of energy, it is important to mention that the work by [35], already briefly presented in the text, includes a connection between the amount of stochastic heating occurring in a laser standing wave and the angle of incidence for the laser light illumination. In our context, we want to achieve complete hitting of the target in the focus and therefore we can assume that the incidence angle is maximum and that we can consider the full period of the laser as the necessary time-in the rest frame-to let the stochastic process occur.