The operator calculation of Helmholtz equation? I am reading the beam propagation method (BPM) in optical imaging paper. I find a paper solve the Helmholtz equation in the inhomogeneous media. The paper is:
Light propagation in graded-index optical fibers, M. D. Feit and J. A. Fleck, Jr
They use the operator form to solve the equation, I never see this kind of operator calculation to solve the classical eletromagnetic proble before.  I can intuitively understand how it is going on. But I want to learn more about this and acquire such method. However, i am not a mathmatical stduent. I don't want to learn a full book of the function analysis or some group theory to just acquire this. I mean is there some good,simple and not very long pages reading material to learn this kind of technique.
Besides, there are some places that i don't understand in the paper.

*

*Question 1:

Like the equation (3)

How to obtain this? Now the the expansion for $(1+x)^{1/2}$ is not hold for the Laplacian operator?


*The second place i don't understand is:


Why it is equivalent to the Helmholtz equation? Intuitively, I know that using the angular spectrum to solve the Helmholtz equation, we can get maybe the top solution. However, is there or mathmatical derivation of this? From top to expression to the equation (10).
 A: Consider $A := \sqrt{\nabla_\perp^2 + k^2I}$ and $B=kI$ where $k$ is a real positive number. Since $[A,B]=0$, it holds $$(A-B)(A+B) = A^2 +AB-BA -B^2 = A^2 +BA-BA -B^2 = A^2-B^2\:.$$ Hence
$$A = \frac{A^2-B^2}{A+B} + B\:.$$
That is your first identity.
Regarding the second issue, the equation you want to solve is
$$\left(\nabla^2_\perp + k^2 + \frac{\partial^2}{\partial z^2}\right)E=0 \tag{1}\:.$$
You can factorize it as
$$\left(\sqrt{\nabla^2_\perp + k^2} + i \frac{\partial }{\partial z}\right) \left(\sqrt{\nabla^2_\perp + k^2} - i \frac{\partial }{\partial z} \right) E=0\:.$$
In summary, if you are able do solve both equations
$$\left(\sqrt{\nabla^2_\perp + k^2} \pm i \frac{\partial }{\partial z}\right) E_\pm=0$$
you are done by using twice the form of the solution.
As the paper declares on its first page, the solutions of the equations above are
$$E_\pm(x,y,z) = e^{\mp i z\sqrt{\nabla^2_\perp + k^2}} E_\pm(x,y,0)\:.$$
Using the initial identity,
$$E_\pm(x,y,z) = e^{\mp i z\left[\frac{\nabla_\perp^2}{\sqrt{\nabla^2_\perp + k^2}+ k}+ k\right]} E_\pm(x,y,0)
=e^{\mp i k z} e^{\mp i z\frac{\nabla_\perp^2}{\sqrt{\nabla^2_\perp + k^2}+k}} E_\pm(x,y,0)\:.$$
So, the information stored in
$$e^{\mp i z\frac{\nabla_\perp^2}{\sqrt{\nabla^2_\perp + k^2}+k}} E_\pm(x,y,0)\:,$$
taking the trivial phase $e^{\mp i k z}$ into account, permits to solve (1).
In this sense the said operation is ``equivalent'' to solve (1). I omitted a huge number of mathematical subtleties which do not affect the physical substance here.
COMMENT. Reading quite superficially the paper I realized that what I indicated by $k$ is actually a function of $x,y$ and not a constant. If it is really true my derivation does not hold, but also the your first identity  declared in the paper is false. It is of crucial relevance that $\nabla_\perp^2$ commutes with $k$.
A: For the first part of the question
we pose : $x=\Delta_{\perp}\;\;, y=\frac{\omega}{c}n$
the second member of the equation becomes:
$\frac{x^{2}}{\sqrt{x^{2}+y^{2}}+y}+y=\frac{x^{2}(\sqrt{x^{2}+y^{2}}-y)}{(\sqrt{x^{2}+y^{2}}+y)(\sqrt{x^{2}+y^{2}}-y)}+y=\frac{x^{2}(\sqrt{x^{2}+y^{2}}-y)}{x^{2}}+y=(x^{2}+y^{2})^{\frac{1}{2}}$
