In this post we will give an example of the Riemann zeta function at negative arguments that naturally comes up in physics. As it has been mentioned in the comment section, the expression for the Casimir force in $d$ spatial dimensions is related to $\zeta(-d)$; more precisely, the force per unit area is given by
\begin{equation}
\frac{F_\mathrm{Casimir}}{L^{d-1}}=-\frac{C(d)}{a^{d+1}} \tag{1}
\end{equation}
where $L^{d-1}$ is the hyper-area of the plates, and $a$ is the separation between them. Here, the constant $C(d)$ is given by
$$
\begin{aligned}
C(d)&=-2^{-d}\pi^{d/2} \Gamma\left(1-\frac{d}{2}\right)\color{red}{\zeta(-d)}\\
&=\frac{2^{-d}}{1+d}\ \xi(-d)
\end{aligned}\tag{2}
$$
(Remark: in the case of fields with spin, one has to multiply the expression above by the number of degrees of freedom; $2s+1$ for massive particles with spin $s$, and $2$ for non-scalar massless particle).
We will next outline the calculation of $F(d)$ in the zeta-regularisation scheme, though, as usual, one may prove that the result is scheme-independent (the proof is left to the reader).
We consider two plates of length $L$ (and hence, hyper-area $L^{d-1}$), separated by a distance $a$, such that the enclosed volume is $L^{d-1}a$. The free-field energy inside the plates, as given by the quantisation of a free massless scalar field, is
$$
E_0=\frac12 L^{d-1} a\int\frac{\mathrm d\boldsymbol p}{(2\pi)^d}\ \omega_{\boldsymbol p}+\text{const.}\tag{3}
$$
where $\boldsymbol p\in\mathbb R^d$, and $\omega_{\boldsymbol p}\equiv+|\boldsymbol p|$. Moreover, the $\text{const.}$ term is a counter-term which will play no role in zeta-regularisation (even though it is important in other schemes), because it happens to vanish. We will omit it in what follows.
If we take $a\ll L$ to be a finite distance, then the momentum along this direction is quantised. For example, imposing periodic boundary conditions,
$$
\boldsymbol p=(\boldsymbol p_\perp,p_n) \qquad \text{where}\qquad p_n=\frac{n\pi}{a}\tag{4}
$$
With this, the energy becomes
$$
E_0=\frac12 L^{d-1}a\frac{\color{blue}{I(d)}}{(2\pi)^d} \tag{5}
$$
where we have defined
$$
\color{blue}{I(d)}\equiv \frac{2\pi}{a}\int\mathrm d^{d-1}\boldsymbol p_\perp\sum_{n=1}^\infty\sqrt{\left(\frac{n\pi}{a}\right)^2+\boldsymbol p_\perp^2}\tag{6}
$$
The integral above has rotational symmetry, so we change into hyper-spherical coordinates:
$$
I(d)=\frac{4}{a}\frac{\pi^{\frac{d+1}{2}}}{\Gamma\left(\frac{d-1}{2}\right)}\int_0^\infty\mathrm dp_\perp\ p^{d-2}_\perp\sum_{n=1}^\infty\sqrt{\left(\frac{n\pi}{a}\right)^2+p_\perp^2}\tag{7}
$$
where
$$
\frac{2\pi^{\frac{d-1}{2}}}{\Gamma\left(\frac{d-1}{2}\right)}\tag{8}
$$
is the area of the $(d-2)$-sphere.
We next introduce a zeta-like regulator,
$$
\sqrt{\left(\frac{n\pi}{a}\right)^2+p_\perp^2}\to \left[\left(\frac{n\pi}{a}\right)^2+p_\perp^2\right]^{x/2}\tag{9}
$$
where the physical limit is $x\to 1$. (Remark: we should also introduce a mass scale to keep the units consistent; as usual, the mass scale only affects the counter-terms, so we will omit it for simplicity).
After introducing the regulator, we can explicitly evaluate the integral (which is just the beta function):
$$
\int_0^\infty \mathrm dp_\perp\ p_\perp^{d-2}\left[\left(\frac{n\pi}{a}\right)^2+ p^2_\perp\right]^{x/2}=\left(\frac{n\pi}{a}\right)^{d-1+x}\frac{\Gamma\left(\frac{d-1}{2}\right)\Gamma\left(\frac{1-d-x}{2}\right)}{2\Gamma\left(-\frac x2\right)}\tag{10}
$$
The integral above does only converge if $\text{re}(d+x)<1$, so we must regard this result as the continuation to complex $x$. Plugging this into $I$, we get
\begin{equation}
I(d)=\frac{2\pi^\frac{d+1}{2}}{a}\frac{\Gamma\left(\frac{1-d-x}{2}\right)}{\Gamma\left(-\frac x2\right)}\color{green}{\sum_{n=1}^\infty \left(\frac{n\pi}{a}\right)^{d-1+x}}\tag{11}
\end{equation}
where again, the sum is a zeta function:
\begin{equation}
\color{green}{\sum_{n=1}^\infty \left(\frac{n\pi}{a}\right)^{d-1+x}}=\pi^{d-1+x}a^{1-d-x}\zeta(1-d-x)\tag{12}
\end{equation}
and therefore
\begin{equation}
I(d)=-\frac{\pi^\frac{3d}{2}}{a}\Gamma\left(-\frac{d}{2}\right)a^{-d}\zeta(-d)\tag{13}
\end{equation}
where we have taken $x\to 1$.
After all this, the vacuum energy is
\begin{equation}
E_0=-2^{-(d+1)}\pi^{d/2} L^{d-1}\Gamma\left(-\frac{d}{2}\right)a^{-d}\zeta(-d)\tag{14}
\end{equation}
and so the force per unit area is
\begin{equation}
\frac{F}{L^{d-1}}=-2^{-(d+1)}\pi^{d/2} d\Gamma\left(-\frac{d}{2}\right)a^{-(d+1)}\color{red}{\zeta(-d)}\tag{15}
\end{equation}
which, after a trivial simplification, coincides with the expression for the Casimir force we claimed at the beginning of this post.
We leave it to the reader to try to generalise this into other geometries; for example, we could consider a discretisation of modes along several spatial dimensions, in which case one would presumably encounter a sum of the form
$$
\sum_{\boldsymbol n} (\boldsymbol n^2)^{x/2}\tag{16}
$$
whose properties in the complex-$x$ plane seem to be rather non-trivial. The interested reader can find some useful information in this Math.SE post. We will not pursue this any further.
Best Answer
A way to do this is using regularization by substracting a continuous integral, ,with the help of the Euler-MacLaurin formula:
You can write :
$$ \sum_{Regularized} =(\sum_{n=0}^{+\infty}f(n) - \int_0^{+\infty} f(t) \,dt) = \frac{1}{2}(f(\infty) + f(0)) + \sum_{k=1}^{+\infty} \frac{B_k}{k!} (f^{(k - 1)} (\infty) - f^{(k - 1)} (0))$$ where $B_k$ are the Bernoulli numbers.
With the function $f(t) = te^{-\epsilon t}$, with $\epsilon > 0$, you have $f^{(k)}(\infty) = 0$ and $f(0) = 0$, so with the limit $\epsilon \rightarrow 0$, you will find : $$\sum_{Regularized} = - \frac{B_1}{1!} f (0) - \frac{B_2}{2!} f' (0) = - \frac{1}{12}$$
because $f(0) = 0$ and $B_2 = \frac{1}{6}$