Past exam of the mathematics course of the University of Cambridge / 2018 / iii / Paper 335 (Direct and Inverse Scattering of Waves)

Use the time-harmonic wave convention $\psi e^{-i\omega t}$ and write $\mathbf{z}=(y,z)$, $E(0,\mathbf{z})=E_0(\mathbf{z})$, with a deterministic incident wave envelope. The printed speed ratio is inconsistent with $k=k_{0}n$: the refractive index used below is $n=c_0/c$. Substituting $\psi =e^{i{k}_{0}x}E$ in the Helmholtz equation givesThe paraxial approximation discards ${\partial }_x^{2}E$. The usual weak-fluctuation model also linearizes $n^2-1=2\mu W+O({\mu }^2{W}^2)$, givingDropping the quadratic contrast is an additional weak-fluctuation assumption, not a consequence of small propagation angles. Although the linearized random potential generates attenuation of order ${\mu }^2$, it does not retain every effect of that order in the literal finite-correlation index: the discarded quadratic contrast can also produce a mean wave phase shift. This linearization must precede a Gaussian white noise limit: the square of ideal white noise has no ordinary pointwise meaning.

The split-step Fourier method alternates free-space diffraction, $U_h=e^{hD}$, with a random phase screen,For jointly Gaussian random fields, Gaussian phase averaging gives $\mathbb{E}{e}^{i{\varphi }_h}=e^{-\mathrm{Var}({\varphi }_h)/2}$. More generally, a product of $p$ fields and $q$ conjugate fields picks up $\exp (i{\sum }_a{ϵ}_a{\varphi }_h({\mathbf{z}}_a))$, with signs ${ϵ}_a=+1$ or $-1$. Its screen average is determined entirely by the wave phase covariance matrix. The deterministic diffraction step acts on each coordinate, with opposite signs on factors formed by complex conjugation. This is the basis of the field-moment equations.

There is an important closure qualification. A stationary Gaussian random field need not have independent longitudinal increments. For finite-correlation fluctuations, the unlinearized parabolic wave equation instead givesIn the weak-fluctuation model, the exact first-moment equation isThe last term cannot in general be replaced by a constant times $m$. An exact generic solution of the linearized model iswhere $\mathcal{T}$ denotes time ordering of the propagation operators. For the unlinearized model, add $i{k}_0{\mu }^2{W}^2/2$ inside the propagation generator. The covariance function of the medium is needed to evaluate this expression; One-point Gaussian distributions alone would not even determine the joint wave phase statistics.

For example, omit diffraction and take a longitudinal autocovariance function $R(s)=\mathbb{E}[W(s,\mathbf{z})W(0,\mathbf{z})]$. Direct Gaussian phase averaging givesEven here, unit-variance stationary Gaussian random fields with different $R$ give different answers. In the general problem the diffraction and multiplication operators do not commute, so this scalar attenuation cannot simply be multiplied by $U_x$ without an additional approximation.

The standard closed answer uses the Markov approximation for a random medium, made explicit in part ii. Let the longitudinally integrated covariance kernel beReplace the medium by longitudinal Gaussian white noise with this strength. A screen of thickness $h$ then has $\mathbb{E}[{\varphi }_h({\mathbf{z}}_a){\varphi }_h({\mathbf{z}}_b)]=k_0^2{\mu }^{2}hB({\mathbf{z}}_a-{\mathbf{z}}_b)$ and is independent of the incoming field. For the moment $M=\mathbb{E}{\prod }_a{E}_a^{({ϵ}_a)}$, expanding both steps to order $h$ givesIn particular, the coherent attenuation in a white-noise random medium and its Fresnel propagator solution areEquivalently, write $Wdx=d{\beta }_x$, where ${\beta }_x(\mathbf{z})$ is Brownian motion in $x$ with transverse covariance kernel $\mathbb{E}[d{\beta }_x(\mathbf{z})d{\beta }_x({\mathbf{z}}^{\prime })]=B(\mathbf{z}-{\mathbf{z}}^{\prime })dx$. The Stratonovich integral formulation is $dE=DEdx+i{k}_0\mu E\circ d{\beta }_x$. Its Itô integral form isThe mean of the Itô integral vanishes, independently confirming the attenuation drift. For $x>0$ in two transverse dimensions,Thus a unit plane wave has $m=e^{-\gamma x}$. For a general incident wave envelope, the Fresnel propagator supplies its spreading. The attenuation is redistribution between the coherent and diffuse wave fields, rather than wave absorption: for a field and its conjugate at the same point, the screen contribution in the second-moment equation cancels.

Interpret the longitudinal Dirac delta function as the Markov approximation for a random medium. The two arguments on the left of the printed covariance must be $x_1$ and $x_2$. The transverse covariance kernel is $B(|{\mathbf{z}}_1-{\mathbf{z}}_2|)$. Here statistical isotropy means transverse isotropy: a medium with a distinguished longitudinal white noise direction and a smooth transverse covariance is not literally isotropic in all three directions. Also, ideal Gaussian white noise replaces the original finite-variance field; it cannot simultaneously satisfy a pointwise normalization $\mathbb{E}{W}^2=1$.

Write $Q(r)={\mu }^{2}B(r)$ for the covariance strength of the refractive index fluctuations. Define the transverse power spectrum using the unnormalized forward Fourier transform printed in equation 3:This is the spectrum of the fluctuations, excluding the deterministic mean index. Equivalently it is the zero-longitudinal-frequency slice of the three-dimensional fluctuation spectrum before the Markov approximation. The angular integral is $2\pi J_0(\nu r)$, where $J_0$ is the Bessel function of the first kind. Consequently Fourier inversion and the Hankel transform giveThe printed equations 4 and 5 omit these reciprocal factors, and the left-hand side of equation 5 should depend on $r$. They are instead a consistent order-zero Hankel transform pair if their $F$ denotes $S_H=S_F/(2\pi )$.

Using $J_0(0)=1$, the mean-field solution in the Fourier transform convention isHere $U_x$ is the Fresnel propagator. Finiteness of the integral ensures a finite screen wave phase variance. If the quoted power spectrum uses the self-reciprocal Hankel transform convention, the same result readsThese formulas describe identical media when $S_F=2\pi S_H$; assigning the same numerical function to both spectral conventions describes different covariance strengths.

For the Fourier transform convention of equation 3, insert the given power spectrum as $S_F$. With $t=\nu L/2$,The Fourier transform of a Gaussian, or the order-zero Hankel transform, also givesThus the coherent attenuation in a white-noise random medium isFor a unit incident plane wave, the Fresnel propagator leaves the initial envelope unchanged, so the answer reduces to $\boxed{m=e^{-k_0^2{\mu }^{2}Lx/(2\pi )}}$. If instead the given numerical function is interpreted as the self-reciprocal Hankel transform spectrum suggested by the printed equations 4 and 5, then $Q(r)=2{\mu }^{2}L{e}^{-r^2/L^2}$ and the alternative isThe factor $2\pi$ ambiguity is present in the original PDF. No unique numerical attenuation coefficient follows until the power spectrum convention is fixed. The covariance strength $Q(0)$ has units of length in this longitudinal white noise model, making both attenuation exponents dimensionless.

In time reversal acoustics, the array records the incoming signal, reverses each recorded time trace, and re-emits it through the same medium. This is phase conjugation in the frequency domain. For the time-harmonic wave convention $e^{-i\omega t}$, reversing a real time trace replaces its positive-frequency wave amplitude by its complex conjugate. The medium must remain unchanged between recording and re-emission.

Let $f(\mathbf{z})={\Psi }_0(0,\mathbf{z})$ and $G_L(\mathbf{\zeta },\mathbf{z})=G(L,\mathbf{\zeta };0,\mathbf{z})$, with the chosen source and array normalizations incorporated into the Green function. Let $a(\mathbf{\zeta })\ge 0$ be the array's aperture weight, equal to the indicator of its receiving region for an ideal uniform array. The recorded field isBy wave reciprocity, back-propagation has the same Green function with the source and receiver exchanged. Therefore the physically re-emitted, back-propagated wave amplitude isIf $Hf=R$ and $P_A$ multiplies by $a$, then ${\Psi }_B=H^T{P}_A\overline{Hf}$, where $T$ denotes the transpose without conjugation. Taking a final complex conjugate instead defines the adjoint reconstruction $\overline{{\Psi }_B}=H^{\ast }{P}_{A}Hf$. This distinction prevents an erroneous conjugation in the time reversal operator.

For a localized Gaussian beam or acoustic point source in a homogeneous medium, a finite aperture admits a limited range of angles. The focal width is of order $\lambda L/A$ when $A$ denotes the aperture diameter. In a random medium, multiple scattering creates paths with a larger angular spread. Each reversed path retraces its route, and the paths interfere constructively at the source. This can produce a larger effective aperture in time reversal and a narrower focus, even though the unreversed field has a complicated speckle pattern.

This comparison concerns a homogeneous reference medium; a deterministic heterogeneous medium can also provide useful multipath propagation. Suitable scale limits or frequency and spatial averaging can make refocusing self-averaging. Such self-averaging is not automatic for every monochromatic source and every random realization. Wave absorption, changing medium parameters, unresolved paths or poor array coverage can spoil refocusing. With complete capture of the propagating modes and a lossless unitary operator $H$, ideal adjoint reconstruction is already exact in either medium. Random scattering can improve finite-aperture wave focusing through angular diversity. See the regime-dependent analysis in Statistical stability in time reversal.

Fix a receiver position $\mathbf{\zeta }$ and define $g_{\mathbf{\zeta }}(\mathbf{z})=G_L(\mathbf{\zeta },\mathbf{z})$. Choose the two scalar fields in the Wigner distribution to be $u=v=g_{\mathbf{\zeta }}$. Fourier inversion in the separation variable givesIndeed, substituting the defining Wigner distribution makes the $\mathbf{p}$ integral a $(2\pi {)}^2$ Dirac delta function, canceling its normalization and setting the separation equal to $\mathbf{z}-{\mathbf{z}}^{\prime }$.

Average over the receiver plane, including the aperture:Then the receiver-averaged Wigner distribution represents the back-propagated field asThe fast decay of $f$ justifies its use as a test function. For oscillatory Green functions that are not integrable, the Wigner distribution and its inversion can be understood as tempered distributions, or derived with smooth cutoffs before taking their limits.

The printed $W_A$ expression can also be used literally, but its displayed ${\mathbf{z}}^{\prime }$ integral is over source-plane coordinates, not receiver positions. The missing receiver average and source weight can be encoded by choosing vector fields whose channels are the receivers. LetContract the receiver channels in the product of the vector fields, so thatWith exactly the printed definitionthe same Fourier inversion now yields the requested compact formFor a finite array, the receiver integral in the contraction is a weighted sum. For a continuous array, these are fields valued in the receiver Hilbert space. An uncontracted vector outer product instead gives a matrix-valued Wigner distribution, whose receiver trace must be taken. Thus the printed formula is usable with these choices and this contraction convention, but its phrase “over the plane of the receiver” does not describe the displayed integral.

Let $H$ be the measured wave scattering response> linear operator, with direct background propagation removed if present. Its time reversal operator is the positive operator $T=H^{\ast }H$. For a normalized emitted signal $h$,The maximizing signal is a right singular vector, equivalently an eigenvector of $T$ with its largest eigenvalue. This is a precise intensity statement independent of a scatterer model.

Under the Born approximation for scalar wave scattering, neglect interactions between point scatterers and writewhere $b_j$ is the receiver response to scatterer $j$, $c_j^{\ast }h$ is its illumination by the source array, and ${\alpha }_j$ is its scattering wave amplitude. Well-resolved point scatterers have approximately orthogonal vectors $b_j$ and $c_j$. In the ideal orthogonal limit,Thus the most reflective, geometrically weighted scatterer gives the largest eigenvalue, and its normalized steering vector $c_j/\Vert c_j\Vert$ gives the corresponding wave focusing signal. If illumination and reception factors are equal for all scatterers, this is precisely the scatterer with largest $|{\alpha }_j{|}^2$. Geometric size by itself is not the quantity being ranked.

Repeated adjoint time reversal acoustics is power iteration on $T$. If its top eigenvalue is simple and the initial signal has a nonzero component along its eigenvector, thenA tied largest eigenvalue leaves a combination in the leading eigenspace, and an initial signal orthogonal to that eigenspace cannot excite it. This is the basis of the DORT method. Mere physical separation is insufficient if the array cannot resolve the scatterers; coherent steering-vector overlap or multiple scattering can mix the modes. Also, if $H$ includes unrestricted homogeneous transmission, its largest eigenvalue need not identify any individual scatterer. The correspondence concerns the resolved wave scattering response>. The ideal correspondence and selective wave focusing are analyzed by Prada and Fink.

Use the physical refractive index convention $k=k_{0}n$, with $k_0$ the background wavenumber and time dependence $e^{-i\omega t}$. Treat $B$ as a bounded region in ${\mathbb{R}}^3$; the printed $B\in {\mathbb{R}}^3$ should be $B\subset {\mathbb{R}}^3$. Define the contrast $q=n^2-1=2n_{ϵ}+n_{ϵ}^2$ and the outgoing Green functionThe Sommerfeld radiation condition selects this sign of the outgoing wave phase. Since $(\Delta +k_0^2)\psi =-k_0^{2}q\psi$, the Lippmann-Schwinger equation isThe positive sign in this integral follows from the minus sign in the defining Green function equation. One may instead define the scattering potential with the opposite sign, provided both equations change consistently.

If $\Vert T_q\Vert <1$ on the chosen field space inside $B$, its Neumann series converges. The Born series for the scattered field isFor observation points outside $B$, restrict the intermediate factors to $B$ and use the same outgoing integral for the final factor. For example, the first two terms areThe $m$th term describes $m$ successive scattering interactions. The first Born approximation for scalar wave scattering keeps one interaction and replaces the field inside the medium by the incident field. Higher terms describe multiple scattering and the resulting feedback on the internal field.

A concrete sufficient condition follows on $L^{\infty }(B)$. If $d_B$ is the diameter of $B$, then for $\mathbf{r}\in B$ the region lies in the ball of radius $d_B$ about $\mathbf{r}$, andThus $\boxed{\rho \ll 1}$ is a conservative sufficient condition for the first Born approximation for scalar wave scattering. When $\rho <1$, the omitted internal-field terms satisfyThe bound ignores cancellation in the oscillatory Green function, so it is sufficient, not necessary. A common physical small-contrast criterion for an extended weak medium is small accumulated extra wave phase, $k_0\ell \Vert n_{ϵ}{\Vert }_{\infty }\ll 1$, together with weak scattering and no resonant internal enhancement. Small local contrast alone is not enough for an arbitrarily large or resonant object.

To first order in $n_{ϵ}$, replace $q$ by $2n_{ϵ}$ and use the known incident field in the outgoing integral. The Born approximation for scalar wave scattering givesBoth the $n_{ϵ}^2$ contrast term and the internal-field correction are second order for fixed geometry in the perturbative regime established in part i.

Specify the data space before taking an adjoint operator. For example, let $X=L^2(B)$ and $Y=L^2(\Gamma )$, where $\Gamma$ is a bounded measurement surface separated from $B$. Define the linear operatorFor bounded incident field this is a Hilbert-Schmidt operator, hence a compact operator. The equation is $A{n}_{ϵ}=y$, with $y={\psi }_s{|}_{\Gamma }$. Other sampling geometries give corresponding data spaces and weights; the paper does not specify one. Using the usual complex $L^2$ inner products, its adjoint operator isThe complex conjugations are required by the adjoint operator identity, not by wave reciprocity alone.

The Landweber iteration starts from $f_0=0$ and applies gradient descent to $\frac{1}{2}\Vert Af-y{\Vert }^2$:Each step back-propagates the data residual. The Landweber relaxation parameter controls stability, and early stopping of Landweber iteration prevents small singular values from amplifying measurement errors. For an explicitly real-valued index contrast, use the real Hilbert space structure and replace $A^{\ast }$ by $\mathrm{Re}{A}^{\ast }$ in this update. Additional sign or support constraints require corresponding projections; none are assumed here.

Write $y={\psi }_s$ and retain the Landweber relaxation parameter $\tau$. The printed $\Vert A\Vert <2$ allows the fixed choice $\tau =1/2$, since then $\tau \Vert A{\Vert }^2<2$. It does not justify a unit step: the scalar linear operator $A=3/2$ has norm below $2$, but unit-step error is multiplied by $1-|A{|}^2=-5/4$ and diverges. We use $0<\tau <2/\Vert A{\Vert }^2$ throughout.

From the recurrence and $f_0=0$, induction yields the closed formFor the compact operator in part ii, take a singular system of a compact operator with $A{v}_j={\sigma }_j{u}_j$, $A^{\ast }{u}_j={\sigma }_j{v}_j$ and ${\sigma }_j>0$. Thus ${\sigma }_j^2$ are the eigenvalues of $A^{\ast }A$; this fixes the notation implicit in the printed hint. In the singular component $j$, the recurrence isSumming the geometric series gives the Landweber spectral filterThere is no contribution from data in $\ker A^{\ast }$ or from $\ker A$ in the reconstruction. Starting at zero is what selects the minimum-norm least-squares solution.

For fixed ${\sigma }_j>0$, $|1-\tau {\sigma }_j^2|<1$, so the numerator tends to $1$. The assumption $y\in D(A^{†})$ is the Picard criterionwith an arbitrary additional component in $\ker A^{\ast }$. The squared reconstruction error isEach term tends to zero and is bounded by a summable term from the Picard criterion. The dominated convergence theorem therefore proves $R_{n}y\to A^{†}y$, where $A^{†}$ is the Moore–Penrose inverse of an operator.

For each finite $n$, the regularization of an inverse problem is stable. Set $t=\tau {\sigma }^2\in (0,2)$. The geometric series and $|1-t|\le 1$ implyConsequently the Landweber noise amplification bound in this general step-size convention is $\Vert R_n\Vert \le \sqrt{2\tau n}$. If additionally $\tau \Vert A{\Vert }^2\le 1$, the sharper numerator bound $\min (nt,1)$ gives $\Vert R_n\Vert \le \sqrt{\tau n}$.

To see the regularization parameter directly, put $\alpha =1/(\tau n)$. Modes with ${\sigma }^2\ll \alpha$ have numerator approximately $n\tau {\sigma }^2$, so their inverse coefficient is approximately $\sigma /\alpha$ rather than $1/\sigma$. Each fixed nonzero mode is eventually restored as $\alpha \downarrow 0$. ThusThe finite iterates suppress unstable small singular values; infinitely many iterations remove this suppression. For noisy data $\Vert y^{\delta }-y\Vert \le \delta$, the noise-bias decomposition for linear regularization givesChoosing $n(\delta )\to \infty$ and ${\delta }^{2}n(\delta )\to 0$ proves noisy-data convergence. The reciprocal iteration index is a regularization parameter, and early stopping controls noise amplification.