Suppose you want to compute the sum $\sum_{n=0}^{\infty} a_n $
You could consider the expression $f(x) = \sum_{n=0}^{\infty} e^{a_n x}$ and try to compute the coefficient of an $x^1$ term in the asymptotic series expansion of $f(x)$.
To see this in practice we can consider $\sum_{n=0}^{\infty} n $. We then consider the series $\sum_{n=0}^{\infty} e^{nx}$ which is equal to $-\frac{1}{x} + \frac{1}{2} – \frac{1}{12}x \dots $
We then make the jump and state $0 + 1 + 2 + 3 + 4 … = -\frac{1}{12}$
We can also try a divergent geometric series that everyone knows how to sum. Suppose we wish to compute $1 + 2 + 4 + 8 + 16 … = \sum_{n=0}^{\infty}{2^n}$
Some calculations here by user 2734364041 tell us that
$$ \sum_{n=0}^{\infty} e^{-\frac{2^n}{x}} = \frac{\log x-\gamma}{\log 2}+\frac{1}{2}+\frac{2}{\log 2}\sum_{j=1}^{\infty}\Re\Big(x^{\frac{2\pi ij}{\log 2}}\Gamma\Big(\frac{2\pi i j}{\log 2}\Big)\Big)-\sum_{j=1}^{\infty}\frac{(-1)^{j}}{j!}\frac{x^{-j}}{2^j-1} $$
It then follows that:
$$ \sum_{n=0}^{\infty} e^{2^n x } = \frac{-\log -x-\gamma}{\log 2}+\frac{1}{2}+\frac{2}{\log 2}\sum_{j=1}^{\infty}\Re\Big(x^{\frac{2\pi ij}{\log 2}}\Gamma\Big(\frac{2\pi i j}{\log 2}\Big)\Big)-\sum_{j=1}^{\infty}\frac{1}{j!}\frac{x^{j}}{2^j-1} $$
It's easy now to see that the coefficient on the "x" term is $-1$ so we conclude that $1 + 2 + 4 + 8 + 16 … = -1$ which is consistent with the geometric series result.
The nice thing about this trick is it seems to be able to squeeze sums that are not analytically tractable in the usual sense. Suppose we wish to sum $\sum_{n=0}^{\infty} 2^{2^n}$
We want to extract the $x$ coefficient of $\sum_{n=0}^{\infty} e^{(2^{2^n}x)} $ this amounts to basically first trying to generate an asymptotic expansion of the expression
$$ \sum_{n=0}^{\infty} e^{-\frac{2^{2^n}}{x}} $$
But this is an extremely slowly growing function so we can turn our attention instead to:
$$ \sum_{n=0}^{\infty} e^{-\frac{2^{2^n}}{2^{2^x}}} = \sum_{n=0}^{\infty} e^{-2^{(2^n – 2^x )}} $$
Per Desmos this expression appears to grow like the step function $\left \lfloor x \right \rfloor $ In fact they differ by only a constant for at least 15 decimal digits (before I ended up breaking Desmos's arithmetic abilities).
So while it'll take quite a bit of effort to wrap this up I want to pre-emptively ask, is there any reason this technique cannot be used for summation of currently unsummable divergent series that consist of monotonically increasing terms? This technique appears to agree with our usual analytic continuation and it looks like it has a shot of working on any series of entirely positive terms whose terms grow sufficiently fast.
More generally if you want to compute $\sum_{n=0}^{\infty} a_n$ then if you can find meromorphic functions $f$ and $g$ such that the first order $x^1$ term in the asymptotic expansion of $f$ isn't 0 and such that $\sum_{n=0}^{\infty} f(a_n g(x))$ has an asymptotic expansion itself, then you can assign the coefficient of the $g(x)^1$ term of the sum divided by the coefficient of the $x^1$ term in $f$ as the "divergent sum" of the series. In our case we just focus on using $f = e^x, g = x$ and looking at the $x$ term $\sum_{n=0}^{\infty} e^{a_n x}$ but as we show in the question it's usually more convenient to consider $f = e^x, g = -\frac{1}{x}$ and to look at the negative of the $\frac{1}{x}$ term of this expansion.
Best Answer
This summation method gives answers that are close to, but do not always match, traditional divergent summation methods. For instance, for constants $a,b>0$, the divergent sum $$ \sum_{n=1}^\infty (an+b) $$ evaluates to $-\frac{a}{12}-\frac{b}{2}$ by most standard divergent summation methods (such as zeta function regularization), but from the asymptotic $$ \sum_{n=1}^\infty e^{-(an+b)x} = \frac{e^{-bx}}{e^{ax}-1} = \frac{1}{ax} - \frac{b}{a} - \frac{1}{2} + (\frac{b^2}{2a} + \frac{b}{2} + \frac{a}{12}) x + \dots$$ we see that this alternate summation method gives the different answer of $-\frac{b^2}{2a} - \frac{b}{2} - \frac{a}{12}$. As observed by OP, the two answers match when $b=0$, but not in general.
This example also illustrates that this method is not a linear method of summation, in contrast to the usual methods. This seems to diminish the potential utility of this summation method for applications.
EDIT: The situation can be clarified by the Euler-Maclaurin formula. For $P(t)$ a non-constant polynomial with positive coefficients, this formula gives (after some additional work) $$ \sum_{n=1}^\infty e^{-P(n) x} = \int_0^\infty e^{-P(t) x}\ dt - \frac{1}{2} + \sum_{k=1}^\infty \frac{B_k}{k!} P^{(k-1)}(0) x + O(x^2)$$ for small positive $x$, where $B_k$ is the $k^{th}$ Bernoulli number (with the $B_1 = 1/2$ sign convention). Note that the expression $-\sum_{k=1}^\infty \frac{B_k}{k!} P^{(k-1)}(0)$ is the usual zeta function regularized value for $\sum_{n=1}^\infty P(n)$. So the discrepancy between the two methods here comes from a possible linear (in $x$) term in the integral $\int_0^\infty e^{-P(t) x}\ dt$. In the monomial case $P(t) = t^d$, this integral is proportional to $x^{-1/d}$ and so does not contribute any linear term, but for more general polynomials it could generate some (presumably unwanted) linear terms.