= Solution
Use $X=\epsilon x$ and the <WKB approximation>
$$
p=e^{-i\Theta(X)/\epsilon}\bigl[P_0(y,X)+\epsilon P_1(y,X)+\cdots\bigr],\qquad k(X)=\Theta_X(X).
$$
At leading order, the divergence-form acoustic <Helmholtz equation> becomes $P_{0,yy}+[k_0(X)^2-k(X)^2]P_0=0$. The sloping-wall <Neumann boundary condition> becomes $P_{0,y}=0$ at $y=\pm R(X)$. In the even transverse family this has <normal modes> $P_0=A(X)\cos(n\pi y/R(X))$, where $n=0,1,\ldots$, and
$$
\boxed{k(X)=\sqrt{k_0(X)^2-\left(\frac{n\pi}{R(X)}\right)^2},\quad
p_0=A(X)\cos\!\left(\frac{n\pi y}{R(X)}\right)
\exp\!\left[-\frac i\epsilon\int^Xk(\xi)d\xi\right].}
$$
For the right-going propagating branch take real positive $k$. The displayed local expression $e^{-ikx}$ in the question must be read as a frozen-coefficient shorthand: the actual slowly varying phase is $\exp[-i\int^x k(\epsilon s)ds]$. Differentiating $e^{-ik(X)x}$ would instead give the erroneous local wavenumber $k+Xk_X$.
The symmetric duct also has odd transverse <normal modes> $\sin[(n+1/2)\pi y/R]$. The integer-cosine form selects the even family requested here, rather than representing every duct mode. The approximation requires smooth slow variation and a positive $k$ separated from cutoff; it fails in a cutoff transition.
Back to article page
