Transparent boundary condition I am interested in the finite-difference beam propagation method and its applications. I try to solve the Helmholtz equation. At first, i would like to solve numerically it for the easiest case, without nonlinearities. Just to make sure I'm on the right way. But i really don't understand how to wright the boundary condition. I chose the transparent boundary condition and i need to write it properly to solve numerically the equation.
So, for a linear, homogeneous and instantaneous medium the Helmholtz equation is writen (in 3D case, z is the propagation direction)
$$
\frac{\partial^{2} E(x,y,z)}{\partial x^{2}} + \frac{\partial^{2} E(x,y,z)}{\partial y^{2}} + \frac{\partial^{2} E(x,y,z)}{\partial z^{2}} = - (k_{0} n )^{2} E(x,y,z)
$$
It can be solved if the initial condition is known, $E(x,y,0)$.
Introducing operator $\hat{S}$
$$
\hat{S} = \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + (k_{0} n )^{2}
$$ 
The equation can be written in the following form
$$
\frac{\partial^{2} E(x,y,z)}{\partial z^{2}} = -\hat{S} \ E(x,y,z)
$$
The solution of this equation is
$$
E(x,y,z) = \exp \left [ - i \sqrt{\hat{S}} z \right ] E^{+}(x,y,0) + \exp \left [  i \sqrt{\hat{S}} z \right ] E^{-}(x,y,0)
$$
Considering only the forward propagating component and introducing the propagation operator $\hat{P}^{+}$ the electric field at $z=\Delta z$ can be written through the value of the field at $z=0$ (initial condition written earlier) and so on.
$$
E(x,y,\Delta z) = \hat{P}^{+}(\Delta z) \ E(x,y,0)
$$
where
$$
\hat{P}^{+}(\Delta z) = \sum \limits_{n=0}^{\infty} \frac{1}{n!}\left[- i \sqrt{\hat{S}} \right]^{n} \Delta z^{n}
$$
Obtained expression can be adopted to the Crank-Nicholson scheme. But it is also necessary to write the boundary condition. How to write the boundary condition if the medium is confined in the transparent walls ? 
 A: Since you seem to be using finite differences, you should look at the paper of Hadley, titled 'Transparent Boundary Condition for the Beam Propagation Method' - without any treatment of the boundary values you're automatically assuming a Dirichlet boundary condition.
You can incorporate the boundary conditions into your square-root differential operator $\sqrt{\hat{S}}$ by approximating it via a Padé approximant, in short: quotients of polynomials in $\hat{S}$. With this you can arrive at equations containing $\hat{S}$ and $E$, and satisfy your boundary conditions for your $E$ field by choosing the right 'boundary entries' for $\hat{S}$'s matrix representation. That's the only way I managed to find and implement, as it is similar to doing the same for the paraxial approximation that reduces the Helmholtz-equation to a form that would only contain $S$ instead of $\sqrt{\hat{S}}$ (I'm not a clever man).
As you'll see from the paper of Hadley, the simplest form of the TBC (transparent boundary condition) involves treating the field at the boundary as a lateral plane wave and making sure that no reflected plane wave gets back into the computation window by modifying the fields at the boundary.
When I started to look into solving the Helmholtz-Equation I came around a dissertation by Filippo Pigozzo that was valuable for getting a very good overview of the matter.
