$\def\L{\operatorname L} \def\M{\operatorname M} $ Our goal is inverting $f(x)=\tanh(x)+ax$ via Laplace inversion. The Laplace transform of $h(x)=f^{-1}(x)$ is:
$$\L_s(h(t))=\int_0^\infty e^{-st}h(t)dt\mathop=^{h(t)\to t}\int_0^\infty t e^{-s f(t)} df(t)=\frac1s\int_0^\infty e^{-s(\tanh(t)+at)}dt\mathop=^{\tanh(t)\to t}\frac1s\int_0^1e^{-st}(1-t)^{\frac{as}2-1}(1+t)^{-\frac{as}2-1}dt \tag1$$
This integral matches the Whittaker M function integral representation:
$$\frac{\Gamma\left(u-k+\frac12\right) \Gamma\left(u+k+\frac12\right)}{(2u)!z^{u+\frac12}2^{-2u}}\M_{k,u}(z)=\int_{-1}^1 e^\frac{zt}2 (1+t)^{u-k-\frac12}(1-t)^{u+k-\frac12}$$
but the bounds are different. A bilateral Laplace transform gives $\int_{-1}^1$ bounds after using the same substitutions, but the resulting integral may diverge. The similar $\int_0^{\frac{\pi}{2}} e^{-x\tan t+\alpha t} \;dt $ almost certainly does not have a closed form with the confluent hypergeometric function, of which Whittaker M function is one. Humbert series:
$$\Phi_1(a,b,c;x,y)=\frac{\Gamma(c)}{\Gamma(a)\Gamma(c-a)}\int_0^1 e^{yt}t^{a-1}(1-t)^{c-a-1}(1-xt)^{-b}dt$$
would work, but it is a last resort. Also, as implied by “Conditions to apply the inverse Laplace transform”, $\L_s(h(t))$ is finite, so the inverse Laplace transform converges. Thus we likely have an integral representation for $f^{-1}(x)$ if $(1)$ is evaluated.
Does $(1)$ have a simpler closed form than with the Humbert series?
Hopefully someone will find a simpler function than the Humbert series, but expanding with $e^y$’s Maclaurin series, the Gauss hypergeometric function, and the Pochhammer symbol gives: $$\frac1s\int_0^1e^{-st}(1-t)^{\frac{as}2-1}(1+t)^{-\frac{as}2-1} dt=\frac2{as^2}\sum_{n=0}^\infty\frac{(1)_n(-s)^n\,_2\text F_1\left(n+1,\frac{as}2+1;n+\frac{as}2+1;-1\right)}{\left(\frac{as}2+1\right)_nn!}$$
whereas
The $\Phi_1$ integral representation, from the question, or single sum definition analytically continues the double sum outside $|x|,|y|<1$. Therefore:
$$\boxed{f(x)=\tanh(x)+ax\implies f^{-1}(x)=\frac1{\pi ia}\int_{c-i\infty}^{c+i\infty}\frac{e^{st}}{s^2}\Phi_1\left(1,\frac{as}2+1,\frac{as}2+1;-1,-s\right)ds}$$
The Mathematica code
InverseLaplaceTransform[Sum[((-s)^n Gamma[(a* s)/2] Hypergeometric2F1[1 + n, 1 + (a *s)/2, 1 + n + (a*s)/2, -1])/(s Gamma[1/2 (2 + 2 n + a *s)]),{n,0,b}],s,x]for $b\to\infty$ and certain $a,x$ confirms this: