Optimized computation of tight focusing of short pulses using mapping to periodic space

When a pulsed, few-cycle electromagnetic wave is focused by an optic with f-number smaller than two, the frequency components it contains are focused to different regions of space, building up a complex electromagnetic field structure. Accurate numerical computation of this structure is essential for many applications such as the analysis, diagnostics and control of high-intensity laser-matter interactions. However, straightforward use of finite-difference methods can impose unacceptably high demands on computational resources, owing to the necessity of resolving far-field and near-field zones at sufficiently high resolution to overcome numerical dispersion effects. Here we present a procedure for fast computation of tight focusing by mapping a spherically curved far-field region to periodic space, where the field can be advanced by a dispersion-free spectral solver. We provide an open-source C++ implementation with Python bindings and demonstrate its use for a desktop machine, where it can handle up to 100 harmonics.


I. INTRODUCTION
There are well-established routines for computing the evolution of an electromagnetic field in many areas of computational physics.Nevertheless, in some applications, the specific properties of the problem allow for tailored approaches that can reduce computational costs.Tight focusing of short electromagnetic pulses is one of such problems.It is of growing interest in many research areas, including attosecond physics [1], high-intensity laser-plasma interactions [2,3], as well as laser-based studies of fundamental quantum systems [4].Apart from being central for reaching high field strengths for a given radiation power [5], the use of tight focusing of few-cycle pulses with specific phase-and polarization properties enables opportunities for controlling laser-matter interaction scenarios [6,7].
While analytical models of a tightly focused electromagnetic field (e.g.paraxial [8][9][10][11][12] or beyond paraxial [13][14][15][16][17]) are widely used in many theoretical studies, numerical computation provides a direct approach that is of particular interest, as it avoids intrinsic inaccuracies associated with the underlying assumptions [18], and is more flexible in terms of initial conditions.Paraxial models, in particular, require that the diffraction angle be small and therefore become increasingly inaccurate as the focusing becomes stronger.Numerical computations can therefore be necessary when tight focusing is used to achieve strong fields at a target, and the interaction process is sensitive to the field structure in the focal region.This is the case in many cases, such as vacuum electron acceleration [19,20] and the generation of radiation using laser-driven individual [21] or collective [22][23][24] dynamics of electrons [19][20][21][22][23][24].Numerical computations can be also crucial when the laser pulse is only few cycles long [25] and, therefore, cannot be treated analytically as monochromatic.
Besides the use in simulations [26,27], numerical computations can also support experimental efforts to achieve high intensity, e.g. by controlling adaptive optics [28] or by retrieving information from the measured output based on the solution of inverse problem [29].Numerical computations are also of clear interest for designing experiments [30][31][32][33][34][35][36][37][38][39] at the next generation large-scale laser facilities [40], where strong electromagnetic fields are likely to be reached by the combination of several, tightly focused laser pulses [41,42].
The numerical computation of electromagnetic field in the focal region for a given field in far-field zone can be performed via direct integration using the Green's function.In the most general case, when the phase and/or polarization properties vary across the transverse direction of the pulse to be focused, the computational costs are proportional to N M , where N and M are the number of nodes of the grids that sample the field in the farand near-field zones respectively.Note that the lowest frequency components of the pulse determine the distance from and, thus, the size of the far-field zone, whereas the highest frequency components determine the necessary resolution.This leads to an extensive growth of computational costs in the case of tight focusing (i.e., large diffraction angle) and short pulse duration, which is of interest in the context of the problems discussed above.The use of finite-difference time-domain (FDTD) methods for evolving the field from the far-field to the near-field zone leads to computational costs that are proportional to N T ∼ N 4/3 , where T ∼ N 1/3 is the number of time steps to be performed.Such computations are subject to numerical dispersion, which cannot be removed along all the directions simultaneously and may affect significantly high-frequency components of the pulse [43].
Another attractive option is to use a spectral field solver, since it is free of numerical dispersion and it is possible to advance the field over large time intervals in a single step, leading to the numerical costs proportional to N ln N , assuming the use of a fast Fourier transform (FFT).In this case, however, the computations can still be hampered by the extensive size of the grid used for sampling the pulse in the far-field zone.The use of supercomputers, on the other hand, is complicated by the nonlocal data processing patterns required by the FFT routine.
In this article, we present a computational scheme that makes use of the fact that a short pulse in the far-field zone is located within a thin spherical layer and thus occupies only a small fraction of the entire simulation domain.Taking advantage of the periodicity of spectral solvers, we map the far-field zone onto a smaller grid, achieving a significant reduction of computational costs.The implementation of the spectral solver and the proposed scheme was performed as part of the hi-χ (Hi-Chi, High-Intensity Collisions and Interactions) project [44] (see also [45]).This open-source software package is a collection of Python-controlled tools for performing simulations and data analysis in the research area of strong-field particle and plasma physics.All the components are being developed in C++ and optimized for high performance CPUs, which allows for efficient performance of resource-intensive calculations.
As part of the Hi-Chi project, we have implemented an easy-to-use and high-performance tool for modeling tightly focused radiation that can be used for computations on desktops (the supercomputer version is under development).In this work we validate the developed computational module and use it to characterize the properties of focused radiation, namely, the dependency of peak amplitude and the location of the attained peak on the f -number.

II. SIMULATION METHOD A. Spectral solvers
An essential part of the scheme we present is the use of a spectral (Fourier) method to advance the electromagnetic field in a region with periodic boundary conditions.In this subsection we briefly outline the main equations and practical aspects of this approach.
The central idea of spectral methods for time-dependent problems is to make use of the fact that, for linear homogeneous equations with a well-determined non-degenerate spectrum, the state can be represented via eigenstates that can be independently advanced over an arbitrary large time interval without computational errors.When eigenstates correspond to plane waves, the transformation from the representation on a spatial grid to the eigenbasis representation (and vice versa) can be done with a FFT.Source terms (i.e.inhomogeneous part) and/or nonlinear terms can be accounted for by the step-splitting technique (see, for example, [46,47]).Since the propagators for these terms do not necessarily commute with that for the homogeneous equation, this technique introduces errors that grow with, and thus restrict, the size of the time step.Nevertheless, since we consider the propagation of electromagnetic radiation in vacuum, the source term (governed by the current) is zero and we thus can do arbitrarily large time steps.
The application of split-step spectral method to the Maxwell equations has been carried out and used in many numerical studies (see, for example, [48][49][50][51][52][53] or chapter 15.9 (b) in [47]).We denote the electromagnetic field and current density in Fourier space by Ê, B and ĵ respectively.We have also that k is the normalized wave vector k, i is the imaginary unit, c is the speed of light, C = cos(ck∆t), and S = sin(ck∆t).Using superscripts n and n + 1 to denote consecutive numerical states, we can express the Fourier method in the following form (see [52], CGS units are used): It is worth noting that Poisson's equation k • Ê = 0 is not satisfied here automatically but maintained approximately.This may lead to the appearance of static electric fields that arise from a spurious charge associated with either the initial field state or accumulated computational errors (in case of computations with non-zero current).Applying the related adjustment we obtain (for more details see [52] or section 6.2 in [53]): As our implementation follows that referred to as PSATD in Ref. [52], we adopt the same terminology for our spectral solver.

B. Problem statement
We assume that, at the start of the simulation, the pulse that is to be advanced to focus is located within a spherical sector, defined by a polar angle θ around the negative x axis that can be as large as π/2.The opening angle θ is related to the f -number, f , as The distance to the far-field zone is determined by the lowest characteristic frequency component of the pulse.Since we are interested in the treatment of few-cycle pulses, we can take the central wavelength λ (the one that corresponds to the mean of the frequency band) as the typical scale for the lowest frequency component.We then must have that the distance R 0 from the initial location of the pulse to the geometrical center is much larger than λ.In practice we set R 0 ≈ 16λ but this can vary depending on the bandwidth of the spectrum and on the required level of accuracy.
The pulse can have arbitrary profiles along the longitudinal and transverse directions, i.e. perpendicular to, and along, the spherical surface respectively.The pulse can have arbitrary variation of shape, polarization and phase along the transverse directions, but the variations should not be too rapid, so that we still can assume that locally the pulse propagates predominantly towards the geometrical center of the sphere.The permitted magnitude of variations depends on the particular problem and on the required accuracy level.This is a matter of numerical verification.
Although the problem statement is not restricted to any specific form, we consider a specific case of a short pulse under tight focusing.This facilitates the description of the method and provides a good demonstration of its capabilities.
To mimic a laser beam being uniformly amplified within its gain medium, we consider a pulse with flat-top transverse profile.In this case we need to define the transverse shape u ts so that the pulse amplitude is non-zero only within the angle θ.In order to provide a smooth transition near this angular limit, we introduce an edge smoothing angle ε and use cos 2 shape: where we use the rectangular function We consider a longitudinal shape of the form: where L is the pulse length, α is the angle between the x axis and the position vector Finally, the spherical pulse can be defined as where P 0 is the input power.We consider linear polarization and set the electromagnetic field as follows: where s 0 and s 1 are, respectively, the normalized vectors p × p × R and p × R and p is the polarisation vector.
An example of simulation of the tight focusing problem by the method described in Sec.II A is demonstrated in Fig. 1.

C. Mapping to and from the computational subregion
The stated problem of advancing electromagnetic field can be solved numerically by sampling and iterating the electromagnetic field in the region Thus only a small fraction of the computational grid actually contains the pulse.If we were using a FDTD method, we could reduce computational costs by updating the grid nodes in the non-empty regions only.However, as we wish to avoid numerical dispersion, we use a Fourier-based solver and thus do not have such a straightforward opportunity.We now describe how we can effectively exclude the empty regions and achieve a significant reduction of computational costs in this case.
Let us consider a periodic field structure containing a series of repetitions of the pulse, spaced apart by a distance D along the x axis (see fig. 2a).If D is large enough, the pulses overlap neither initially nor during their propagation towards the focal region (see fig. 2b).
When each pulse reaches focus, a sequence of identical structures of focused pulses is formed (see fig. 2c).These structures are spaced apart by the distance D and have only a minor overlap due to diffractive spreading of the signal from the fringes of the initial structure.
One can see that we can simulate this process by computing the evolution of the field in a region that has periodic boundary conditions and size D along x axis.The volume of this region is 4DR 2 0 4R 3 0 (we assume that D ∼ L R 0 ).Periodic boundary conditions are effectively provided by the Fourier solver by default.In this way we can obtain the numerical solution of the stated problem at reduced computational cost, with respect to both run time and memory.Let us describe particular details of our implementation.
The choice of parameter D can be restricted by requiring that neighbouring pulses do not overlap.From the geometrical consideration shown in fig.3, one can obtain the following restriction: where fringes and the required accuracy.There is a trade-off between increasing D, which reduces the errors, but increases the computational costs of the method.
Although this choice of periodic domain is arbitrary with respect to translation along x, it is useful to describe a particular implementation in detail.We suggest to carry out calculations in a layer x ∈ [−R min − D, −R min ] with periodic boundary conditions.In order to set the initial conditions we need to perform the following mapping: the field state of each point (x, ỹ, z) of the simulated subregion is assigned to that of point (x, y, z) of the original unlimited space, for which the problem is stated: Here {•} is the fractional part of an argument.
Unfortunately, there is no reverse mapping.However, by considering the geometrical propagation of the pulse, we can "cut out" the primary pulse from the periodic space.We FIG. 3. Schematic illustration of the proposed method in 3D space (a) and the cross-section in the x-y plane at z = 0 (b).The restriction given by Eq. 12 is derived from requiring that the initial (red) and replicated (blue) pulses do not overlap.The shaded region is the subregion where the simulation is performed.
do so in the following way: for each point R of the original space the field state u = (E, B) is assigned using the field state ũ (T (R)) within the simulated region according to: 1, in other cases ( 15) The described method for cutting out a single pulse from the periodic sequence is shown in fig. 2 using black contours.
This method does not always provide an ideal cut-off of secondary pulses, which can still fall into the region S(t) denoted by Eq. 15-16, if D is not large enough.However, in practice, it is usually sufficient to take the parameter D equal to 3L − 4L in order to avoid notable presence of the secondary pulses.However, depending on the shape of the transverse component of the spherical wave, the fringes may show up in the obtained results and indicate the necessity of choosing larger value of D. The choice of the parameter D is described in more detail in Sec.III.

III. VERIFICATION AND ACCURACY DETERMINATION
The mapping described in Sec.II C allows us to reduce the time and memory costs significantly.However, the question arises as to how accurate the computational results are.
As compared to the direct computation, the only possible inaccuracy is the one caused by the overlap of the pulse in the periodic space.The accuracy therefore can be controlled by the ratio D/L, which defines the relative size of margins used for avoiding overlaps.For a given required accuracy, the acceptable value of D/L depends on the problem parameters, primarily on the f-number and on the sharpness of the fringes mentioned above.To provide a general idea about the achievable accuracy and the corresponding values of D/L, we compare the results of the described method to the results obtained via direct computation.
In order to verify and examine the capabilities of the method in the most extreme conditions we take the field configuration described in Sec.II B and consider a large value of the opening angle θ = 1 (f-number ≈ 0.3) and short pulse duration L = 2λ.We take the edge smoothing angle ε = 0.1.We assume that the pulse is focused along the x-axis and is initially located within the spherical layer in the region x < 0 at a sufficiently large distance R 0 = 16λ from the origin.For this case the direct computation can be performed was calculated as the maximal difference of the field measured between the cases of direct computation and the computation with the considered method.
The computations for D = 7 = 3.5L show that the error introduced by the presented method is 0.12% relative to the peak amplitude obtained in the direct computation.The time required for the computation is reduced by a factor of almost 6.The memory costs for the presented method were 55 GB RAM, instead of 300 GB RAM required for the direct computation.
The dependence of the relative error on D for a grid with the resolution of 12 points per wavelength along x axis is shown in Fig. 4. As one can see, even at D min /L ≈ 1.9 the error does not exceed 0.9%, and it quickly falls below 0.1% at D/L ≈ 3.5.Fig. 2 shows that the small error values at D/L > 3.5 can be attributed to the diffraction of the fringes that do not focus with the remainder of the pulse.

IV. EXAMPLES
In this section we provide two illustrative examples.In the first example we compute focusing of Gaussian beams for different values of f-number and compare the focused beam field with that obtained analytically within the paraxial approximation.This example both provides independent validation based on the known analytical result and indicates the importance of using the complete computation for small values of f-number.In the second example we again consider the described case of a short pulse with a circular flat-top transverse profile and compute the peak field strength as a function of f-number.The obtained data provides an opportunity to assess the potential benefit of using tight focusing for increasing peak field strength at current and future laser systems.
A. Focusing of a Gaussian laser pulse In this subsection we compare our simulation results with theoretical predictions for a focused Gaussian beam.In this case, a simulation is initialized with a linearly polarized, spherical pulse with transverse electric field ).Here n is the number of wavelengths corresponding to the pulse duration (full width at half maximum) and the diffraction angle θ 0 is related to the f-number f by tan θ 0 = 1/(2f ).The amplitude E 0 is determined by requiring that the instantaneous power transmitted through the spherical surface r = R 0 (the pulse peak) is equal to a fixed value P 0 , i.e., E 0 = P 0 /[πε 0 cR 2 0 F (θ 0 )] where F (θ 0 ) = π 0 sin(θ) exp(−2θ 2 /θ 2 0 ) dθ.For θ 0 < 2 (equivalent to f > 0.2), this auxiliary function is well approximated as F (θ 0 ) θ 2 0 /4 − θ 4 0 /48 + θ 6 0 /960.We set the wavelength λ = 1.0 micron, the number of cycles n = 4, the initial position of the pulse R 0 = 20λ, the peak power P 0 = 1 W, and vary f from 0.25 to 2. The simulation covers the domain −30λ < z < 10λ and −20λ < x, y < 20λ, which is represented by 512 grid points in each dimension.
The peak intensity at focus, as obtained from simulations, is shown in Fig. 5.If f < 1, there are substantial deviations from the theoretical prediction for a paraxial Gaussian beam , where P 0 is the instantaneous power passing through the focal plane z = 0 and the waist w 0 = 2f λ/π (see [12]).As the latter is calculated to lowest order in the diffraction angle, or for large f , this is to be expected.A more detailed comparison for the specific case f = 1 is shown in Fig. 5(b, c, d).While the transverse field along the laser axis is reasonably well predicted by theory, the longitudinal component is not.This supports the necessity of going beyond the paraxial approximation when considering tightly focused laser pulses.However, we also confirm that our simulations correctly reproduce the expected theoretical result in the limit that the focusing becomes weak.

B. Focusing of a laser pulse with a circular flat-top transverse profile
In this subsection we apply the developed method for computing tight focusing of laser beams with more realistic structure.To account for the uniform amplification within the gain medium and the importance of not exceeding the breakdown intensity threshold, we consider the flat-top model of the pulse described in Sec.II B with the field configuration parameters presented in Sec.III.We then apply the described method to compute the peak field strength, as well as when and where it is reached.Fig. 6 shows the dependency of the peak field amplitude and the coordinate at which it is reached on the f-number.The distance D between pulses was taken to be 2L, so that the error of the scheme does not exceed 1% of the peak amplitude (Sec.III).The size of the computational grid was 128×1024×1024 nodes, and the instant of reaching peak power was determined with an accuracy of ∆t = 0.01.The complete code by which these results were obtained is given in the hi-χ repository [44].
Note that for f-number 0.6 the peak field strength is achieved almost at the geometrical center of the initial spherical phase front.Nevertheless, for larger values of f-number the place of reaching the strongest field moves away from the geometrical centre with increase of f-number.This is the result of competition between geometrical propagation and diffraction.
As the f -number increases, the pulse must be further from the origin before its wavefronts become approximately spherical: in the near-field zone, they are planar.Thus a pulse initialized with a spherical wavefront has a smaller radius of curvature and focuses 'early', at an x coordinate that tends to negative infinity.

V. CONCLUSION
We have presented an open-source implementation (see [44]) of a method for optimized numerical computation of tight focusing of short electromagnetic pulses with arbitrary longitudinal and transverse variation of polarization, phase and amplitude.The method is based on mapping 3D space to and from a periodic domain, which can have a much smaller size than that required for the straightforward computation.For example, for the case of a two-cycle pulse focused with f-number = 0.3, the computational costs for both run time and memory can be reduced by a factor of almost 6, while the introduced error with respect to the field strength does not exceed 0.1% of the peak value.
We use a spectral solver to avoid numerical dispersion and to make it possible to ad-vance the field state over an arbitrary time interval within one iteration.This provides an opportunity to perform computations of practical interest within minutes on a personal computer.
The method and its implementation has been validated by the comparison with the direct computation and also with the results of analytical computations based on the paraxial approximation for Gaussian beams.We have also considered tight focusing of a laser pulse with a circular flat-top transverse profile.We apply the method to compute the peak achievable intensity as a function of f-number, indicating the potential benefit of using tight focusing for reaching high field strength with current and future laser facilities.

VI. DATA AVAILABILITY
The scripts required to reproduce the numerical results in this paper are included in the hi-χ github repository [44].The code is freely available and may be downloaded from [44].

2
FIG. 2. Illustration of the idea of the proposed method.The pictures show the evolution of the electromagnetic field composed by the equidistant replication of the initial pulse along the x axis: the field intensity is shown for (a) t = −R 0 /c, (b) t = −R 0 /2c, (c) t = 0 and (d) t = R 0 /c.
FIG. 4. Relative error of the peak amplitude, comparing the computations in a subregion to that in the full domain, as a function of D, the subregion size.Here f-number = 0.3.

7 FIG. 5 .
FIG. 5. Comparison of numerical and analytical results for a tightly focused pulse with Gaussian temporal and angular profiles and fixed peak power P = 1 W: (a) The peak intensity at focus.The (b) transverse and (c) longitudinal fields in the plane y = 0, at focus, from simulations with an f -number of 1.(d) Comparison of theoretical (solid colours) and simulation (black, dashed) results for the transverse field along x = y = 0 (blue) and the longitudinal field along x = λ and y = 0, at focus, for f = 1.

FIG. 6 .
FIG. 6.The dependence of (a) the peak value of cycle-averaged intensity on f-number and (b)the coordinate x at which the peak is reached, at a fixed peak power P 0 = 1 W and wavelength λ = 1 micron.
which has volume equal to approximately 4R 3 0 (assuming that R 0 L).However, geometric arguments indicate that the pulse should remain, approximately, in a thin spherically curved layer R 0