Linear scattering theory of short-pulse reflectometry

The linear scattering theory of short pulse reflectometry is presented. An expression for the scattering signal is obtained by applying the perturbation theory approach for solving the Helmholtz equation. Based on these analytical results, a method for measuring the radial wavenumber spectrum of the turbulence is suggested. Analytical results are validated against full-wave numerical modelling.


Introduction
Microwave diagnostics are widely used in fusion research for both routine measurements of electron density profiles [1] and the study of plasma turbulence [2][3][4][5]. While they benefit from simplicity of use, interpretation of the measurements is often challenging due to the complex propagation of electromagnetic waves in plasmas. In particular, the radar pulse reflectometry technique was proposed and employed for density profile measurements [6]. The diagnostic has the benefit of operating in the time-domain, making it possible to directly measure the delay of the reflected pulses and to separate them from parasitic reflections and scattering away from the cut-off.
Considerable experimental [7][8][9] and theoretical [10,11] work was carried out to develop the diagnostic, most of it focused on density profile measurements rather than turbulence studies. Recently, a short-pulse reflectometry * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. (SPR) system has been developed and employed on the TCV tokamak [12,13]. It is capable of using very short (<ns) microwave probing pulses and digitally recording the envelope of the reflected pulse.
So far, the interpretation of the corresponding experimental measurements has not taken into account the presence of electron density fluctuations induced by the plasma turbulence. These density fluctuations could potentially influence the measured delays, affecting the reconstructed density profile. Moreover, if the properties of the pulse can be directly linked to those of the turbulence, as was theorized for the radar enhanced scattering technique [14], the diagnostic could potentially be used for turbulence studies.
The goal of this work is an investigation of the possibility of using the SPR pulse shape to obtain the turbulence properties. The scattering SPR pulse is analysed theoretically using perturbation theory applied to the Helmholtz equation, similar to analysis performed previously for other reflectometry diagnostics [15,16]. The novelty of the current paper lies in considering pulse reflectometry specifically and obtaining expressions for the scattering signal in the temporal domain. Temporal domain was considered in [14], however the analysis was done in the presence of the upper hybrid resonance rather than a cut-off, qualitatively changing expressions for the scattered wave.
To validate the results of the theoretical analysis, full-wave numerical modelling is employed. This paper is organised as follows. Details of the derivation and final expressions are presented in section 2 (with additional details in appendix). Section 3 covers full-wave numerical modelling carried out to validate the linear theory results. A discussion on the model's limitations, experimental relevance and future prospects is presented in section 4. Lastly, the main points of the paper are summarized in the final section.

Theoretical analysis
The application of perturbation theory to the Helmholtz equation is a well-established approach for analyzing microwave scattering [15][16][17]. Similar to other works, we will consider a slab geometry and ordinary (O) polarization of the probing wave. A Cartesian coordinate system (illustrated in figure 1) will be employed. The background density profile will be assumed to be linear along the radial coordinate x and uniform along the poloidal coordinate y. Variations along the z-direction, aligned with the applied magnetic field, are neglected in this work.
The value L(k y ) in figure 1 designates the radial (x) coordinate of the turning point corresponding to the probing beam's Fourier component with poloidal wavenumber k y . We will consider a monostatic antenna configuration, which is a conventional reflectometry setup.

Basic equations
The Helmholtz equation for a single frequency O-mode can be written as [18]: where ω is the probing frequency, n(x, y) = n bg (x) + δn(x, y) is the plasma electron density and n c = m e ω 2 /(4 πe 2 ) is the critical cut-off density (in cgs units). We are assuming that the background density profile has linear dependence n bg (x) = n c x Lω , where L ω is the density profile gradient scale length at the cut-off (as well as the cut-off position). It should be noted that the value of L ω depends quadratically on the probing frequency, which is why we use the subscript ω.
In the case of a linear density profile, the analytical solution of the unperturbed equation (1) for each poloidal harmonic is expressed in terms of the Airy function Ai(x): (2) Here, α = (L ω c 2 /ω 2 ) 1/3 is the Airy scale and k y is the poloidal wavenumber corresponding to a specific harmonic of the probing wave. Due to quadratic dependence of L ω on the probing frequency, α does not depend on the probing frequency. Next, L ω (k y ) designates the radial position of the turning point for a specific poloidal wavenumber k y . In the case of normal probing, k y = 0, the turning point position coincides with the plasma cut-off L ω (0) = L ω . This solution is normalized to have amplitude A in (ω) of its incident component at the probing antenna (x = 0) and the coefficient in front of the integral is determined by the asymptotic form of the Airy function. We will consider a Gaussian probing pulse, which means the factor A in (ω) takes the following form: In equation (5) ω 0 corresponds to the central probing frequency. t 0 is the pulse launching time and t p is the pulse halfwaist. To further simplify the following derivation, we will assume that the probing antenna generates a Gaussian pattern: The value K corresponds to the central poloidal wavenumber of the antenna pattern and ρ is the half-waist of the beam. In our model, K is related to antenna tilt angle θ (figure 1) and is defined as K ω = ω c sin θ. Using equation (2), we can directly obtain the 0th order approximation corresponding to unperturbed reflected pulse. The details of the derivation are given in appendix. After employing additional assumptions t 2 p >> ρ 2 /c 2 >> L ω0 /(ω 0 c), which seem reasonable for experimental parameters [12], we arrive at the final result for the main beam harmonic k y = K: Here, we introduced a designation L 0 = L ω0 for the cutoff position corresponding to the central probing frequency. Equation (8) for the delay is in agreement with a simple geometrical optics estimate and is appropriate for comparison with the next section.

Linear scattering signal
Within this chapter, we will follow the approach introduced in [15] and employ the same notations as in its previous applications [16,19]. To obtain the first order (linear with regard to density fluctuations) solution of the Helmholtz equation for monostatic configuration in ordinary (O) polarization, the reciprocity theorem [20] can be employed in the following form: Here, A s represents the amplitude of the scattering signal recorded by the receiver antenna. P corresponds to the power of the probing beam and δn Ω is the density perturbation possessing characteristic frequency Ω. This characteristic frequency for realistic parameters is much smaller than the width of the frequency envelope (5) and can be neglected.
This integral for a linear background density profile has been well studied in the literature, and in particular was shown in Ref [15] to have the following form (obtained by utilizing equation (2) and Fourier representation of the Airy function): The integral over k y can be estimated by applying the saddle-point method [16,19]: ] .
Up to this point, our results mostly repeat those available in the literature. Now, adapting the theory to the case of pulse reflectometry, we will return to the temporal domain by utilizing Fourier transform.
we will assume ∆ω << ω 0 and approximate ω ≈ ω 0 in the denominator: If we take a look at the exponential in the integral, particularly the part that depends on ω, we can see that all but one term are the same as in the case of simple reflection. The only new term iL ω (q/2 ) κ depends on the radial wavenumber κ of the turbulence. This term is proportional to the square of the frequency and thus provides a change both in the time delay of the pulse and its width: As for the other phase terms, the expansion into Taylor series over ∆ω is the same as in equation (A.3).
Assuming that the beam is much wider than the poloidal correlation length characterizing the k y spectrum of the turbulence ρ >> l cy , we can once again assume that the integral over q is dominated by the antenna pattern: Terms that do not depend on κ are the same as the ones obtained in 0th order of perturbation theory described in appendix and not connected to the scattering. The term 4i L0 c cos ϑ∆ω is responsible for the delay in equation (8) and corresponds to the time it takes the probing pulse to return after reflecting off the cut-off. The term 2i L0 cω0 ( cos ϑ + 1 cos ϑ ) ∆ω 2 describes dispersive broadening, which can be neglected when the condition t 2 p >> L 0 /(ω 0 c) is fulfilled. The phase terms proportional to κ in the exponent provide a change both in the delay time of the pulse and its width. Carrying out the integration over frequency, we arrive to the final expression for the amplitude of the scattering signal: Here, same as in appendix, we have neglected all the pulse broadening terms (given by (A.5)) not connected to scattering, according to the condition t 2 p >> ρ 2 /c 2 >> L 0 /(ω 0 c). Due to the fact that that |κ| < 2ω/c (for Bragg resonance to be possible), the broadening term depending on κ can also be neglected.
This expression is impossible to interpret further due to the randomness of the density perturbation δn(κ, 2K 0 ). One way to move forward is to consider the statistical average of the previous formula, i.e. by considering the statistical average of the scattering signal power. If we assume that density fluctuations are random and statistically homogeneous so that harmonics with different wavenumbers are independent: There are three distinct parts in the above integral: the radial wavenumber spectrum of the turbulence |δn(κ)| 2 , the timedependent exponential in the numerator and the denominator.
If the denominator can be neglected, the final pulse shape will be determined by the 'broader' of the two remaining terms. This implies that for the case of 4 √ 2L 0 >> l cx w 0 t p , the time dependent exponential in (19) essentially corresponds to a Dirac δ distribution centered at κ = (t d −t)ω0 2L0 and the scattering pulse will reproduce the shape of the turbulence radial wavenumber spectrum: The pulse broadening would be observable even for a less strict condition 4 √ 2L 0 > l cx w 0 t p , but in that case the exact shape of the radial wavenumber spectrum (assuming it is not Gaussian) may not be reproduced perfectly. Such a case will be presented in the next chapter.
However, in most cases, the denominator in the integrand cannot be neglected. This denominator is commonly attributed to poorly localised forward scattering, which provides a disproportionately high contribution to the scattering signal [16]. This term, possessing a singularity, will result in a particular in which case the resulting integral can be roughly evaluated by its value at κ = κ fs : In this scenario, the scattered pulse shape coincides with the reflected pulse, but the delay is slightly different. From the experimental standpoint, equation (20) corresponds to the case of larger probing angles, when small-angle scattering is naturally suppressed, while equation (21) corresponds to smaller probing angles. That includes the case of normal probing, for which κ fs = 0 and equation (21) predicts that the scattered pulse will coincide with the reflected one (equations (A.5) and (7)).

Numerical validation
To confirm the analytical results, we employed full-wave numerical modelling with the GPU-enabled code CUWA [21].

Computational setup
The computations were carried out on a 2D grid, using frozen density perturbations and a linear background density profile in slab (x, y) geometry as shown in figure 1. The source term was adapted for pulsed operation and the scattering signal was recorded throughout all timesteps to obtain the returning pulses.
Similar to previous works [17,22], a Gaussian turbulence spectrum was used to produce random density perturbations: Each sample was obtained by performing the Fourier transform of the spectrum, including a random phase ∆ϕ(κ, q) sampled uniformly over [0, 2π]. A separate full-wave computation was carried out for each sample and the resulting average scattering pulse power was obtained by averaging over 1000 random samples. This amount of samples was found to be sufficient for the average pulse shapes to converge (e.g. the difference between average pulses computed with 500 and 1000 samples was about 2%). An example of a CUWA computation is given in figure 2.
The emitting antenna was simulated at the computation boundary (x = 0) with a Gaussian pattern. The emitted power was modulated to produce a Gaussian pulse described by equation (4) with t 0 = 3 t p . This offset was used to avoid fast variations of the emitted power that would generate additional frequencies in the probing field's frequency spectrum. To obtain the signal from the computation a receiving antenna was simulated by integrating the electric field at the grid boundary (corresponding to antenna position) with antenna pattern. Due to the fact that we are interested specifically in the scattering pulse, the reflected pulse (obtained from simulation with no fluctuations) was numerically removed from the computations by subtracting it from the recorded complex signals. Accounting for the reflected signal in this way was important for the computations at low probing angles presented later in this section. Finally, to ensure that the linear scattering theory was applicable, the amplitude of density perturbations was kept at r.m.s.(δn) = 2 × 10 −4 n c , which is expected to be low enough to not enter the nonlinear scattering regime  [17,23]. Additionally, computations at different turbulence amplitudes were carried out to confirm the linearity of the scattering signal. An example of such computations will be presented in section 4.

Pulse broadening
First, to verify the concept of the measurement of the radial wavenumber of the turbulence, we selected a set of parameters which would satisfy both the condition of pulse broadening 4 √ 2L 0 > l cx w 0 t p and the condition for small-angle scattering suppression given in [16]. The parameters selected were f = 50 GHz, L 0 = 40 cm, ρ = 3 cm, ϑ = 35 • , t p = 0.8 ns, l cx = 0.5 cm and l cy = 0.6 cm. In this case, according to equation (21) we would expect the average (over all samples) scattering pulse to reproduce the radial wavenumber spectrum of the turbulence. The result of the full-wave computation is given in figure 3.
The full-wave computation demonstrates good agreement with analytical expectations (given by equation (20)) and reproduces the Gaussian turbulence spectrum. Small differences in delay (1.6%) and pulse waist (2%) are likely caused by small-angle scattering. Despite being mostly suppressed, small angle scattering will still make radial wavenumbers closer to κ fs provide slightly larger contribution into the scattering signal resulting in a shift. This effect is much more pronounced when small-angle scattering dominates the scattering signal. Overall, the full-wave modelling confirms the possibility to determine the turbulence radial wavenumber spectrum from the scattered pulse. Next, a computation for the same set of parameters but a different radial wavenumber spectrum was carried out. The spectrum used in this case has two distinct parts-constant level for large scales and power law for smaller scales: In this case, the pulse does not perfectly reproduce the shape of the radial wavenumber spectrum. This is explained by the fact that for this type of spectrum, small-angle scattering plays a significant role. Due to the slower decay of the spectrum over radial wavenumbers, the criterion introduced in [16] is no longer sufficient to eliminate the influence of the small-angle scattering.
The pulse presented in figure 4 corresponds to an intermediate case with regards to the small-angle scattering, which has enough of an effect to cause asymmetry of the pulse in the vicinity of κ fs = − , but is not strong enough to dominate the scattering signal completely and results in the scattering signal simply reproducing the shape of the input pulse [given by equation (21)].

Normal probing
To study the case of small-angle scattering dominated signal, modelling was carried out for the case of normal probing. This case should be well described by equation (21). In the case of normal probing (standard SPR measurement), we can expect the scattering signal to reproduce the probing pulse and to have the same delay as the reflected pulse. The parameters selected for this computation were f = 30 GHz, L = 40 cm, ρ = 8 cm, ϑ = 0 • , t p = 0.9 ns, l cx = 0.1 cm and l cy = 2 cm. In this case, the following radial wavenumber spectrum was selected: The small value of l cx was selected to guarantee that the broadening condition 4 √ 2L 0 >> l cx w 0 t p was fulfilled. That way, the absence of broadening could be attributed directly to the small-angle scattering. The resulting computation is presented in figure 5.
From this figure, it is clear that the significant broadening predicted by equation (20) (red curve) is not realised and the scattering signal reproduces the probing pulse given by equation (21) (yellow curve) instead. It should also be noted that in this case the pulse delay is exactly 4L0 c , the same as for the reflected pulse. This means that in the case of standard SPR configuration, linear scattering would not affect the shape or delay of the pulses, making the diagnostic robust to linear scattering effects.
Further considering normal probing, a computation with a radial wavenumber spectrum designed to compensate for the effects of small-angle scattering (all other parameters unchanged) was carried out to reinforce the conclusion about the influence of small-angle scattering on the pulse shape. The spectrum had the following shape: This spectrum was expected to counteract the denominator in equation (19) (which in the case of normal probing takes the simple form 1 |κ| ), resulting in a very broad scattering pulse, similar to the red curve in figure 5. The results of the full-wave computation for this case are presented in figure 6.
While the computed pulse indeed demonstrates significant broadening, there is a noticeable gap at its centre (corresponding to the scattering signal for κ = 0). This gap is explained by the fact that the scattering efficiency, given by the equation (11), is only an approximation and does not correctly describe the special point κ fs = 0 in the case of normal probing. The saturation of the scattering efficiency for large scales [15,24] is not accounted for in our model, which results in the spectrum of equation (25) 'over-correcting' and suppressing the scattering on the turbulence with the largest scale. Nevertheless, figure 6 demonstrates that small-angle scattering is the effect behind the absence of the pulse broadening in figure 5.

Model limitations
Before discussing the obtained results and their implications, it would be useful to outline a few factors limiting the applicability of the results.
First of all, O-mode was considered within the presented paper. In the case of X-mode, the wave vector has a more complicated dependence on the plasma electron density. It is, however, possible to account for the difference by introducing an additional factor in equation (11), similar to an approach presented in [23]. Such an approach would make it possible to generalize the method to X-mode polarization.
A second, the simplification we used is the assumption t 2 p >> ρ 2 /c 2 >> L 0 /(ω 0 c), which allowed us to neglect the broadening of the pulse not associated with the scattering [25]. This broadening was partially described in equation (A.5), but is neglected in the rest of the analytical derivation. A potential analytical generalization is possible, but it was not explored within this paper since the condition t 2 p >> ρ 2 /c 2 >> L 0 /(ω 0 c) is reasonable for the relevant experimental parameters [12].
One more unrealistic assumption of the model is the linear density profile, which made it possible to use the Airy function as the solution for the unperturbed Helmholtz equation. However, while the pulse delay will change when a realistic profile is used, the shape of the pulse and the delay variations are determined by the region in immediate vicinity of the turning point. Thus, the results can easily be generalized to the case of an arbitrary density profile by interpreting L 0 as the density gradient scale length at the turning point.
Another simplification is the use of the slab geometry and thus discarding the curvature of magnetic surfaces. This approach is justified when considering large magnetic confinement devices (for which the curvature radius of magnetic surface is larger than L 0 ), but in the case of smaller machines the geometry effects might play a significant role. Nevertheless, previous linear studies in cylindrical geometry [26] have shown no qualitative difference compared to slab geometry. Thus, one can expect the technique to still be applicable in smaller devices.
Finally, the most important simplification of the model is the use of linear scattering theory. To validate it, one has to justify both the exclusion of the zeroth order reflected part of the signal and higher order nonlinear scattering effects.
In the case of SPR, separating reflected and scattering signals is a non-trivial task due to the fact that measurements do not necessarily contain information about the frequency composition of the scattering signal. However, as shown in the previous sections, the turbulence radial wavenumber spectrum measurements can only be carried out at oblique incidence of the probing beam with respect to magnetic surface, which naturally eliminates the reflected part of the signal from the measurements in monostatic configuration.
When it comes to nonlinear scattering effects, the turbulence amplitudes at which they can be neglected are well-defined in the literature [17,22,23,27]. Unfortunately, benchmarking of experimental Doppler reflectometry measurements with numerical modelling results [28,29] suggests that it is quite common for the back-scattering diagnostics to operate in the nonlinear scattering regime. While our theoretical model does not cover the nonlinear regime, we have carried out full-wave computations, which make no assumptions with regards to density fluctuations' amplitude. Thus, our computations are able provide information regarding nonlinear scattering regime. In the following figures we show the results of full-wave modelling for different r.m.s. values of density perturbation (the same fluctuation samples multiplied by a constant value were used). The parameters for this case were chosen to be computationally inexpensive (due to smaller grid): f = 50 GHz, L = 20 cm, ρ = 8 cm, ϑ = 0 • , t p = 0.9 ns, l cx = 0.5 cm and l cy = 2 cm. For these parameters the computation was carried out for a set of density perturbation amplitudes, and the average scattered pulse amplitude (corresponding to square root of pulse power) is presented in figure 7.
In the linear scattering regime, the pulse amplitude should depend linearly on the amplitude of the input density perturbations. This was verified to be the case for all the results presented in the previous section. In the case plotted in figure 7 however, only the first five points demonstrate a linear dependence. The following points related to higher fluctuation amplitude demonstrate a lower order dependence, which corresponds to a saturated nonlinear regime [23], which is expected for normal probing [17]. While in the linear regime there is no change to the pulse shape or delay in the case of normal probing, that is not necessarily true in the nonlinear regime. To investigate this point, the average signal powers (including both reflected and scattered parts) corresponding to different amplitudes are presented in figure 8.
It can be seen from the figure that as the scattering regime becomes more nonlinear, the pulse delay decreases and the pulse width increases. The decrease of delay is expected, as in the nonlinear regime the probing pulse power is partially exhausted before it reaches cut-off. As a result the front slope of the scattered pulse will have higher amplitude than the back one, effectively decreasing the delay. Additionally, in a strongly nonlinear regime, the cut-off position will on average be shifted outwards due to density perturbations.
As for the pulse broadening, there is no clear explanation for this effect, but it was found to be consistent (in the case of normal probing) among all the computations in the nonlinear regime. Thus, the change of the pulse shape with respect to the probing pulse in the case of normal probing can be a marker of the nonlinear scattering regime in SPR measurements.

Discussion
Overall, the full-wave computations are in agreement with the conclusions drawn in section 2. It seems that the SPR diagnostic in its 'Doppler' (non-normal probing) configuration could be used to measure the radial wavenumber spectrum of the turbulence. However, several conditions need to be fulfilled for that [16]: Pulse broadening ρ >> l cy ; Antenna pattern selecting precisely q = 2K 0 While the condition on the beam waist usually holds true in experiment, the other two conditions are less straightforward. The pulse broadening condition requires the plasma to have a relatively flat density profile and thus could only be used in particular discharges and at particular radial positions.
Finally, the small-angle scattering suppression condition requires a particular set of probing parameters, but with enough flexibility of the probing angle and narrow enough beam, it would be possible to fulfill.
All three conditions require one to have a prior estimate of the correlation lengths l cx/cy , which is not necessarily possible. One advantage of this method is that if the optimal conditions are not fulfilled, it would be immediately obvious from the fact that the shape of the measured scattered pulse would coincide with that of the probing pulse.
We can use the TCV SPR system [12,13] to give a quantitative example. The possible probing beam parameters for the TCV SPR are f = 50 GHz, ρ ≈ 4 cm and t p =0.7 ns. We will assume that the radial and poloidal correlation lengths of the turbulence are of the order of 1 cm. This means that the broadening condition in this case requires us to have L 0 > 39 cm. The L-mode profile given in [13] has values of L 0 varying between 10 and 20 cm, making it not suitable for such measurements. The H-mode profile on the other hand can reach L 0 of 35 cm, which is close to fulfilling the broadening condition. Using an even shorter pulse t p = 0.5 ns (which is possible with the current hardware), would lead to the broadening condition L 0 > 28, making the measurement in the H-mode pedestal top viable.
The condition ρ >> l cy is fulfilled without significantly limiting the measurements.
Finally, for L 0 = 35 cm, the small-angle suppression criterion gives us the condition (sin ϑ) 2 > 0.12 resulting in a required angle of 21 • , which is well within the capabilities of the diagnostic.
As can be seen from this example, very specific experimental parameters, mainly dictated by the pulse broadening condition, are necessary for the method to be viable. Nevertheless, if using lower probing frequencies and shorter probing pulses becomes possible, the method could be used in a wider range of density profiles.
As far as the available alternatives to SPR go, the most similar diagnostics seem to be radial correlation reflectometry [2,3] and radial correlation Doppler reflectometry (RCDR) [30]. These diagnostics also suffer from the influence of the small-angle scattering, but a mathematical method was suggested to account for it [19,31]. A similar method could be derived for the SPR pulse analysis, as part of the future theoretical development. While the potential Doppler SPR shares some disadvantages of RCDR, its main appeal is the simplicity of the setup. Instead of having to probe plasma with multiple continuous frequencies (possibly over multiple discharges), the measurement is done with a single probing frequency within one discharge. Moreover, with a modern setup, such as the one used at the TCV tokamak [12], it is possible to perform the measurements at several probing frequencies, allowing one to compare the radial wavenumber spectra corresponding to several radial positions.
Next, let us still point out that the numerical validation performed within the paper only confirms the basic feasibility of the method without taking into account the full complexity of the experimental setup. This work is meant to introduce the concept of the Doppler SPR approach to turbulence studies. Further validation for the case of realistic geometry and physical turbulence spectrum will be the subject of future studies, including potential experimental implementation.
Another subject for future work is the possible use of additional information provided by the SPR for studying plasma turbulence, such as the statistical properties of pulse delays and their correlation. While this paper focuses specifically on the pulse shape, other characteristics could also provide insight into amplitude and correlation lengths of turbulent density profile perturbations.

Conclusions
Within this paper, perturbation theory was applied to the Helmholtz equation to obtain the linear scattering signal for the conditions of short pulse reflectometry. The particular interest lay in investigating the shape of the scattering pulse and its delay. Expressions describing both have been derived for an arbitrary probing angle.
The potential of the diagnostic to measure radial wavenumber spectrum of the turbulence was demonstrated and the necessary conditions, in particular the use of a 'Doppler' configuration, were outlined. Additionally, robustness of the classic SPR to scattering effects was demonstrated.
Analytical results have been validated with full-wave modelling. The modelling was also employed to illustrate and quantify nonlinear scattering effects and their impact on the SPR measurements. A potential experimental indicator for having reached the nonlinear scattering regime was suggested. Finally, possible directions of future studies were outlined.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix. Delay of the reflected pulse
The complex signal recorded by the antenna, following the reciprocity theorem, can be expressed as the integral of the electric field with the one that would be generated by the receiver antenna: To obtain this expression, we have used the asymptotic limit of the Airy function. Separating the sine function into two exponential terms (Euler's formula), one immediately notices that one of them corresponds to the original probing pulse. The second term corresponds to the returning pulse and the integral for it takes the form: dω.

(A.2)
Accounting for the quadratic dependence of L ω over probing frequency as well as the fact that the pulse is centred around the carrying frequency (∆ω << ω 0 ), we can expand the last phase term into Taylor series: Here, we designated by L 0 the position of the turning point corresponding to the central probing frequency ω 0 . To estimate the integral over k y , we can assume that the antenna pattern dominates the integral, i.e. ρ 2 >> l 0 c/ω, corresponding to the near field of the probing antenna. Under such conditions, the integral over poloidal wavenumbers is determined entirely by the antenna patterns f(k y ) f(−k y ), which select k y = 0 and we get: The opposite case of antenna far field can be considered by employing the paraxial approximation (as in [16]). It will not result in a qualitatively different conclusion and we will not consider it within this work. Integrating over frequencies and multiplying by the complex conjugate to obtain power, we get the final expression for the reflected pulse and its delay within our model: (A.5) Here, we have neglected smaller exponential terms according to the condition ρ 2 >> l 0 c/ω. Given the typical experimental parameters [12], we can also assume t p > ρ/c and neglect the modification of the reflected pulse width. A detailed description of the reflected pulse width modification can be found in [25].
When probing under an oblique angle, the reflected signal is significantly attenuated by the antenna pattern, which corresponds to the factor K 2 ρ 2 in the exponential of equation (A.5). If the receiving antenna was instead aligned to receive the main harmonic (k y = K) of the probing beam, the resulting pulse would instead have the form: The conditions t 2 p >> ρ 2 /c 2 >> L 0 /(ω 0 c) that we used are not strictly necessary to obtain the expression for the reflected or the linear scattering pulses. However, they allow us to neglect all the pulse broadening not associated with scattering and they generally hold true in experiment. Without using these conditions one would have to use a cumbersome expression similar to equation (A.5) and account for the modification of the reflected pulse compared to the probing one [25].