What parameter can you expand QFT in for a convergent series? It is known that the expansion in terms of Feynman diagrams with a series in terms of the coupling constant is an asymptotic series (i.e. it starts off a good approximation but eventually diverges). (After renormalisation is taken into account)
That said, is there a way to expand the formula in terms of another parameter?
e.g. take the integral:
$$I(a,b)=\int\limits_{-\infty}^\infty e^{-a x^2 - b x^4 } dx\tag{1}$$
Expanding this in the parameter $b$ gives an asymptotic series. But we could also expand this in the parameter $a$ and surely this would converge as $x^4$ will always grow bigger than $x^2$ even if it has a small coefficient.
In terms of QFT an expression such as (Wick-rotated $\phi^4$):
$$f(x,y;m,\lambda)=\int \phi(x)\phi(y)e^{\int(-\phi(z)(\partial^2+m^2)\phi(z) - \lambda \phi(z)^4) d^4z} D\phi\tag{2}$$
could perhaps be expanding using $\partial^2+m^2$ as the parameter?
Is this possible? And if so are there disadvantages in this, i.e. would the series converge too slowly as to be useless?
Thus is it true to say that for a small coupling constant, the asymptotic series is the best method? Has no one worked out a way to have best of both worlds; a series that converges quickly but is non-asymptotic?
 A: The glib answer would be to suggest you expand in powers of $1$ so the path integral $Z$ admits expansion $Z = Z$ in powers of $1$. Indeed, this expansion not only converges, but is given in closed form!
The point I am making here is that it's not enough to just find "another parameter." You also need to find a parameter such that the expansion coefficients are calculable themselves. The reason we don't simply expand in powers of $1$ is obvious: computing the coefficients would require us to solve a problem we don't know how to solve.
With all that said, let me point out a couple of things about OP's statements. Firstly, it does not make sense to expand in power of $\partial^2+m^2$ as it is not small. This becomes clear if we were to transform to momentum space where the operator in question is $-k^2+m^2$ and we are integrating over $k$. So not only is $k$ an integration parameter (and hence does not make sense outside the integral of the action), but it's also arbitrarily large...not a very good candidate for an expansion parameter.
To put this statement in terms of the simple example OP gives, which rightly points out (so far as I'm aware) that the expansion of the 1-dimensional integral in terms of $a$ has a non-zero radius of curvature, you would need to generalize the 1-dimensional integral to an $N$-dimensional integral. There the coefficient $a$ would need to be a matric $A$, though we can always diagonalize $A$ and write this as a sum over squares. But now your "expansion in $a$" is an "expansion in the eigenvalues of $A$." Such an expansion obviously only makes sense if the eigenvalues of $A$ are small. What I pointed out about the operator $\partial^2+m^2$ is that it's spectrum of eigenvalues is unbounded, and hence there are eigenvalues in the spectrum which make no sense to expand around.
Next, I will point out that there's a subtle difference between the expansion in powers of coupling and loop order. Note that all perturbative calculations are done in terms of loop order rather than coupling order (the precise statement is that depending on the theory, there are bounding on the coupling order at a fixed loop order, so the two aren't that different).
The parameter $\hbar$ in the exponential of the path integral serves as a useful order-tracking parameter for loop orders (which is why the loop order expansion is sometimes called the $\hbar$ expansion as in WKB). You can find the details of this counting argument in various places, but in particular, I know it's detailed in the "External Field Methods" chapter of Weinberg's QFT volume 2.
The expansion about $\hbar$ small is equivalent to steepest descent approximation, which is a generic integral approximation method well-known to produce asymptotic series.
The only other general approximation scheme I'm aware of (though I don't know anything about convergence properties myself) would be lattice methods. These are obviously numerical rather than analytic, but they are also valid in strong-coupling regimes, which is why they are used for QCD calculations.
A: *

*For what its worth, by change of the integration variable $x$, one may show that OP's 1st integral
$$I(a,b)~=~\frac{I(1,\frac{b}{a^2})}{\sqrt{a}}~=~\frac{I(\frac{a}{\sqrt{b}},1)}{\sqrt[4]{b}}. $$
Depending of whether one is interested in the weak coupling limit $|\frac{b}{a^2}|\ll 1$ or the strong coupling limit $|\frac{a}{\sqrt{b}}|\ll 1$, one may want to choose different perturbative expansions.


*Similarly for OP's 2nd integral. Of course a great advantage of Gaussian integrals is that they can easily be evaluated exactly. That's why  perturbation theory is usually performed around a free theory.
