$\newcommand{\R}{\mathbb R}$Sorry for being too sketchy in the following answer, time permitting, I'll try to expand.
Step 1. First, consider the Poisson problem in $D \times \R$, with boundary data given by $f : \partial D \times \R \to \R$ (let us denote the boundary data by $f$ rather than $\varphi$, which we will need for the eigenfunctions). The solution is given by what is called the harmonic measure, which, due to translation invariance of the problem, is translation invariant itself: $$ u(x, y) = \int_{\partial D \times \R} f(\xi, \eta) P_{D \times \R}(x, y, d\xi d\eta) $$ for an appropriate measure $P_{D \times \R}(x, y, d\xi d\eta)$. If $D$ is nice enough — say $C^{1,1}$ — then $P_{D \times \R}(x, y, d\xi d\eta)$ has a density function $P_{D \times \R}(x, y, \xi, \eta)$ with respect to the surface measure $\sigma(d\xi) d\eta$, and this density is called the Poisson kernel. Thus, $$ u(x, y) = \int_{\R} \int_{\partial D} f(\xi, \eta) P_{D \times \R}(x, y, \xi, \eta) \sigma(d\xi) d\eta . $$ Translation invariance means that $P_{D \times \R}(x, y, \xi, \eta) = P_{D \times \R}(x, 0, \xi, \eta - y)$.
Step 2. How fast does $P_{D \times \R}(x, 0, \xi, \eta)$ decay with $|\eta|$? In my comment to the question, I sketched a probabilistic argument which shows exponential decay. Here is a more analytic (but still potential-theoretic) version of the same argument.
Let $p_t^D(x, \xi)$ be the heat kernel in $D$, and $p_t^{D \times \R}(x, y, \xi, \eta)$ be the heat kernel in $D \times \R$. Thus, $$ p_t^{D \times \R}(x, y, \xi, \eta) = p_t^D(x, \xi) (4 \pi t)^{-1/2} e^{-(\eta - y)^2 / (4t)} . $$ Again if $D$ is nice enough ($C^{1,1}$ is more than enough, Lipschitz is already fine), than $p_t^D$ is known to be intrinsically ultracontractive. In particular, $$ p_t^D(x, \xi) \approx C e^{-\lambda_1 t} \varphi_1(x) \varphi_1(\xi) $$ for $t > 1$. Here $\approx$ means that the ratio is bounded from above and below by positive constants.
Some standard estimates (for $t < 1$) and direct integration give a bound for the Green function when $|\eta|$ is large enough: $$\begin{aligned} G_{D \times \R}(x, 0, \xi, \eta) & = \int_0^\infty p_t(x, 0, \xi, \eta) dt \\ & \approx \varphi_1(x) \varphi_1(\xi) \int_0^\infty e^{-\lambda_1 t} t^{-1/2} e^{-\eta^2 / (4 t)} dt \\ & \approx \varphi_1(x) \varphi_1(\xi) e^{-\sqrt{\lambda_1} |\eta|} . \end{aligned} $$ Now the Poisson kernel is the normal derivative of the Green function. Thus, if $D$ is a $C^{1,1}$ set, $$ P_{D \times \R}(x, 0, \xi, \eta) \approx \varphi_1(x) e^{-\sqrt{\lambda_1} |\eta|} . $$
Step 3. Now it remains to translate this into a result on $D \times (-L, L)$ with zero Neumann boundary condition on the bases. This, however, is pretty standard: if $u$ is the solution of the problem on $D \times (-L, L)$, then the function $v$ given by $$ v(x, y + 4 n L) = u(x, y) , \qquad v(x, y + 2 n L) = u(x, -y) $$ whenever $x \in D$< $y \in (-L, L)$ and $n \in \mathbb Z$, is a solution of the corresponding Poisson problem in $D \times \R$. Using the Poisson representation for $v$, we find that $$ u(x, y) = \int_{(-L, L)} \int_{\partial D} f(\xi, \eta) \sum_{n = -\infty}^\infty (P_{D \times \R}(x, y, \xi, \eta + 4 n L) + P_{D \times \R}(x, y, \xi, -\eta + 2 n L)) \sigma(d\xi) d\eta . $$ By using the estimate for the Poisson kernel found above, we easily see that again $$ u(x, 0) \approx \varphi_1(x) \int_{(-L, L)} \int_{\partial D} f(\xi, \eta) e^{-\sqrt{\lambda_1} \eta} \sigma(d\xi) d\eta $$ uniformly in $L$ large enough and $f$.