Computation of a series.
First, we will transform $S_n$ to a form easier to manipulate.
Let $C \subset \mathbb{C}$ be a circle of radius $r \ll 1$ centered at $0$. For any $n \in \mathbb{N}$, $m \in \mathbb{N}^n$, we can single out those $m \in \mathfrak{M}_n$ with help of contour integrals of the form:
$$\delta_n(m) \stackrel{def}{=} \frac{1}{2\pi i} \oint_{C} s^{\sum_{k=1}^n k m_k} \frac{ds}{s^{n+1}} = \begin{cases}1, & m \in \mathfrak{M}_n\\ 0, & \text{ otherwise }\end{cases}$$ Together with following integral representation of factorial:
$$n! = \Gamma(n+1) = \int_0^\infty t^n e^{-t}dt$$ We have
$$\begin{align} S_n &= \sum_{m \in \mathbb{N}^n} \delta_n(m) \int_0^\infty \prod_{k=1}^n \frac{1}{m_k!} \left(\frac{t}{k+1}\right)^{m_k} t^n e^{-t} dt\\ &= \frac{1}{2\pi i}\sum_{m \in \mathbb{N}^n} \oint_C \left[\int_0^\infty \prod_{k=1}^n \frac{1}{m_k!} \left(\frac{ts^k}{k+1}\right)^{m_k} \left(\frac{t}{s}\right)^ne^{-t} dt \right] \frac{ds}{s}\\ &= \frac{1}{2\pi i}\oint_C \left[\int_0^\infty \prod_{k=1}^n \left(\sum_{m_k=0}^\infty \frac{1}{m_k!} \left(\frac{ts^k}{k+1}\right)^{m_k}\right) \left(\frac{t}{s}\right)^ne^{-t} dt \right] \frac{ds}{s}\\ &= \frac{1}{2\pi i}\oint_C \left[\int_0^\infty \prod_{k=1}^n \exp\left(\frac{ts^k}{k+1}\right) \left(\frac{t}{s}\right)^ne^{-t} dt \right] \frac{ds}{s}\\ &= \frac{1}{2\pi i}\oint_C \left[\int_0^\infty \exp\left(\sum_{k=1}^n\frac{ts^k}{k+1}\right) \left(\frac{t}{s}\right)^ne^{-t} dt \right] \frac{ds}{s}\\ &\stackrel{\color{blue}{[1]}}{=} \frac{1}{2\pi i}\oint_C \left[\int_0^\infty \exp\left(\sum_{k=1}^\infty\frac{ts^k}{k+1}\right) \left(\frac{t}{s}\right)^ne^{-t} dt \right] \frac{ds}{s}\\ &= \frac{1}{2\pi i}\oint_C \left[\int_0^\infty \exp\left[-t\left(\frac{\log(1 - s)}{s} + 2\right)\right]\left(\frac{t}{s}\right)^n dt \right] \frac{ds}{s}\tag{*1}\\ \end{align} $$ Next, let
- $S(x) \stackrel{def}{=} \sum_{n=0}^\infty S_n \frac{x^n}{n!}$ be the EGF (exponential generating function) for $S_n$.
- $\Delta(x) = \sum_{n=0}^\infty S_n \frac{x^{n+1}}{(n+1)!}$ be the series we want to study its convergence.
They are related by the relation $\Delta(x) = \int_0^x S(t) dt$.
For any $x$ with $|x| \ll r$, $(*1)$ implies
$$\begin{align} S(x) &= \frac{1}{2\pi i}\oint_C \left[\int_0^\infty \exp\left[-t\left(\frac{\log(1 - s)}{s} + 2 - \frac{x}{s}\right)\right] dt \right] \frac{ds}{s}\\ &= \frac{1}{2\pi i} \oint_C \frac{ds}{\log(1-s) + 2s - x} \end{align} $$ Change variable to $y = -\log(1-s) \iff s = 1 - e^{-y}$. When $r$ is small, the image of $C$ in $y$-space is close to circle $C$. We can deform the contour back to $C$ without changing the integral. This leads to
$$S(x) = \frac{1}{2\pi i}\oint_C \frac{dy}{P(x,y)} \quad\text{ where }\quad P(x,y) = (2-x-y)e^y - 2 $$ Under the condition $|x| < |y| = r \ll 1$, we have
$$P(x,y) \approx (2 - x - y) (1 + y + O(r^2)) - 2 \approx y - x + O(r^2)$$
This means for fixed $x$ and as a function in $y$, $P(x,y)$ has only one root inside $C$. Furthermore, the root in $y$ is close to $x$. Let $\eta$ be that root, we have
$$\begin{align} P(x,\eta) = 0 &\iff (2-x-\eta)e^\eta - 2 = 0 \iff (\eta + x - 2)e^{\eta + x - 2} = -2e^{x-2}\\ & \implies 2 - x - \eta = -W(-2e^{x-2}) \end{align} $$ where $W(z)$ is a branch of the Lambert-W function. In terms of $\eta$, we have
$$\begin{align} S(x) &= \text{Res}_{y=\eta}\left(\frac{1}{P(x,y)}\right) = \left.\frac{1}{\frac{\partial}{\partial y}P(x,y)}\right|_{y=\eta} = \frac{1}{(1 - x - \eta)e^\eta}\\ &= \frac{2-x-\eta}{2(1-x-\eta)} = \frac{W(-2e^{x-2})}{2(1+W(-2e^{x-2}))} \end{align} $$
Since $S(0) = 1$, we need to choose a branch of Lambert W function with $W(-2e^{-2}) = -2$. The correct branch is the "lower branch" described in above wiki link. It is usually denoted as $W_{-1}(\cdot)$. In terms of it, we find
$$S(x) = \frac{W_{-1}(-2e^{x-2})}{2(1+W_{-1}(-2e^{x-2}))}$$
Notice the branches of Lambert W function satisfies ODE $$z\frac{d}{dz}W(z) = \frac{W(z)}{1+W(z)}\tag{*2}$$ We can integrate $(*2)$ and deduce a closed form expression for $\Delta(x)$: $$\Delta(x) = \frac12 \int_0^x \left[ z\frac{dW_{-1}(z)}{dz} \right]_{z=-2e^{t-2}} dt = 1 + \frac12 W_{-1}(-2e^{x-2})\tag{*3}$$
$W_{-1}(z)$ has two branch cuts, one terminated at $z = -\frac1e$, the other at $z = 0$. The closest singularity of $\Delta(x)$ to origin is located at $x = 1 - \log(2)$. As a result, $r_0$, the radius of convergence of the power series expansion of $\Delta(x) = \sum\limits_{n=0}^\infty S_n\frac{x^{n+1}}{(n+1)!}$, $r_0$ equals to $1 - \log(2)$. A corollary of this is$\color{blue}{{}^{[2]}}$
$$\frac{S_n}{(n+1)!} \sim o(\rho^n)\quad\text{ for any }\; \rho > \frac{1}{1-\log(2)} \approx 3.258891353270929$$
As a double check, we evaluate the power series expansion of $\Delta(x)$
using following command Series[1+1/2*LambertW[-1,-2*Exp[x-2]],{x,0,8}]
on WA (wolfram alpha). WA returns
$$\begin{align} \Delta(x) = & x+\frac{{x}^{2}}{2}+\frac{5\,{x}^{3}}{6}+\frac{41\,{x}^{4}}{24}+\frac{469\,{x}^{5}}{120}+\frac{6889\,{x}^{6}}{720}\\ & +\frac{24721\,{x}^{7}}{1008}+\frac{2620169\,{x}^{8}}{40320}+\frac{64074901\,{x}^{9}}{362880} + \cdots \end{align}$$ Translate back to $S_n$, this is equivalent to
$$( S_0,S_1,\ldots ) = (1, 1, 5, 41, 469, 6889, 123605, 2620169, 64074901,\ldots )$$ For $n \le 5$, I have checked by hand this is indeed the correct value.
An OEIS search return the sequence OEIS A032188. Up to $n = 18$, I've verified the $S_n$ extract from expansion of $(*3)$ matches the numbers on OEIS. Look at references there and see whether there is anything useful for your purposes.
Notes
$\color{blue}{[1]}$ - As a function of $s$, $$\exp\left(\sum_{k=1}^n\frac{ts^k}{k+1}\right) \left(\frac{t}{s}\right)^n\frac{e^{-t}}{s} = \frac{1}{s^{n+1}}A(s)\quad\text{ and }\quad \exp\left(\sum_{k=n+1}^\infty\frac{ts^k}{k+1}\right) = 1 + s^{n+1}B(s) $$ where $A(s), B(s)$ are analytic over the disc bounded by $C$. This implies $$\exp\left(\sum_{k=1}^\infty\frac{ts^k}{k+1}\right) \left(\frac{t}{s}\right)^n\frac{e^{-t}}{s} = \exp\left(\sum_{k=1}^n\frac{ts^k}{k+1}\right) \left(\frac{t}{s}\right)^n\frac{e^{-t}}{s} + A(s)B(s)$$ Changing the upper bound in the sum within the exponent from $n$ to $\infty$ modifies the integrand by a function analytic over the disc bounded by $C$. The value of the contour integral over $C$ remains the same.
$\color{blue}{[2]}$ - A more detailed analysis suggests for large $n$, $S_n$ has following approximation: $$S_n \approx \frac{(2n)!}{\sqrt{8r_0}n!(4r_0)^n}\left( 1 - \frac{r_0}{6(2n-1)} + \cdots \right)\quad\text{ where }\quad r_0 = 1 - \log(2)$$ For $n$ as small as $4$, this formula gives a relative error below $10^{-4}$ (checked against numbers from OEIS). The leading behavior of coefficients of $\Delta(x)$ should be: $$\frac{S_n}{(n+1)!} \sim O\left(\frac{r_0^{-n}}{\sqrt{8\pi r_0} n^{3/2}}\right)$$