This problem is intimately bound up with tunnelling (wholly analogous to quantum tunnelling) and evanescent waves. The results you get from analyses like those below depend critically on how one assumes the hole interacts with the incident field.

Let us first that the screen is perfectly reflecting, or infinitely absorbing, so that, even though the screen may be infinitely thin, on the output side of the hole the boundary conditions are (we use polar co-ordinates in the transverse plane):

$$E_x(r,\,\phi) = \left\{\begin{array}{lcl}E_0\,\hat{X}&&r < \frac{d}{2}\\0&&\text{elsewhere}\end{array}\right.$$

i.e. we simply assume a straight screening off of the field and there is negligible reaction field (e.g. generated in a conductive screen).

What happens here if you analyse the problem is detail (see my appendix) is that the form of the farfield variation with position is almost exactly what Huygens's principle would tell us: the small hole looks exactly like a point source (with a dipole radiation field shape). What Huygens's principle will not tell us, however, is the scaling constant for this variation. The amplitude of the farfield (if I've done a stationary phase integral right) is

$$\psi(R,\,\theta) = i\,E_0\,\frac{k\,d^2}{8\,R}\, \cos\theta \,e^{i\,k\,R}$$

where $R$ is the distance from the hole, $\theta$ the colatitude (spherical $\theta$ co-ordinate) and $\lambda$ the wavelength.

Huygens's principle would simply assume that all of the power incident on the hole would reach the farfield as propagating waves, and that we have a point source. Without Kirchoff's famous obliquity factor, therefore, Huygens's principle would foretell the following scalar field:

$$\psi(R,\,\theta) = E_0\,\frac{d}{\sqrt{8}\,R}\,e^{i\,k\,R}$$

the above gotten simply by equating the incident power (equal to $\mathcal{Z}_0\,E_0^2\, \pi\,d^2/4$, where $\mathcal{Z}_0$ is the characteristic impedance of freespace) and the total farfield power (equal to $\mathcal{Z}_0\,E_{far}^2\, 2\,\pi\,R^2$) and solving for $E_{far}$. Even if we put in the dipolar dependence (the $\cos\theta$ term) ad-hoc, Huygens's principle does not give the right wavelength dependence. What's going on?

Let's first look at how valid Huygens's principle is. The rigorous grounding of this principle is the Kirchoff diffraction integral (Born and Wolf, "Principles of Optics", Sixth Edition, §8.3)

$$ \psi\left(\mathbf{r}\right) = \frac{1}{4\pi} \int_{\partial V} \left(\psi\left(\mathbf{u}\right) \nabla \frac{e^{\pm i\,k\,\left|\mathbf{u} -\mathbf{r}\right|}}{\left|\mathbf{u} -\mathbf{r}\right|} - \frac{e^{\pm i\,k\,\left|\mathbf{u} -\mathbf{r}\right|}}{\left|\mathbf{u} -\mathbf{r}\right|}\nabla \psi\left(\mathbf{u}\right) \right) \cdot \hat{\mathbf{n}}\left(\mathbf{u}\right) \mathrm{d}^2 u $$

for the scalar field $\psi\left(\mathbf{r}\right)$ within a volume $V$ in terms of its values and normal derivatives $\nabla \psi\left(\mathbf{u}\right) \cdot \hat{\mathbf{n}}\left(\mathbf{u}\right)$ on the boundary $\partial\,V$ of this volume. We need to assume reasonable values for the boundary: Huygens's principle can be shown (in Born and Wolf §8.3 for example) to result when we assume the values of the field and its normal derivatives at holes in a screen are equal to the values that would prevail without the screen there. If the holes are of a reasonable size, this is a good approximation. However, amongst other things, the introduction of the screen distorts the values of the normal derivatives severely at the edges of holes. This is not a problem for big holes: the error term is dwarfed by the contribution from the hole itself. However, this assumption most certainly goes awry for subwavelength holes.

So Huygens, whilst telling is the general character of the field passing through the hole, cannot tell us the fine details. Indeed, not all the power incident on the hole passes through it. In a scalar theory, the reason for this is the evanescent field. A full Maxwell theory for e.g. a perfectly conducting screen, also shows the presence of the evanescent field, but furthermore that this evanescent field is set up as a result of a back reflexion of some of the power incident on the hole. A nice Maxwellian analysis to give intuition for what is going on here is of a perfectly conducting screen, with a hole of a nonzero thickness. In this case, the hole is like a cylindrical waveguide with conducting boundaries. Such a waveguide can cut off all its modes and make them evanescent. So, propagating power is output from the other side of the screen by dint of tunnelling: enough of the evanescent mode reaches the other side to become a propagating wave again. The modes of a cylindrical waveguide vary like:

$$\psi_k(r,\,z) = J_0(\omega_n\,\frac{2\,r}{d})\,\exp\left(i\,\sqrt{k^2-4\frac{\omega_n^2}{d^2}}\,z\right) = J_0(\omega_n\,\frac{2\,r}{d})\,\exp\left(-\sqrt{4\frac{\omega_n^2}{d^2}-k^2}\,z\right)\approx J_0(\omega_n\,\frac{2\,r}{d})\,\exp\left(-2\,\frac{\omega_n\,z}{d}\right)$$

since $d<<\lambda$. Here $\omega_n$ is the $n^{th}$ zero of the Bessel function $J_0$: the slowest dwindling field is the one corresponding to $\omega_1 \approx 2.405$. So now we have an expression for how much power gets through a small hole as a function of how thick the conductive screen is: it is:

$$P\propto \exp\left(-2\,\sqrt{4\frac{\omega_0^2}{d^2}-k^2}\,t\right)=\exp\left(-4\,\frac{\omega_0\,z}{d}\right)$$

where $t$ is the screen's thickness. This formula lets us see how thick we need to make a mesh in a microwave oven, for example, so as to quell the microwaves within the oven to a safe level for the hungry animal outside gawking in there for his dinner to be ready (and in so doing, exposing the most microwave-vulnerable part of his body - the eyes - to the show).

#Evanescent Waves in Detail

What are these weird creatures, these evanescent waves? Let's go back to the perfectly absorbing, zero thickness screen.

The Cartesian components of our EM field fulfil the Helmholtz equation, therefore, if the field comprises only plane waves in the positive $z$ direction, our algorithm for analysis of propagation is:

$$\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 \sqrt{k^2 - k_x^2-k_y^2}\, 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}\tag{1}$$

This looks fearsome, but is simply the use of the Fourier transform to break the field up into plane waves, propagate those plane waves (as known solutions to the Helmholtz equation) and then re-assemble the propagated plane waves at any $z$-plane we want. I.e. we are simply using the linear superposition principle for the linear Helmholtz equation. In more detail:

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(k-\sqrt{k^2 - k_x^2-k_y^2}\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.

In our simple, circularly symmetric problem, we can use zeroth order Hankel transforms instead: the "prototype" radially symmetric solution to Helmholtz's equation is:

$$\psi_k(r,\,z) = J_0(k_\perp\,r)\,\exp\left(i\,\sqrt{k^2-k_\perp^2}\,z\right)\tag{2}$$

where $k_\perp$ is the transverse component of the wavevector (this "prototype" solution is actually the sum of all the plane waves with different $k_x,\,k_y$ in (1) whose transverse wavevector is $k_\perp=\sqrt{k_x^2+k_y^2}$). So (1) becomes:

$$\begin{array}{lcl}\psi(r,z) &=& \int_0^\infty k_\perp\,J_0(k_\perp\,r)\,\Psi(k_\perp)\,\exp\left(i\,\sqrt{k^2-k_\perp^2}\,z\right){\rm d} k_\perp\\ \Psi(k_\perp)&=&\int_0^\infty r\,J_0(k_\perp\,r)\,\psi(r,\,0)\,{\rm d}r \\&=&\int_0^{\frac{d}{2}} r\,J_0(k_\perp\,r)\,{\rm d}r\\ &=&\frac{d}{2\, k_\perp}\,J_1\left(\frac{d\, k_\perp}{2}\right)\end{array}\tag{3}$$

Now of course the weight function $\Psi(k_\perp) = \frac{d}{2\, k_\perp}\,J_1\left(\frac{d\, k_\perp}{2}\right)$ does not have a compact support: it is nonzero for values of $k_\perp>k$. The superposition components for which hese of course lead to a $z$ variation in our prototype Helmholtz equation solution (2) and in the first equation of (3) given by $\exp\left(i\,\sqrt{k^2-k_\perp^2}\,z\right) = \exp(-\sqrt{k_\perp^2-k}\,z)$ which are thus swiftly attenuated with increasing $z$. These are the evanescent waves, they represent inductive and capacitive nonpropagating stores of energy that stay "bound" to the screen. They give rise to the so-called "near field" (warning: take heed that I use this term to mean strictly the evanescent field: this term is used differently in connexion with Fresnel refraction, which deals wholly with propagating, non-evanescent fields).

However, there is still a propagating, free photon component of the field left over when $z\to\infty$ and the evanescent fields have vanished:

$$\begin{array}{lcl}\psi(r,z) &=& \frac{d}{2}\,\int_0^k J_0(k_\perp\,r)\,J_1\left(\frac{d\, k_\perp}{2}\right)\,\exp\left(i\,\sqrt{k^2-k_\perp^2}\,z\right){\rm d} k_\perp\\&\approx& \frac{d^2}{16}\,\int_0^k\,k_\perp\, J_0(k_\perp\,r)\,\exp\left(i\,\sqrt{k^2-k_\perp^2}\,z\right){\rm d} k_\perp\\ &=& \frac{d^2}{32\,\pi}\,\int_0^k\,k_\perp\, \int_0^{2\pi}\exp\left(i\,k_\perp\,r\,\cos\phi+\,i\,\sqrt{k^2-k_\perp^2}\,z\right){\rm d}\phi\,{\rm d} k_\perp\end{array}\tag{4}$$

where, in the last line, I have converted the Bessel function by its integral definition $J_0(z) = \frac{1}{2\pi}\int_0^{2\pi}\exp(i\,z\,\cos\phi){\rm d}\phi$, so that I can bring the method of stationary phase to bear (we recall that the Bessel function is is a highly oscillatory function of $r$). The integrand has stationary phase when $(k_\perp,\,\phi) = (k\, r/\sqrt{r^2+z^2},\,0)$ and $(k_\perp,\,\phi) = (-k\, r/\sqrt{r^2+z^2},\,\pi)$ and which lets us approximate:

$$\psi(r,z)\approx i\,\frac{k\,d^2}{8\,R}\, \cos\theta \,e^{i\,k\,R}\tag{5}$$

In essence, the form of this variation is what we would get from Huygens's principle, but Huygens's principle does not give us the amplitude of the farfield.