Lenses and wave optics? We all have studied lenses in the framework of geometrical optics, but how do they work within wave optics? I figure that the topic is quite broad, but I would appreciate any hints, like which equations can be borrowed etc
 A: The speed of light is slower inside the glass. If you have a convex lens and some point source at one of its foci, this compensates for how much sooner each wave front arrives at the center than at the edges. "The light arrives sooner at the center, but has to take a longer path through the slow medium".
As a result, behind the lens you have not a spherical wave anymore but a plane wave.
That's actually pretty much it. If you want to do this for some general lenses, you mainly need to consider how much the light is "delayed" by the lens, depending on where it passes through it.
A: The thin lens formula carries over with a few approximations directly to wave optics. This formula, together with an appropriate analysis of diffraction, will let you get a first approximation to the behaviour of many systems.
Leftaroundabout's answer gives the overall picture. A good way to translate this into a wave picture is with the following assumptions:

*

*The paraxial approximation: i.e. all fields are a superposition of plane waves that are propagating at small angles to the optical axis;


*Gaussian scalar fields. These are fields that, on a transverse plane, vary like $\exp\left(-((x-x_0)^2+(y-y_0)^2)\left(\frac{1}{2\,\sigma^2} + \frac{i\,k,\kappa}{2}\right)\right)$ where $\sigma$ is the spotsize and $\kappa$ the wavefront curvature.
The "algorithms" are as follows.
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{align}\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(k-\sqrt{k^2 - k_x^2-k_y^2}\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{align}\tag{1}$$
In words:

*

*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)$;

*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(k-\sqrt{k^2 - k_x^2-k_y^2}\right) z\right)$;

*Propagate each such plane wave from the $z=0$ plane to the general $z$ plane using the plane wave solution noted in step 2;

*Inverse Fourier transform the propagated waves to reassemble the field at the general $z$ plane.

Now we make the paraxial approximation to the propagation relationship in step 2 above, i.e. we assume that the plane waves aren't skewed at too steep angles relative to the $z$ axis so that $k_x^2+k_y^2 \ll k^2$. Then our two propagation equations above become the Fresnel propagation integral:
$$\begin{align}\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 \frac{k_z^2+k_y^2}{2\,k} 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{align}\tag{2}$$
Now witness that a beam $\psi(x,y,0)$ with Gaussian dependence on $x$ and $y$ becomes, under the Fresnel diffraction integral, $\psi(x,y,z)$ with Gaussian dependence on $x$ and $y$. The Fourier transform of a Gaussian is a Gaussian, a fact which does not change when we multiply by the $\exp\left(i \frac{k_z^2+k_y^2}{2\,k} z\right)$ kernel in the Fresnel diffraction integral, and of course a Gaussian is recovered by the inverse Fourier transform.
Furthermore, thin lenses can be thought of simply as phase masks, i.e. a field whose transverse variation is $\psi(x,y)$ passing through them is transformed by:
$$\psi(x,y)\to \exp\left(-i\,k\,\frac{(x-x_0)^2+(y-y_0)^2}{2\,f}\right)\psi(x,y)\tag{3}$$
So the Gaussian form is preserved by all the diffraction and lensing operations.
Best of all, for Gaussian beams the diffraction integrals split up into a product of separate $x$ and $y$ dependences and these separate dependences are also operated on independently by the phase mask (3), so propagation analysis can be done as a product of two decoupled one-transverse dimensional diffraction problems.  Therefore the propagation through a system comprising homogeneous, diffractive mediums and thin lenses can be done wholly in closed form expressions (use Mathematica though!) and you end up with something slightly more general than the thin lens formula that becomes the thin lens formula when the axial distances involved become longer than the Rayleigh diffraction length (of the order of wavelengths).
This method assumes perfectly unaberrated waves. However, witness that any lens whose surface sag (height) near its vertex can be described by an analytic function of distance $r$ from the optical axis (axis of symmetry) will impart a phase mask well described by the above in the paraxial limit, i.e. by limiting the numerical aperture (maximum skew angle) in any field such that the support of any field in the transverse plane $z=0$ is small enough. As long as the width of the support needed to validate the analysis above stays big compared with a wavelength, the analysis above will be valid in the paraxial limit. This means it works in the paraxial limit for all lenses of practical curvatures in the neighbourhood of the optical axis.
