When I check the matlab function kummerU, it is said that for kummerU function: $$ U(a,b,z)=\frac{1}{\Gamma(a)}\int\limits_{0}^{\infty}e^{-z t}t^{a-1}(1+t)^{b-a-1}dt $$ But only for $\operatorname{Re}(a)>0$ and $\operatorname{Re}(z)>0$.
However, I want to find the integration formula for negative integer $a$ and imaginary $z$: $$ a\in \mathbb{Z}^-\quad \text{ and } \quad z=xi,\; x\in \mathbb R. $$ But what I found is just a summation formula (EDITED, thanks to Dr. Hintze): $$U(n,m,x)=(-1)^n(m)_{-n}\sum_{j=0}^{-n}\dfrac{(n)_j}{(m)_j(1)_j}x^j\quad $$ Any suggestion is highly appreciated!
§1 Drawing on literature
EDIT
Unfortunately, previously I have misunderstood the Wikipedia entry so that the polynomial formula was wrong. This is corrected here using another source.
The behaviour of the function defined for $\Re(a) \gt 0, \Re(z) \gt 0$ by the integral
$$u(a,b,z) = \frac{1}{\Gamma (a)}\int_{0}^{\infty}t^{a-1}e^{- t z} (t+1)^{-a+b-1}\tag{1}$$
can be desribed in the context of the OP as follows:
For negative integer $a$ the function reduces to a polynomial of degreen $n$ in the variables $b$ and $z$.
Indeed, from https://en.wikipedia.org/wiki/Laguerre_polynomials#Generalized_Laguerre_polynomials we find that $u$ for non-positive integer $a = -n$ is a multiple of a generalized Laguerre polynomia, explicitly
$$\boxed{u(-n,b,z) = (-1)^n n! \sum _{k=0}^n (-1)^k \binom{b+n-1}{n-k}\frac{z^k }{k!}}\tag{2}$$
The first few polynomials in the in the format {$n$,$u$} are
$$\begin{array}{l} \{0,1\} \\ \{1,z-b\} \\ \left\{2,b^2-2 b z+b+z^2-2 z\right\} \\ \left\{3,-b^3+3 b^2 z-3 b^2-3 b z^2+9 b z-2 b+z^3-6 z^2+6 z\right\} \\ \end{array}$$
Notice that $(2)$ is also a correction to the formula given in the OP.
§2 Analytic continuation of the integral to negative values of $a$
Although the problem of the OP was solved in §1 it might be interesting to see how an analytic continuation to negative $a$ can be performed.
To this end we continue the integral expression analytically with respect to the parameter $a$ still assuming $\Re(z) \gt 0$.
As noticed in the OP the integral $(1)$ is defined only for $\Re(a) \gt 0$ because of the divergence at $t=0$ and $\Re(z)\gt 0$ because of the divergence at $t=\infty$.
As suggested in my comment, we can perform an anaytic continuation to negative values of $a$ in a certain interval by integrating the integral $u$ by parts
$$\int_{r}^{s}U' V \;dt= U V |^s_r- \int_{r}^{s} U V'\;dt$$
letting
$$\begin {align} U=\int t^{a-1}\;dt = t^a/a\\V=\frac{1}{\Gamma (a)}e^{- t z} (t+1)^{-a+b-1}\end {align} \tag{4}$$
$U'V$ gives back the integrand.
Now the integrated part in the region $(a \gt 0 , z \gt 0)$ vanishes and we find for the remaining integral $\int_{r}^{s} U V'\;dt$ the following expression
$$u(a,b,z) = - \frac{1}{\Gamma (a)}\int_{0}^{\infty}t^a e^{-t z} \left((-a+b-1) (t+1)^{-a+b-2}\\-z (t+1)^{-a+b-1}\right)\tag{5}$$
Now the integral is convergent down to $a>-1$. Writing this as a sum of two integrals and using the definition $(1)$ of $u$ we have found that
$$u(a,b,z)=(a-b+1) u(a+1,b,z)+z\; u(a+1,b+1,z)\tag{6}$$
where $u$ stands for the defining integral.
Repeating this process using $(6)$ we obtain integral representations valid in aribtrary negative intervals of $a$.
Specifically, for negative integer $a$, we find from $(6)$, assuming for a moment that $u(0,b,z)=1$ recursively the following first expressions
$$\begin{align} u(a=-1,b,z) &=-b+z\\ u(a=-2,b,z) &= (-b-1) u(-1,b,z)+z\; u(-1,b+1,z) \\ &= (-b-1)(-b+z) + z(-b-1+z)\\ &= b^2- b z + b - z-b z - z + z^2 \\ &= b^2-2 b z + b -2 z + z^2\\ &= (b-z)^2+ b -2 z \\ \end{align}$$
and so on, and it can be easily verified that the sum $(2)$ satisfies the recursion relation $(6)$, which, in turn, proves the validity of $(2)$
Finally, let us fill the last gap which is to prove that $u(0,b,z)=1$. Indeed, consider the integral $(1)$ and notice that for $a \to 0$ the term $\frac{1}{\Gamma(a)}$ gives a factor $a$ in front of the integral which in turn diverges. Hence for $u$ we have an expression of the form $0 \times \infty$. Now noticing that for $a\to 0$ the main contribution to the integral comes from the term $t^{a-1}$ in the region $t\simeq 0$ we set $t=0$ in all other terms of the integrand so that the integral gives trivially $\frac{1}{a}$. Hence $u(a \to 0,b,z) \simeq \frac{1}{a \Gamma(a)} = \frac{1}{ \Gamma(a+1)}\to \frac{1}{ \Gamma(1)}= 1$. QED.
Summarizing we have solved the problem of the OP in a self-contained manner based only on the integral $(1)$ without using unproved formulas from the literature.
§3 Miscellaneous
(i) $z=0$
Changeing variables $t\to \frac{s}{1-s}$ , $dt \to \frac{ds}{(1-s)^2} $ transforms $(1)$ into
$$u_s(a,b,z) = \frac{1}{\Gamma (a)}\int_0^1 s^{a-1} (1-s)^{-b} e^{-\frac{s z}{1-s}} \, ds\tag{7}$$
Here we can safely set $z=0$ which leads to
$$u_s(a,b,0) = \frac{1}{\Gamma (a)}\int_0^1 s^{a-1} (1-s)^{-b} \, ds= \frac {\Gamma(1-b)}{\Gamma(a+1-b)}\tag{8}$$
where we have used the definition of the Beta function and the relation
$$B(r,s) = \int_{0}^{1}x^{r-1}(1-x)^{s-1}\;dx = \frac{\Gamma(r) \Gamma(s)}{\Gamma(r+s)}$$
(ii) Special case $b=1+a$
Here the integral reduces to
$$\begin{align} u(a,1+a,z)=\frac{1}{\Gamma(a)}\int_{0}^{\infty} t^{a-1} e^{-z t}\;dt \\\overset{zt\to u}=z^{-a}\frac{1}{\Gamma(a)}\int_{0}^{\infty} u^{a-1} e^{-u}\;du=z^{-a} \end{align}\tag{9} $$
Except for integer $a$ this function has a branch cut from $z=0$ to $z=-\infty$.
Here is a plot of real and imaginary part of $u$ as a function of complex $z$ for fairly general parameters $a$ and $b$
The branch cut is clearly visible.