Fourier propagation is propagating backwards

I am trying to write a code which propagates a beam a certain distance through Fourier propagation. Here is the relevant excerpt of the code:

import numpy as np
import scipy
from matplotlib import pyplot as plt
from numpy.fft import fft, ifft, fft2, ifft2, fftfreq, fftshift, ifftshift

exp=np.exp
pi=np.pi
sqrt=np.sqrt
tg=np.tan
arctg=np.arctan2
angle=np.angle
array=np.array
plt.rcParams["figure.autolayout"] = True

class beam:
  def __init__(self, z_0: float, compOnda: float):
    self.z_0 = z_0 #distância de Rayleigh
    self.compOnda = compOnda #comprimento de onda

    self.k = 2*pi / self.compOnda #número de onda
    self.W_0 = sqrt(self.compOnda * self.z_0 / pi)
    self.teta_0 = self.W_0 / self.z_0
    self.A_0 = sqrt(2 / (pi * self.W_0)) #amplitude em 0; definimos ela de forma a normalizar o pulso (gaussiano - preciso verificar se normaliza o laguerre gaussiano também)
    self.I_0 = self.A_0 * self.A_0.conjugate() #intensidade em 0

  def W(self,z: np.array):
    return self.W_0 * sqrt(1 + (z / self.z_0)**2)

  def R(self,z: np.array):
    z_=np.where(z==0.0,1e-20,z) #evitando divisão por zero
    R=z_ * (1 + (self.z_0 / z_)**2)
    R=np.where(z==0.0,np.inf,R)
    return R

  def zeta(self,z: np.array):
    return arctg(z, self.z_0)

class GaussianBeam(beam):
  def __init__(self, z_0: float, compOnda: float):
    super().__init__(z_0,compOnda)

  def Amplitude(self,z: array = 0.0, rho: array = 0.0, phi: array = 0.0) -> array:
    z=np.asarray(z,dtype=float) 

    z_0=self.z_0
    A_0=self.A_0
    k=self.k
    W_0=self.W_0

    W=self.W(z)
    R=self.R(z)
    zeta=self.zeta(z)

    U=A_0*(W_0/W)*exp(-rho**2/W**2 + (-k*z-k*(rho**2)/(2*R)+zeta)*1.0j)
    return U

wavelength=800*(10**-9) #nm
z_0=pi*(10**(-3))/2 #m ->approx. 0.0015. W_0=20.10⁻6
dz=0.006 #propagation distance (m)

def Propag(U,dz,gr,d):
  #frequency coordinates
  fx = fftfreq(gr, d=d)
  fy = fftfreq(gr, d=d)
  fx, fy = np.meshgrid(fx, fy)
  kx = 2*pi*fx
  ky = 2*pi*fy
  kz = sqrt(k**2 - kx**2 - ky**2 + 0j) #propagation constant in z

  #radial complex amplitude
  Fourier_U = fft2(U) #angular spectrum
  Propag_FU = Fourier_U * exp(1j * kz * dz)
  U_propag = ifft2(Propag_FU)
  return U_propag

fx=GaussianBeam(z_0,wavelength)
W_0=fx.W_0
gr=1000
k=fx.k
x = np.linspace(-4*W_0,4*W_0,gr) #grid size
y = np.linspace(-4*W_0,4*W_0,gr)
d=8*W_0/(gr-1)
X,Y = np.meshgrid(x,y)
rho=sqrt(X**2+Y**2)
phi=arctg(Y,X)
U=fx.Amplitude(z_0,rho,phi) #this is a bidimensional Gaussian

U_propag=Propag(U,dz,gr,d)

def plot(u):
  modulo = abs(u)
  fase = angle(u)

  # Gráfico
  plt.figure(figsize=(12, 5))

  plt.subplot(1, 2, 1)
  plt.pcolor(X, Y, modulo, shading='auto', cmap='viridis')
  plt.title('Módulo')
  plt.xlabel('x')
  plt.ylabel('y')
  plt.colorbar()

  plt.subplot(1, 2, 2)
  plt.pcolor(X, Y, fase, shading='auto', cmap='twilight')
  plt.title('Fase')
  plt.xlabel('x')
  plt.ylabel('y')
  plt.colorbar()

  plt.tight_layout()
  plt.show()
  
plot(U_propag)

The code works, but it propagates the beam backwards - which can be seen by the shrunkage of the profile when dz=z_0, meaning that the beam is propagating towards its waist. I want to fix this and make the beam propagate forward. I read that the Fourier transform in optics has a positive sign, but swapping fft and ifft produced the same results. I then tried to give exp(1j * kz * dz) a negative sign, but that gave me an error:

/tmp/ipython-input-17-2768534464.py:14: RuntimeWarning: overflow encountered in exp
  Propag_FU = Fourier_U * exp(-1j * kz * dz)
/usr/local/lib/python3.11/dist-packages/numpy/fft/_pocketfft.py:94: RuntimeWarning: invalid value encountered in fft
  return ufunc(a, fct, axes=[(axis,), (), (axis,)], out=out)

Why is my beam being propagated backwards, and how can I propagate it forward?

If you want to take the complex conjugate of a numpy array, the function .conj() will do that. If I remember my Gaussian optics correctly (a big if, at this point), that should change the direction of the beam.

If I use .conj() instead of changing 1j to -1j, the calculation works without error. The beam appears to be propagating away from its waist, but I didn't spend a whole lot of time on this. I think you have more of a physics problem here than a programming problem.

In other words, replace this line:

Propag_FU = Fourier_U * exp(1j * kz * dz)

with:

Propag_FU = Fourier_U * exp(1j * kz * dz).conj()

Sorting out the various factors of 1j, -1j, and which version of the fft to use here is a bit outside my current pay grade.