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.
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."