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

Muchstronger: the density function of the exit time from $D$ (of the 2-D Brownian) motion decays exponentially fast, so the harmonic measure in the infinite case also decays exponentially fast. By periodization, the same is true for the reflecting boundary. That is, the "weight" associated with $\partial D \times \{x\}$ is roughly $\exp(-\sqrt{\lambda_1} |x|) / \sqrt{4 \lambda_1}$, where $\lambda_1$ is the smallest eigenvalue of $\Delta$ in $D$ with Dirichlet boundary conditions. $\endgroup$