 Solution

Use $X=ϵx$ and the WKB approximationAt 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,\dots$, andFor 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(ϵs)ds]$. Differentiating $e^{-ik(X)x}$ would instead give the erroneous local wavenumber $k+X{k}_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.