How to do a Inverse Discrete Fourier which contains frequency component?

I need do a integral is shown below

It coulde be writen as bellow which kz is related to kx and ky :

I used FourierTransform. But because I need to calculate multiple iterations of ∆z，time gets super long and it did’t work right. So I want to use Discrete Fourier to reduce time. But I don’t know how to deal kz. My idea is to create a matrix of frequency variables to represent it, but I don’t know its frequency range and how to use matrix to represent the frequency in InverseFourierTransform.

That’s what I have tried. The problem happened on :

w = 1;
initialEi[x_] = E^(-x^2/w^2);
InverseFourierTransform[FourierTransform[initial, x, kx]*
          E^(I*Sqrt[k^2 - kx^2]*Δz), kx, x]

For a function like this, InverseFourierTransform didn’t work.

A simplified version of the entire code is shown below:

Clear["Global`*"]
λ = 1.2398/10; 
k = 2*(Pi/λ); 
w = 1;
initialEi[x_] = E^(-x^2/w^2);

Beam[initial_, Δz_] := 
 InverseFourierTransform [FourierTransform[initial, x, kx]*
           E^(I*Sqrt[k^2 - kx^2]*Δz), kx, x]
BeamPropagate[initial2_, Distance_, Δz_] := 
  NestList[Beam, initial2, Ceiling[Distance/Δz]];
listofBeam = Block[
  {x = 1, Δz = .05, Distance = 1}, 
  BeamPropagate[initialEi[x], Distance, Δz]
  ]

This is actually a beam propagation problem, the result I want to get is a graph as below, given the initial state of E[x,z] at z, we can deduce the state of the light wave at nΔz. This method of deriving beam propagation is called Multi- slice Method. I’m not a theoretical professional, so I have a lot of confusion about the implementation of these formulas. If anyone knows this stuff, I’d really appreciate it if you could give me some advices.

By the way this is the Peltio’s answer I referenced about this problem.

xStart = -50;
xEnd = 50;
Δx = .1;

Nx = (xEnd - xStart)/Δx + 1;
kxStart = -π/Δx;
kxEnd = π/Δx;
Δkx = (2 π)/(xEnd - xStart);

x = Table[i, {i, xStart, xEnd, Δx}];
kx = Table[j, {j, kxStart, kxEnd, Δkx}];