Solution to wave equations with predefined phase I want to construct special solutions to wave equations, specifically to Helmholtz-equation and paraxial wave equation (PWE).

Let us consider the PWE
  $$(\partial_x^2 + \partial_y^2 - 2ik\partial_z)\psi(x,y,z)=0$$
  with the constraint
  $$\psi(x,y,z=0)=A(x,y)\cdot \exp(i\cdot g(x,y))$$
  with a specific predefined $g: (x,y) \to \mathbb{R}$, and arbitrary $A: (x,y) \to \mathbb{R}$, with natural boundary conditions $\lim_{x||y \to \infty} \psi(x,y,z)=0$. How can I find $\psi(x,y,z)$ for general $g(x,y)$?

Example:
Suppose we are in cylindrical coordinates, and I set g(r,$\phi$)=$m\cdot \phi$.
In that case I will get $A(r,\phi)=\exp(-r^2)\cdot r^m$, which are Laguerre-Gaussian solutions.
 A: I'm not sure that your original question is well posed. If only $g$ is specified, then the boundary conditions are not fully-specified and the problem of solving for $\psi$ is ill-posed. I will continue, assuming that both $A$ and $g$ are specified and that you want to solve the problem for $\psi$ given the fully-specified boundary conditions.
We have the PWE and boundary conditions
$$
\left|
\begin{array}{l}
(\partial_x^2 + \partial_y^2 - 2ik\partial_z)\psi(x,y,z)=0, \\
\psi(x,y,0) = A(x,y) \exp(ig(x,y)), \\
\psi(x,y,z)\text{ goes to zero fast enough for large }x\text{ and }y,
\end{array}
\right.
$$
where $A$ and $g$ are given functions, regular and smooth enough for our purposes. You'll notice I handwave the degree of regularity and the rate of decay as $x\rightarrow\infty$ and $y\rightarrow\infty$ because, well, I'm not mathematician.
Suppose $\psi$ has the following Fourier integral expansion
$$
\psi(x,y,z) = \int_{-\infty}^\infty\int_{-\infty}^\infty \tilde\psi(k_x,k_y,z) e^{ik_x x+ik_y y}\;dk_x\;dk_y.
$$
Inserting into the original equation, we have
$$
\int_{-\infty}^\infty\int_{-\infty}^\infty(-k_x^2 - k_y^2 - 2ik\partial_z)\tilde\psi(k_x,k_y,z) e^{ik_x x+ik_y y}\;dk_x\;dk_y=0.
$$
For that to integrate to zero, it is sufficient (and I'm certain it can be proved necessary) for the integrand to identically be zero for all values of $k_x$ and $k_y$,
$$
(-k_x^2 - k_y^2 - 2ik\partial_z)\tilde\psi = 0,
$$
where I've dropped the explicit $k_x$ and $k_y$ notation and cancelled the non-zero term $e^{ik_x x+ik_y y}$. Now noting the $z$ derivative with a prime ($'$), we have
$$
2ik\tilde\psi' = -(k_x^2 + k_y^2)\tilde\psi.
$$
Separating $\tilde\psi$ terms all on one side gives
$$
\frac{\tilde\psi'}{\tilde\psi} = -\frac{k_x^2 + k_y^2}{2ik}.
$$
Now integrating $dz$,
$$
\int \frac{1}{\tilde\psi}\;d\tilde\psi = -\int\frac{k_x^2 + k_y^2}{2ik}\;dz,
$$
$$
\log\tilde\psi = -\frac{k_x^2 + k_y^2}{2ik}z + c,
$$
for some constant $c$ (a function of $k_x$ and $k_y$). Now exponentiating gives
$$
\tilde\psi = C(k_x,k_y)\exp\left(-\frac{k_x^2 + k_y^2}{2ik}z\right),
$$
for some other constant (in $z$, but still a function of $k_x$ and $k_y$) $C=\exp(c)$. Now inserting back into the original integral representation, we have
$$
\psi(x,y,z) = \int_{-\infty}^\infty\int_{-\infty}^\infty C(k_x,k_y)\exp\left(ik_x x+ik_y y-\frac{k_x^2 + k_y^2}{2ik}z\right) \;dk_x\;dk_y.
$$
This is our general solution, and you can verify that is solves the original PWE. To meet the boundary condition you require, we must find the specific $C(k_x,k_y)$ that satisfies the condition. So we plug in the value $z=0$ into the general solution, giving
$$
\psi(x,y,0) = \int_{-\infty}^\infty\int_{-\infty}^\infty C(k_x,k_y)\exp\left(ik_x x+ik_y y\right) \;dk_x\;dk_y = A(x,y)\exp(i g(x,y))
$$
It's pretty clear then that $C(k_x,k_y)$ should be the Fourier transform of $A(x,y)\exp(i g(x,y))$. In the Fourier convention I am using (a pretty non-standard one but I like it), the transform from $(x,y)$ to $(k_x,k_y)$ works out to 
$$
C(k_x,k_y) = \frac{1}{(2\pi)^2}\int_{-\infty}^\infty\int_{-\infty}^\infty A(x,y)\exp(i g(x,y)) \exp(-ik_xx-ik_yy)\;dx\;dy
$$
So, for some particular amplitude/phase screen $A\exp(ig)$, you could in principle calculate the Fourier transform above, and insert it into the general solution. You would then have an integral expression for $\psi$ that solves the PWE and the boundary condition. For some particular choices of $A$ and $g$, closed form expressions would result. However, for arbitrary boundary conditions, I expect that no closed form integrals could be calculated.
