We solve the coupled equations numerically by converting Eq. (4) into Fourier space along the transverse dimensions $x$ and $y$, giving Display Formula
$2\pi ikp(fxkp,x+fykp,y)A\u02dc+kp,zkp\u2202A\u02dc\u2202z=F{\u2212i\kappa e\u2212i\varphi (x,y)B}2\pi ikp(fxkd,x+fykd,y)B\u02dc+kd,zkp\u2202B\u02dc\u2202z=F{\u2212i\kappa ei\varphi (x,y)A},$(6)
where $A\u02dc$ and $B\u02dc$ are the Fourier transforms of $A$ and $B$, respectively, and $fx$ and $fy$ are the spatial frequencies along the $x$ and $y$ axes, respectively. To solve these equations, we split the propagation and energy transfer between the waves into two discrete steps and successively propagate and transfer energy between waves over several small propagation steps. To calculate the propagation of the beam, the right side of Eq. (6) is assumed to be zero. In this case, the Fourier amplitudes will have a solution of the form Display Formula$A\u02dc(fx,fy,z+\Delta z)=A\u02dc(fx,fy,z)e[\u2212i2\pi kp,z(fxkp,x+fykp,y)\Delta z]B\u02dc(fx,fy,z+\Delta z)=B\u02dc(fx,fy,z)e[\u2212i2\pi kd,z(fxkd,x+fykd,y)\Delta z].$(7)
Note that this is only exact in the case where the right side of Eq. (6) truly equals zero, but for small ($\u223c100\u2009\u2009nm$) propagation steps, this is a reasonable approximation. To account for energy transfer, we piecewise integrate the right side of Eq. (5) with the Euler method and add it to the inverse Fourier transform of Eq. (7) Display Formula$A(x,y,z+\Delta z)=F\u22121{A\u02dc(fx,fy,z+\Delta z)}\u2212i\kappa e\u2212i\varphi (x,y)B(x,y,z)\Delta zB(x,y,z+\Delta z)=F\u22121{B\u02dc(fx,fy,z+\Delta z)}\u2212i\kappa ei\varphi (x,y)A(x,y,z)\Delta z.$(8)