A paper of Kneser (1904) strongly suggests that the idea does (indeed) go back to Cauchy, in connection with Sturm-Liouville problems (i.e. ordinary differential operators, as opposed to the Laplacian in the body of your question). Given functions $g,k,l$ and writing $\smash{L=\frac d{dx}\left(k\frac{d}{dx}\,\cdot\right) - l},$ Kneser considers the [for us: “eigenvalue”] problem
$$
LV + rgV=0
$$
with boundary conditions
$$
\left[k\frac{dV}{dx}-hV\right](0)=0,\qquad
\left[k\frac{dV}{dx}+HV\right](X)=0.
\tag1
$$
Addressing the Sturm-Liouville (1837) question whether any $f(x)$ can be expanded into a series of solutions $V_\nu$ belonging to [“eigenvalues”] $r_\nu$:
$$
f(x) = A_1V_1+A_2V_2+\cdots,
\tag3
$$
Kneser writes:
The particular analytic developments I use for this are inspired or drawn from relevant work of Dini (1880), Harnack (1887), Poincaré (1894, 1895) and Stekloff (1901); but the basic idea can be explained as follows.
These recent authors all use a device introduced by Cauchy (1827) in his study of Fourier series: they build a function of a complex variable $r$ containing $x$ as a parameter, having poles at $r=r_\nu$ as its only singularities, and producing as residues the corresponding terms of the series $(3)$. Poincaré, apparently, first pointed out [I guess here: (1894, 1895)] that Cauchy's auxiliary function is the solution of the equation [of the “resolvent” $\smash{(L+rg)^{-1}}$]
$$
LV+rgV+f(x)=0
$$
satisfying conditions $(1)$.
While this early literature is no easy reading, the treatises of Picard (1893, pp. 167-183), Poincaré (1895, pp. 210-223) and Watson (1922, pp. 576-617) have chapter-long expositions of what may be the first three cases historically:
$(g,k,l)=(1,1,0)$ on $[0,\pi]$ with Dirichlet boundary conditions. Then $V_\nu=\sin(\nu x)$, $\smash{r_\nu=\nu^2}$, and (3) is a Fourier sine series. Or Neumann conditions, $V_\nu=\cos(\nu x)$, and cosine series; which Picard (p. 177) and Poincaré (p. 220) attribute to Cauchy (1827, pp. 364-365).
$(g,k,l)=(1,1,0)$ on $[0,1]$ with Fourier's sphere cooling conditions (1822, pp. 340-342): $V(0)=0$ and $V(1)=AV'(1)$ for some $A>1$. Then $V_\nu=\sin(k_\nu x)$ and $\smash{r_\nu=k_\nu^2}$ where the $k_\nu$ are the positive solutions of $\tan(k)=Ak$. So (3) is a “nonharmonic Fourier series,” which Picard (pp. 178-183) and Poincaré (pp. 168-179, 220-223) attribute to Cauchy.
$(g,k,l)=(x,x,a^2/x)$ on $[0,1]$ with $V(1)=0$, and no condition at the singular endpoint $0$. Then $V_\nu=J_a(k_\nu x)$ and $\smash{r_\nu=k_\nu^2}$ where the $k_\nu$ are the roots of the Bessel function $J_a$. So (3) is a “Fourier-Bessel series,” which Watson (pp. 582-591) attributes to Schläfli (1876).
Note added: An earlier (and perhaps clearest) statement by Cauchy occurs in his Application du calcul des résidus à l’intégration des équations différentielles linéaires et à coefficients constants (Exercices de mathématiques 1 (1826) 202-204 = Œuvres (2) 6 (1887) 252-255):
Consider first the task of integrating the differential equation
$$
\frac{d^ny}{dx^n} + a_1
\frac{d^{n-1}y}{dx^{n-1}} + a_2
\frac{d^{n-2}y}{dx^{n-2}} + \ldots + a_{n-1}
\frac{dy}{dx} + a_ny=0,
\tag1
$$
where $a_1, a_2,\dots a_{n-1}, a_n$ denote constant coefficients; and let, for short
$$
F(r) = r^n + a_1r^{n-1} + a_2r^{n-2}+\dots+ a_{n-1}r + a_n.
\tag2
$$
It is clear that, to satisfy the equation (1), it will suffice to take
$$
y = \raise{-1ex}{\huge{\mathcal E}}\,\frac{\varphi(r)\,e^{rx}}{((\,F(r)\,))},
\tag3
$$
where $\varphi(r)$ denotes any function of $r$ which does not become infinite for values of $r$ that verify the formula
$$
F(r)=0.
\tag4
$$
(Of course, this isn’t yet framed in terms of eigenvalues of the differential operator. It becomes so if we replace $a_n$ by $a_n -\lambda$, but Cauchy didn’t call (4) “characteristic equation” until 1839, and names for its roots seem to have come even later — I’m not sure when.)
Also, for symmetric operators (or quadratic forms) on $\mathbf R^n$ this is all in Weierstrass (1859), cf. p. 219.