# [Physics] Plane wave focused by lens to a point

geometric-opticslensesoptics

I want to mathematically recover the situation in the following picture, that is a plane wave which is incident on a (thin) lens, such that the outgoing beam focuses to a finite spot size at a distance equal to the focal length of the lens.

I have my incident plane wave $$e^{ikz}$$ and I know that a thin lens causes a phase delay to the wave-front of $$e^{-\frac{ik}{2f}(x^2+y^2)}$$.

I have plot it on Mathematica, just looking at the $$x$$-axis so $$y = 0$$, I get a constant curvature in the correct direction but never a focus:

(the horizontal axis is $$z$$, the direction of propagation)

Just to be clear here I have plotted $$e^{ikz}\quad \mathrm{for} \quad z<0$$
and $$e^{ikz} \exp\left ({-\frac{k x^2}{\lambda f}} \right) \quad \mathrm{for} \quad z>0.$$

I checked for very large $$z$$ and it never goes to a spot.

Is this because I am using the thin lens approximation?

The equations you write are not those of a focussing wave (I think you're missing an $i\,\pi$ in your exponent for the $x$ variation for $z>0$). There is no way for the wavefront curvature to change with $z$ in your equations, thus no focussing. You need to model the effect of diffraction on your wavefront.

The easiest way to model diffraction is to assume a Gaussian intensity variation in the input, instead of a plane wave as you have done. You simply have to have the spotsize large enough to model the beamwidth you are dealing with. Then you impart the thin lens phase mask so that the field at $z=0$ which is immediately to the right of the lens output has the $x$ variation:

$$E(x,\,0) = \exp\left(-\frac{x^2}{2\,\sigma^2}\right) \,\exp\left(i\frac{k\, x^2}{2\,f}\right)\tag{1}$$

One can model the effect of diffraction on a field variation like this by taking heed of the following formula for the propagation of a generalized Gaussian beam in a homogeneous medium:

$$E(x,\,z) = \frac{1}{\sqrt{z-z_0 + i\,z_R}}\, \exp\left(-i \,k\, \frac{x^2}{2 \,(z-z_0 + i\,z_R)}\right)\tag{2}$$

where $z_R$ is the Rayleigh length for the beam. So your task is to find $z_0$ and $z_R$ in (2) to match (1) and then you can use (2) to propagate the Gaussian beam.

Note that the above is not the exact scalar diffraction operator; it makes a paraxial approximation that the transverse component $k_x$ of the wavevector is small compared to $k$; alternatively, that the beam's numerical aperture is small (less than about 0.3, depending on the accuracy you need).

Otherwise, you need to calculate the exact diffraction integral, which I outline below.

Full Diffraction Calculation

You begin with the Helmholtz equation in a homogeneous medium $(\nabla^2 + k^2)\psi = 0$. If the field comprises only plane waves in the positive $z$ direction then we can represent the diffraction of any scalar field on any transverse (of the form $z=c$) plane by:

$$\begin{array}{lcl}\psi(x,y,z) &=& \frac{1}{2\pi}\int_{\mathbb{R}^2} \left[\exp\left(i \left(k_x x + k_y y\right)\right) \exp\left(i \left(\sqrt{k^2 - k_x^2-k_y^2}-k\right) z\right)\,\Psi(k_x,k_y)\right]{\rm d} k_x {\rm d} k_y\\ \Psi(k_x,k_y)&=&\frac{1}{2\pi}\int_{\mathbb{R}^2} \exp\left(-i \left(k_x u + k_y v\right)\right)\,\psi(x,y,0)\,{\rm d} u\, {\rm d} v\end{array}$$

To understand this, let's put carefully into words the algorithmic steps encoded in these two equations:

1. Take the Fourier transform of the scalar field over a transverse plane to express it as a superposition of scalar plane waves $\psi_{k_x,k_y}(x,y,0) = \exp\left(i \left(k_x x + k_y y\right)\right)$ with superposition weights $\Psi(k_x,k_y)$;
2. Note that plane waves propagating in the $+z$ direction fulfilling the Helmholtz equation vary as $\psi_{k_x,k_y}(x,y,z) = \exp\left(i \left(k_x x + k_y y\right)\right) \exp\left(i \left(\sqrt{k^2 - k_x^2-k_y^2}-k\right) z\right)$;
3. Propagate each such plane wave from the $z=0$ plane to the general $z$ plane using the plane wave solution noted in step 2;
4. Inverse Fourier transform the propagated waves to reassemble the field at the general $z$ plane.

Here is a short, well tested Mathematica code of mine to implement the above. In the following $f$ is a square array of complex values of the input field, $d$ the axial ($z$) distance we wish to diffract the wave, $k$ the wavenumber and $Dx,\,Dy$ are the widths of the simulation domain in the $x$ and $y$ directions. It is easy to modify this code to cope with one transverse direction.

 Diffract[f_, d_, k_, Dx_, Dy_] /; If[MatrixQ[f], True, Message[Diffract::nnarg] False] := Module[{lenX, lenY, phase, mask, jx, jy}, ( lenX = Length[f]; lenY = Length[f[[1]]]; phase = Table[N[(jx - If[jx > lenX/2, lenX, 0])/Dx]^2 + N[(jy - If[jy > lenY/2, lenY, 0])/Dy]^2, {jx, 0, lenX - 1}, {jy, 0, lenY - 1}]; phase = k^2 - (4 Pi^2 phase); mask = Table[If[phase[[jx, jy]] < 0, 0, 1], {jx, 1, lenX}, {jy, 1, lenY}]; phase = Table[If[phase[[jx, jy]] < 0, 0, d (k - Sqrt[phase[[jx, jy]]])], {jx, 1, lenX}, {jy, 1, lenY}]; Return[InverseFourier[Fourier[f] Exp[-I phase] mask]]; );];