$\newcommand{\R}{\mathbb R}$Sorry for being too sketchy in the following answer, time permitting, I'll try to expand.

***

*Step 0.* Some more-or-less classical potential theory. Let $D$ be an open set in $\R^d$ (with $d \geqslant 3$ for simplicity), and assume that $D$ is sufficiently regular (for example, Lipschitz). Let $p_t^D(x, \xi)$ be the *heat kernel* in $D$: the fundamental solution of the heat equation $$\frac{\partial u}{\partial t}(t, x) = \Delta u(t, x) $$ with $u(t, x) = 0$ when $x \in \partial D$. In other words, $$u(t, x) = \int_D p_t^D(x, \xi) f(\xi) d\xi$$ is the (unique) solution of the heat equation with initial condition $u(0, x) = f(x)$.

The *Green function* in $D$ can be defined by the time integral of the heat kernel:
$$G_D(x, \xi) = \int_0^\infty p_t^D(x, \xi) dt .$$
Note that this is always finite (except when $x = \xi$, of course), because $p_t^D(x, \xi) \leqslant p_t^{\R^d}(x, \xi)$, and the integral of $p_t^{\R^d}(x, \xi)$ is the Newtonian potential kernel $c_d |x - \xi|^{2 - d}$. Furthermore, $G_D(x, \xi)$ is zero on the boundary (this is not quite obvious, though), and — formally — we have $$ \Delta_x G_D(x, \xi) = \int_0^\infty \Delta_x p_t^D(x, \xi) dt = \int_0^\infty \frac{\partial p_t}{\partial t}(x, \xi) dt = 0 - \delta_\xi(x) ,$$ so that $G_D(x, \xi)$ is the fundamental solution for the Poisson problem $$ \Delta u(x) = -f(x) $$ in $D$, with $u(x) = 0$ for $x \in \partial D$. And indeed, one can rigorously prove that if $$u(x) = \int_D G_D(x, \xi) f(\xi) d\xi$$ for, say, continuous and bounded $f$, then indeed $\Delta u = -f$ in $D$ and $u = 0$ on $\partial D$.

Finally, if $D$ is regular enough ($C^{1,1}$ is the usual condition), then the Dirichlet problem in $D$: $$\Delta u(x) = 0$$ with Dirichlet boundary condition $u(x) = f(x)$ for $x \in \partial D$ can be solved using the *Poisson kernel*: $$u(x) = \int_{\partial D} f(\xi) P_D(x, \xi) \sigma(d\xi),$$ where $\sigma$ is the surface measure on $\partial D$ and the Poisson kernel $P_D(x, \xi)$ is the boundary derivative of the Green function: $$P_D(x, \xi) = \frac{\partial G_D(x, \cdot)}{\partial \nu}(\xi) = \lim_{s \to 0^+} \frac{G_D(x, \xi + s \nu)}{s} \, ,$$ where $\nu$ is the inward normal vector at $\xi$. This follows relatively easily from the divergence theorem (or Green's identities).

By the way, for a general open set $D$, the solution of the Dirichlet problem is given in terms of the *harmonic measure*: $$ u(x) = \int_{\partial D} f(\xi) P_D(x, d\xi) ,$$
which is again closely related to the Green function, but this is a completely different story.

***

*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; here $f$ is an arbitrary bounded and continuous function). The solution is given by the harmonic measure, which, due to translation invariance of the problem, is translation invariant itself:
$$ \begin{aligned} u(x, y) & = \int_{\partial D \times \R} f(\xi, \eta) P_{D \times \R}(x, y, d\xi d\eta) \\ & = \int_{\partial D \times \R} f(\xi, y + \eta) P_{D \times \R}(x, 0, d\xi d\eta) \end{aligned} $$
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.

Using an estimate $0 \leqslant p_t^D(x, \xi) \leqslant p_t^{\R^2}(x, \xi)$ (where $p_t^{\R^2}$ is the usual Gauss–Weierstrass kernel) for $t < 1$ and intrinsic ultracontractiviety for $t > 1$, by direct integration, we find the following estimate of the Green function, valid when $|\eta|$ is large enough (here I omit the details):
$$\begin{aligned} G_{D \times \R}(x, 0, \xi, \eta) & = \int_0^\infty p_t^{D \times \R}(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|} . $$
This may look as if we "differentiate both sides of an inequality", but it is not the case: since the Green function is zero on the boundary, the normal derivative reduces to a simple limit of $G_D(x, 0, \xi + s \nu, \eta) / s$, where $\nu$ is the inward normal vector at $\xi$.

***

*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$.