Yes there is. For this, it is best to go back to the origin of the Fresnel propagator. A possible starting point is the Helmholtz equation for stationary waves in a scalar theory of light:
$$
\partial_z^2\phi+\Delta_{2D}\phi +\omega^2\phi = 0
$$
In the paraxial approximation, you recover Schrödinger's equation:
$$
i2\omega\partial_z\phi+\Delta_{2D}\phi = 0
$$
The Fresnel propagator is just the Green's function of the equation in the initial value, i.e. it satisfies:
$$
i2\omega\partial_zG+\Delta_{2D}G = 0 \\
G(z=0) = \delta
$$
You typically compute it by going in Fourier space:
$$
\begin{align}
G &= \int e^{-i(k^2+l^2)z/2\omega+i(kx+ly)}\frac{dkdl}{(2\pi)^2} \\
&= \left(\int e^{-ik^2z/2\omega+ikx}\frac{dk}{2\pi}\right)\left(\int e^{-il^2z/2\omega+ily}\frac{dl}{2\pi}\right) \\
&= \frac{\omega}{2\pi iz}e^{ix^2\omega/2z+iy^2\omega/2z}
\end{align}
$$
The unitarity of the transfer function translates as:
$$
\int G(x,y)^*G(x+\xi,y+\eta)dxdy = \delta(\xi)\delta(\eta)
$$
The first step is discretise space. The simplest approach would be to let $x,y\in a\mathbb Z\times b\mathbb Z$ on a lattice and replace the laplacian by a discrete one like:
$$
\Delta_{2D}\phi \to \frac{\phi(x+a,y)+\phi(x-a,y)-2\phi(x,y)}{a^2}+\frac{\phi(x,y+b)+\phi(x,y-b)-2\phi(x,y)}{b^2}
$$
In Fourier space, you therefore get a finite Brillouin zone $k,l\in(\mathbb R/(2\pi/a)\mathbb Z)\times (\mathbb R/(2\pi/b)\mathbb Z)$ and you just need to replace in the exponential:
$$
k^2+l^2 \to \frac{2(1-\cos(k a))}{a^2}+\frac{2(1-\cos(lb))}{b^2}
$$
Note that you do recover the continuous limit when $a,b\to0$. Back in real space, if you let $x_m = am,y = bn$ you get:
$$
\begin{align}
G_{mn} &= \int \exp\left(\frac{iz}{2\omega}\left(\frac{2(1-\cos(k a))}{a^2}+\frac{2(1-\cos(lb))}{b^2}\right)+i(kx_m+ly_n)\right)\frac{dkdl}{(2\pi)^2} \\
&= \left(\int \exp\left(\frac{iz}{\omega}\frac{1-\cos(ka)}{a^2}+ikx_m\right)\frac{dk}{2\pi}\right)\left(\int \exp\left(\frac{iz}{\omega}\frac{1-\cos(lb)}{b^2}+ily_n\right)\frac{dl}{2\pi}\right) \\
&= \frac{(-i)^{m+n}}{ab}\exp\left(\frac{iz}{\omega a^2}+\frac{iz}{\omega b^2}\right)J_m\left(\frac z{\omega a^2}\right)J_n\left(\frac z{\omega b^2}\right)
\end{align}
$$
The Green's function is still unitary because the relevant hermitian product that has a continuum limit (to recover the Dirac deltas) is:
$$
\langle \phi,\psi\rangle = \sum \phi_{mn}^*\psi_{mn} ab
$$
so unitarity translates as:
$$
\sum G_{mn}^*G_{m+m',n+n'}ab = \frac{\delta_{m'}}a\frac{\delta_{n'}}b
$$
If you want a finite spatial domain then the simplest would be to impose periodic boundary conditions $x,y\in(\mathbb R/A\mathbb Z)\times (\mathbb R/B\mathbb Z)$. The equations in real space do not change, but now the wave vectors are discretised: $k,l\in 2\pi/A\mathbb Z\times 2\pi/B\mathbb Z$ and the new Green's function is (essentially the Riemann sum of the integral):
$$
\begin{align}
G^{(A,B)} &= \sum_{k,l\in 2\pi/A\mathbb Z\times 2\pi/B\mathbb Z} e^{i(k^2+l^2)z/2\omega+i(kx+ly)} \frac1{AB}\\
&= \left(\sum_{k\in 2\pi/A\mathbb Z} e^{ik^2z/2\omega+ikx}\frac1A\right)\left(\sum_{l\in 2\pi/B\mathbb Z}e^{il^2z/2\omega+ily}\frac1B\right) \\
&= \frac{\vartheta\left(x,\frac{2\pi z}{\omega A^2}\right)}A\frac{\vartheta\left(y,\frac{2\pi z}{\omega B^2}\right)}B
\end{align}
$$
Finally, if you combine both to $A=aM,B = bN$, you get your desired Green's function by simply periodising it:
$$
G^{(M,N)}_{mn} = \sum_{\mu,\nu\in\mathbb Z^2}G_{m+\mu M,n+\nu N}
$$
or more relevant for numerical applications, by sampling and rescaling the integral:
$$
G^{(M,N)}_{mn} = \left(\sum_{\mu=0}^{M-1} \exp\left(\frac{iz}{\omega}\frac{1-\cos(k_\mu a)}{a^2}+ik_\mu x_m\right)\frac1A\right)\left(\sum_{\nu=0}^{N-1} \exp\left(\frac{iz}{\omega}\frac{1-\cos(l_\nu b)}{b^2}+il_\nu y_n\right)\frac1B\right)
$$
where $k_\mu = 2\pi \mu/A,l_\nu = 2\pi \nu/B$. As always, you need to rescale the hemitian product in the discrete setting by $ab$ and the operator is unitary with respect to it. As you anticipated, it is the DFT of the unitary group generated by the discrete Laplacian with periodic boundary conditions. I do not think these finite sums have nice closed form expressions in terms of standard functions.
At the end of the day, the choice of spatial discretisation and finiteness with the discrete Laplacian is rather arbitrary (but this choice should be the simplest). Depending on these, you can get a corresponding discrete Fresnel propagator, but all should agree in the infinite, continuous limit.
For the 1D version, you just get more simply:
$$
\begin{align}
G^{(M)}_m &= \frac1A\sum_{\mu=0}^{M-1} \exp\left(\frac{iz}{\omega}\frac{1-\cos(2\pi \mu/M)}{a^2}+i\frac{2\pi\mu m}{M}\right) \\
&= \frac{\exp\left(\frac{iz}{\omega a^2}\right)}a\sum_{\mu\in\mathbb Z}(-i)^{m+\mu M}J_{m+\mu M}\left(\frac z{\omega a^2}\right)
\end{align}
$$
so you matrix propagator is explicitly:
$$
T^{(M)}_{mm'} = \frac1{A}\sum_{\mu=0}^{M-1} \exp\left(\frac{iz}{\omega}\frac{1-\cos(2\pi \mu/M)}{a^2}+i\frac{2\pi\mu(m-m')}{M}\right)
$$
which you can check directly is unitary:
$$
\sum_{m'=0}^{M-1} T_{mm'}T_{m'm''} a= \frac{\delta_{mm''}}a
$$
