I know from the angular spectrum method that given a field $U$ with a wavelength $\lambda$, we can decompose it with Fourier transform.
\begin{equation} U(x, y,0) = \iint{\tilde{U}_0(f_x,f_y)} e^{2 \pi i (f_x x + f_y y)} d f_x d f_y \end{equation}
And propagate it along a distance $z$
\begin{equation} U(x, y, z) = \iint{\tilde{U}_0(f_x,f_y)} e^{2 \pi i (f_x x + f_y y + {z \sqrt{\frac{1}{\lambda^2} - f^2_x - f^2_y}})} d f_x d f_y \end{equation}
When I have $\frac{1}{\lambda^2} - f^2_x - f^2_y <0$ I will have evanescent fields
\begin{equation} e^{ - 2 \pi { z \sqrt{ \bigl| \frac{1}{\lambda^2} - f^2_x - f^2_y}\ \bigr|}} \end{equation}
So the resolution is limited to $\lambda/2$, and this can also be applied for free propagation.
But by studying resolution with an imaging lens (and so the point spread function), I know that the resolution is limited to $ \frac{\lambda}{2 NA}$. The numerical aperture (NA) can be greater than 1, so I can have a better resolution than $\lambda/2$.
How can this happen?