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

Use the time-harmonic wave convention $\mathrm{Re}(\psi e^{-i\omega t})$, with a positive background wavenumber $k$. The scalar Helmholtz equation for a constant-density acoustic model isWriting $\psi =e^{ikx}E$ gives the exact envelope equationThe paraxial approximation neglects $E_{xx}$ relative to $2ik{E}_x$, giving the parabolic wave equationThe choice of carrier and the direction of propagation matter: this is a forward, slowly varying envelope approximation. Sufficient scale conditions are $|E_{xx}|\ll 2k|E_x|$, transverse spectral components $|p|\ll k$, and medium/envelope variation on longitudinal scales large compared with $1/k$. With the carrier fixed at the background $k$, small $|n-1|$ makes the refractive phase vary slowly too. Large-angle propagation, appreciable backscattering, or rapid longitudinal variation violates the approximation. Acoustic models with variable mass density can have additional gradient terms, so the scalar Helmholtz equation itself is a model assumption. For the Gaussian beam in part (b), useful initial conditions are $kD\gg 1$ and $D/|F|\ll 1$, with the second condition absent for an initially uncurved beam. A narrow angular spectrum, rather than merely the label “Gaussian”, justifies the paraxial approximation.

PutIn free space, the parabolic wave equation is $E_x=i{E}_{zz}/(2k)$. Under the Fourier transform convention $\hat{E}(p)=\int E(z)e^{-ipz}dz$, it becomesThe Gaussian integral is legitimate because $\mathrm{Re}A=1/D^2>0$. Invert the transform after multiplying by $e^{-i{p}^{2}x/(2k)}$. A second Gaussian integral, or equivalently one-dimensional transverse Fresnel propagation, givesChoose the square-root branch continuously from $q(0{)}^{-1/2}=1$; for real $x\ge 0$ there is no zero of $q$. This gives the correct incident field at $x=0$, and direct differentiation verifies the free parabolic wave equation.

For clarity, the squared envelope magnitude isThe Gaussian beam with one transverse coordinate remains Gaussian, with one-transverse-coordinate amplitude factor $q^{-1/2}$, not the $q^{-1}$ of a beam with two transverse coordinates. The negative initial quadratic phase produces focusing for $F>0$; diffraction prevents a singularity at $x=F$. These expressions describe the paraxial approximation to free propagation, rather than an exact unrestricted Helmholtz equation beam.

No: the first moment generally depends on the transverse coordinate. At entry to the random medium, the initial data are deterministic, sowhich is a nonconstant Gaussian beam profile. Statistical homogeneity of the medium says that shifting both the medium and the incident data gives correspondingly shifted field statistics. It does not make the response to fixed, localized incident data translation invariant.

In particular, Gaussian coherent-field propagation in the Markov approximation, derived in part (d), givesfor the linear weak-index model. A homogeneous attenuation factor multiplies the varying beam profile. Thus spatially homogeneous medium statistics can coexist with a transversely inhomogeneous coherent field. The first moment is not the mean wave intensity, and disappearance of the coherent component is not itself disappearance of the total wave energy.

Write $m=⟨E⟩$ and $L_0=i{\partial }_z^2/(2k)$. Substituting $n=1+\mu W$ into part (a), without prematurely replacing $n^2$ by its linear part, givesTaking the expectation yields the exact mean equation within the paraxial approximation:The correlation term cannot be discarded just because $⟨W⟩=0$: the field depends on the same random medium as $W$. Nor can $⟨W^{2}E⟩$ be factored exactly from the given variance alone. The printed statistics do not specify the autocorrelation function of a random field, so they do not determine a unique local closed equation for $m$.

A useful covariance-dependent equation follows to second order in $\mu$. Defineand letbe the Fresnel propagator kernel for one transverse coordinate. The Duhamel principle applied to the linear fluctuation term giveswhere $E^{(0)}=E_{\mathrm{free}}$. Multiply by $W(x,z)$ and average. Also $⟨W(x,z{)}^{2}E(x,z)⟩=E^{(0)}(x,z)+O(\mu )$. Since $m=E^{(0)}+O({\mu }^2)$, replacing $E^{(0)}$ by $m$ in terms already multiplied by ${\mu }^2$ preserves second-order accuracy. We obtain the weak-fluctuation memory equation for the coherent fieldThis perturbation statement is for fixed propagation intervals with sufficient covariance regularity and moment bounds. The kernel is understood as an oscillatory integral. If the field is jointly Gaussian, odd perturbative contributions vanish and the remainder can be improved under the corresponding expansion hypotheses. Joint Gaussianity itself still leaves the covariance unspecified. The term $ik{\mu }^{2}m/2$ is the quadratic refractive-index shift in the coherent field; it comes from the mean of $W^2$.

For the usual local answer, impose the additional Markov approximation for a random medium. Its longitudinal correlation length must be much shorter than the distance on which the envelope varies, and diffraction across that correlation distance must not resolve substantial transverse covariance variation. Away from the entry layer, replace the slowly varying factors in the memory integral by their local values and extend the longitudinal integral to infinity. Define the integrated covarianceFor an integrable real stationary covariance, $C(s,0)$ is even, and $B(0)\ge 0$ is its zero-longitudinal-frequency spectral density. The conventional leading weak-index model keeps only the linear random potential $ik\mu W$; in that model the local first-moment equation isThe factor $1/2$ can be checked independently in the white-noise model. If ${\mathcal{B}}_x(z)$ is a Brownian field with $⟨d{\mathcal{B}}_x(z)d{\mathcal{B}}_x(z^{\prime })⟩=B(z-z^{\prime })dx$, writeConverting the Stratonovich integral to an Itô integral adds $-k^2{\mu }^{2}B(0)Edx/2$. The remaining stochastic integral has zero mean, giving the boxed equation andIf the finite-correlation model retains $n^2$ consistently through order ${\mu }^2$, the local second-order expansion also has $+ik{\mu }^{2}m/2$, adding the phase factor $e^{ik{\mu }^2(x-x_0)/2}$. This phase is computed before the white-noise idealization: squaring ideal Gaussian white noise is not the finite-variance operation $W^2$.

The need for a correlation assumption can be seen without any closure argument. The stationary Gaussian random field $W(x,z)=Z$, with a single standard normal variable $Z$, has all the printed one-point statistics. In the linear weak-index model it givesIts attenuation is quadratic in propagation distance, whereas the Markov model gives linear-distance exponential attenuation. The local attenuation equation requires the additional correlation/Markov assumption; it does not follow from statistical homogeneity and unit variance alone.

Take the same $e^{-i\omega t}$ convention and an outgoing field. Define the scattering potential with the positive-contrast sign,The outgoing Green function for the three-dimensional Helmholtz equation isThus the Lippmann-Schwinger equation has a plus sign:Replacing the total field inside the integral by the incident field gives the Born approximation for scalar wave scatteringMore generally, iterating the integral equation gives the Born series ${\psi }_s={\sum }_{j\ge 1}{T}^j{\psi }_i$; its first term is the boxed expression. A sufficient convergence condition is $\Vert T\Vert <1$ on the chosen space of fields in $V$, and accurate single-scattering truncation needs the neglected terms to be small. Weak refractive contrast is useful, but for a sufficiently large or resonant scatterer it is not by itself a guarantee of negligible multiple scattering.

For $\mathbf{r}=R\hat{\mathbf{r}}$ with the observation distance much larger than the scatterer, expand $|\mathbf{r}-{\mathbf{r}}^{\prime }|=R-\hat{\mathbf{r}}\cdot {\mathbf{r}}^{\prime }+O(R^{-1})$. The far-field pattern in the Born approximation for scalar wave scattering is consequentlyThis convention includes the factor $1/(4\pi )$ in $f_{\infty }$; moving that factor into the outgoing-wave definition would change the pattern's normalization.

The absolute energy flux depends on what physical amplitude $\psi$ represents. For acoustic pressure in a homogeneous fluid of mass density ${\rho }_0$ and wave speed $c_0$, the momentum equation gives the particle-velocity amplitude $\mathbf{v}=\nabla \psi /(i{\rho }_0\omega )$. The time average of harmonic power and acoustic intensity then giveThe radial derivative of the outgoing factor is $(ik-1/R)e^{ikR}/R$. Its imaginary current is therefore $k|f_{\infty }{|}^2/R^2$ to leading order. The incident unit-amplitude plane wave has imaginary current $k{\hat{\mathbf{r}}}_0$. Hence, with $k=\omega /c_0$,The scattered energy flux points along $\hat{\mathbf{r}}$ at leading order, and the incident energy flux along ${\hat{\mathbf{r}}}_0$. If $\psi$ is a velocity-potential amplitude instead, both leading expressions use the common prefactor ${\rho }_0\omega k/2$ in place of $k/(2{\rho }_0\omega )$. A normalized scalar current uses $k$ as the common prefactor. The paper does not specify this amplitude normalization; the ratio in part (c) is independent of that choice.

Use the same physical normalization for both fields. Part (b) has $E_s=E_i|f_{\infty }^{(1)}{|}^2/R^2+O(R^{-3})$, so the common energy flux prefactor cancels. The differential cross-section for scalar wave scattering isThe far-field pattern has dimensions of length in this normalization, so the differential cross-section for scalar wave scattering has dimensions of area per unit solid angle. Only the scattered field is used in this ratio; interference between the incident and scattered fields is a different contribution to the total field's energy flux.

Let the medium's autocorrelation function of a random field and power spectrum use the explicit conventionSince $n^2-1=2\mu W+{\mu }^2{W}^2$, only the linear term in $\mu$ is needed to obtain the cross-section through second order:Squaring and taking the expectation givesThis is the finite-volume formula in the Born approximation for scalar wave scattering; the zero mean of $W$ does not make its squared scattered amplitude vanish. For a jointly centered Gaussian random field, the third-order contribution is zero, so the next perturbative contribution is fourth order when that expansion is valid.

Introduce the overlap volume $H_V(\mathbf{s})=|V\cap (V-\mathbf{s})|$. Changing variables rewrites the double integral asWhen the linear dimensions of a regular scattering volume are large compared with the correlation length, and $C$ is integrable, $H_V(\mathbf{s})\simeq |V|$ over the significant covariance range. This gives the finite-volume random scattering spectrum resultThere are two distinct approximations here: the weak-scattering expansion in $\mu$ and the negligible-boundary-overlap approximation in the volume-to-correlation-scale ratio. If $S$ instead denotes the power spectrum of the index fluctuation $n-1=\mu W$, then $S={\mu }^2{S}_W$ and the external factor ${\mu }^2$ is absorbed into $S$.

Statistical homogeneity alone allows an anisotropic spectrum, so its argument is generally the vector ${\mathbf{k}}_s$. If statistical isotropy is additionally assumed, $S_W$ depends only onThe scalar notation $S(k_s)$ is valid in that isotropic case. For the alternative convention $S_{\mathrm{norm}}=(2\pi {)}^{-3}{S}_W$, the bulk formula is $2\pi {\mu }^2{k}^4|V|S_{\mathrm{norm}}({\mathbf{k}}_s)$. Changing the sign of the scattering vector does not change the spectrum of a real stationary scalar field; changing the Fourier transform normalization does change the displayed $2\pi$ factors. Both conventions have therefore been stated explicitly.

For noisy data, exact equality $Ax=y$ may be impossible: the data can have a component outside the operator range. A least-squares solution of a linear inverse problem instead minimizesFor any variation $h\in X$, the adjoint operator givesThus statistical homogeneity in every direction is the normal equation for a linear inverse problem, $A^{\ast }Ax=A^{\ast }y$. It is also sufficient for a global minimum: when it holds,The residual is orthogonal to the closure of the range, since it lies in $\ker A^{\ast }$. A normal-equation solution therefore fits the component of the data representable by the forward map, while leaving the perpendicular component unfitted. If there is more than one such solution, select the minimum-norm least-squares solution, which lies in $(\ker A{)}^{\perp }$.

These facts explain the formulation, but do not make the inverse problem well posed. A compact operator of infinite rank typically has nonclosed operator range, so even the least-squares infimum need not be attained. A least-squares solution exists exactly when the orthogonal projection of $y$ onto $\overline{\mathrm{Ran}A}$ belongs to $\mathrm{Ran}A$. In a singular system of a compact operator, this becomes the Picard criterion. Also, the nonzero eigenvalues of $A^{\ast }A$ are the squares of the singular values of $A$: solving the operator normal equation directly can worsen numerical conditioning. The operator normal equation expresses least squares; it does not itself regularize the inverse problem. The subsequent finite Landweber iteration and its stopping rule provide the inverse-problem regularization.

Choose a fixed positive Landweber relaxation parameter $\tau$ and rearrange the normal equation for a linear inverse problem as the fixed point equationSuccessive substitution gives Landweber iterationSet $B=A^{\ast }A$ and $P=I-\tau B$. Induction yields the closed polynomial-operator expressionFor $n=0$ the sum is empty, giving zero. The expression is well defined without an inverse of $A^{\ast }A$, which can have a kernel of a linear map and an unbounded inverse on its range. Each iterate lies in the closure of $\mathrm{Ran}{A}^{\ast }$, so initialization at zero eliminates the arbitrary null-space component of a solution. The unrelaxed choice $\tau =1$ gives $x_{n+1}=(I-A^{\ast }A)x_n+A^{\ast }y$ and the corresponding sum with $\tau =1$; it needs the convergence restriction in part (d).

Use the convention required by the printed expansion:where $u_i\in Y$ and $v_i\in X$ form the positive-singular-value parts of the singular system of a compact operator. Consequently $A^{\ast }A{v}_i={\sigma }_i^2{v}_i$. Apply part (b)'s polynomial to each such eigenvector, using noisy data:The geometric series then gives the full inverse coefficient, or Landweber full inverse filter,In particular, for the unrelaxed iteration $\tau =1$, $g_{\alpha }(\sigma )=[1-(1-{\sigma }^2{)}^{1/\alpha }]/\sigma$. The dimensionless Landweber spectral filter is instead $\sigma g_{\alpha }(\sigma )$; the printed expression uses the full coefficient and therefore needs the denominator $\sigma$.

For each fixed finite $n$, the filter extends continuously at zero by $g_{1/n}(0)=0$, since $g_{1/n}(\sigma )=n\tau \sigma +O({\sigma }^3)$ as $\sigma \to 0$. Small singular values are suppressed rather than immediately divided into noisy data. A component of $y^{(\delta )}$ in $\ker A^{\ast }$ contributes nothing. The parameter here is discrete: $1/\alpha$ is an integer. This matters when $1-\tau {\sigma }^2<0$; arbitrary real powers are not intended. A family for every $\alpha >0$ can use $n=\lfloor 1/\alpha \rfloor$ instead.

For nonzero $A$, a sufficient step-size condition for strong convergence of relaxed Landweber iteration isFor the unrelaxed iteration, this reads $\Vert A{\Vert }^2<2$. The data must also admit a minimum-norm least-squares solution $x^{†}$; for the exact problem, $y\in \mathrm{Ran}A$ is sufficient. More generally, the projected-data condition in part (a), or the Picard criterion ${\sum }_i|⟨y,u_i⟩{|}^2/{\sigma }_i^2<\infty$, ensures the required solution.

Every positive singular value satisfies $|1-\tau {\sigma }_i^2|<1$, so for fixed ${\sigma }_i>0$,This is the correct inverse coefficient. For exact data with a least-squares solution, the norm error isThe summable majorant is precisely the Picard criterion, so the dominated convergence theorem proves strong convergence, despite the absence of a uniform contraction constant when the singular values approach zero. Initialization at zero selects the solution orthogonal to $\ker A$. At $\tau =2/\Vert A{\Vert }^2$, a top singular mode of a nonzero compact operator has multiplier $-1$ and can oscillate indefinitely; the strict step bound avoids this failure.

Pointwise convergence of filters alone is not a full noisy-data inverse-problem regularization proof. For each finite $n$, put $t=\tau {\sigma }^2$; the step restriction gives $|1-t|\le 1$, andThis proves the uniform noise bound for relaxed Landweber iterationTogether with the exact-data convergence,whenever $n=n(\delta )\to \infty$ and $\delta \sqrt{n(\delta )}\to 0$. Equivalently, the a priori regularization parameter choice satisfies $\alpha (\delta )\to 0$ and $\delta /\sqrt{\alpha (\delta )}\to 0$. For example, $n(\delta )=\lfloor {\delta }^{-1}\rfloor$ works as $\delta \to 0$. The sharper familiar bound $\delta \sqrt{n\tau }$ is available when $\tau \Vert A{\Vert }^2\le 1$.

Thus finite iteration followed by an appropriate stopping rule is a convergent regularization of an inverse problem. Taking $n\to \infty$ for fixed noisy data is not the same operation and can amplify the small-singular-value noise without bound. If $A=0$, the zero starting iterate stays zero and is the minimum-norm least-squares solution for every datum; the step restriction is needed only for nonzero $A$.