Perturbation Theory and Thermodynamic Limit
As far as I know, in many-body condensed-matter theory, a simplistic finite-order Taylor expansion is (almost) never used, unless for pedagogical reasons. In fact, one always needs to “organize the series expansion more cleverly” [as higgsss commented] in form of infinite perturbative resummations or asymptotic expansions.
A simple Taylor-expansion, as in $ e^{-x} = 1 - x + \frac{1}{2}x^2 + \mathcal{O}(x^3) $, works only for some single-particle problems in physics.
In words of Ref. [1] (§5.1, p. 195):
The moral to be taken ... is that perturbative expansions should not be confused with rigorous Taylor expansions. Rather they represent asymptotic expansions, in the sense that, for weaker and weaker coupling, a partial resummation of the perturbation series leads to an ever more precise approximation to the exact result. For weak enough coupling the distinction between Taylor expansion and asymptotic expansion becomes academic (at least for physicists). However, for intermediate or strong coupling theories, the asymptotic character of perturbation theory must be kept in mind.
Consult §5.1 of Ref. [1] for an extensive discussion.
In some contexts, eg. in the perturbative renormalization group, one actually “re-exponentiates” the finite-order perturbative expansion; for instance,
$$ Z = \sum_{\text{states}} e^{-\beta (E_0 + g \Delta E)} = Z_0 \, \langle e^{-\beta \, g\Delta E} \rangle_0 \approx Z_0 \, \langle 1 - \beta g \Delta E \rangle_0 \approx Z_0 \, e^{- \beta g \langle \Delta E \rangle} ~, $$
so that the problem observed in the question does not appear at all; see, for instance, §8.2 of Ref. [1].
[1] Altland, A., and B. D. Simons. “Condensed matter field theory”. CUP (2010) [wcat].
As you noticed, the first thing to do for the thermodynamic limit is to consider quantities which have a chance of making sense in this limit. The partition function $Z$ is not good for that, but the limit of $\frac{1}{N}\log Z$ (related to the free energy density or the pressure) is indeed the right quantity. Correlation functions also work. Using your notations, one has $$ Z=Z_0 \langle e^{-g\beta\Delta H}\rangle_0 $$ so $$ \log Z=\log Z_0+\log \langle e^{-g\beta\Delta H}\rangle_0 =\log Z_0+\sum_{m=1}^{\infty}\frac{(-g\beta)^m}{m!} \langle \Delta H,\ldots,\Delta H\rangle_0^{\rm T} $$ where $\langle \Delta H,\ldots,\Delta H\rangle_0^{\rm T}$ is the $m$-th joint cumulant of the random variable $\Delta H$ with itself (or simply the $m$-th cumulant of $\Delta H$). The underlying probability measure corresponds to expectations $\langle\cdots\rangle_0$. If you expanded the logarithm as one of the comments suggested, you would have seen, e.g., that the quadratic term would involves the variance $$ \langle \Delta H,\Delta H\rangle_0^{\rm T}=\langle (\Delta H)^2\rangle_0- \langle \Delta H\rangle_0^2\ . $$ One may naively think that, since $\Delta H$ is intensive, the cumulant $$ \langle \Delta H,\ldots,\Delta H\rangle_0^{\rm T} $$ scales like $N^m$, but this is not true. It scales like $N$ for any $m$. There are several ingredients coming together for this property. The first is that a joint cumulant $\langle X_1,\ldots,X_m\rangle^{\rm T}$ is always zero if the random variables $X_1,\ldots,X_m$ can be arranged in two groups which are statistically independent. This forces a "connectedness" property. The second is that a decent $\Delta H$ should have a built-in decay property in the size and spatial extent of the support of the interaction. More concretely, take the Ising model with possibly long-range interactions. Then $H_0=0$ (if you want a less trivial $H_0$ turn on a magnetic field), and $$ \Delta H=\sum_{x,y}J_{x,y} \sigma_x\sigma_y $$ for some translation invariant couplings that satisfy a decay condition when $x,y$ get far apart. In particular, one wants a condition like $$ \sum_y |J_{x,y}|<\infty $$ for $x$ fixed, otherwise $\Delta H$ would not be "extensive". Inserting the definition of $\Delta H$ in the cumulant, you see that you have to sum over a configuration of $m$ edges thrown into the lattice. They must form a connected graph. One basically has to perform a sum $$ \sum_{x_1,y_1}\cdots\sum_{x_m,y_m} J_{x_1,y_1}\cdots J_{x_m,y_m} \langle \sigma_{x_1}\sigma_{y_1},\ldots,\sigma_{x_m}\sigma_{y_m} \rangle_0^{\rm T}\ . $$ This results in an overall factor of $N$ to pick the location of $x_1$ say.Then by translation invariance you can assume that this point is pinned at the origin say. But you have to sum over the $2m-1$ remaining points. This is where the connectedness and summability condition on the $J$'s will save you.
This is a simple case of a very general method called a cluster expansion.
To learn more about this, have a look a the notes I wrote
- "Notes on the cluster expansion for the polymer gas, a.k.a., the Mayer expansion."
- "Notes on the Brydges-Kennedy-Abdesselam-Rivasseau forest interpolation formula" needed for Lemma 1 in the previous note and for more sophisticated versions of the cluster expansion.
This was from a graduate course I taught a while ago.
I think the problem arises because you assume that $\langle \Delta H^2\rangle_0$ is in the same order as $\langle \Delta H\rangle_0^2$, i.e in $O(N^2)$. In fact, for any canonical ensemble the energy fluctuation, $$\frac {\langle H^2\rangle}{\langle H\rangle^2}\propto \frac{1}{N},$$ (see, for example in here). Thus, following your perturbation case, $\langle \Delta H^2\rangle_0$ should be in $O(N)$; you can check this at the higher orders (notice that this has some connections to Abdelmalek's answer).
Therefore for any $\left|\frac{g\Delta H}{H_0}\right|\ll1$ the partition function should converge .