Fresnel approximation

In choosing a theory, one should pay attention to simplicity in hypotheses only. Simplicity in computation can be of no weight in the balance of probabilities. Nature is not embarrassed by difficulties of analysis. She avoids complication only in means. Nature seems to have proposed to do much with little: it is a principle that the development of physics constantly supports by new evidence.

–Augustin-Jean Fresnel

The study of focused sources motivates the Fresnel ("paraxial") approximation. The Fresnel approximation is not limited to focused sources, however, and it allows for analytical ease in the study of general diffraction phenomena. Nor is the Fresnel approximation limited in $z$: indeed, a far-field approximation of the Fresnel approximation can be taken.

Contents:

• Focused soruces
• Fresnel approximation
• Far field of Fresnel approximation
• Gaussian beams
• Circular piston
• Rectangular piston
• Derivation of parabolic equation

← Return to home

Focused sources

Consider a spherical wave converging at $z=d>0$. The pressure field is then $$\begin{array}{rl}{p}_{\omega }(x,y,z)& =A\frac{e^{-ikR}}{R}\\ & =A\frac{e^{ik\sqrt{x^2+y^2+(d-z{)}^2}}}{\sqrt{x^2+y^2+(d-z{)}^2}},\end{array}$$ where the dimensions of $A$ are pressure $\times$ distance. At $z=0$, the pressure field is $$\begin{array}{}\text{(1)}& p_{\omega }(x,y,0)& =A\frac{e^{-ik\sqrt{x^2+y^2+d^2}}}{\sqrt{x^2+y^2+d^2}}.\end{array}$$ In the absence of diffraction (ray theory), any finite source at $z=0$ with the above phasing will focus at $(x,y,z)=(0,0,d)$:

Let us assume that $p_0(x,y)$ or $u_0(x,y)\sim 0$ for $x^2+y^2\ge a^2$, where $a$ is the characteristic source radius, and that $a^2\ll d^2$. This is to say that the so-called $f$-number: $N=d/2a\sim 1/2\alpha$ so $N^2\gg 1$.

The phase in Eq. $\text{(1)}$ can be approximated by $$\begin{array}{rl}-ik\sqrt{x^2+y^2+d^2}& =-ikd(1+\frac{x^2+y^2}{d^2}{)}^{1/2}\\ & =-ikd+ik\frac{x^2+y^2}{2d}+\mathcal{O}(k{a}^2/d^3).\end{array}$$ Meanwhile, the amplitude is approximated as $$\begin{array}{r}\frac{A}{\sqrt{x^2+y^2+d^2}}\simeq \frac{A}{d},\end{array}$$ which has units of pressure (pressure $\times$ distance $÷$ distance $=$ pressure). So the source condition, Eq. $\text{(1)}$, becomes $$\begin{array}{rl}{p}_{\omega }(x,y,0)& =A\frac{e^{-ikd}}{d}{e}^{-ik(x^2+y^2)/2d}.\end{array}$$ Thus if $$\begin{array}{r}{u}_0(x,y),p_0(x,y)=\text{ unfocused source distribution,}\end{array}$$ focusing is achieved by multiplying by $e^{-ik(x^2+y^2)/2d}$, or in axisymmetric form, $e^{-ik(x^2+y^2)/2d}$. The factor $A{e}^{-ikd}/d$ can be neglected because phase is a relative quantity, and because $A/d$ is simply a pressure amplitude that can be included in the pressure source amplitude $p_0$ or the velocity source amplitude $u_0=p_0/{\rho }_0{c}_0$. Reintroducing the time dependence renders focusing as time-advancing: $$\begin{array}{r}{e}^{-i\omega t}{e}^{-ik{\rho }^2/2d}=e^{-i\omega (t+\rho /2c_{0}d)}.\end{array}$$ That is to say, focusing is achieved in the time domain by the transformation $$\begin{array}{r}t\mapsto t+\frac{{\rho }^2}{2c_{0}d},\end{array}$$ i.e., the further from the origin, the earlier the waveform must be launched.

Field in the focal plane of a focused source.

In the focal plane $z=d$, consider a focused velocity source $$\begin{array}{rl}u(x,y,0,t)& =u_0(x,y)e^{-ik(x^2+y^2)/2d}{e}^{-i\omega t}.\end{array}$$ The Rayleigh integral evaluated at the focal plane is $$\begin{array}{}\text{(R)}& p_{\omega }(x,y,d)& =-\frac{ik{\rho }_0{c}_0}{2\pi }{\iint }_{-\mathrm{\infty }}^{\mathrm{\infty }}{u}_0(x_0,y_0)e^{-ik(x_0^2+y_0^2)/2d}\frac{e^{ikR}}{R}d{x}_{0}d{y}_0,\end{array}$$ where $$\begin{array}{rl}ikR& =ik\sqrt{(x-x_0{)}^2+(y-y_0{)}^2+d^2}\\ & =ikd(1+\frac{x^2+y^2}{d^2}-2\frac{x{x}_0+y{y}_0}{d^2}+\frac{x_0^2+y_0^2}{d^2}{)}^{1/2}\\ & =ikd+ik\frac{x^2+y^2}{2d}-ik\frac{x{x}_0+y{y}_0}{d}+ik\frac{x_0^2+y_0^2}{2d}+\mathcal{O}(k{a}^2/d^3).\end{array}$$ In the amplitude, let $1/R\simeq 1/d$. Equation $\text{(R)}$ becomes $$\begin{array}{}\text{(2)}& p_{\omega }(x,y,d)& =-\frac{ik{\rho }_0{c}_0}{2\pi d}{e}^{ikd+ik(x^2+y^2)/2d}{\iint }_{-\mathrm{\infty }}^{\mathrm{\infty }}{u}_0(x_0,y_0)e^{-ikx{x}_0/d-iky{y}_0/d}d{x}_{0}d{y}_0,\end{array}$$ the magnitude of which is $$\begin{array}{}\text{(3)}& \boxed{{\displaystyle |p_{\omega }(x,y,d)|=\frac{k{\rho }_0{c}_0}{2\pi d}|U_0(kx/d,ky/d)|,}}\end{array}$$ or in axisymmetric form, since $U_0(\kappa )=2\pi U_{0H}(\kappa )$, $$\begin{array}{}\text{(3')}& \boxed{{\displaystyle |p_{\omega }(\rho ,d)|=\frac{k{\rho }_0{c}_0}{d}|U_{0H}(k\rho /d)|,}}\end{array}$$ That is to say, the field in the focal plane is but a spatial mapping of the far field into the focal plane: both are given in terms of the 2D spatial Fourier transform of the source condition.

Field at the the focal point of a focused source.

Further insight can be gained from assessing the field at the focal point $(x,y,z)=(0,0,d)$, Eq. $\text{(3)}$ becomes $$\begin{array}{r}|p_{\omega }(0,0,d)|=\frac{k{\rho }_0{c}_0}{2\pi d}Q\end{array}$$ Let $Q={\overline{u}}_{0}S$, where $S$ is the surface area of the source, and where the mean source particle velocity is $$\begin{array}{r}{\overline{u}}_0=\frac{Q}{S}=\frac{U_0(0,0)}{S}=\frac{1}{S}{\int }_S{u}_0(x,y)dS.\end{array}$$ Then, since $|p_{\omega }(0,0,d)|=\frac{k{\rho }_0{c}_0}{2\pi d}Q$, the ratio of the pressure at the focal point to that in a plane wave with mean source particle velocity is $$\begin{array}{r}\frac{|p_{\omega }(0,0,d)|}{{\rho }_0{c}_0{\overline{u}}_0}=\frac{kS}{2\pi d}=\frac{S}{\lambda d}.\end{array}$$ Define $$\begin{array}{r}\boxed{{\displaystyle G=\frac{S}{\lambda d}=\text{ focusing gain},}}\end{array}$$ which for a uniform circular piston of surface area $S=\pi a^2$ reduces to $$\begin{array}{r}G=\frac{k{a}^2}{2d}=\frac{z_0}{d}.\end{array}$$ Note: the geometric focus at $(x,y,z)=(0,0,d)$ is typically further beyond the location of the maximum axial amplitude.

The so-called ''spot size'' for a source of characteristic size $a$ corresponds to the $k$-space ''radius.'' From Eqs. $\text{(2)}$ or $\text{(3')}$, in the focal plane, this corresponds to $$\begin{array}{r}\frac{k\rho }{d}\sim \frac{1}{a}\end{array}$$ or $$\begin{array}{r}\rho \sim \frac{d}{ka}=a\frac{d}{k{a}^2}\sim \frac{a}{G}.\end{array}$$ Thus beam radius is reduced by $\sim 1/G$ in the focal plane.

Focal point in the time domain.

Recall Eq. $\text{(3)}$, copied below for convenience $$\begin{array}{r}|p_{\omega }(x,y,d)|=\frac{k{\rho }_0{c}_0}{2\pi d}|U_0(kx/d,ky/d)|.\end{array}$$ In the time domain, the field at the focal point is $$\begin{array}{rl}|p_{\omega }(0,0,0,t)|& =-i\omega \frac{{\rho }_{0}Q}{2\pi d}{e}^{i(kd-\omega t)}\\ & =\frac{{\rho }_0\dot{Q}(t-d/c_0)}{2\pi d}.\end{array}$$ The expansion of $ikR$ led to linear order led to the discarding of a term that was $\mathcal{O}(k{a}^4/d^3)$. So $$\begin{array}{rl}\frac{k{a}^4}{8d^3}& =\frac{1}{4}\frac{k{a}^2}{2d}\frac{a^2}{d^2}=\frac{G}{4}(\frac{a}{d}{)}^2.\end{array}$$ For diagnostic medical ultrasound, $G=4$, $a/d\simeq \frac{1}{8}$. So $$\begin{array}{r}\frac{G}{4}(\frac{a}{d}{)}^2\sim {10}^{-2}\ll 2\pi \end{array}$$ In Lucas and Muir [JASA 72, 1289 (1982)], $G\sim 40$ and $\frac{a}{d}\simeq \frac{1}{4}$, so $$\frac{G}{4}(\frac{a}{d})\sim 0.6\ll 2\pi .$$

Fresnel approximation

Begin with the Rayleigh integral $$\begin{array}{}\text{(1)}& p_{\omega }(x,y,z)=-\frac{ik{\rho }_0{c}_0}{2\pi }{\iint }_{-\mathrm{\infty }}^{\mathrm{\infty }}{u}_0(x_0,y_0)\frac{e^{ikR}}{R}d{x}_{0}d{y}_0.\end{array}$$ Now expand $R$ in powers of $1/z$ (rather than $1/r$, as was done in the Fraunhofer approximation): $$\begin{array}{rl}kR& =k\sqrt{(x-x_0{)}^2+(y-y_0{)}^2+z^2}\\ & =kz[1+\frac{(x-x_0{)}^2}{z^2}+\frac{(y-y_0{)}^2}{z^2}{]}^{1/2}\\ & =kz+\frac{k}{2z}[(x-x_0{)}^2+(y-y_0{)}^2]+\text{ H.O.T.}.\end{array}$$ Terminating the above at $\mathcal{O}(1/z)$ for the phase is less restrictive than the Fraunhofer approximation because we have retained the term $$\begin{array}{r}\frac{k(x_0^2+y_0^2)}{2z},\end{array}$$ i.e., it is not required that $k{a}^2/2z\ll 1$ as in the Fraunhofer approximation. Now the restriction appears to be $$\begin{array}{r}\frac{k}{8z^3}[(x-x_0{)}^2+(y-y_0{)}^2{]}^2\sim \frac{k{a}^4}{8z^3}\ll 2\pi \end{array}$$ or $$\begin{array}{r}\frac{z}{a}\gg \frac{ka}{16\pi }\end{array}$$ or $$\begin{array}{r}\boxed{{\displaystyle \frac{z}{a}\gtrsim (ka{)}^{1/3},}}\end{array}$$ though this restriction can be weaker as the main contribution to the integral can come from points $(x_0,y_0)\sim (x,y)$ due to phase variations. Substituting the approximation of $kR$ into $\text{(1)}$ gives $$\begin{array}{}\text{(2)}& \boxed{{\displaystyle p_{\omega }(x,y,z)=-\frac{ik{\rho }_0{c}_0}{2\pi }\frac{e^{ikz}}{z}{\iint }_{-\mathrm{\infty }}^{\mathrm{\infty }}{u}_0(x_0,y_0)e^{\frac{ik}{2z}[(x-x_0{)}^2+(y-y_0{)}^2]}d{x}_{0}d{y}_0.}}\end{array}$$ For axisymmetric sources, $u_0(x,y)=u_0(\rho )$, so $$\begin{array}{rl}(x-x_0{)}^2+(y-y_0{)}^2& =(\mathit{\rho }-{\mathit{\rho }}_0)\cdot (\mathit{\rho }-{\mathit{\rho }}_0)\\ & ={\rho }^2+{\rho }_0^2-2\rho {\rho }_0\cos {\varphi }_0\end{array}$$ for $\mathit{\rho }=\rho {\mathbf{e}}_x$. Thus $\varphi =0$, so Eq. $\text{(2)}$ becomes $$\begin{array}{rl}{p}_{\omega }(\rho ,z)& =-\frac{ik{\rho }_0{c}_0}{2\pi }\frac{e^{ikz}}{z}{e}^{ik{\rho }^2/2z}{\int }_0^{\mathrm{\infty }}{u}_0({\rho }_0)e^{ik{\rho }_0^2/2z}{\rho }_{0}d{\rho }_0{\int }_0^{2\pi }{e}^{-i(k\rho {\rho }_0/z)\cos {\varphi }_0}d{\varphi }_0.\end{array}$$ Taking the angular integral results in a Bessel function $2\pi J_0(k\rho {\rho }_0/z)$: $$\begin{array}{}\text{(3)}& p_{\omega }(\rho ,z)& =-ik{\rho }_0{c}_0\frac{e^{ikz}}{z}{e}^{ik{\rho }^2/2z}{\int }_0^{\mathrm{\infty }}{u}_0({\rho }_0)J_0(k\rho {\rho }_0/z)e^{ik{\rho }_0^2/2z}{\rho }_{0}d{\rho }_0.\end{array}$$ In general, $$\begin{array}{rl}{p}_{\omega }(x,y,z)& =q_{\omega }(x,y,z)e^{ikz},\end{array}$$ where $q_{\omega }$ varies slowly on the scale of a wavelength. In other words, if $p_{\omega }$ is a plane wave, then $q_{\omega }$ is a constant (the most slowly varying function).

Far field of Fresnel approximation

From Eq. $\text{(2)}$, for $z\gg z_0=\frac{1}{2}k{a}^2$, one obtains $$\begin{array}{rl}{p}_{\omega }(x,y,z)& =-\frac{ik{\rho }_0{c}_0}{2\pi }\frac{e^{ikz}}{z}{e}^{ik(x^2+y^2)/2z}{\iint }_{-\mathrm{\infty }}^{\mathrm{\infty }}{u}_0(x_0,y_0)e^{-i(kx{x}_0/z_0+ky{y}_0/z_0)}d{x}_{0}d{y}_0.\end{array}$$ Thus the far field of the Fresnel approximation is $$\begin{array}{}\text{(4)}& p_{\omega }(x,y,z)& =-\frac{ik{\rho }_0{c}_0}{2\pi }\frac{e^{ikz}}{z}{e}^{ik(x^2+y^2)/2z}{U}_0(kx/z,ky/z).\end{array}$$ For an axisymmetric beam, the far-field approximation of the Fresnel approximation reads $$\begin{array}{r}{p}_{\omega }(z,\theta )=-ik{\rho }_0{c}_0\frac{e^{ikz}}{z}{e}^{i\frac{1}{2}kz{\tan }^2\theta }{U}_H(k\tan \theta ).\end{array}$$ Compare Eq. $\text{(4)}$ with the Fraunhofer approximation in spherical coordinates, $$\begin{array}{}\text{(5)}& p_{\omega }(r,\theta ,\varphi )& =-\frac{ik{\rho }_0{c}_0}{2\pi }\frac{e^{ikr}}{r}{U}_0(k\alpha ,k\beta ),\end{array}$$ where $$\begin{array}{rl}\alpha & =x/r=\sin \theta \cos \varphi ,\\ \beta & =y/r=\sin \theta \sin \varphi ,\end{array}$$ where Eq. $\text{(5)}$ is exact for all $(\theta ,\varphi )$. Some important distinctions between Eqs. $\text{(4)}$ and $\text{(5)}$ are made:

• Wavefront curvature. From Eq. $\text{(5)}$, $$e^{ikr}=e^{ik\sqrt{x^2+y^2+z^2}}=e^{ikz[1+(x^2+y^2)/z^2]}$$ Only near the $z$ axis can we expand and write $$e^{ikr}=e^{ikz}{e}^{ik(x^2+y^2)/2z}.$$
• Angular dependence. From Eq. $\text{(5)}$, $$p_{\omega }(\theta ,\varphi )\propto |U_0(k\sin \theta \cos \varphi ,k\sin \theta \sin \varphi )|.$$ Only for small $\theta$ can the above equation be written as $$\begin{array}{rl}{p}_{\omega }(\theta ,\varphi )& \propto |U_0(k\tan \theta \cos \varphi ,k\tan \theta \sin \varphi )|\\ & =|U_0(kx/z,ky/z)|.\end{array}$$ Thus in the paraxial region, $$\begin{array}{r}\sin \theta \sim \tan \theta \sim \theta ,(\theta \lesssim {20}^{\circ }).\end{array}$$
• Range of validity in $z$. The Fresnel approximation is thus appropriate for sound beams for which $ka\gg 1$ and $z/a\gtrsim (ka{)}^{1/3}$, which is much closer to the source than in the Fraunhofer approximation.

Gaussian beams

The properties of Gaussian beams, which were introduced in the above examples, are now studied in more depth.

Unfocused Gaussian beam

As derived above, the unfocused Gaussian beam is given in cylindrical coordinates by $$\begin{array}{r}{p}_{\omega }(\rho ,z)=\frac{p_0{e}^{ikz}}{1+iz/z_0}\exp \{-\frac{{\rho }^2/a^2}{1+iz/z_0}\},\end{array}$$ which is $p_0(\rho )$ for $z=0$. Note that the solution is in the form $q_{\omega }(\rho ,z)=p_{\omega }(\rho ,z)e^{ikz}$, where $q_{\omega }(\rho ,z)$ is slowly varying over the scale of a wavelength. The magnitude of the Gaussian beam is $$\begin{array}{rl}|p_{\omega }(\rho ,z)|& =\frac{p_0{e}^{ikz}}{\sqrt{1+z/z_0}}\exp \{-\frac{{\rho }^2/a^2}{1+(z/z_0{)}^2}\}\\ & \simeq p_0{e}^{-{\rho }^2/a^2},z\ll z_0\\ & \simeq \frac{z_0}{z}{p}_0{e}^{-(z_0/z{)}^2{\rho }^2/a^2},z\gg z_0\\ & =\frac{z_0}{z}{p}_0{e}^{-(ka\tan \theta {)}^2/4},\tan \theta =\rho /z,\end{array}$$ and the amplitude profile as a function of $z$ is sketched below on a dB scale:

Focused Gaussian beam

Inclusion of focusing modifies the source condition to

$$\begin{array}{rl}{p}_0(\rho )& =p_0{e}^{-{\rho }^2/a^2}{e}^{ik{\rho }^2/2d}\\ & =p_0{e}^{-(1+ik{a}^2/2d){\rho }^2/a^2}\\ & =p_0{e}^{-(1+iG){\rho }^2/a^2},G=k{a}^2/2d\\ & =p_0{e}^{-{\rho }^2/{\tilde{a}}^2},{\tilde{a}}^2=\frac{a^2}{1+iG}.\end{array}$$ To obtain the focused Gaussian beam solution, simply replace $a^2$ by ${\tilde{a}}^2$ in the unfocused solution. The quantities ${\rho }^2/a^2$ and $iz/z_0$ become $$\begin{array}{rl}{\rho }^2/a^2& \mapsto p^2/a^2=(1+iG){\rho }^2/a^2\\ \frac{iz}{z_0}& \mapsto \frac{iz{z}_0}{k{\tilde{a}}^2}=i(1+iG)\frac{z}{z_0}=\frac{i-G}{G}\frac{z}{d}=\frac{iz}{Gd}-\frac{z}{d}.\end{array}$$ The pressure field is therefore $$\begin{array}{}\text{(G)}& p_{\omega }(\rho ,z)& =\frac{p_0{e}^{ikz}}{1-z/d+iz/Gd}\exp \{-\frac{(1+iG){\rho }^2/a^2}{1-z/d+iz/Gd}\}.\end{array}$$ The field in the focal plane is found by setting $z=d$ in Eq. $\text{(G)}$: $$\begin{array}{rl}{p}_{\omega }(\rho ,d)& =-iG{p}_0{e}^{-G^2{\rho }^2/a^2}{e}^{i(kd+G{\rho }^2/a^2)}⟹\frac{|p(\rho ,d)|}{p_0}=G{e}^{-G^2{\rho }^2/a^2}=G{e}^{-{\rho }^2/(a/G{)}^2}.\end{array}$$ The axial field is found by setting $z=0$ in Eq. $\text{(G)}$: $$\begin{array}{rl}{p}_{\omega }(0,z)& =\frac{p_0{e}^{i\varphi (z)}}{\sqrt{(1-z/d{)}^2+(z/d{)}^2/G^2}},\\ \text{where }\varphi (z)& =kz-\arctan \frac{z/d}{G(1-z/d)}.\end{array}$$ Sketches of the axial magnitude and phase of the focused Gaussian beam are shown below:

Circular piston

The uniform circular piston is now considered. It will be seen that the discontinuous edges of the piston introduce complications in the Fresnel approximation:

The velocity source condition for the uniform circular piston is $$\begin{array}{r}{u}_0(\rho )=\{\begin{array}{ll}{u}_0,& \rho \le a\\ 0,& \rho >a\end{array}.\end{array}$$ In this case, the Fresnel diffraction integral reads $$\begin{array}{rl}{p}_{\omega }(\rho ,z)& =-ik{\rho }_0{c}_0{u}_0\frac{e^{ikz}}{z}{e}^{ik{\rho }^2/2z}{\int }_0^{\mathrm{\infty }}{J}_0(k\rho {\rho }_0/z)e^{ik{\rho }_0^2/2z}{\rho }_{0}d{\rho }_0.\end{array}$$ For the axial field, $\rho =0$, so the Bessel function $\to 1$. Meanwhile, let $x\equiv ik{\rho }_0^2/2z$ and ${\rho }_{0}d{\rho }_0=(z/ik)dy$. The integral becomes $$\begin{array}{}\text{(1)}& p_{\omega }(0,z)& =-{\rho }_0{c}_0{u}_0{e}^{ikz}{\int }_0^{ik{a}^2/2z}{e}^{x}dx={\rho }_0{c}_0{u}_0{e}^{ikz}(1-e^{ik{a}^2/2z}).\end{array}$$ Compare Eq. $\text{(1)}$ to the exact solution of the Helmholtz equation derived by Rayleigh for the (exact) axial field of a uniform circular piston: $$\begin{array}{rl}{p}_{\omega }(0,z)& ={\rho }_0{c}_0{u}_0(e^{ikz}-e^{ik\sqrt{z^2+a^2}})\\ \text{(2)}& & ={\rho }_0{c}_0{u}_0{e}^{ikz}[1-e^{ik(\sqrt{z^2+a^2}-z)}]\end{array}$$ Now expand Eq. $\text{(2)}$ in powers of $1/z$: $$\begin{array}{rl}k(\sqrt{z^2+a^2}-z)& =kz[(1+\frac{a}{z}{)}^{1/2}-1]\\ & =kz[(1+\frac{a^2}{2z^2}-\frac{a^4}{8z^4}+\dots )-1]\\ & =\frac{k{a}^2}{2z}-\frac{k{a}^4}{8z^3}+\mathcal{O}(1/z^5).\end{array}$$ Thus Eq. $\text{(1)}$ $\simeq$ $\text{(2)}$ for $$\begin{array}{rl}\frac{k{a}^4}{8z^3}& \ll 2\pi \\ (z/a{)}^3& \gg ka/16\pi \\ z/a& \gtrsim (ka{)}^{1/3},\end{array}$$ which is the very condition of the Fresnel approximation. This suggests that the Fresnel approximation can be taken a priori (before any calculations have been made made), or a posteriori (after an exact result has been found).

Note that Eq. $\text{(2)}$ is singular $z=0$, because the complex exponential as its argument approaches $\mathrm{\infty }$ executes circles in the complex plane: $$\begin{array}{rl}\frac{p_{\omega }(0,z)}{{\rho }_0{c}_0{u}_0}& =e^{ikz}{e}^{i{z}_0/z}(e^{-i{z}_0/2z}-e^{i{z}_0/2z}),z_0=k{a}^2/2\\ & =-i2e^{i(kz+z_0/2z)}\sin (z_0/2z).\end{array}$$ So $$\begin{array}{rl}\frac{p_{\omega }(0,z)}{{\rho }_0{c}_0{u}_0}& =2|\sin (z_0/2z)|=z_0/z,z\gg z_0.\end{array}$$ The far field starts at $z/z_0=1/\pi$.

Meanwhile, the far field of the Fresnel approximation is given by $$\begin{array}{rl}{p}_{\omega }(z,\theta )& =-ik{\rho }_0{c}_0\frac{e^{ikz}}{z}{e}^{ik{\rho }^2/2z}{U}_H(k\tan \theta ).\end{array}$$ In this case, for the the uniform circular piston $$\begin{array}{r}{U}_H(\kappa )=\mathcal{H}\{u_0(\rho )\}=u_0\mathcal{H}\{\mathrm{circ}(\rho /a)\}=u_0\frac{a^2}{2}\frac{2J_1(\kappa a)}{\kappa a}.\end{array}$$ The far field of the Fresnel approximation for a uniform circular piston is $$\begin{array}{rl}\frac{p_{\omega }(z,\theta )}{{\rho }_0{c}_0{u}_0}& =-i\frac{z_0}{z}{e}^{ikz}{e}^{i(kz/2){\tan }^2\theta }\frac{2J_1(ka\tan \theta )}{ka\tan \theta }.\end{array}$$ Meanwhile, for the sake of comparison, note that the Fraunhofer approximation of the uniform circular piston is given in spherical coordinates by $$\begin{array}{rl}\frac{p_{\omega }(r,\theta )}{{\rho }_0{c}_0{u}_0}& =-i\frac{z_0}{r}{e}^{ikr}\frac{2J_1(ka\sin \theta )}{ka\sin \theta }\end{array}$$

See here for Dr. Hamilton's handwritten notes on the axial pressure of the focused circular piston. This was not covered formally in class, but is uploaded here for reference.

Rectangular piston

The field due to a rectangular piston is now calculated in the paraxial approximation.

Since the source pressure in the paraxial approximation is related to particle velocity through the plane wave impedance relation, $p={\rho }_0{c}_{0}u$, the following source condition is considered: $$\begin{array}{r}{p}_0(x,y)=p_0,x\in [-a,a],y\in [-b,b].\end{array}$$ The field is given in Cartesian coordinates by $$\begin{array}{rl}{p}_{\omega }(x,y,z)& =-\frac{ik{p}_0}{2\pi }\frac{e^{ikz}}{z}{I}_x{I}_y\end{array}$$ where $$\begin{array}{rl}{I}_x& ={\int }_{-a}^a{e}^{ik(x-x_0{)}^2/2z}d{x}_0,\\ I_y& ={\int }_{-b}^b{e}^{ik(y-y_0{)}^2/2z}d{y}_0.\end{array}$$ To take the integrals above, let $$\begin{array}{rl}{t}^2& =-\frac{ik}{2z}(x-x_0{)}^2=\frac{k}{i2z}(x-x_0{)}^2\\ ⟹x-x_0& =\sqrt{\frac{i2z}{k}}t,d{x}_0=-\sqrt{\frac{i2z}{k}}dt.\end{array}$$ The $x$ integral becomes $$\begin{array}{rl}{I}_x& =-\frac{i2z}{k}{\int }_{\sqrt{k/i2z}(x+a)}^{\sqrt{k/i2z}(x-a)}{e}^{-t^2}dt.\end{array}$$ The integral evaluates to $$\begin{array}{rl}{\int }_{\alpha }^{\beta }{e}^{-t^2}dt& =\{{\int }_0^{\beta }-{\int }_0^{\alpha }\}{e}^{-t^2}dt\\ & =\frac{\sqrt{\pi }}{2}(\mathrm{erf}\beta -\mathrm{erf}\alpha ),\end{array}$$ where $$\begin{array}{rl}\mathrm{erf}z& =\frac{2}{\sqrt{\pi }}{\int }_0^z{e}^{-t^2}dt=\text{error function}\\ & \to 1,z\to \mathrm{\infty }.\end{array}$$ Substituting in the appropriate limits of integration, $I_x$ becomes $$\begin{array}{rl}{I}_x& =-\sqrt{\frac{i2z}{k}}\frac{\sqrt{\pi }}{2}\{\mathrm{erf}[\sqrt{\frac{k}{i2z}}(x-a)]-\mathrm{erf}[\sqrt{\frac{k}{i2z}}(x+a)]\}.\end{array}$$ Since the error function is odd, i.e., $$\begin{array}{}\text{(}\star \text{)}& \mathrm{erf}(-z)=-\mathrm{erf}(z),\end{array}$$ $I_x$ can be written as $$\begin{array}{rl}{I}_x& =\frac{1}{2}\sqrt{\frac{i2\pi z}{k}}\{\mathrm{erf}[\sqrt{\frac{z_a}{iz}}(1+x/a)]+\mathrm{erf}[\sqrt{\frac{z_a}{iz}}(1-x/a)]\},\end{array}$$ where $z_a=k{a}^2/2$ and $z_b=k{b}^2/2$. The expression for $I_y$ is very similar.

The pressure field is therefore $$\begin{array}{rl}{p}_{\omega }(x,y,z)& =\frac{1}{4}{p}_0{e}^{ikz}\{\mathrm{erf}[\sqrt{\frac{z_a}{iz}}(1+x/a)]+\mathrm{erf}[\sqrt{\frac{z_a}{iz}}(1-x/a)]\}\\ \text{(1)}& & \times \{\mathrm{erf}[\sqrt{\frac{z_b}{iz}}(1+y/b)]+\mathrm{erf}[\sqrt{\frac{z_b}{iz}}(1-y/b)]\}.\end{array}$$ The axial field is found by setting $x=y=0$: $$\begin{array}{}\text{(2)}& p_{\omega }(0,0,z)& =p_0{e}^{ikz}\mathrm{erf}\sqrt{\frac{z_a}{iz}}\mathrm{erf}\sqrt{\frac{z_b}{iz}}.\end{array}$$ Note that $$\begin{array}{r}\sqrt{\frac{z_a}{iz}}=\pm (1-i)\sqrt{\frac{z_a}{2z}},\end{array}$$ but either sign gives the same result in Eqs. $\text{(1)}$ and $\text{(2)}$ by Eq. $\text{(}\star \text{)}$.

Derivation of parabolic equation (paraxial equation)

The Fresnel approximation is now considered from the perspective of partial differential equations. Consider the Helmholtz equation: $$\begin{array}{}\text{(1)}& {\mathrm{\nabla }}^2{p}_{\omega }+k^2{p}_{\omega }=0.\end{array}$$ Let $$\begin{array}{}\text{(2)}& q(x,y,z)=q_{\omega }(x,y,z)e^{ikz}.\end{array}$$ $q$ is slowly varying in $z$ relative to the wavelength. Equation $\text{(2)}$ is now inserted into Eq. $\text{(1)}$. The appropriate derivatives are first taken: $$\begin{array}{rl}\frac{\partial p_{\omega }}{\partial z}& =(\frac{\partial q_{\omega }}{\partial z}+ik{q}_{\omega })e^{ikz}\\ \frac{{\partial }^2{p}_{\omega }}{\partial z^2}& =(\frac{{\partial }^2{q}_{\omega }}{\partial z^2}+i2k\frac{\partial q_{\omega }}{\partial z}-k^2{q}_{\omega })e^{ikz}.\end{array}$$ Equation $\text{(1)}$ becomes $$\begin{array}{}\text{(1')}& \frac{{\partial }^2{q}_{\omega }}{\partial z^2}+i2k\frac{\partial q_{\omega }}{\partial z}+{\mathrm{\nabla }}_{\perp }^2{q}_{\omega }=0,\end{array}$$ where ${\mathrm{\nabla }}_{\perp }^2$ is the Laplacian in the transverse direction, e.g., ${\partial }^2/\partial x^2+{\partial }^2/\partial y^2$ in Cartesian coordinates. Note that $$\begin{array}{r}\frac{{\partial }^2{q}_{\omega }/\partial z^2}{i2k\partial q_{\omega }/\partial z}\sim \frac{q_{\omega }/z_0^2}{k{q}_{\omega }/z_0}\sim \frac{1}{k{z}_0}\sim \frac{1}{(ka{)}^2}\ll 1.\end{array}$$

Thus Eq. $\text{(1')}$ is approximated by $$\begin{array}{}\text{(3)}& \boxed{{\displaystyle i2k\frac{\partial q_{\omega }}{\partial z}+{\mathrm{\nabla }}_{\perp }^2{q}_{\omega }=0.}}\end{array}$$ Equation $\text{(3)}$ is called the parabolic approximation of the Helmholtz equation. The parabolic equation is first order in $z$, reducing the elliptic equation (Helmholtz equation) to a parabolic equation.

Note that in the time domain, the parabolic equation is $$\begin{array}{r}\frac{{\partial }^{2}p}{\partial z\partial \tau }=\frac{c_0}{2}{\mathrm{\nabla }}_{\perp }^{2}p,\end{array}$$ where $p=p(x,y,z,\tau )$, and $\tau =t-z/c_0$.

To solve Eq. $\text{(3)}$, one can use the standard Fourier acoustics procedure. See Dr. Hamilton's notes here in which the Fresnel diffraction integral is recovered.