Which values of the Riemann zeta function at negative arguments come up in physics? For my bachelor's thesis, I am investigating Divergent Series. Apart from the mathematical theory behind them (which I find fascinating), I am also interested in their applications in physics. 
Currently, I am studying the divergent series that arise when considering the Riemann zeta function at negative arguments. The Riemann zeta function can be analytically continued. By doing this, finite constants can be assigned to the divergent series. For $n \geq 1$, we have the formula: 
$$ \zeta(-n) = - \frac{B_{n+1}}{n+1} . $$
This formula can be used to find: 


*

*$\zeta(-1) = \sum_{n=1}^{\infty} n = - \frac{1}{12} . $ This formula is used in Bosonic String Theory to find the so-called "critical dimension" $d = 26$. For more info, one can consult the relevant wikipedia page. 

*$\zeta(-3) = \sum_{n=1}^{\infty} n^3 = - \frac{1}{120} $ . This identity is used in the calculation of the energy per area between metallic plates that arises in the Casimir Effect. 


My first question is: do more of these values of the Riemann zeta function at negative arguments arise in physics? If so: which ones, and in what context? 
Furthermore, I consider summing powers of the Riemann zeta function at negative arguments. I try to do this by means of 
Faulhaber's formula. Let's say, for example, we want to compute the sum of $$p = \Big( \sum_{k=1}^{\infty} k \Big)^3  . $$ 
If we set $a = 1 + 2 + 3 + \dots + n = \frac{n(n+1)}{2} $, then from Faulhaber's formula we find that $$\frac{4a^3 - a^2}{3} = 1^5 + 2^5 + 3^5 + \dots + n^5 , $$ from which we can deduce that $$ p = a^3 = \frac{ 3 \cdot \sum_{k=1}^{\infty} k^5 + a^2 }{4} .$$
Since we can also sum the divergent series arising from the Riemann zeta function at negative arguments by means of Ramanujan Summation (which produces that same results as analytic continuation) and the Ramanujan Summation method is linear, we find that the Ramanujan ($R$) or regularised sum of $p$ amounts to $$R(p) = R(a^3) = \frac{3}{4} R\Big(\sum_{k=1}^{\infty} k^5\Big) + \frac{1}{4} R(a^2) . $$
Again, we know from Faulhaber's Formula that $a^2 = \sum_{k=1}^{\infty} k^3 $ , so $R(a^2) = R(\zeta(-3)) = - \frac{1}{120} $, so $$R(p) = \frac{3}{4} \Big(- \frac{1}{252} \Big) + \frac{ ( \frac{1}{120} )} {4} = - \frac{1}{1120} .  $$
My second (bunch of) question(s) is: Do powers of these zeta values at negative arguments arise in physics? If so, how? Are they summed in a manner similar to process I just described, or in a different manner? Of the latter is the case, which other summation method is used? Do powers of divergent series arise in physics in general? If so: which ones, and in what context?
My third and last (bunch of) question(s) is: which other divergent series arise in physics (not just considering (powers of) the Riemann zeta function at negative arguments) ?
I know there are whole books on renormalisation and/or regularisation in physics. However, for the sake of my bachelor's thesis I would like to know some concrete examples of divergent series that arise in physics which I can study. It would also be nice if you could mention some divergent series which have defied summation by any summation method that physicists (or mathematicians) currently employ. Please also indicate as to how these divergent series arise in physics. 
I also posted a somewhat more general and improved version of this question on MO.
 A: 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.
