In physics, when we calculate the renormalized sum of $S=\sum_{n=1}^\infty n$, we usually use an exponential regulator and instead first calculate
$$S_\epsilon = \sum_{n=1}^\infty ne^{-\epsilon n} = \frac{e^\epsilon}{(e^\epsilon -1)^2}. $$
There's nothing controversial about this, but an example by Lubos Motl can be found here.
Now, if we expand this around small $\epsilon>0$ we obtain that
\begin{align*} S_\epsilon &= \frac{1}{\epsilon^2} - \frac{1}{12} + {\cal O}(\epsilon^2)\\ &=\zeta(-1) + \frac{1}{\epsilon^2} + {\cal O}(\epsilon^2). \end{align*}
So that when we take $\epsilon \to 0$ (in which limit $S_\epsilon \to S)$, we obtain that
$$S \to \zeta(-1) + \text{divergent piece.}$$
Now, I know that this "decomposition", if you will, into a convergent and divergent part is independent of the sum regulator, and hence always yields the same result. Furthermore, the domain of validity of the series representation of the $\zeta$ function is only when $\mathbb{R}(s) >0$. Therefore, we ought to expect it to be of no use to us when naively applying it in regions where $\mathbb{R}(s)<1$.
However, it is useful. Namely, its renormalized sum gives us the correct value for the analytic continuation of the $\zeta$ function. Furthermore, since the analytic continuation of a complex function is unique I feel there must be a connection between these two things.
Question: How does the series representation "know" about its analytic continuation? That is, what exactly is the connection between the naive summation and its analytic continuation?