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

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.

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 1961). 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 2022a; Ono et al. 2022b) 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 2022) 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 1961; Orlov & Tropkin 1980; Maj, Pereverzev & Poli 2009; Maj, Balakin & Poli 2010; Lopez, Kur & Strozzi 2023) 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 2023). 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. 2010):

(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 1980; Maj et al. 2009; Lopez et al. 2023) 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 1980; Maj et al. 2009; Lopez et al. 2023). 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 2024).

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. (2023) 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. (2023), since these phenomena are only present with a finite beam waist (McGuirk & Carniglia 1977). Hence, one should use these expressions rather than those of Lopez et al. (2023) when the beam waist is sufficiently large.

4.3 Example: softened critical focusing with finite beam waist

As discussed in Lopez (2023), the peak intensity of a critically focused wave ((4.9) being satisfied) can exceed the standard Airy intensity peak by orders of magnitude.¹ 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. 2023); 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 (1980), Maj et al. (2009) and Lopez et al. (2023), 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) as²

(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 be³

(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. 1993). 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 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 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. (2023) by removing the Goos–Hänchen shifts that cannot arise without a finite beam waist. That said, the general prediction of Lopez et al. (2023) 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. 2009, 2010) 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 2020, 2021, 2022; Lopez, Højlund & Senstius 2024; Højlund Marholt, Senstius & Nielsen 2024) 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. 2021), 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. 2016) 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.