Deriving the non-paraxial form of Rayleigh criterion

Background

The Rayleigh criterion of imaging resolution says that two incoherent point sources are barely resolved by a diffraction-limited system with a circular aperture where the center of the Airy intensity pattern generated by one point source falls exactly on the first zero of the Airy pattern generated by the second. Usually the Airy pattern of a circular aperture is derived as the Fraunhofer diffraction of the aperture, with the following result:

$I(r) = left(frac{A}{lambda z}right)^2 left[2 frac{J_1(kwr/z)}{kwr/z} right]^2$

where $A$ is the area of the aperture, $lambda$ and $k$ are wavelength and wavenumber respectively, $w$ is the radius of the aperture, $r$ is the radius coordinate of the observation plane, and $z$ is the distance between the aperture and the observation plane. From this expression, one can derive that the Rayleigh resolution is

$delta = 0.61 frac{lambda z}{w}$

Question

Fraunhofer diffraction or Fourier optics in general is based on the paraxial approximation of the Kirchhoff diffraction formula. So you cannot trust the expression above to be valid in the non-paraxial situation often encountered in microscopy. In fact, “Introduction to Fourier Optics” by Goodman says (pg 158) the corresponding result in the nonparaxial case can be shown to be

$delta = 0.61 frac{lambda}{sintheta} = 0.61 frac{lambda}{NA}$

where $NA$ is the numerical aperture. But the textbook does not provide any reference for this derivation. If you look at Wikipedia article on Rayleigh criterion, you will see that $z/w$ is related to F-number, and $text{F-number}approx 1/(2NA)$, but again that relation only holds up on the paraxial case. I would like to see a full derivation of the non-paraxial expression for the Rayleigh resolution.

Personal attempt at solution

Often, a high NA objective is paired with a long focal length tube lens. So when we measure the point spread function (PSF) of a point source with this (objective + tube lens) system, we are measuring the magnified version of the true PSF. You can say that in most cases, the paraxial approximation is valid on the tube lens side. That is, the PSF we measure is described by Rayleigh resolution of $0.61 frac{lambda}{NA_{text{imaging space}}}$. If we can say that magnification of the imaging system is $M = frac{NA_{text{objective space}}}{NA_{text{imaging space}}}$, then we can say that the Rayleigh resolution of the high NA objective is indeed $0.61 frac{lambda}{NA}$. The fact that the magnification is described as a ratio of two numerical apertures is described by the Abbe sine condition, which must be obeyed by diffraction-limited imaging system.

I think the above argument is sound, but it relies on the fact that you use a long focal lens tube lens paired with an infinity-corrected objective. I would like to see a more general derivation.