The Rayleigh Expansions and Diffraction in Periodic Media Motivation/Background: The Rayleigh Expansions are well known in the electromagnetics community. They were derived by Lord Rayleigh and give a way of representing solutions to the Helmholtz equation with $\alpha$-quasiperiodic boundary conditions. As discussed in the first chapter of Electromagnetic Theory of Gratings by Roger Petit, they also provide a way of representing outgoing solutions as "propagating" (which are $p\in\mathcal{U}$ in the notation below) and incoming solutions as "evanescent" (i.e. exponential decay and are $p\not\in\mathcal{U}$ in the notation below).
Problem: Suppose we are given the following boundary value problem
\begin{align}
\Delta u + k^2 u &=0,\quad z>g(x),\tag{1a}\\
u(x+d,z)&=e^{i\alpha d}u(x,z),\tag{1b}\\
\end{align}
where $(1a)$ is the Helmholtz equation and $(1b)$ is known as $\alpha$-quasiperiodicity. The quantity $g(x)$ is known as the surface where we consider solutions in the upper layer and use the notation $z$ (the height) > $g(x)$ (the surface). Lord Rayleigh deduced that solutions of this boundary value problem take the form
\begin{align}
u(x,z)=\sum_{p=-\infty}^{\infty}a_pe^{i\alpha_px + i\gamma_p z}+ \sum_{p=-\infty}^{\infty}b_pe^{i\alpha_px - i\gamma_p z}, \tag{2}   
\end{align}
where
\begin{align}
\alpha_p := \alpha + \left(\frac{2\pi}{d}\right)p, \quad \gamma_p:= \begin{cases} 
      \sqrt{k^2-\alpha_p^2}, & p\in\mathcal{U},\tag{4} \\
      i\sqrt{\alpha_p^2-k^2}, & p\not\in\mathcal{U},
   \end{cases}
\end{align}
and
\begin{align}
\mathcal{U}:=\{p\in\mathbb Z ~|~ \alpha_p^2 < k^2\}\tag{5}.
\end{align}
The condition $(5)$ is also known as "propagating modes."
My Attempts:
We split $(1\text{a})$ into a set of ordinary differential equations by considering
$$u(x,z)=X(x)Z(z).$$
Substituting this product into the Helmholtz equation, we obtain
$$X''Z + XZ'' + k^2XZ=0.$$
Dividing by $u = XZ$ and rearranging terms, we get
\begin{align}
\frac{X''}{X}= -k^2 - \frac{Z''}{Z}.\tag{6}
\end{align}
Equation $(6)$ exhibits one Separation of Variables. The left-hand side is a function of $x$ alone, whereas the right-hand side depends only on $z$ and not on $x$. But $x$ and $z$ are independent coordinates. The equality of both sides depending on different variables means that the behavior of $x$ as an independent variable is not determined by $z$. Therefore, each side must be equal to a constant which is known as the constant of separation. We choose
\begin{align}
  \frac{X''}{X}&=\lambda^2,\tag{7} \\ 
  -k^2 - \frac{Z''}{Z}&=\lambda^2.\tag{8}
\end{align}
Now, turning our attention to $(8)$, we obtain
\begin{align}
 \frac{Z''}{Z} = -\left(\lambda^2 +k^2\right).\tag{9}
\end{align}
Our goal is to solve the ODEs $(7)$ and $(9)$ through our boundary conditions. Rearranging $(7)$ as $X''-\lambda^2 X=0$ and solving the auxiliary equation gives
\begin{align}
\quad X_n(x)=a_ne^{\lambda x} + b_ne^{-\lambda x}.\tag{10}
\end{align}
Similarly, we write $(9)$ as $Z'' + Z\left(\lambda^2+k^2\right)$ and solve the auxiliary equation to find
\begin{align}
\quad Z_n(z)=c_ne^{i\sqrt{\lambda^2+k^2}z} + d_ne^{-i\sqrt{\lambda^2+k^2}z}.\tag{11}
\end{align}
We first analyze $(10)$ and focus on the case where $x=0$. The boundary condition $(1\text{b})$ for quasiperiodic $u(x+d,z)=e^{i\alpha d}u(x,z)$ becomes $X(d)=e^{i\alpha d}X(0)$ which, after some simplification, delivers
$$ \lambda = i\alpha.$$
Evaluating $(10)$ when $\lambda=0$ yields
$$\alpha = \left(\frac{2\pi n}{d}\right),\quad n\in\mathbb  Z.$$
Therefore, the auxiliary equation $(10)$ becomes
\begin{align}
\quad X_n(x)=a_ne^{i\alpha x} + b_ne^{-i\alpha x}.\tag{12}
\end{align}
Next, we analyze the auxiliary equations when $z=0.$ The boundary condition $(1\text{b})$ becomes $X(x+d)=e^{i\alpha d}X(x)$ which doesn't simplify $(11)$ or $(12)$. Inserting $\lambda=i\alpha$ into $(11)$ forms
\begin{align}
\quad Z_n(z)=c_ne^{i\sqrt{k^2-\alpha^2}z} + d_ne^{-i\sqrt{k^2 -\alpha^2}z}.\tag{13}
\end{align}
By $(4)$, we see that $(13)$ is equivalent to (this is clearly wrong since $\alpha \neq \alpha_p$. However, I don't see what I did wrong.)
\begin{align}
\quad Z_n(z)=c_ne^{i\gamma_p z} + d_ne^{-i\gamma_p z}.\tag{14}
\end{align}
By $(12),(14),$ the superposition principle, and replacing the integer $n$ by $p$, the general solution is
\begin{align*}
u(x,z)&=\sum_{p=-\infty}^{\infty}\Big( a_pe^{i\alpha x} + b_pe^{-i\alpha x}\Big)\Big(c_pe^{i\gamma_p z} + d_pe^{-i\gamma_p z}\Big)\\&=
\sum_{p=-\infty}^{\infty}e_pe^{i\alpha x + i\gamma_p z} + \sum_{p=-\infty}^{\infty}f_pe^{i\alpha x - i\gamma_p z}.
\end{align*}
My Concerns: In step $(13)$ I should have obtained
\begin{align}
\quad Z_n(z)=c_ne^{i\sqrt{k^2-\alpha_p^2}z} + d_ne^{-i\sqrt{k^2 -\alpha_p^2}z},\tag{15}
\end{align}
where (as defined by $(4)$)
$$\alpha_p := \alpha + \left(\frac{2\pi}{d}\right)p.$$
Then the general solution would be
\begin{align*}
u(x,z)&=\sum_{p=-\infty}^{\infty}\Big( a_pe^{i\alpha_p x} + b_pe^{-i\alpha_p x}\Big)\Big(c_pe^{i\gamma_p z} + d_pe^{-i\gamma_p z}\Big)\\&=
\sum_{p=-\infty}^{\infty}e_pe^{i\alpha_p x + i\gamma_p z} + \sum_{p=-\infty}^{\infty}f_pe^{i\alpha_p x - i\gamma_p z},
\end{align*}
and would match what Lord Rayleigh found. I don't see how to obtain this result by Separation of Variables with the given boundary condition.
 A: The mistake made in the question above happens after equation $(11)$. I forgot that logarithms are multivalued in nature. As I am used to working in the field of real numbers, I incorrectly deduced that when $d\neq 0$,  $e^{\lambda d}=e^{i\alpha d}$ implies $\lambda = i\alpha.$ To correct my analysis, I first recall the definition of the complex logarithm.
Let $z\in\mathbb C$ be a nonzero complex number. We would like to solve for $w$ in the equation
$$e^w=z\tag{1}.$$
If $\theta=\text{Arg}(z)$ with $-\pi < \theta \le \pi$, then $z$ and $w$ can be written as follows
$$z=re^{i\theta}\quad \text{and} \quad  w=u+iv.$$
Then equation $(1)$ becomes
$$e^ue^{iv}=re^{i\theta}.$$
Thus, we have
$$e^u = r \quad \text{and} \quad  v=\theta+2n\pi, \quad n\in\mathbb Z.$$
Since $e^u=r$ is equivalent to the expression $u=\ln r$, it follows that equation $(1)$ is satisfied if and only if $w$ has one of the values
$$w=\ln r+i(\theta+2n\pi),\quad n\in\mathbb Z.$$
Therefore, the (multivalued) logarithmic function of a nonzero complex variable $z=re^{i\theta}$ is defined by the formula
$$\log z=\ln r+i(\theta+2n\pi),\quad n\in\mathbb Z.\tag{2}$$
Hence to solve for $\lambda$ in $e^{\lambda d}=e^{i\alpha d}$ when $d\neq 0$, we may observe for $z=e^{\lambda d},r=1,$ and $\theta=\alpha d$, we may take the logarithm of both sides to deduce
$$\lambda d = i(\alpha d + 2n\pi),\quad n\in\mathbb Z.$$
Since $d\neq 0$, the above expression is equivalent to
$$\lambda = i\alpha + i\left(\frac{2\pi}{d}\right)n,\quad n\in\mathbb Z.\tag{3}$$
To be complete, I also fix a minor mistake made in the analysis in the question.

Problem: Suppose we are given the following boundary value problem
\begin{align}
\Delta u + k^2 u &=0,\quad z>g(x),\tag{4a}\\
u(x+d,z)&=e^{i\alpha d}u(x,z),\tag{4b}\\
\end{align}
Show that solutions take the form
\begin{align}
u(x,z)=\sum_{p=-\infty}^{\infty}a_pe^{i\alpha_px + i\gamma_p z}+ \sum_{p=-\infty}^{\infty}b_pe^{i\alpha_px - i\gamma_p z}, \tag{5}   
\end{align}
where
\begin{align}
\alpha_p := \alpha + \left(\frac{2\pi}{d}\right)p, \quad \gamma_p:= \begin{cases} 
      \sqrt{k^2-\alpha_p^2}, & p\in\mathcal{U},\tag{6} \\
      i\sqrt{\alpha_p^2-k^2}, & p\not\in\mathcal{U},
   \end{cases}
\end{align}
and
\begin{align}
\mathcal{U}:=\{p\in\mathbb Z ~|~ \alpha_p^2 < k^2\}\tag{7}.
\end{align}

Solution: We split $(4\text{a})$ into a set of ordinary differential equations by considering
$$u(x,z)=X(x)Z(z).$$
Substituting this product into the Helmholtz equation, we obtain
$$X''Z + XZ'' + k^2XZ=0.$$
Dividing by $u = XZ$ and rearranging terms, we get
\begin{align}
\frac{X''}{X}= -k^2 - \frac{Z''}{Z}.\tag{8}
\end{align}
Both sides of $(8)$ must be equal to a constant where we choose
\begin{align}
  \frac{X''}{X}&=\lambda^2, \tag{9}\\ 
  -k^2 - \frac{Z''}{Z}&=\lambda^2. \tag{10}
\end{align}
Now, turning our attention to $(10)$, we obtain
\begin{align}
 \frac{Z''}{Z} = -\left(\lambda^2 +k^2\right).  \tag{11}
\end{align}
Our goal is to solve the ODEs $(9)$ and $(11)$ through our boundary conditions. Rearranging $(9)$ as $X''-\lambda^2 X=0$ and solving the auxiliary equation gives
\begin{align}
\quad X_n(x)=a_ne^{\lambda x} + b_ne^{-\lambda x}.\tag{12}
\end{align}
Similarly, we write $(11)$ as $Z'' + Z\left(\lambda^2+k^2\right)$ and solve the auxiliary equation to find
\begin{align}
\quad Z_n(z)=c_ne^{i\sqrt{\lambda^2+k^2}z} + d_ne^{-i\sqrt{\lambda^2+k^2}z}.\tag{13}
\end{align}
We first analyze $(12)$ in combination with the boundary condition $(4\text{b})$ for quasiperiodic $u(x+d,z)=e^{i\alpha d}u(x,z)$. The case $d=0$ gives the trivial solution so we will assume $d\neq 0$. Now we analyze $(12)$ and focus on the case where $x=0$ and $d\neq 0$. The boundary condition $(4\text{b})$ becomes $X(d)=e^{i\alpha d}X(0)$ which, upon invoking the definition of the complex logarithm and simplifying as in $(3)$, delivers
$$ \lambda = i\alpha+i\left(\frac{2\pi }{d}\right)p,\quad p\in\mathbb  Z, ~d\neq 0.$$
It is clear from the left-hand side of $(6)$ that we may represent $\alpha_p$ by
$$ \lambda = i\alpha+i\left(\frac{2\pi }{d}\right)p=i\alpha_p,\quad p\in\mathbb  Z, ~d\neq 0.$$
Next, we analyze the auxiliary equations when $z=0.$ The boundary condition $(4\text{b})$ becomes $X(x+d)=e^{i\alpha d}X(x)$ which, after some elementary computations, gives the same representation for $\lambda$. Therefore no further simplifications are necessary for $(12)$ and $(13)$. In light of this, the auxiliary equation $(12)$ becomes
\begin{align}
\quad X_n(x)=a_ne^{i\alpha_p x} + b_ne^{-i\alpha_p x}.\tag{14}
\end{align}
Inserting $\lambda=i\alpha_p$ into $(13)$ forms
\begin{align}
\quad Z_n(z)=c_ne^{i\sqrt{k^2-\alpha_p^2}z} + d_ne^{-i\sqrt{k^2 -\alpha_p^2}z}.\tag{15}
\end{align}
By the right-hand side of $(6)$, we see that $(15)$ is equivalent to
\begin{align}
\quad Z_n(z)=c_ne^{i\gamma_p z} + d_ne^{-i\gamma_p z}.\tag{16}
\end{align}
By $(14),(16),$ and the superposition principle, the general solution is
\begin{align*}
u(x,z)&=\sum_{p=-\infty}^{\infty}\Big( a_pe^{i\alpha_p x} + b_pe^{-i\alpha_p x}\Big)\Big(c_pe^{i\gamma_p z} + d_pe^{-i\gamma_p z}\Big)\\&=
\sum_{p=-\infty}^{\infty}e_pe^{i\alpha_p x + i\gamma_p z} + \sum_{p=-\infty}^{\infty}f_pe^{i\alpha_p x - i\gamma_p z}.\tag{17}
\end{align*}
A similar calculation was made by Lord Rayleigh whose initial focus was investigating the electromagnetic problem of the diffraction of a plane wave at normal incidence in a periodic structure. If you are familiar with plane waves then you may recognize that equation $(17)$ is often called the  Rayleigh expansion. In fact,  every term  of  this  expansion  represents  a  propagating  plane  wave  only if $p\in\mathcal{U}$ where $\mathcal{U}$ is defined in $(7)$. If $p\not\in \mathcal{U}$ then the associated term represents an evanescent wave that is exponentially damped. To enforce the requirement that solutions to $(4)$ are outgoing and bounded, we require $f_p\equiv 0.$ If $f_p\not\equiv 0,$ then solutions will be inward propagating for $p\in \mathcal{U}$, and unbounded for $p\not \in \mathcal{U}$. In classical optics, practitioners are mainly interested in the propagating waves. This is the reason for the definitions $(6)$ and $(7)$: for $p\in\mathcal{U}$, solutions are outgoing and the modes are "propagating." In contrast, if $p\not\in\mathcal{U}$, then solutions decay exponentially and the modes are known as "evanescent."
