Next Article in Journal
Numerical Study on Temperature Distribution of Steel Truss Aqueducts under Solar Radiation
Next Article in Special Issue
Adjustment of Subwavelength Rippled Structures on Titanium by Two-Step Fabrication Using Femtosecond Laser Pulses
Previous Article in Journal
Advanced Chirp Transform Spectrometer with Novel Digital Pulse Compression Method for Spectrum Detection
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Optimized Computation of Tight Focusing of Short Pulses Using Mapping to Periodic Space

1
Department of Mathematical Software and Supercomputing Technologies, Lobachevsky State University of Nizhni Novgorod, 603950 Nizhny Novgorod, Russia
2
Mathematical Center, Lobachevsky State University of Nizhni Novgorod, 603950 Nizhny Novgorod, Russia
3
Institute of Applied Physics, Russian Academy of Sciences, 603950 Nizhny Novgorod, Russia
4
Department of Physics, University of Gothenburg, 41296 Gothenburg, Sweden
5
Department of Physics, Umeå University, 90187 Umeå, Sweden
*
Authors to whom correspondence should be addressed.
Appl. Sci. 2021, 11(3), 956; https://doi.org/10.3390/app11030956
Submission received: 28 December 2020 / Revised: 15 January 2021 / Accepted: 18 January 2021 / Published: 21 January 2021
(This article belongs to the Special Issue Ultra-Short Laser Pulses and its Application in Physics)

Abstract

:
When a pulsed, few-cycle electromagnetic wave is focused by optics 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. In many cases of interest, the mapping reduces both run time and memory requirements by a factor of order 10, making it possible to carry out simulations on a desktop machine or a single node of a supercomputer. We provide an open-source C++ implementation with Python bindings and demonstrate its use for a desktop machine, where the routine provides the opportunity to use the resolution sufficient for handling the pulses with spectra spanning over several octaves. The described approach can facilitate the stability analysis of theoretical proposals, the studies based on statistical inferences, as well as the overall development and analysis of experiments with tightly-focused short laser pulses.

1. 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. Tight focusing 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 situations, 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].
Tight focusing of few-cycle pulses is an emerging topic with a broad field of applications [43]. Recently, the few-cycle laser technology has reached the state to produce higher energies much beyond the mJ level, which in combination with tight focusing provides strongly relativistic peak intensities [25] for various promising sources of isolated attosecond XUV pulses [23] and relativistic electron bunches [24]. However, this tight focusing is still very challenging [5,28], especially with the almost octave spanning spectrum of few cycle laser pulses, and it will profit from further support by quick numerical algorithms.
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 numbers of nodes of the grids that sample the field in the far- and near-field zones respectively. Note that the lowest spatial frequency components of the pulse determine the distance to and, thus, the size of the far-field zone, whereas the highest spatial 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 NTN 4 / 3 , where TN 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 [44].
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 [45] (see also [46]). 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.
The presented method and the developed routine can be of interest for both theoretical and experimental studies that involve tightly focused pulses, especially if the pulses have few-cycle durations. We would like to outline three specific potential applications that involve the calculation of the tight focus a large number of times: (1) stability/robustness analysis of the outcome of theoretical proposals with respect to spatial laser-pulse imperfections; (2) iterative determination of the input parameters (beam profile and wavefront) from other experimental observables, e.g., electron spatial distribution in a specific experiment, which can even be performed on many single-shot measurements; (3) generation of large sets of simulations for training machine learning methods employed for diagnosis and improvement of tight focusing by means of adaptive optics.

2. Materials and Methods

2.1. 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 nondegenerate 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, [47,48]). 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, [49,50,51,52,53,54] or chapter 15.9 (b) in [48]). We denote the electromagnetic field and current density in Fourier space by E ^ , B ^ and j ^ respectively. Furthermore, k ^ is the normalized wave vector k , i is the imaginary unit, c is the speed of light, and C = cos ( c k Δ t ) , and S = sin ( c k Δ t ) . Using superscripts n and n + 1 to denote consecutive numerical states, we can express the Fourier method in the following form (see [53], CGS units are used):
E ^ n + 1 = C E ^ n + i S k ^ × B ^ n 4 π S k c j ^ n + 1 / 2 + ( 1 C ) k ^ ( k ^ · E ^ n ) + 4 π k ^ k ^ · j ^ n + 1 / 2 S k c Δ t ,
B ^ n + 1 = C B ^ n i S k ^ × E ^ n + 4 π i 1 C k c k ^ × j ^ n + 1 / 2 .
It is worth noting that Poisson’s equation k ^ · E ^ = 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 nonzero current). To correct for this, we need to subtract the longitudinal component of electric field k ^ ( k ^ · E ^ ) from the initial electric field E ^ once, before the simulation starts (for more details see [53] or Section 6.2 in [54]). This will ensure that Poisson’s equation is satisfied throughout the simulation and that the corresponding term in Equation (1) vanishes. However, to make our solver capable of handling continuous sources, and to avoid the accumulation of numerical errors due to imperfect satisfaction of the continuity equation in nonvacuum simulations (beyond the scope of this article), we carry out this adjustment at each iteration:
E ^ n + 1 = C ( E ^ n k ^ ( k ^ · E ^ n ) ) + i S k ^ × B ^ n
B ^ n + 1 = C B ^ n i S k ^ × E ^ n .
As our implementation follows that referred to as PSATD in Ref. [53], we adopt the same terminology for our spectral solver.

2.2. 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 θ = arctan 1 / ( 2 f ) .
The distance to the far-field zone is determined by the lowest characteristic spatial 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 spatial frequency component. Then we require that the distance R 0 from the initial location of the pulse to the geometrical center (focus position) is much larger than λ . In practice we set R 0 16 λ but this can vary depending on the spectral bandwidth 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 profile, polarization and phase (wavefront) 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 particular 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 the amplifying medium, we consider a pulse with flat-top transverse profile. In this case we need to define the transverse shape u t s so that the pulse amplitude is nonzero 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:
u t s ( α ) = ω α , 0 , θ ε 2 + cos π α θ + ε 2 2 ε 2 · ω α , θ ε 2 , θ + ε 2 ,
where we use the rectangular function
ω x , x m i n , x m a x = 1 for x x m i n , x m a x 0 for x x m i n , x m a x .
We consider a longitudinal shape of the form:
u l ( x ) = sin 2 π x λ · cos 2 π x L · ω x , L 2 , L 2 ,
where L is the pulse length, α is the angle between the x axis and the position vector R = ( x , y , z ) .
Finally, the spherical pulse can be defined as
u x , y , z = u R = A R · u l R R 0 · u t s arcsin y 2 + z 2 R
A = 4 P 0 c ( 1 cos θ ) ,
where P 0 is the input power. We consider linear polarization and set the electromagnetic field as follows:
E ( R ) = u ( R ) s 0
B ( R ) = u ( R ) s 1 ,
where s 0 and s 1 are, respectively, the normalized vectors d × d × R and d × R , and d is the polarization vector.
An example of simulation of the tight focusing problem by the method described in Section 2.1 is demonstrated in Figure 1.

2.3. 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 x R 0 L / 2 , L / 2 y R 0 , R 0 z R 0 , R 0 , which has volume equal to approximately 4 R 0 3 (assuming that R 0 L ). However, geometric arguments indicate that the pulse should remain, approximately, in a thin spherically curved layer R 0 L / 2 c t < R < R 0 + L / 2 c t . Thus only a small fraction of the computational grid actually contains the pulse. If we were using the Finite Difference Time Domain (FDTD) method, we could reduce computational costs by updating the grid nodes in the nonempty 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 Figure 2a). If D is large enough, the pulses overlap neither initially nor during their propagation towards the focal region (see Figure 2b). When each pulse reaches its focus, a sequence of identical structures of focused pulses is formed (see Figure 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. 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 4 D R 0 2 4 R 0 3 (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 neighboring pulses do not overlap. From the geometrical consideration shown in Figure 3, one can obtain the following restriction:
D D m i n = R m i n cos θ + R m a x 2 R m i n 2 sin 2 θ ,
where R m a x = R 0 + L / 2 and R m i n = R 0 L / 2 are the boundaries of the initial pulse in spherical coordinates. Note that one can avoid this restriction by considering a superposition of the overlapping replicated pulses. In this case we can still obtain the separated structure of a single pulse at the focus if D > L due to the linearity of the Maxwell’s equations. However, due to overlap, in this case we are not able to obtain (in a straightforward way) the intermediate states of the pulse on the way to the focus. The choice of D can also depend on the extent of fringes and the required accuracy. There is a trade-off: increasing D 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 m i n D , R m i n with periodic boundary conditions. In order to set the initial conditions we need to perform the following mapping T R : the field state of each point ( x ˜ , y ˜ , 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:
x ˜ , y ˜ , z ˜ = T R = x m i n + D x x m i n D , y , z .
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 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 u ˜ T R within the simulated region according to:
u ( R , t ) = u ˜ ( T R ) · I S ( t ) ( R )
I S ( t ) ( R ) = 0 , r ( t ) < 0 and R l ( t ) or R < r ( t ) or ( x > 0 ) 0 , l ( t ) 0 and r ( t ) 0 and x < 0 and R > L 0 , l ( t ) 0 and r ( t ) 0 and x 0 and R > l ( t ) 0 , l ( t ) > 0 and R l ( t ) or R > r ( t ) or ( x < 0 ) 1 , in other cases
l ( t ) = R m a x + c t , r ( t ) = R m i n + c t , R = ( x , y , z ) , R = | | R | |
The described method for cutting out a single pulse from the periodic sequence is shown in Figure 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 Equations (15) and (16), if D is not large enough. However, in practice, it is usually sufficient to take the parameter D equal to 3 L 4 L 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 Section 3.1.
We would like to conclude this section by making a few general remarks. Firstly, the choice of R 0 in many cases is determined by the transition from the geometrical to wave optics. If the deviations from the spherical shape of the phase front are sufficiently minor, one can combine the proposed method with the computation of the earlier field evolution within geometrical optics. This implies advancing the pulse along the rays pointing to the geometrical center in combination with the adjustment of the field amplitude inverse proportionally to the distance to the center. Note that the effect of imperfect phase front can be also estimated and taken into consideration. Anyway, in many cases it should be possible to compute the pulse at the entrance to the near-field zone, which our method is designed for. Secondly, our method makes use of the solver periodicity along the main direction of propagation (x axis) and effectively reduces the computation in a cube to the computation in a thin layer that spans to a large distance along the transverse directions (y and z axes) only. Using this as a basis, one can further duplicate the pulse several times along the transverse directions, making use of the solver periodicity in these directions and of the linearity of the Maxwell’s equations. In the considered case the effective speedup scales as R 0 , whereas the described follow-up development would yield a speedup that scales as R 0 3 .
Finally, we estimate the expected benefit of using the method. According to the discussion above, the method reduces both computational time and memory demands by a factor of R 0 / D (neglecting the logarithmic factor in the computational complexity of the FFT routine). While the choice of R 0 depends on the required level of accuracy, an estimate appropriate in many cases is 10 times the wavelength of the radiation. For long pulses this would mean R 0 10 λ . However, since the method is designed for very short pulses (with duration equivalent to a few cycles) we assume an octave-spanning spectrum and account for the shape of the pulses by using the longest wavelength component L, i.e., R 0 10 L . As we will see further (see Section 3.1) a relative error in the peak field amplitude of less than 1% is reached at D 2 L . Thus we estimate that the described method can reduce the computational time and memory demands by a factor of order 10 (Section 3.1), depending on the level of required accuracy.

3. Results and Discussion

3.1. Verification and Accuracy Determination

The mapping described in Section 2.3 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 Section 2.2 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 in the region x , y , z [ 20 λ , 20 λ ] . The proposed method can be applied using the region with reduced size D along the x-axis. Because of the large value of θ , the peak intensity is reached almost exactly at the geometrical focus, i.e., at time t = R 0 / c since the start of the simulation. Due to the properties of the spectral solver (Section 2.1), we can do only one iteration of the method with a time step of Δ t = R 0 / c to obtain the electromagnetic field at the instant of focusing. For the simulation, a computational grid with the resolution of 48 points per wavelength along x-axis and 1536 × 1536 points along y- and z-axis ( y , z [ 20 λ , 20 λ ] ) was used. The error 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.5 L 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 Figure 4a. As one can see, even at D m i n / L 1.9 the error does not exceed 0.9%, and it quickly falls below 0.1% at D / L 3.5 . Figure 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.
All of the results in this and further sections were obtained with the open source Hi-Chi framework [45]. The computation time of a typical simulation for different values of parameter D is presented in Figure 4b. The experiments were performed on the MVS-10P supercomputer at the Joint Supercomputing Center of the Russian Academy of Sciences. Each computational node of the supercomputer consists of 2×18-core Intel Xeon Gold 6154 CPU (3.0 GHz) and 192 GB of RAM. The code was built with Intel C++ Compiler 17.0.4, and the high-performance multithreaded FFT implementation from the Intel Math Kernel Library 2017 (MKL) was employed. For double precision, we find that the computation time and memory cost per grid cell are 230 ns (in single-thread mode) and 48 bytes, respectively. Therefore, while the simulation in entire region ( D / L = 20 ) for a grid with the resolution of 28 points per wavelength takes 256 s and 40 GB of memory, the simulation using the presented scheme with D / L = 2 takes only 18 s and 4 GB of memory. Thus, we are able to reduce the calculation time and memory costs by a factor of 10, while the error introduced by the method did not exceed 0.85%.

3.2. Examples

In this section we provide three 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. Finally, in the third example, we apply the method to illustrate the degradation of the intensity profile at focus when using real laser pulses. This shows the necessity to describe accurately the laser pulse for obtaining realistic computational predictions and for interpreting the experimental data.

3.2.1. 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 E = E 0 R 0 sin [ k ( r R 0 ) ] u l ( r R 0 ) u t s ( θ ) / r , where the longitudinal profile u l ( ζ ) = exp [ ζ 2 ln 4 / ( n λ ) 2 ] and the transverse profile u t s ( θ ) = exp ( θ 2 / θ 0 2 ) . 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 / ( 2 f ) . 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 c R 0 2 F ( θ 0 ) ] where F ( θ 0 ) = 0 π sin ( θ ) exp ( 2 θ 2 / θ 0 2 ) d θ . For θ 0 < 2 (equivalent to f > 0.2 ), this auxiliary function is well approximated as F ( θ 0 ) θ 0 2 / 4 θ 0 4 / 48 + θ 0 6 / 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 λ < x < 10 λ and 20 λ < y , z < 20 λ , which is represented by 512 grid points in each dimension.
The peak intensity at focus, as obtained from simulations, is shown in Figure 5. If f < 1 , there are substantial deviations from the theoretical prediction for a paraxial Gaussian beam I 0 = 2 P 0 / ( π w 0 2 ) , where P 0 is the instantaneous power passing through the focal plane z = 0 and the waist w 0 = 2 f λ / π (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 Figure 5b–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.

3.2.2. 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 Section 2.2 with the field configuration parameters presented in Section 3.1. We then apply the described method to compute the peak field strength, as well as when and where it is reached. Figure 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 2 L , so that the error of the scheme does not exceed 1% of the peak amplitude (Section 3.1). 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 s. The complete code by which these results were obtained is given in the Supplementary Materials.
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. One can understand the nature of this mismatch through considering the limit of large f-numbers. When approaching to this limit, the pulse is limited by an increasingly small angle θ and the phase front is increasingly close to plane under this limitation. In this case the pulse almost immediately starts transverse spreading due to diffraction, i.e., the field amplitude starts to decrease. That is why the peak field is achieved far away from the geometrical center of the sphere that defines the phase front. (If we fix the transverse size of the pulse, the center itself tends to be infinitely far from the initial location of the pulse in the limit of large f-numbers.) From this explanation we see that this mismatch is originated from the diffraction of the pulse in the initial state. This means that for larger values of f-number we would need to take a larger value of R 0 to have the initial pulse in the zone of geometrical optics, i.e., so that it is well-described by a spherical wavefront. A simple estimate is to scale R 0 with the Rayleigh length, i.e., R 0 10 λ f 2 . In brief, for large values of f-number, one should either account for the outlined shift (if approximate calculation is sufficient) or adjust the value of R 0 accordingly (preferred).

3.2.3. Focusing of Realistic Laser Pulses

We now modify the flat top profile used previously to introduce perturbation of the phase that could stand for the imperfections of an experimental wavefront, in the case of a f-number f = 1 . The perturbation is described by a sum of Gaussians randomly located in the initial transverse profile, whose individual amplitude varies randomly between ± ε ϕ , but with fixed width σ = 0.3 θ , with θ the opening angle. The number of individual perturbation scales is ( 2 θ / σ ) 2 = 44 . Figure 7 shows how the maximum intensity at focus varies with the standard deviation of the phase σ ( ϕ ) . This clearly indicates that initial perturbations of the wavefront can significantly influence the intensity profile at focus, with a steady decrease of the maximum intensity as the perturbation increases. Thus, one has to give some thoughts to the quality of the laser system when designing experiments, as applications of tight focusing pulses are highly dependent on the maximum intensity.
Another useful feature is that the described method allows for a fast computation of the laser focal spot given the initial beam profile and wavefront measured in an experiment. In Figure 8a the beam profile of the red spectral range (700–1000 nm) of the Light Wave Synthesizer 20 (LWS-20) at Umeå is plotted, while (b) shows the wavefront on the focusing optics, which was originally corrected at the end of the laser and collected extra aberrations during transport in the vacuum compressor and beam line including about 10 mirrors. The calculated LWS-20 focal spot with f = 3 without further corrections is highlighted in Figure 8c. If applied directly in an experiment such a spot is applicable for simulation purposes. Note that other methods can also be used to retrieve beam profiles using experimental inputs, such as the Gerchberg–Saxton algorithm [55] which was, for example, used for laser wakefield simulations with realistic laser profiles [56]. Alternatively, the experimental focus can be optimized with the help of our method coupled with machine learning techniques and using the measured beam profile before focusing and focal spot. This compensates for all added aberrations from the beam transport as well as the focusing optics.

4. Conclusions

We have presented an open-source implementation (see [45]) 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 advance 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 have been validated by the comparison with the straightforward computations (i.e., without using the described method) 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. Two realistic applications are also demonstrated including the focusing of an ideal and a real laser beam profile with imperfect wavefronts. These extend the potential utilization of our technique with performing simulations with real experimental data and an enhanced correction of experimental focus quality even for low repetition rate lasers that will have a high impact on future experiments.

Supplementary Materials

The Hi-Chi project is available online at https://github.com/hi-chi/pyHiChi. The scripts required to reproduce the numerical results are located in the “doc/Tight Focusing” folder.

Author Contributions

Conceptualization, A.G.; methodology, A.G.; software, E.P. and V.V.; validation, E.P., T.B. and A.G.; formal analysis, E.E. and A.G.; investigation, E.P., T.B., J.F., A.M., A.D.A.G., P.F. and L.V.; resources, A.G., I.M. and L.V.; data curation, E.P., T.B., J.F., A.G. and L.V.; writing—original draft preparation, A.G., E.P., T.B., J.F., I.M. and L.V.; writing—review and editing, E.E., T.B., M.M., L.V. and I.M.; visualization, E.P., T.B., J.F. and L.V.; supervision, I.M. and A.G.; project administration, I.M. and A.G.; funding acquisition, A.G. and I.M. All authors have read and agreed to the published version of the manuscript.

Funding

The work is supported by the Lobachevsky University academic excellence program (Project 5-100). A.G. acknowledges the support of the Swedish Research Council (Grant No. 2017-05148), and L.V. also acknowledges the support of the Swedish Research Council (2016-05409 and 2019-02376).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The scripts required to reproduce the numerical results in this paper are included in the Supplemental Material. The code is freely available and may be downloaded from https://github.com/hi-chi/pyHiChi.

Acknowledgments

The authors acknowledge the use of computational resources provided by the Joint Supercomputer Center of the Russian Academy of Sciences and by the Swedish National Infrastructure for Computing (SNIC).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Krausz, F.; Ivanov, M. Attosecond physics. Rev. Mod. Phys. 2009, 81, 163–234. [Google Scholar] [CrossRef] [Green Version]
  2. Mourou, G.A.; Tajima, T.; Bulanov, S.V. Optics in the relativistic regime. Rev. Mod. Phys. 2006, 78, 309–371. [Google Scholar] [CrossRef] [Green Version]
  3. Marklund, M.; Shukla, P.K. Nonlinear collective effects in photon-photon and photon-plasma interactions. Rev. Mod. Phys. 2006, 78, 591–640. [Google Scholar] [CrossRef] [Green Version]
  4. Di Piazza, A.; Müller, C.; Hatsagortsyan, K.Z.; Keitel, C.H. Extremely high-intensity laser interactions with fundamental quantum systems. Rev. Mod. Phys. 2012, 84, 1177–1228. [Google Scholar] [CrossRef] [Green Version]
  5. Yanovsky, V.; Chvykov, V.; Kalinchenko, G.; Rousseau, P.; Planchon, T.; Matsuoka, T.; Maksimchuk, A.; Nees, J.; Cheriaux, G.; Mourou, G.; et al. Ultra-high intensity- 300-TW laser at 0.1 Hz repetition rate. Opt. Express 2008, 16, 2109. [Google Scholar] [CrossRef] [PubMed]
  6. Chatziathanasiou, S.; Kahaly, S.; Skantzakis, E.; Sansone, G.; Lopez-Martens, R.; Haessler, S.; Varju, K.; Tsakiris, G.; Charalambidis, D.; Tzallas, P. Generation of Attosecond Light Pulses from Gas and Solid State Media. Photonics 2017, 4, 26. [Google Scholar] [CrossRef] [Green Version]
  7. Harvey, C.N.; Gonoskov, A.; Ilderton, A.; Marklund, M. Quantum Quenching of Radiation Losses in Short Laser Pulses. Phys. Rev. Lett. 2017, 118, 105004. [Google Scholar] [CrossRef] [Green Version]
  8. Davis, L.W. Theory of electromagnetic beams. Phys. Rev. A 1979, 19, 1177–1179. [Google Scholar] [CrossRef]
  9. Barton, J.P.; Alexander, D.R. Fifth-order corrected electromagnetic field components for a fundamental Gaussian beam. J. Appl. Phys. 1989, 66, 2800–2802. [Google Scholar] [CrossRef]
  10. Sheppard, C.J.R.; Saghafi, S. Electromagnetic Gaussian beams beyond the paraxial approximation. J. Opt. Soc. Am. A 1999, 16, 1381. [Google Scholar] [CrossRef]
  11. Sepke, S.M.; Umstadter, D.P. Analytical solutions for the electromagnetic fields of tightly focused laser beams of arbitrary pulse length. Opt. Lett. 2006, 31, 2589. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Salamin, Y. Fields of a Gaussian beam beyond the paraxial approximation. Appl. Phys. B 2006, 86, 319–326. [Google Scholar] [CrossRef]
  13. Couture, M.; Belanger, P.A. From Gaussian beam to complex-source-point spherical wave. Phys. Rev. A 1981, 24, 355–359. [Google Scholar] [CrossRef] [Green Version]
  14. Narozhny, N.B.; Fofanov, M.S. Scattering of relativistic electrons by a focused laser pulse. J. Exp. Theor. Phys. 2000, 90, 753–768. [Google Scholar] [CrossRef]
  15. Lin, Q.; Zheng, J.; Becker, W. Subcycle Pulsed Focused Vector Beams. Phys. Rev. Lett. 2006, 97, 253902. [Google Scholar] [CrossRef]
  16. Fedotov, A.M.; Korolev, K.Y.; Legkov, M.V. Exact analytical expression for the electromagnetic field in a focused laser beam or pulse. In Physics of Intense and Superintense Laser Fields; Attosecond Pulses; Quantum and Atomic Optics; and Engineering of Quantum Information, Proceedings of the International Conference on Coherent and Nonlinear Optics, Minsk, Belarus, 28 May–1 June 2007; Fedorov, M.V., Sandner, W., Giacobino, E., Kilin, S., Kulik, S., Sergienko, A., Bandrauk, A., Sergeev, A.M., Eds.; SPIE: Bellingham, WA, USA, 2007. [Google Scholar] [CrossRef] [Green Version]
  17. Sapozhnikov, O.A. An exact solution to the Helmholtz equation for a quasi-Gaussian beam in the form of a superposition of two sources and sinks with complex coordinates. Acoust. Phys. 2012, 58, 41–47. [Google Scholar] [CrossRef]
  18. Yu, B.; Lin, Z.; Pu, J. Comparative study on the paraxial approximation errors of tightly focused fundamental Gaussian beams. Opt. Eng. 2019, 58, 1. [Google Scholar] [CrossRef]
  19. Popov, K.I.; Bychenkov, V.Y.; Rozmus, W.; Sydora, R.D. Electron vacuum acceleration by a tightly focused laser pulse. Phys. Plasmas 2008, 15, 013108. [Google Scholar] [CrossRef]
  20. Bochkarev, S.G.; Bychenkov, V.Y. Acceleration of electrons by tightly focused femtosecond laser pulses. Quantum Electron. 2007, 37, 273–284. [Google Scholar] [CrossRef]
  21. Harvey, C.; Marklund, M.; Holkundkar, A.R. Focusing effects in laser-electron Thomson scattering. Phys. Rev. Accel. Beams 2016, 19, 094701. [Google Scholar] [CrossRef]
  22. Gonoskov, A.A.; Korzhimanov, A.V.; Kim, A.V.; Marklund, M.; Sergeev, A.M. Ultrarelativistic nanoplasmonics as a route towards extreme-intensity attosecond pulses. Phys. Rev. E 2011, 84, 046403. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Kormin, D.; Borot, A.; Ma, G.; Dallari, W.; Bergues, B.; Aladi, M.; Földes, I.B.; Veisz, L. Spectral interferometry with waveform-dependent relativistic high-order harmonics from plasma surfaces. Nat. Commun. 2018, 9, 4992. [Google Scholar] [CrossRef] [PubMed]
  24. Cardenas, D.E.; Ostermayr, T.M.; Lucchio, L.D.; Hofmann, L.; Kling, M.F.; Gibbon, P.; Schreiber, J.; Veisz, L. Sub-cycle dynamics in relativistic nanoplasma acceleration. Sci. Rep. 2019, 9, 7321. [Google Scholar] [CrossRef] [Green Version]
  25. Rivas, D.E.; Borot, A.; Cardenas, D.E.; Marcus, G.; Gu, X.; Herrmann, D.; Xu, J.; Tan, J.; Kormin, D.; Ma, G.; et al. Next Generation Driver for Attosecond and Laser-plasma Physics. Sci. Rep. 2017, 7, 5224. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Thiele, I.; Skupin, S.; Nuter, R. Boundary conditions for arbitrarily shaped and tightly focused laser pulses in electromagnetic codes. J. Comput. Phys. 2016, 321, 1110–1119. [Google Scholar] [CrossRef] [Green Version]
  27. Pérez, F.; Grech, M. Oblique-incidence, arbitrary-profile wave injection for electromagnetic simulations. Phys. Rev. E 2019, 99, 033307. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Bahk, S.W.; Rousseau, P.; Planchon, T.A.; Chvykov, V.; Kalintchenko, G.; Maksimchuk, A.; Mourou, G.A.; Yanovsky, V. Characterization of focal field formed by a large numerical aperture paraboloidal mirror and generation of ultra-high intensity (1022 W/cm2). Appl. Phys. B 2005, 80, 823–832. [Google Scholar] [CrossRef] [Green Version]
  29. Gonoskov, A.; Wallin, E.; Polovinkin, A.; Meyerov, I. Employing machine learning for theory validation and identification of experimental conditions in laser-plasma physics. Sci. Rep. 2019, 9, 7043. [Google Scholar] [CrossRef]
  30. Gonoskov, A.; Gonoskov, I.; Harvey, C.; Ilderton, A.; Kim, A.; Marklund, M.; Mourou, G.; Sergeev, A. Probing Nonperturbative QED with Optimally Focused Laser Pulses. Phys. Rev. Lett. 2013, 111, 060404. [Google Scholar] [CrossRef] [Green Version]
  31. Gonoskov, A.; Bashinov, A.; Gonoskov, I.; Harvey, C.; Ilderton, A.; Kim, A.; Marklund, M.; Mourou, G.; Sergeev, A. Anomalous Radiative Trapping in Laser Fields of Extreme Intensity. Phys. Rev. Lett. 2014, 113, 014801. [Google Scholar] [CrossRef] [Green Version]
  32. Gelfer, E.G.; Mironov, A.A.; Fedotov, A.M.; Bashmakov, V.F.; Nerush, E.N.; Kostyukov, I.Y.; Narozhny, N.B. Optimized multibeam configuration for observation of QED cascades. Phys. Rev. A 2015, 92, 022113. [Google Scholar] [CrossRef] [Green Version]
  33. Gonoskov, A.; Bashinov, A.; Bastrakov, S.; Efimenko, E.; Ilderton, A.; Kim, A.; Marklund, M.; Meyerov, I.; Muraviev, A.; Sergeev, A. Ultrabright GeV Photon Source via Controlled Electromagnetic Cascades in Laser-Dipole Waves. Phys. Rev. X 2017, 7, 041003. [Google Scholar] [CrossRef] [Green Version]
  34. Vranic, M.; Grismayer, T.; Fonseca, R.A.; Silva, L.O. Electron–positron cascades in multiple-laser optical traps. Plasma Phys. Control. Fusion 2016, 59, 014040. [Google Scholar] [CrossRef] [Green Version]
  35. Gong, Z.; Hu, R.H.; Shou, Y.R.; Qiao, B.; Chen, C.E.; He, X.T.; Bulanov, S.S.; Esirkepov, T.Z.; Bulanov, S.V.; Yan, X.Q. High-efficiency γ-ray flash generation via multiple-laser scattering in ponderomotive potential well. Phys. Rev. E 2017, 95, 013210. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Efimenko, E.S.; Bashinov, A.V.; Bastrakov, S.I.; Gonoskov, A.A.; Muraviev, A.A.; Meyerov, I.B.; Kim, A.V.; Sergeev, A.M. Extreme plasma states in laser-governed vacuum breakdown. Sci. Rep. 2018, 8, 2329. [Google Scholar] [CrossRef] [Green Version]
  37. Efimenko, E.S.; Bashinov, A.V.; Gonoskov, A.A.; Bastrakov, S.I.; Muraviev, A.A.; Meyerov, I.B.; Kim, A.V.; Sergeev, A.M. Laser-driven plasma pinching in ee+ cascade. Phys. Rev. E 2019, 99, 031201. [Google Scholar] [CrossRef] [Green Version]
  38. Magnusson, J.; Gonoskov, A.; Marklund, M.; Esirkepov, T.Z.; Koga, J.K.; Kondo, K.; Kando, M.; Bulanov, S.V.; Korn, G.; Bulanov, S.S. Laser-Particle Collider for Multi-GeV Photon Production. Phys. Rev. Lett. 2019, 122, 254801. [Google Scholar] [CrossRef] [Green Version]
  39. Magnusson, J.; Gonoskov, A.; Marklund, M.; Esirkepov, T.Z.; Koga, J.K.; Kondo, K.; Kando, M.; Bulanov, S.V.; Korn, G.; Geddes, C.G.R.; et al. Multiple colliding laser pulses as a basis for studying high-field high-energy physics. Phys. Rev. A 2019, 100, 063404. [Google Scholar] [CrossRef] [Green Version]
  40. Danson, C.N.; Haefner, C.; Bromage, J.; Butcher, T.; Chanteloup, J.C.F.; Chowdhury, E.A.; Galvanauskas, A.; Gizzi, L.A.; Hein, J.; Hillier, D.I.; et al. Petawatt and exawatt class lasers worldwide. High Power Laser Sci. Eng. 2019, 7. [Google Scholar] [CrossRef]
  41. Bulanov, S.S.; Esirkepov, T.Z.; Thomas, A.G.R.; Koga, J.K.; Bulanov, S.V. Schwinger Limit Attainability with Extreme Power Lasers. Phys. Rev. Lett. 2010, 105, 220407. [Google Scholar] [CrossRef] [Green Version]
  42. Gonoskov, I.; Aiello, A.; Heugel, S.; Leuchs, G. Dipole pulse theory: Maximizing the field amplitude from 4π focused laser pulses. Phys. Rev. A 2012, 86, 053836. [Google Scholar] [CrossRef]
  43. Naumova, N.M.; Nees, J.A.; Mourou, G.A. Relativistic attosecond physics. Phys. Plasmas 2005, 12, 056707. [Google Scholar] [CrossRef] [Green Version]
  44. Blinne, A.; Schinkel, D.; Kuschel, S.; Elkina, N.; Rykovanov, S.G.; Zepf, M. A systematic approach to numerical dispersion in Maxwell solvers. Comput. Phys. Commun. 2018, 224, 273–281. [Google Scholar] [CrossRef] [Green Version]
  45. hi-χ Project. Available online: https://github.com/hi-chi (accessed on 12 December 2020).
  46. Muraviev, A.; Bashinov, A.; Efimenko, E.; Volokitin, V.; Meyerov, I.; Gonoskov, A. Strategies for particle resampling in PIC simulations. Comput. Phys. Commun. 2021, 107826. [Google Scholar] [CrossRef]
  47. Weideman, J.A.C.; Herbst, B.M. Split-Step Methods for the Solution of the Nonlinear Schrödinger Equation. SIAM J. Numer. Anal. 1986, 23, 485–507. [Google Scholar] [CrossRef]
  48. Birdsall, C.K.; Langdon, A.B. Plasma Physics via Computer Simulation; IOP: Bristol, UK, 1991. [Google Scholar]
  49. Haber, I.; Lee, R.; Klein, H.; Boris, J. Advances in electromagnetic simulation techniques. In Proceedings of the Sixth Conference Numerical Simulation Plasmas, Berkeley, CA, USA, 18 July 1973; pp. 46–48. [Google Scholar]
  50. Lin, A.T. Application of electromagnetic particle simulation to the generation of electromagnetic radiation. Phys. Fluids 1974, 17, 1995. [Google Scholar] [CrossRef]
  51. Buneman, O.; Barnes, C.; Green, J.; Nielsen, D. Principles and capabilities of 3-D, E-M particle simulations. J. Comput. Phys. 1980, 38, 1–44. [Google Scholar] [CrossRef]
  52. Gustafsson, B.; Kreiss, H.O.; Oliger, J. Time Dependent Problems and Difference Methods; John Wiley & Sons: Hoboken, NJ, USA, 1995; Volume 24. [Google Scholar]
  53. Vay, J.L.; Haber, I.; Godfrey, B.B. A domain decomposition method for pseudo-spectral electromagnetic simulations of plasmas. J. Comput. Phys. 2013, 243, 260–268. [Google Scholar] [CrossRef]
  54. Gonoskov, A. Ultra-Intense Laser-Plasma Interaction for Applied and Fundamental Physics. Ph.D. Thesis, Umeå University, Umeå, Sweden, 2013. [Google Scholar]
  55. Gerchberg, R.W.; Saxton, W. A practical algorithm for the determination of phase from image and diffraction plane pictures. Optik 1972, 35, 237–246. [Google Scholar]
  56. Ferri, J.; Davoine, X.; Fourmaux, S.; Kieffer, J.; Corde, S.; Phuoc, K.T.; Lifschitz, A. Effect of experimental laser imperfections on laser wakefield acceleration and betatron source. Sci. Rep. 2016, 6, 27846. [Google Scholar] [CrossRef]
Figure 1. An example of simulation of the tight focusing problem. The electric field intensity is shown for (a) t = R 0 / c , (b) t = R 0 / 2 c , (c) t = 0 and (d) t = R 0 / 2 c with f-number = 0.3 ( θ 1 rad), λ = 1 μm, R 0 = 16 λ , L = 2 λ , ε = 0.1 rad, P 0 = 1 W. The pulse propagates towards positive x direction, the transverse directions are spanned by y and z axes (see Equation (8)).
Figure 1. An example of simulation of the tight focusing problem. The electric field intensity is shown for (a) t = R 0 / c , (b) t = R 0 / 2 c , (c) t = 0 and (d) t = R 0 / 2 c with f-number = 0.3 ( θ 1 rad), λ = 1 μm, R 0 = 16 λ , L = 2 λ , ε = 0.1 rad, P 0 = 1 W. The pulse propagates towards positive x direction, the transverse directions are spanned by y and z axes (see Equation (8)).
Applsci 11 00956 g001
Figure 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 / 2 c , (c) t = 0 and (d) t = R 0 / 2 c with f-number = 0.3 ( θ 1 rad), λ = 1 μm, R 0 = 16 λ , L = 2 λ , ε = 0.1 rad, P 0 = 1 W.
Figure 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 / 2 c , (c) t = 0 and (d) t = R 0 / 2 c with f-number = 0.3 ( θ 1 rad), λ = 1 μm, R 0 = 16 λ , L = 2 λ , ε = 0.1 rad, P 0 = 1 W.
Applsci 11 00956 g002
Figure 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 Equation (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.
Figure 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 Equation (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.
Applsci 11 00956 g003
Figure 4. (a) 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 . (b) Computational time of one iteration of the method (grid initialisation, forward FFT, spectral solver execution, and backward FFT), as a function of the parameter D / L . All the simulations were performed using Hi-Chi. Computation time of sequential (all lines except the brown one) and multithreaded (brown line, 36 physical cores) versions of the code are presented. Size of the computational grid varies with the parameter D and is equal to ( 56 D / L ) × 896 × 896.
Figure 4. (a) 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 . (b) Computational time of one iteration of the method (grid initialisation, forward FFT, spectral solver execution, and backward FFT), as a function of the parameter D / L . All the simulations were performed using Hi-Chi. Computation time of sequential (all lines except the brown one) and multithreaded (brown line, 36 physical cores) versions of the code are presented. Size of the computational grid varies with the parameter D and is equal to ( 56 D / L ) × 896 × 896.
Applsci 11 00956 g004
Figure 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 z = 0 , at focus, from simulations with an f-number of 1. (d) Comparison of theoretical (solid colors) and simulation (black, dashed) results for the transverse field along y = z = 0 (blue) and the longitudinal field along y = λ and z = 0 , at focus, for f = 1 .
Figure 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 z = 0 , at focus, from simulations with an f-number of 1. (d) Comparison of theoretical (solid colors) and simulation (black, dashed) results for the transverse field along y = z = 0 (blue) and the longitudinal field along y = λ and z = 0 , at focus, for f = 1 .
Applsci 11 00956 g005
Figure 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 μm (see the explanation in the text).
Figure 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 μm (see the explanation in the text).
Applsci 11 00956 g006
Figure 7. Evolution of the peak amplitude a 0 at focus as a function of the standard deviation of the initial phase ( σ ( ϕ ) in radian, solid line) and standard deviation of a 0 over 20 runs (dashed lines). Inlets exhibit examples of the normalized transverse intensity profile obtained at focus for two different values of σ ( ϕ ) . Transverse directions are normalized to λ .
Figure 7. Evolution of the peak amplitude a 0 at focus as a function of the standard deviation of the initial phase ( σ ( ϕ ) in radian, solid line) and standard deviation of a 0 over 20 runs (dashed lines). Inlets exhibit examples of the normalized transverse intensity profile obtained at focus for two different values of σ ( ϕ ) . Transverse directions are normalized to λ .
Applsci 11 00956 g007
Figure 8. Example of focal spot calculated from realistic profile: (a) typical transverse intensity profile (a. u.) and (b) wavefront (in radians) without correction of extra aberrations during beam transport from the end of the laser to the experiment and (c) normalized focal intensity obtained applying our method with a f-number f = 3 .
Figure 8. Example of focal spot calculated from realistic profile: (a) typical transverse intensity profile (a. u.) and (b) wavefront (in radians) without correction of extra aberrations during beam transport from the end of the laser to the experiment and (c) normalized focal intensity obtained applying our method with a f-number f = 3 .
Applsci 11 00956 g008
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Panova, E.; Volokitin, V.; Efimenko, E.; Ferri, J.; Blackburn, T.; Marklund, M.; Muschet, A.; De Andres Gonzalez, A.; Fischer, P.; Veisz, L.; et al. Optimized Computation of Tight Focusing of Short Pulses Using Mapping to Periodic Space. Appl. Sci. 2021, 11, 956. https://doi.org/10.3390/app11030956

AMA Style

Panova E, Volokitin V, Efimenko E, Ferri J, Blackburn T, Marklund M, Muschet A, De Andres Gonzalez A, Fischer P, Veisz L, et al. Optimized Computation of Tight Focusing of Short Pulses Using Mapping to Periodic Space. Applied Sciences. 2021; 11(3):956. https://doi.org/10.3390/app11030956

Chicago/Turabian Style

Panova, Elena, Valentin Volokitin, Evgeny Efimenko, Julien Ferri, Thomas Blackburn, Mattias Marklund, Alexander Muschet, Aitor De Andres Gonzalez, Peter Fischer, Laszlo Veisz, and et al. 2021. "Optimized Computation of Tight Focusing of Short Pulses Using Mapping to Periodic Space" Applied Sciences 11, no. 3: 956. https://doi.org/10.3390/app11030956

APA Style

Panova, E., Volokitin, V., Efimenko, E., Ferri, J., Blackburn, T., Marklund, M., Muschet, A., De Andres Gonzalez, A., Fischer, P., Veisz, L., Meyerov, I., & Gonoskov, A. (2021). Optimized Computation of Tight Focusing of Short Pulses Using Mapping to Periodic Space. Applied Sciences, 11(3), 956. https://doi.org/10.3390/app11030956

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop