Hostname: page-component-74656f9884-2dmjx Total loading time: 0 Render date: 2026-11-10T17:13:08.292Z Has data issue: false hasContentIssue false
  • English
  • Français

Exact boundary-value solution for an electromagnetic wave propagating in a linearly varying index of refraction

Part of: Invited Contributions from the 20th European Fusion Theory Conference

Published online by Cambridge University Press:  12 November 2024

N.A. Lopez Open the ORCID record for N.A. Lopez [Opens in a new window]
N.A. Lopez*
Affiliation:
Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU, UK
*
†Email address for correspondence: nicolas.lopez@tokamakenergy.com
Rights & Permissions [Opens in a new window]

Abstract

The propagation of electromagnetic waves in a linearly varying index of refraction is a fundamental problem in wave physics, being relevant in fusion science for describing certain wave-based heating and diagnostic schemes. Here, an exact solution is obtained for a given incoming wavefield specified on the boundary transverse to the direction of inhomogeneity by performing a spectral, rather than asymptotic, matching. Two case studies are then presented: a Gaussian beam at oblique incidence and a speckled wavefield at normal incidence. For the Gaussian beam, it is shown that when the waist $W$ is sufficiently large, oblique incidence manifests simply as rigid translation and focal shift of the corresponding diffraction pattern at normal incidence. The destruction of the hyperbolic umbilic caustic (corresponding to a critically focused beam) as $W$ is reduced is then demonstrated. The caustic disappears once $W \lesssim \delta _a \sqrt {L}$ ($L$ being the medium length scale normalized by the Airy skin depth $\delta _a$), at which point the wave behaviour is increasingly described by Airy functions, but experiences less focusing as a result. To maximize the intensity of a launched Gaussian beam at a turning point, one should therefore minimize the imaginary part of the launched complex beam parameter while having the real part satisfy critical focusing. For the speckled wavefield, it is shown that the transverse speckle pattern only couples to the Airy longitudinal pattern when the coupling parameter $\eta = \sqrt {L}/f_{\#}$ is large, with $f_{\#}$ being the f-number of the launching aperture. When $\eta \ll 1$, a reduced description of the total wavefield can be obtained by simply multiplying the incoming speckle pattern with the Airy swelling.

Keywords

fusion plasmaplasma waves

Information

Type
Research Article
Information
Journal of Plasma Physics , Volume 90 , Issue 6 , December 2024 , 995900602
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
Copyright © The Author(s), 2024. Published by Cambridge University Press

1 Introduction

The description of an electromagnetic wave propagating in a linearly varying medium is a classic problem in wave physics, sometimes referred to as the linear-layer problem (Ginzburg Reference Ginzburg1961). Its relevance to plasma physics and fusion research is predominantly as a local description near a turning point for wave-based heating and diagnostics applications. For example, a recently developed scheme to perform fundamental X-mode heating and current drive in startup plasmas (Ono, Bertelli & Shevchenko Reference Ono, Bertelli and Shevchenko2022a; Ono et al. Reference Ono, Bertelli, Shevchenko, Idei and Hanada2022b) involves an obliquely launched X-mode resonating with electrons near the low-density cut-off; the efficiency of the heating process will therefore depend on the wavefield intensity behaviour near the cut-off. Similarly, the Doppler backscattering diagnostic (Hall-Chen, Parra & Hillesheim Reference Hall-Chen, Parra and Hillesheim2022) relies on the swelling of a wavefield near a turning point to probe turbulence spectra via nonlinear scattering of the diagnostic beam; the intensity details near the turning point therefore determine the diagnostic sensitivity and localization.

Many authors (Ginzburg Reference Ginzburg1961; Orlov & Tropkin Reference Orlov and Tropkin1980; Maj, Pereverzev & Poli Reference Maj, Pereverzev and Poli2009; Maj, Balakin & Poli Reference Maj, Balakin and Poli2010; Lopez, Kur & Strozzi Reference Lopez, Kur and Strozzi2023) have obtained solutions to the linear-layer problem in terms of a prescribed field $\psi$ along the boundary $z = 0$, where $z$ is the direction of medium inhomogeneity. However, the full field $\psi$ is made up of both the incoming and reflected components, and often only the incoming component is known in practice. To remedy this apparent shortcoming in the obtained solutions, most of the aforementioned authors performed an asymptotic matching to isolate only the incoming component. This resolves the issue with the boundary condition, but at the cost of introducing an asymptotic validity criterion into the analysis: the resulting formulas are not valid when the incoming field is specified too close to the turning point, as might occur when beams are launched obliquely. A universally valid formula that matches to a prescribed incoming field would be more desirable for use in applications.

Here we obtain such a formula by performing the matching spectrally instead of asymptotically. No new asymptotic validity criteria are introduced into the problem, and it is shown explicitly that the obtained solution exactly reproduces the incoming wavefield at the boundary, regardless of the distance between the boundary and the turning point. The solution involves an integral whose kernel contains the standard Airy function, which is expected to arise in this problem, and also the related Scorer function, which is less well known. To demonstrate the flexibility of the new solution, two special cases are considered: an obliquely launched Gaussian beam and a normally incident speckled wavefield produced by a bilevel random phase plate (RPP). For this first example, it is shown how the hyperbolic umbilic caustic created at critical focusing gets degraded as the beam waist becomes smaller for a variety of injection angles, including angles at which the asymptotic validity criterion for the previous solutions is violated. For the second example, an explicit coupling parameter is derived and demonstrated that governs whether the speckles influence the behaviour of the wavefield near the turning point. Both these examples can serve as starting points to developing reduced models of waves near turning points with more comprehensive physics content, which would be useful for the fusion applications mentioned in the first paragraph, among other applications.

This paper is organized as follows. In § 2 the linear-layer problem is introduced. In § 3 a spectral matching is performed to allow the solution to the linear-layer problem to be expressed only in terms of the incoming field at the boundary. This is the main result of this paper. In § 4 the special case of an incoming Gaussian beam is studied, with particular emphasis on its behaviour near critical focusing. In § 5 the special case of an incoming speckled wavefield is studied, with focus on characterizing the coupling between speckles and the Airy pattern. Finally, § 6 summarizes the main conclusions. Auxiliary calculations are presented in appendices.

2 Background

Let us consider an electromagnetic wave propagating in $N + 1$ spatial dimensions in a medium whose index of refraction varies as a linear function. We take one dimension, denoted $z$, to be aligned with the direction of inhomogeneity, and label the remaining $N$ dimensions by the vector ${\boldsymbol {x}}$. Assuming time-harmonic modes with a single angular frequency $\omega$, the wavefield amplitude can be shown to satisfy the Helmholtz equation:

(2.1)\begin{equation} \partial_{{\boldsymbol{x}}}^2 \psi({\boldsymbol{x}},z) + \partial_{z}^2 \psi({\boldsymbol{x}},z) + \frac{\ell - z}{\delta_a^3} \psi({\boldsymbol{x}},z) = 0, \end{equation}

where $\delta _a$ is a constant with units of length sometimes called the ‘Airy skin depth’ (Michel Reference Michel2023). In terms of the angular frequency $\omega$, the medium length scale $\ell$ and the speed of light in vacuum $c$, $\delta _a$ is given as

(2.2)\begin{equation} \delta_a = \sqrt[3]{ \frac{\ell c^2 }{\omega^2} }. \end{equation}

We shall maintain $N$ unspecified in the following analysis, but note that practical calculations will have either $N = 1$ or $N = 2$ for two- or three-dimensional propagation, respectively. We also only seek solutions that are stable, such that they become evanescent and decay to zero as $z \to +\infty$.

Let us now introduce normalized spatial coordinates

(2.3a)\begin{equation} {\boldsymbol{x}} = \delta_a {\boldsymbol{X}} , \quad z = \delta_a Z , \quad \ell = \delta_a L , \end{equation}

such that (2.1) becomes

(2.4)\begin{equation} \partial_{{\boldsymbol{X}}}^2 \psi({\boldsymbol{X}},Z) + \partial_{Z}^2 \psi({\boldsymbol{X}},Z) + (L - Z) \psi({\boldsymbol{X}},Z) = 0 . \end{equation}

Let us also adopt the following convention for the Fourier transform (FT):

(2.5a)\begin{gather} \smash{\widetilde{\psi}}({\boldsymbol{K}}_x,K_z) = \int \frac{\mathrm{d} {\boldsymbol{X}} \, \mathrm{d} Z}{(2{\rm \pi})^{N + 1}} \psi({\boldsymbol{X}},Z) \exp({- {\rm i} {\boldsymbol{K}}_x \boldsymbol{\cdot} {\boldsymbol{X}} - {\rm i} K_z Z}), \end{gather}
(2.5b)\begin{gather}\psi({\boldsymbol{X}},Z) = \int \,\mathrm{d} {\boldsymbol{K}}_x \, \mathrm{d} K_z \, \smash{\widetilde{\psi}}({\boldsymbol{K}}_x,K_z) \exp({{\rm i} {\boldsymbol{K}}_x \boldsymbol{\cdot} {\boldsymbol{X}} + {\rm i} K_z Z}) . \end{gather}

Taking the FT of (2.4) then gives

(2.6)\begin{equation} {\rm i} \partial_{K_z} \smash{\widetilde{\psi}}({\boldsymbol{K}}_x,K_z) = \left( L - |{\boldsymbol{K}}_x|^2 - K_z^2 \right) \smash{\widetilde{\psi}}({\boldsymbol{K}}_x,K_z) , \end{equation}

with solution given by

(2.7)\begin{equation} \smash{\widetilde{\psi}}({\boldsymbol{K}}_x,K_z) = \smash{\widetilde{\psi}}_0({\boldsymbol{K}}_x) \exp\left[ {\rm i} \frac{K_z^3}{3} + {\rm i} \left( |{\boldsymbol{K}}_x|^2 - L \right) K_z \right]\!. \end{equation}

Here, $\smash {\widetilde {\psi }}_0({\boldsymbol {K}}_x) \equiv \smash {\widetilde {\psi }}(0,{\boldsymbol {K}}_x)$ is an arbitrary function that can eventually be matched to boundary conditions, as we show in the next section.

The general solution to (2.4) is then obtained by taking an inverse FT of (2.7):

(2.8)\begin{equation} \psi({\boldsymbol{X}},Z) = \int \,\mathrm{d} {\boldsymbol{K}}_x \, \mathrm{d} K_z \, \smash{\widetilde{\psi}}_0({\boldsymbol{K}}_x) \exp\left[ {\rm i} \frac{K_z^3}{3} + {\rm i} \left( |{\boldsymbol{K}}_x|^2 + Z - L \right) K_z + {\rm i} {\boldsymbol{K}}_x \boldsymbol{\cdot} {\boldsymbol{X}} \right]\!. \end{equation}

Note that the decaying boundary condition has been tacitly imposed through our use of an FT to obtain the solution (2.8).

3 Boundary-value solution for prescribed incoming wavefield

3.1 Isolating the incoming contribution

Although (2.8) constitutes the general boundary-value solution for a prescribed $\smash {\widetilde {\psi }}_0({\boldsymbol {K}}_x)$, this boundary-value solution is formulated in the spectral domain, which is not useful for most applications. Here we instead obtain the boundary-value solution in the coordinate domain by prescribing an incoming wavefield, since this is often what is known in practice.

To do this, consider the field on the $Z = 0$ plane, i.e. $\psi ({\boldsymbol {X}},0)$. We can split this into ‘incoming’ and ‘outgoing’ components by using a spectral filter based on the sign of $K_z$ as follows:

(3.1a)\begin{gather} \psi_\text{in}({\boldsymbol{X}}) = \int \,\mathrm{d} {\boldsymbol{K}}_x \, \smash{\widetilde{\psi}}_0({\boldsymbol{K}}_x) {\rm e}^{{\rm i} {\boldsymbol{K}}_x \boldsymbol{\cdot} {\boldsymbol{X}}} \int_0^\infty \,\mathrm{d} K_z \exp\left[ {\rm i} \frac{K_z^3}{3} + {\rm i} \left( |{\boldsymbol{K}}_x|^2 - L \right) K_z \right]\!, \end{gather}
(3.1b)\begin{gather}\psi_\text{out}({\boldsymbol{X}}) = \int \,\mathrm{d} {\boldsymbol{K}}_x \, \smash{\widetilde{\psi}}_0({\boldsymbol{K}}_x) {\rm e}^{{\rm i} {\boldsymbol{K}}_x \boldsymbol{\cdot} {\boldsymbol{X}}} \int_{-\infty}^0 \,\mathrm{d} K_z \exp\left[ {\rm i} \frac{K_z^3}{3} + {\rm i} \left( |{\boldsymbol{K}}_x|^2 - L \right) K_z \right]\!. \end{gather}

One then has the exact decomposition

(3.2)\begin{equation} \psi({\boldsymbol{X}},0) = \psi_\text{in}({\boldsymbol{X}}) + \psi_\text{out}({\boldsymbol{X}}) . \end{equation}

To proceed further, one is required to compute the integral

(3.3)\begin{align} I(\zeta) & = \int_0^\infty \,\mathrm{d} K_z \exp\left( {\rm i} \frac{K_z^3}{3} + {\rm i} \zeta K_z \right)\nonumber\\ & \equiv \int_0^\infty \,\mathrm{d} K_z \cos\left( \frac{K_z^3}{3} + \zeta K_z \right) + {\rm i} \int_0^\infty \,\mathrm{d} K_z \sin\left( \frac{K_z^3}{3} + \zeta K_z \right)\!, \end{align}

where $\zeta$ is a real-valued parameter. Note that the other relevant integral is

(3.4)\begin{equation} \int_{-\infty}^0 \,\mathrm{d} K_z \exp\left[ {\rm i} \frac{K_z^3}{3} + {\rm i} \zeta K_z \right] = \left[ I(\zeta) \right]^*\!. \end{equation}

Both the integrals involved in the real and imaginary parts of (3.3) can be solved in terms of Airy-related functions (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010):

(3.5a)\begin{equation} \int_0^\infty \,\mathrm{d} K_z \cos \left( \frac{K_z^3}{3} + \zeta K_z \right) = {\rm \pi}\operatorname{Ai}(\zeta) , \quad \int_0^\infty \,\mathrm{d} K_z \sin \left( \frac{K_z^3}{3} + \zeta K_z \right) = {\rm \pi}\operatorname{Gi}(\zeta), \end{equation}

where $\operatorname {Ai}$ denotes the Airy function and $\operatorname {Gi}$ denotes the Scorer function. These functions are plotted in figure 1 for reference.

Figure 1.

The Airy function $\operatorname {Ai}(x)$ and the Scorer function $\operatorname {Gi}(x)$ plotted against their argument. Both functions behave qualitatively similar, exhibiting exponential decay for $x > 0$ and oscillatory behaviour (with relative phase shift) for $x < 0$.

One therefore obtains

(3.6)\begin{equation} I(\zeta) = {\rm \pi}\left[ \operatorname{Ai}\left( \zeta \right) + {\rm i} \operatorname{Gi}\left( \zeta \right) \right] . \end{equation}

This implies that the incoming and outgoing wavefields can be expressed as

(3.7)\begin{equation} \psi_{\text{in}, \text{out}}({\boldsymbol{X}}) = {\rm \pi}\int \,\mathrm{d} {\boldsymbol{K}}_x \, \smash{\widetilde{\psi}}_0({\boldsymbol{K}}_x) {\rm e}^{{\rm i} {\boldsymbol{K}}_x \boldsymbol{\cdot} {\boldsymbol{X}}} \left[ \operatorname{Ai}\left( |{\boldsymbol{K}}_x|^2 - L \right) \pm {\rm i} \operatorname{Gi}\left( |{\boldsymbol{K}}_x|^2 - L \right) \right] , \end{equation}

where the top $(+)$ sign corresponds to the incoming wavefield and the bottom $(-)$ sign corresponds to the outgoing wavefield. Again, we emphasize that these are exact relationships.

3.2 Fourier inversion to obtain general solution

Let us now obtain the desired boundary-value solution (in coordinate space) using (3.7). To do so, let us perform an FT with respect to only the transverse coordinates ${\boldsymbol {X}}$. Analogous to (2.5), this transverse FT takes the form

(3.8a)\begin{gather} \hat{\psi}_\text{in}({\boldsymbol{K}}_x) = \int \frac{\mathrm{d} {\boldsymbol{X}}}{(2{\rm \pi})^N } \psi_\text{in}({\boldsymbol{X}}) {\rm e}^{- {\rm i} {\boldsymbol{K}}_x \boldsymbol{\cdot} {\boldsymbol{X}}}, \end{gather}
(3.8b)\begin{gather}\psi_\text{in}({\boldsymbol{X}}) = \int \,\mathrm{d} {\boldsymbol{K}}_x \, \hat{\psi}_\text{in}({\boldsymbol{K}}_x) {\rm e}^{{\rm i} {\boldsymbol{K}}_x \boldsymbol{\cdot} {\boldsymbol{X}}} . \end{gather}

Then, we can clearly identify from (3.7) the relationship between the transverse FT and the standard FT images of $\psi$:

(3.9)\begin{equation} \hat{\psi}_\text{in}({\boldsymbol{K}}_x) = {\rm \pi}\smash{\widetilde{\psi}}_0({\boldsymbol{K}}_x) \left[ \operatorname{Ai}\left( |{\boldsymbol{K}}_x|^2 - L \right) + {\rm i} \operatorname{Gi}\left( |{\boldsymbol{K}}_x|^2 - L \right) \right] . \end{equation}

Since the quantity $\operatorname {Ai}( |{\boldsymbol {K}}_x|^2 - L ) + {\rm i} \operatorname {Gi}( |{\boldsymbol {K}}_x|^2 - L )$ is always non-zero (see Appendix A), we can then invert the relationship (3.9) to obtain

(3.10)\begin{equation} \smash{\widetilde{\psi}}_0({\boldsymbol{K}}_x) = \frac{1}{\rm \pi} \frac{ \hat{\psi}_\text{in}({\boldsymbol{K}}_x) }{ \operatorname{Ai}\left( |{\boldsymbol{K}}_x|^2 - L \right) + {\rm i} \operatorname{Gi}\left( |{\boldsymbol{K}}_x|^2 - L \right) } . \end{equation}

The general solution can then be written as

(3.11)\begin{equation} \psi({\boldsymbol{X}},Z) = \int \,\mathrm{d} {\boldsymbol{K}}_x \, \frac{ 2 \operatorname{Ai}\left( |{\boldsymbol{K}}_x|^2 + Z - L \right) \hat{\psi}_\text{in}({\boldsymbol{K}}_x) }{ \operatorname{Ai}\left( |{\boldsymbol{K}}_x|^2 - L \right) + {\rm i} \operatorname{Gi}\left( |{\boldsymbol{K}}_x|^2 - L \right) } {\rm e}^{{\rm i} {\boldsymbol{K}}_x \boldsymbol{\cdot} {\boldsymbol{X}}} , \end{equation}

where $\hat {\psi }_\text {in}({\boldsymbol {K}}_x)$ is the spectrum of the incoming beam, related to the prescribed boundary value via (3.8a). In essence, we have determined the arbitrary function $\smash {\widetilde {\psi }}_0({\boldsymbol {K}}_x)$ in (2.8) that matches a prescribed boundary condition $\psi _\text {in}({\boldsymbol {X}})$ exactly, without appealing to asymptotic approximations (Orlov & Tropkin Reference Orlov and Tropkin1980; Maj et al. Reference Maj, Pereverzev and Poli2009; Lopez et al. Reference Lopez, Kur and Strozzi2023) that necessarily restrict the validity of the resulting expressions.

As a sanity check, one can confirm that (3.11) reproduces the known result when $\psi _\text {in}({\boldsymbol {X}})$ is a constant; in this case, $\hat {\psi }_\text {in}({\boldsymbol {K}}_x) \propto \delta ({\boldsymbol {K}}_x)$ such that subsequent integration gives $\psi ({\boldsymbol {X}}, Z) \propto \operatorname {Ai}(Z - L)$ as desired. Also, note that the incoming and outgoing components to (3.11) at $Z = 0$ can be identified by performing the rearrangement

(3.12)\begin{equation} \frac{ 2 \operatorname{Ai}\left( |{\boldsymbol{K}}_x|^2 - L \right) }{ \operatorname{Ai}\left( |{\boldsymbol{K}}_x|^2 - L \right) + {\rm i} \operatorname{Gi}\left( |{\boldsymbol{K}}_x|^2 - L \right) } = 1 + \frac{ \operatorname{Ai}\left( |{\boldsymbol{K}}_x|^2 - L \right) - {\rm i} \operatorname{Gi}\left( |{\boldsymbol{K}}_x|^2 - L \right) }{ \operatorname{Ai}\left( |{\boldsymbol{K}}_x|^2 - L \right) + {\rm i} \operatorname{Gi}\left( |{\boldsymbol{K}}_x|^2 - L \right) } . \end{equation}

The incoming or outgoing component corresponds to the subsequent integration of the first or second factor, respectively. In particular, one recovers (3.8b) exactly.

4 Special case: incident Gaussian beam in two dimensions

Let us now consider the case when the incoming wavefield is a Gaussian beam, which is of practical importance. We also specialize to only consider two-dimensional propagation ($N = 1$), since this will facilitate comparisons with other published formulas in the literature (Orlov & Tropkin Reference Orlov and Tropkin1980; Maj et al. Reference Maj, Pereverzev and Poli2009; Lopez et al. Reference Lopez, Kur and Strozzi2023). Specifically, we take

(4.1)\begin{equation} \psi_\text{in}(X) = E_0 \exp\left( {\rm i} \sqrt{L} \, X \sin \theta - {\rm i} \frac{X^2 \cos^2 \theta}{2 \sqrt{L}\, q_c} \right)\!, \end{equation}

where $E_0$ is a constant and $q_c$ is the complex beam parameter normalized by the plasma length scale $\ell$, with $\text {Im}(q_c) \ge 0$. Note that we have made use of the relationship

(4.2)\begin{equation} \frac{\omega \delta_a}{c} = \sqrt{L} . \end{equation}

Note also that one has

(4.3)\begin{equation} - \frac{{\rm i}}{q_c} ={-} {\rm i} \frac{\text{Re}(q_c)}{|q_c|^2} - \frac{ \text{Im}(q_c)}{|q_c|^2} . \end{equation}

Hence, one can identify $|q_c|^2/\text {Re}(q_c)$ as the radius of curvature and $\sqrt {2 L^{-3/2} |q_c|^2/ \text {Im}(q_c)}$ as the beam waist (both normalized by $\ell$). Focusing occurs when $\text {Re}(q_c) > 0$. The linear phase term in (4.1) simply rotates the phase fronts according to the angle of incidence $\theta$ (with $\theta = 0$ being normal incidence), while the additional factor of $\cos ^2 \theta$ in the quadratic phase term in (4.1) accounts for the stretching that occurs for oblique incidence. It is also worth mentioning that when $\theta \neq 0$, (4.1) corresponds to a well-collimated beam (long Rayleigh range) such that the variation of $q_c$ along the incident boundary $z = 0$ can be neglected (Belyaev, Banks & Chapman Reference Belyaev, Banks and Chapman2024).

4.1 Exact solution

The transverse spectrum of the incoming wavefield is computed to be

(4.4)\begin{equation} \hat{\psi}_\text{in}(K_x) = \mathcal{E} \exp\left[ \frac{{\rm i}}{2} \sqrt{L} \, q_c \frac{(K_x - \sqrt{L} \sin\theta)^2}{\cos^2 \theta} \right]\!, \end{equation}

where $\mathcal {E} \doteq E_0\sqrt {{\sqrt {L}\, q_c}/{2{\rm \pi} {\rm i} \cos ^2 \theta }}$ is the overall constant. Equation (3.11) therefore takes the form

(4.5)\begin{align} \psi(X,Z) & = \mathcal{E} \int \,\mathrm{d} K_x \frac{ 2 \operatorname{Ai}\left( K_x^2 + Z - L \right) }{ \operatorname{Ai}\left( K_x^2 - L \right) + {\rm i} \operatorname{Gi}\left( K_x^2 - L \right) } {\rm e}^{{\rm i} K_x X}\nonumber\\ & \quad \times \exp\left[ \frac{{\rm i}}{2} \sqrt{L} \, q_c \frac{(K_x - \sqrt{L} \sin\theta)^2}{\cos^2 \theta} \right]\!. \end{align}

Equation (4.5) depends on four free parameters, three that characterize the boundary value of the incident beam and one that characterizes the medium. They are: (i) $L$, the medium length scale normalized by $\delta _a$ defined in (2.2), related to normalization by the vacuum wavelength via (4.2); (ii) $\text {Re}(q_c)$, which parametrizes the incident radius of curvature normalized by the medium length scale; (iii) $\text {Im}(q_c) \ge 0$, which parametrizes the incident beam waist normalized by the medium length scale; and (iv) $\theta \in [0, {\rm \pi}/2)$, the angle of incidence with respect to the direction of inhomogeneity $Z$.

4.2 Oblique injection as rigid translation and focal shift

Although the injection angle is nominally a free parameter, when the complex beam parameter is purely real, $\text {Im}(q_c) = 0$, then the injection angle $\theta$ can also be removed by a coordinate transformation, up to an overall phase and focal length dilation (focal shift). Indeed, by completing the square, the initial condition (4.1) can be rewritten as

(4.6)\begin{equation} \psi_\text{in}(X) = \mathcal{G} \left( X -L q_c \frac{\tan \theta}{\cos \theta}, L, \frac{q_c}{\cos^2 \theta} \right) \exp\left( \frac{{\rm i}}{2} L^{3/2} q_c \tan^2 \theta \right)\!, \end{equation}

where $\mathcal {G}$ is the Gaussian profile of the incoming beam at normal incidence:

(4.7)\begin{equation} \mathcal{G}(X;L,q_c) = E_0 \exp\left( - {\rm i} \frac{X^2}{2 \sqrt{L}\, q_c} \right)\!. \end{equation}

Clearly, when $\text {Im}(q_c) = 0$, (4.6) corresponds to rigid translation $\varDelta$ in the transverse $X$ direction and a transformed focal length $q_{c,\text {eff}}$ given by

(4.8a)\begin{equation} \varDelta = L q_c \frac{\tan \theta}{\cos \theta} , \quad q_{c,\text{eff}} = \frac{q_c}{\cos^2 \theta} . \end{equation}

In particular, it was shown in Lopez et al. (Reference Lopez, Kur and Strozzi2023) that critical focusing occurs at normal incidence when $q_c = 2$. Hence, for oblique incidence the critical focusing occurs when

(4.9)\begin{equation} q_{c,{\rm crit}} = 2\cos^2 \theta . \end{equation}

Note that expressions (4.8a,b) and (4.9) do not contain the additional Goos–Hänchen and focal shifts contained in the analogous formula from Lopez et al. (Reference Lopez, Kur and Strozzi2023), since these phenomena are only present with a finite beam waist (McGuirk & Carniglia Reference McGuirk and Carniglia1977). Hence, one should use these expressions rather than those of Lopez et al. (Reference Lopez, Kur and Strozzi2023) when the beam waist is sufficiently large.

4.3 Example: softened critical focusing with finite beam waist

As discussed in Lopez (Reference Lopez2023), the peak intensity of a critically focused wave ((4.9) being satisfied) can exceed the standard Airy intensity peak by orders of magnitude.Footnote 1 However, this only occurs when the beam waist is formally infinite. At normal incidence, the characteristic ‘detuning width’ to still achieve critical focusing is $\varDelta _q = 1/\sqrt {L}$ (Lopez et al. Reference Lopez, Kur and Strozzi2023); if this is entirely accounted for by a finite beam waist ($\text {Im}(q_c) \neq 0$), then (4.8a,b) implies that for critical focusing to occur, one must have

(4.10)\begin{equation} \text{Im}(q_c) \lesssim \frac{\cos^2 \theta}{\sqrt{L}} . \end{equation}

If one further assumes that $L \gg 1$, which is equivalent to assuming that $q_{c,{\rm crit}} \gg \text {Im}(q_c)$ with $q_{c,{\rm crit}}$ given by (4.9) and $\text {Im}(q_c)$ given by (4.10), then one can use the relationship between the beam waist $W$ and $q_c$ introduced following (4.3), namely

(4.11)\begin{equation} W = \delta_a L^{1/4} \sqrt{ \frac{ 2 |q_c|^2 }{ \text{Im}(q_c) } } \approx \delta_a L^{1/4} \sqrt{ \frac{ 2 q_{c,{\rm crit}}^2 }{ \text{Im}(q_c) } } , \end{equation}

to obtain a more intuitive condition for maintaining critical focusing compared with (4.10):

(4.12)\begin{equation} \frac{W}{\delta_a} \gtrsim \sqrt{L} \cos \theta . \end{equation}

Figures 2–4 demonstrate how the hyperbolic umbilic diffraction pattern corresponding to critical focusing gets destroyed by a finite beam waist at various values of $\theta$, confirming the prediction of (4.12). (Note that they also confirm the formulas (4.8a,b) and (4.9) relating oblique injection with translations and focal shifts for infinitely wide beams.) Also shown are cases when the beam waist is minimized at fixed $\text {Re}(q_c)$ (which occurs for $\text {Im}(q_c) = \text {Re}(q_c)$), and when $\text {Im}(q_c) \gg \text {Re}(q_c)$. By softening the hyperbolic umbilic caustic, the peak intensity decreases with increasing $\text {Im}(q_c)$; this has clear consequences for applications that require strong focusing near the turning point (such as those discussed in the introduction) and suggests that one should minimize $\text {Im}(q_c)$ as much as possible.

Figure 2.

Morphology of the hyperbolic umbilic caustic that occurs at critical focusing (4.9) as $\text {Im}(q_c)$ is increased from zero at normal incidence ($\theta = 0$) with $L = 10$. The plots are obtained from numerically integrating the exact solution (4.5). From left to right, the panels show the diffraction pattern at critical focusing, at the expected detuning value of $\text {Im}(q_c)$ (4.10), at the $\text {Im}(q_c)$ that minimizes the beam waist and at $\text {Im}(q_c) \gg \text {Re}(q_c)$. Note that the colourbar axis changes for each plot.

Figure 3.

Same as figure 2 but for $\theta = 30^\circ$.

Figure 4.

Same as figure 2 but for $\theta = 60^\circ$.

4.4 Comparison with existing asymptotic formulas

For comparison purposes, let us also list the existing asymptotic formulas provided by Orlov & Tropkin (Reference Orlov and Tropkin1980), Maj et al. (Reference Maj, Pereverzev and Poli2009) and Lopez et al. (Reference Lopez, Kur and Strozzi2023), which we refer to respectively as the O80, the M09 and the L23 formulas. (We reiterate that to the best of our knowledge, only asymptotic solutions to the linear-layer problem with prescribed boundary value have appeared in the literature.) These formulas are given for the initial condition (4.1) asFootnote 2

(4.13)\begin{align} \psi_\text{O80}(X,Z) & = \mathcal{E} \int \,\mathrm{d} K_x \frac{ 2 \operatorname{Ai}\left( K_x^2 + Z - L \right) }{ \operatorname{Ai}\left( K_x^2 - L \right) - {\rm i} \operatorname{Ai}'\left( K_x^2 - L \right)/\sqrt{L - K_x^2 } } {\rm e}^{{\rm i} K_x X} \nonumber\\ & \quad \times \exp\left[ \frac{{\rm i}}{2} \sqrt{L} \, q_c \frac{(K_x - \sqrt{L} \sin\theta)^2}{\cos^2 \theta} \right]\!, \end{align}
(4.14)\begin{align} \psi_\text{M09}(X,Z) & = \mathcal{E} \frac{2{\rm \pi}}{\sqrt{{\rm \pi} {\rm i}}} \int_{-\sqrt{L}}^{\sqrt{L}} \,\mathrm{d} K_x \left(L - K_x^2 \right)^{1/4} \operatorname{Ai}\left( K_x^2 + Z - L \right) {\rm e}^{{\rm i} K_x X} \nonumber\\ & \quad \times \exp\left[ \frac{{\rm i}}{2} \sqrt{L} \, q_c \frac{(K_x - \sqrt{L} \sin\theta)^2}{\cos^2 \theta} + {\rm i} \frac{2}{3} \left( L - K_x^2 \right)^{3/2} \right]\!, \end{align}
(4.15)\begin{align} \psi_\text{L23}(X,Z) & = \mathcal{E}_{L} \int \,\mathrm{d} K_x \, \operatorname{Ai}\left( K_x^2 + Z - L \right) {\rm e}^{{\rm i} K_x X}\nonumber\\ & \quad \times \exp\left[ \frac{{\rm i}}{2} \sqrt{L} \frac{ q_c - 2 \cos 2\theta \cos \theta }{ \cos^2 \theta } \left( K_x - \sqrt{L} \sin \theta \frac{ q_c + \sin 2\theta \sin \theta }{ q_c - 2 \cos 2\theta \cos \theta } \right)^2 \right]\!, \end{align}

where we have introduced

(4.16)\begin{equation} \mathcal{E}_{L} = \mathcal{E} 2{\rm \pi} \sqrt{ \frac{\sqrt{L} \cos \theta}{{\rm \pi} {\rm i}} } \exp\left( \frac{{\rm i}}{6} L^{3/2} \cos^3 \theta \frac{ 4 q_c - 7 \cos \theta - \cos 3\theta }{ q_c - \cos\theta - \cos 3 \theta } \right)\!. \end{equation}

As discussed in Appendix B, each of the formulas (4.13)–(4.15) are asymptotically equivalent to the exact solution (4.5) when $L$ is much greater than all wavevectors contained within the incoming spectrum. The region of validity for this condition is plotted in figure 5. That said, for finite $L - K_x^2$, each of the listed formulas have some shortcomings: (4.13) contains a singularity at $L = K_x^2$ that is clearly unphysical; (4.14) truncates the integral for $K_x^2 \ge L$ and thereby neglects evanescent contributions to the total field; and (4.15) performs a subsidiary Taylor expansion in $K_x^2$ to highlight the caustic structure, with additional loss of accuracy expected as a result. This includes the incorrect predictions for the transverse shift $\varDelta$ and the focal shift for a wide beam at oblique incidence discussed in the previous section.

Figure 5.

Validity boundary for the asymptotic solutions (4.13)–(4.14) as determined by (B1), either as a function of $q_c$ and $\theta$ for specified $L$ (top row) or as a function of $q_c$ and $L$ for specified $\theta$ (bottom row). Within the red hatched region, only the exact solution (4.5) remains valid. No valid asymptotic region exists for $L \le 1$ and for $\theta = {\rm \pi}/2$. The white vertical rectangles in the panels for $\theta = 0^\circ$, $\theta = 30^\circ$ and $\theta = 60^\circ$ correspond to the range of special cases considered in figures 2–4 respectively.

5 Special case: speckled plane wave at normal incidence

Let us now consider an incoming field obtained via paraxial propagation from a lens aperture with focal length $d$ in two dimensions. If the field is evaluated at best focus (i.e. following a propagation distance $d$), then the (far-field) Fraunhofer diffraction formula gives the solution to beFootnote 3

(5.1)\begin{equation} \psi_\text{in}(X) = E_0 \int \,\mathrm{d} Y \, \psi_0(Y) \exp\left( - {\rm i} \frac{\sqrt{L}}{D} X Y \right)\!, \end{equation}

where $D = d/\delta _a$ and $\psi _0$ is the field profile that illuminates the lens aperture. Let us choose $\psi _0$ to be a uniformly illuminated RPP array that consists of $M$ identical elements of width $W/M$ (so the total width is $W$, normalized by $\delta _a$ as usual):

(5.2)\begin{equation} \psi_0(Y) = \sum_{m={-}(({M - 1})/{2})}^{({M - 1})/{2}} {\rm e}^{{\rm i} \phi_{m}} \, \text{rect}\left( \frac{M}{W} Y - m \right)\!, \end{equation}

where $\phi _{m}$ is the corresponding phase shift, set to be either $0$ or ${\rm \pi}$ (Dixit et al. Reference Dixit, Thomas, Woods, Morgan, Henesian, Wegner and Powell1993). We take $M$ to be odd for simplicity. Also, $\text {rect}(x)$ is the unit rectangular function that is non-zero only within the interval $-1/2 \le x \le 1/2$. One then has

(5.3)\begin{equation} \hat{\psi}_\text{in}(K_x) = E_0 \frac{D}{\sqrt{L}} \sum_{m=1}^{M} {\rm e}^{{\rm i} \phi_{m}} \, \text{rect}\left( \frac{M}{\eta} K_x + m \right)\!, \end{equation}

where we have introduced the coupling coefficient $\eta$ as

(5.4a)\begin{equation} \eta = \frac{\sqrt{L}}{f_{\#}} , \quad f_{\#} \doteq \frac{D}{W} . \end{equation}

Note that $f_{\#}$ is the f-number of the launching aperture. Consequently, (3.11) reads

(5.5)\begin{align} \psi(X,Z) & = E_0 \frac{D}{\sqrt{L}} \sum_{m=1}^{M} \exp\left( {\rm i} \phi_{m} + {\rm i} \eta X \frac{m}{M} \right) \nonumber\\ & \quad \times \int_{- ({\eta}/{2 M})}^{{\eta}/{2 M}} \,\mathrm{d} K_x \frac{ 2 \operatorname{Ai}\left[ \left(K_x + \eta \dfrac{m }{M}\right)^2 + Z - L \right] }{ \operatorname{Ai}\left[ \left(K_x + \eta \dfrac{m }{M}\right)^2 - L \right] + {\rm i} \operatorname{Gi}\left[ \left(K_x + \eta \dfrac{m }{M}\right)^2 - L \right] } {\rm e}^{{\rm i} K_x X} . \end{align}

Equation (5.5) is the exact solution for the incoming speckled wavefield given in (5.1). A complete statistical study of this solution is beyond the scope of the present paper, but one property can be seen immediately. Let us take the number of RPP elements to be sufficiently large, $M \gg \eta$, such that the integral in (5.5) can be approximated by its central value. This yields

(5.6)\begin{equation} \psi(X,Z) \approx E_0 \frac{D}{\sqrt{L}} \sum_{m=1}^{M} \exp\left( {\rm i} \phi_{m} + {\rm i} \eta X \dfrac{m}{M} \right) \dfrac{ 2 \operatorname{Ai}\left[ \left(\eta \dfrac{m }{M}\right)^2 + Z - L \right] }{ \operatorname{Ai}\left[ \left(\eta \dfrac{m }{M}\right)^2 - L \right] + {\rm i} \operatorname{Gi}\left[ \left(\eta \dfrac{m }{M}\right)^2 - L \right] } . \end{equation}

Hence when the coupling is small,

(5.7)\begin{equation} \eta \ll 1 , \end{equation}

the transverse speckle pattern and the longitudinal Airy pattern are decoupled from each other. In this limit one needs to simply multiply a given incoming speckle pattern with a single Airy profile along $Z$ to obtain the full solution.

As shown by the definition of $\eta$ in (5.4a,b), decoupling occurs when a large f-number beam is launched into a plasma with relatively short length scale. Conversely, a small f-number beam launched into a plasma with shallow gradient will exhibit a significantly more complicated longitudinal profile. This is demonstrated in figures 6 and 7, which compare the solution (5.5) for a larger and smaller coupling coefficient. The decoupled speckles exhibit a regular periodicity along the propagation direction, while the coupled speckles exhibit a considerably more complicated interference pattern. The longitudinal swelling of the coupled speckles also differs significantly from the standard Airy swelling, with the maximum intensity no longer necessarily occurring near the turning point.

Figure 6.

Exact solution (5.5) for an incoming speckled wavefield at various values of the coupling coefficient $\eta$ and $L = 10$ and $M = 100$ RPP elements. The transverse direction is normalized by the nominal speckle width (Michel Reference Michel2023), which in the normalized coordinates is given by $x_s \sim 2{\rm \pi} /\eta$. For each pair of panels, the top panel shows the beam profile over the nominal envelope width $M x_s$, while the bottom panel shows the region within the white dashed lines of the top panel.

Figure 7.

Lineouts along $Z$ (top) of the exact solution presented in figure 6 at the location of the brightest speckle, along with the solution averaged over the entire $X$ range (bottom). Also shown are Airy function fits obtained by matching the first peak value and location for each respective plot.

6 Conclusion

In this work, the classic problem of an electromagnetic wave propagating in a linearly varying index of refraction (the ‘linear-layer problem’) is revisited. Specifically, we consider when the wavefield throughout the domain is given by a prescribed incoming field at the boundary transverse to the direction of inhomogeneity. Previous studies have only obtained asymptotic solutions to this problem set-up; here we obtain the exact result by using a spectral matching scheme. The resulting solution (3.11) is expressed as an integral involving the FT of the prescribed incoming boundary field multiplied by a kernel involving the familiar Airy function and also the related Scorer function. The kernel is never singular and therefore is uniformly valid regardless the relative sizes of the launched wavelength and the medium length scale.

An incident Gaussian beam is then studied as a test case, with corresponding exact solution given by (4.5). It is shown that when the beam is sufficiently wide (i.e. a Gaussian-focused plane wave), oblique angle of incidence results in a simple shift of the focal length and a rigid translation of the diffraction pattern. Explicit expressions are provided by (4.8a,b) and (4.9), which modify the analogous expressions presented in Lopez et al. (Reference Lopez, Kur and Strozzi2023) by removing the Goos–Hänchen shifts that cannot arise without a finite beam waist. That said, the general prediction of Lopez et al. (Reference Lopez, Kur and Strozzi2023) that the hyperbolic umbilic caustic is structurally stable with respect to angle of incidence still holds true.

It is then demonstrated how the hyperbolic umbilic caustic corresponding to critical focusing of a wide beam gets softened and ultimately disappears as the beam waist is reduced. This softening of the caustic is accompanied by a reduction in the peak intensity obtained by the Gaussian beam at the turning point, which suggests that for applications reliant on the intensity of a wavefield near a turning point, one should make the imaginary part of the complex beam parameter for the launched beam as small as possible. More quantitatively, the hyperbolic umbilic caustic deteriorates once the beam waist $W$ becomes smaller than $W/\delta _a \lesssim \sqrt {L}$, where $\delta _a$ is the Airy skin depth (2.2) and $L$ is the medium length scale normalized by $\delta _a$. This means that beam-tracing solutions (Maj et al. Reference Maj, Pereverzev and Poli2009, Reference Maj, Balakin and Poli2010) are fundamentally incapable of describing the hyperbolic umbilic caustic, since the validity of beam tracing requires $W/\delta _a \sim \sqrt [4]{L}$ (with $L$ also being large). Advanced ray-tracing methods (Lopez & Dodin Reference Lopez and Dodin2020, Reference Lopez and Dodin2021, Reference Lopez and Dodin2022; Lopez, Højlund & Senstius Reference Lopez, Højlund and Senstius2024; Højlund Marholt, Senstius & Nielsen Reference Højlund Marholt, Senstius and Nielsen2024) may be able to describe it, however, as they have no such restriction on the beam waist. Conversely, this analysis suggests that any anomalous focusing observed in beam-tracing solutions is not due to the hyperbolic umbilic caustic, but instead due to the relatively simpler Airy (fold) caustic.

Finally, an incident speckled wavefield is also studied. It had been observed previously in numerical parameter scans that large f-number beams in steep gradients experience decoupled speckle and Airy behaviour (Lopez et al. Reference Lopez, Kur, Chapman, Strozzi and Michel2021), but no concise coupling parameter had been identified. Here we derive the coupling parameter to be $\eta = \sqrt {L}/f_{\#}$, and show that speckles only influence the longitudinal profile of the total wavefield when $\eta \gtrsim 1$, or equivalently, when $\ell \gtrsim \lambda f_{\#}^3/2{\rm \pi}$. This suggests that a reduced model of laser beams near turning points is obtained in the low-coupling regime by simply multiplying the transverse speckle pattern with the Airy swelling factor. In the strong-coupling regime, however, the interference patterns are seen to be considerably more complicated, and more work remains to develop reduced models in this limit. To put this finding in context, in terms of the NIF laser (Spaeth et al. Reference Spaeth, Manes, Kalantar, Miller, Heebner, Bliss, Spec, Parham, Whitman and Wegner2016) the transition to strong coupling occurs for a $351$ nm f/22 beam when the plasma length scale exceeds $3.7$ mm, or when it exceeds $0.2$ mm for an f/8 quad beam.

Acknowledgements

The author thanks Dr J. Ruiz Ruiz of Oxford University for useful conversations.

Editor P. Ricci thanks the referees for their advice in evaluating this article.

Funding

This work was supported by STFC (grant number ST/W000903/1).

Declaration of interests

The author reports no conflict of interest.

Appendix A. Interlacing of Airy and Scorer zeros

It is important to note that

(A1)\begin{equation} \operatorname{Ai}(x) \pm {\rm i} \operatorname{Gi}(x) \neq 0 \end{equation}

because the zeros of $\operatorname {Ai}$ and $\operatorname {Gi}$ never overlap. Since this fact of $\operatorname {Gi}$ is difficult to find in the literature,Footnote 4 we present a simple argument here for completeness.

By definition, one has (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010)

(A2)\begin{equation} \operatorname{Gi}(x) = \operatorname{Bi}(x) \int_x^\infty \operatorname{Ai}(z) \,\mathrm{d} z + \operatorname{Ai}(x) \int_0^x \operatorname{Bi}(z) \,\mathrm{d} z . \end{equation}

If one evaluates $\operatorname {Gi}$ at a zero of $\operatorname {Ai}$, denoted $x_0$ (and one notes that $x_0 < 0$), one obtains

(A3)\begin{equation} \operatorname{Gi}(x_0) = \operatorname{Bi}(x_0) \int_{x_0}^\infty \operatorname{Ai}(z) \,\mathrm{d} z . \end{equation}

Since the zeros of $\operatorname {Ai}$ and $\operatorname {Bi}$ are interlaced, one has $\operatorname {Bi}(x_0) \neq 0$. Also, numerical investigation (figure 8) shows that

(A4)\begin{equation} \int_{x}^\infty \operatorname{Ai}(z) \,\mathrm{d} z > 0, \end{equation}

for all values of $x$. Hence $\operatorname {Gi}(x_0) \neq 0$.

Figure 8.

Evaluating the integral $\int _{x}^\infty \operatorname {Ai}(z) \,\mathrm {d} z$ for different values of the lower limit $x$. Note that since $\operatorname {Ai}(\zeta ) > 0$ for $\zeta > 0$, the integral is manifestly positive for $x > 0$ and asymptotes to zero as $x \to +\infty$.

Appendix B. Asymptotic equivalence of the exact solution with other published expressions

Here we show that the exact solution (4.5) is asymptotically equivalent to the other expressions listed in (4.13)–(4.15) in the large $L \to \infty$ limit. Specifically, we take $L$ to be much larger than all $K_x^2$ that appear in the incoming spectrum, requiring

(B1)\begin{equation} L \left[ 1 - \left( \sin \theta + \frac{\cos \theta}{L^{3/4} \sqrt{|q_c|}} \right)^2 \right] \gg 1 . \end{equation}

Note that to derive (B1), we have estimated the mean wavevector to be $\sqrt {L} \sin \theta$ and the spectral width to be $\cos \theta / |q_c|^{1/2} L^{1/4}$, per (4.4). When (B1) is satisfied, we can take $\zeta \doteq L - K_x^2$ to be sufficiently large, i.e. $\zeta \gg 1$, such that the following asymptotic approximations hold:

(B2a)\begin{gather} \operatorname{Ai}(-\zeta) \sim \frac{ \cos \left( \dfrac{2}{3} \zeta^{3/2} - \dfrac{\rm \pi}{4}\right) }{\sqrt{\rm \pi} \zeta^{1/4}}, \end{gather}
(B2b)\begin{gather}\operatorname{Ai}'(-\zeta) \sim \zeta^{1/4} \frac{ \sin \left( \dfrac{2}{3} \zeta^{3/2} - \dfrac{\rm \pi}{4}\right) }{\sqrt{\rm \pi} }, \end{gather}
(B2c)\begin{gather}\operatorname{Gi}(-\zeta) \sim{-} \frac{ \sin \left( \dfrac{2}{3} \zeta^{3/2} - \dfrac{\rm \pi}{4}\right) }{\sqrt{\rm \pi} \zeta^{1/4}} . \end{gather}

B.1 Asymptotic equivalence with O80 formula

The exact solution (4.5) and the O80 formula (4.13) differ only via the denominator in the integrand: (4.5) has $\operatorname {Ai}(-\zeta ) + {\rm i} \operatorname {Gi}(-\zeta )$ while (4.13) has $\operatorname {Ai}(-\zeta ) - {\rm i} \operatorname {Ai}'(-\zeta )/\sqrt {\zeta }$. When $\zeta \gg 1$, the asymptotic relations (B2) then give

(B3a)\begin{gather} \operatorname{Ai}(-\zeta) + {\rm i} \operatorname{Gi}(-\zeta) \sim \frac{ \exp \left( - {\rm i} \dfrac{2}{3} \zeta^{3/2} + {\rm i} \dfrac{\rm \pi}{4}\right) }{\sqrt{\rm \pi} \zeta^{1/4}} , \end{gather}
(B3b)\begin{gather}\operatorname{Ai}(-\zeta) - {\rm i} \frac{\operatorname{Ai}'(-\zeta)}{\sqrt{\zeta}} \sim \frac{ \exp \left( - {\rm i} \dfrac{2}{3} \zeta^{3/2} + {\rm i} \dfrac{\rm \pi}{4}\right) }{\sqrt{\rm \pi} \zeta^{1/4}} . \end{gather}

The integrands to both formulas have the same asymptotic limit when $\zeta \gg 1$; hence the two formulas should agree when (B1) is satisfied. That said, it is clear that the integrand of (4.13) has a singularity at $\zeta = 0$, so it will not hold uniformly in $\zeta$.

B.2 Asymptotic equivalence with M09 formula

Let us again assume that $\zeta \gg 1$. Applying the asymptotic formula (B3a) to (4.5) then gives the integrand of (4.14). However, (4.14) also restricts the integration bounds from the entire real line to the interval $K_x \in [-\sqrt {L}, \sqrt {L}]$. This ensures that $\zeta$ remains positive over the integration, but neglects evanescent modes (with long decay lengths) that may contribute to the exact solution. These modes should be absent when (B1) is well satisfied, but may be present when it is only marginally true.

B.3 Asymptotic equivalence with L23 formula

Again, let us take $\zeta \gg 1$ and apply (B3a) to (4.5) to obtain the integrand of (4.14). We then perform a subsidiary expansion in the limit of small spectral width $\varDelta _k \doteq K_x - K_0$ about the mean wavevector $K_0 = \sqrt {L} \sin \theta$. This condition reads $L \cos ^2 \theta \gg |(\varDelta _k + 2 K_0)\varDelta _k|$, which is equivalent to (B1) when the same estimate for $\varDelta _k$ is used. The lowest-order approximation for the amplitude gives

(B4)\begin{equation} \zeta^{1/4} = L^{1/4} \sqrt{\cos \theta} + O(\varDelta_k), \end{equation}

while a higher-order approximation for the phase gives

(B5)\begin{equation} \frac{2}{3} \zeta^{3/2} = \frac{2}{3} L^{3/2} \cos^3 \theta - \epsilon L \sin 2 \theta - \epsilon^2 \sqrt{L} \, \frac{\cos 2 \theta}{\cos \theta} + O(\varDelta_k^3) . \end{equation}

Following some algebra, it can then be shown that inserting (B4) and (B5) into (4.5) recovers (4.15). Therefore, (4.5) and (4.15) are asymptotically equivalent when (B1) is well satisfied.

Footnotes

1 The ratio of the two intensities formally diverges in the short-wavelength asymptotic limit.

2 Note that we have corrected an error in Maj et al. (Reference Maj, Pereverzev and Poli2009); their equation (17) is missing an extra factor of $\omega \ell /c$ in the exponent.

3 See, for example, the general discussion in Lopez (Reference Lopez2022), and specifically the cascaded system given by the matrix product of their equations (3.137) and (3.146).

4 For example, Gil, Segura & Temme (Reference Gil, Segura and Temme2003) only compares the zeros of $\operatorname {Gi}$ and $\operatorname {Bi}$, not $\operatorname {Gi}$ and $\operatorname {Ai}$ as we require here.

References

Belyaev, M.A., Banks, J. & Chapman, T. 2024 Exact wave solver for nonparaxial laser beam propagation. Phys. Plasmas 31 (5), 053901.CrossRefGoogle Scholar
Dixit, S.N., Thomas, I.M., Woods, B.W., Morgan, A.J., Henesian, M.A., Wegner, P.J. & Powell, H.T. 1993 Random phase plates for beam smoothing on the nova laser. Appl. Opt. 32 (14), 2543.CrossRefGoogle ScholarPubMed
Gil, A., Segura, J. & Temme, N.M. 2003 On the zeros of the scorer functions. J. Approx. Theory 120 (2), 253.CrossRefGoogle Scholar
Ginzburg, V.L. 1961 Propagation of Electromagnetic Waves in Plasma. Gordon and Breach.Google Scholar
Hall-Chen, V.H., Parra, F.I. & Hillesheim, J.C. 2022 Beam model of doppler backscattering. Plasma Phys. Control. Fusion 64 (9), 095002.CrossRefGoogle Scholar
Højlund Marholt, R., Senstius, M.G. & Nielsen, S.K. 2024 Demonstration of metaplectic geometrical optics for reduced modeling of plasma waves. Phys. Rev. E 110 (2), 025208.CrossRefGoogle ScholarPubMed
Lopez, N.A. 2022 Metaplectic geometrical optics. PhD thesis, Princeton University.Google Scholar
Lopez, N.A. 2023 New solution to Airy's equation for modeling beams near turning points. In 49th EPS Conference on Plasma Physics, EPS 2023, p. P2.011. arXiv:2309.15108.Google Scholar
Lopez, N.A. & Dodin, I.Y. 2020 Restoring geometrical optics near caustics using sequenced metaplectic transforms. New J. Phys. 22 (8), 083078.CrossRefGoogle Scholar
Lopez, N.A. & Dodin, I.Y. 2021 Metaplectic geometrical optics for modeling caustics in uniform and non-uniform media. J. Opt. 23 (2), 025601.CrossRefGoogle Scholar
Lopez, N.A. & Dodin, I.Y. 2022 Metaplectic geometrical optics for ray-based modeling of caustics: theory and algorithms. Phys. Plasmas 29 (5), 052111.CrossRefGoogle Scholar
Lopez, N.A., Højlund, R. & Senstius, M.G. 2024 Regarding the extension of metaplectic geometrical optics to modelling evanescent waves in ray-tracing codes. Phys. Plasmas 31 (8), 083901.CrossRefGoogle Scholar
Lopez, N.A., Kur, E., Chapman, T.D., Strozzi, D.J. & Michel, P.A. 2021 Influence of speckles on laser intensity profile near turning points. Bull. Am. Phys. Soc. 66 (13), Abstract NP11.00063.Google Scholar
Lopez, N.A., Kur, E. & Strozzi, D.J. 2023 Intensity of focused waves near turning points. Phys. Rev. E 107 (5), 055204.CrossRefGoogle ScholarPubMed
Maj, O., Balakin, A.A. & Poli, E. 2010 Effects of aberration on paraxial wave beams: beam tracing versus quasi-optical solutions. Plasma Phys. Control. Fusion 52 (8), 085006.CrossRefGoogle Scholar
Maj, O., Pereverzev, G.V. & Poli, E. 2009 Validation of the paraxial beam-tracing method in critical cases. Phys. Plasmas 16 (6), 062105.CrossRefGoogle Scholar
McGuirk, M. & Carniglia, C.K. 1977 An angular spectrum representation approach to the Goos-Hanchen shift. J. Opt. Soc. Am. 67 (1), 103.CrossRefGoogle Scholar
Michel, P. 2023 Introduction to Laser-Plasma Interactions. Springer.CrossRefGoogle Scholar
Olver, F.W.J., Lozier, D.W., Boisvert, R.F. & Clark, C.W. 2010 NIST Handbook of Mathematical Functions. Cambridge University Press.Google Scholar
Ono, M., Bertelli, N. & Shevchenko, V. 2022 a Multi-harmonic electron cyclotron heating and current drive scenarios for non-inductive start-up and ramp-up in high field ST-40 spherical tokamak. Nucl. Fusion 62 (10), 106035.CrossRefGoogle Scholar
Ono, M., Bertelli, N., Shevchenko, V., Idei, H. & Hanada, K. 2022 b Efficient electron cyclotron current drive regime for plasma current start-up in fusion reactors. Phys. Rev. E 106 (2), L023201.CrossRefGoogle ScholarPubMed
Orlov, Y.I. & Tropkin, S.K. 1980 Field of a planar aperture antenna in an inhomogeneous medium. Radiophys. Quantum Electron. 23 (12), 979.CrossRefGoogle Scholar
Spaeth, M.L., Manes, K.R., Kalantar, D.G., Miller, P.E., Heebner, J.E., Bliss, E.S., Spec, D.R., Parham, T.G., Whitman, P.K., Wegner, P.J., et al. 2016 Description of the NIF laser. Fusion Sci. Technol. 69 (1), 25.CrossRefGoogle Scholar
Figure 0

Figure 1. The Airy function $\operatorname {Ai}(x)$ and the Scorer function $\operatorname {Gi}(x)$ plotted against their argument. Both functions behave qualitatively similar, exhibiting exponential decay for $x > 0$ and oscillatory behaviour (with relative phase shift) for $x < 0$.

Figure 1

Figure 2. Morphology of the hyperbolic umbilic caustic that occurs at critical focusing (4.9) as $\text {Im}(q_c)$ is increased from zero at normal incidence ($\theta = 0$) with $L = 10$. The plots are obtained from numerically integrating the exact solution (4.5). From left to right, the panels show the diffraction pattern at critical focusing, at the expected detuning value of $\text {Im}(q_c)$ (4.10), at the $\text {Im}(q_c)$ that minimizes the beam waist and at $\text {Im}(q_c) \gg \text {Re}(q_c)$. Note that the colourbar axis changes for each plot.

Figure 2

Figure 3. Same as figure 2 but for $\theta = 30^\circ$.

Figure 3

Figure 4. Same as figure 2 but for $\theta = 60^\circ$.

Figure 4

Figure 5. Validity boundary for the asymptotic solutions (4.13)–(4.14) as determined by (B1), either as a function of $q_c$ and $\theta$ for specified $L$ (top row) or as a function of $q_c$ and $L$ for specified $\theta$ (bottom row). Within the red hatched region, only the exact solution (4.5) remains valid. No valid asymptotic region exists for $L \le 1$ and for $\theta = {\rm \pi}/2$. The white vertical rectangles in the panels for $\theta = 0^\circ$, $\theta = 30^\circ$ and $\theta = 60^\circ$ correspond to the range of special cases considered in figures 2–4 respectively.

Figure 5

Figure 6. Exact solution (5.5) for an incoming speckled wavefield at various values of the coupling coefficient $\eta$ and $L = 10$ and $M = 100$ RPP elements. The transverse direction is normalized by the nominal speckle width (Michel 2023), which in the normalized coordinates is given by $x_s \sim 2{\rm \pi} /\eta$. For each pair of panels, the top panel shows the beam profile over the nominal envelope width $M x_s$, while the bottom panel shows the region within the white dashed lines of the top panel.

Figure 6

Figure 7. Lineouts along $Z$ (top) of the exact solution presented in figure 6 at the location of the brightest speckle, along with the solution averaged over the entire $X$ range (bottom). Also shown are Airy function fits obtained by matching the first peak value and location for each respective plot.

Figure 7

Figure 8. Evaluating the integral $\int _{x}^\infty \operatorname {Ai}(z) \,\mathrm {d} z$ for different values of the lower limit $x$. Note that since $\operatorname {Ai}(\zeta ) > 0$ for $\zeta > 0$, the integral is manifestly positive for $x > 0$ and asymptotes to zero as $x \to +\infty$.