An "analytic continuation" of power series coefficients Cauchy residue theorem tells us that for a function
$$f(z) = \sum_{k \in \mathbb{Z}} a(k) z^k,$$
the coefficient $a(k)$ can be extracted by an integral formula
$$a(k) = \frac{1}{2\pi i}\oint f(z) z^{-k-1},$$
with a contour around zero. Now, there is nothing to prevent us from thinking of $a\colon \mathbb{Z} \to \mathbb{C}$ as a function over the complex domain, defined by the above integral. In this way, we have naturally "continued" the function $a(k)$ over the integers to one over the complex numbers. Is there any meaning to this beyond just something amusing?
For example, consider $f(z)=\exp(z)$. In this case, $a(k)=1/k!$ which one would think would be continued to $1/\Gamma(1+k)$. But at least if you attempt to plot what happens above, this is not the case. We get a function with Bessel-like oscillations (blue) which nevertheless coincides with $1/\Gamma(1+k)$ (orange) on integral points! 

One can experiment with other choices of $f(z)$. When $f(z)=(z+1)^n$, it seems like the continuation actually coincides with what one would naturally think (the Beta function). Which is actually quite strange. For example for $n=5$ and $f(z) = 1 + 5 z + 10 z^2 + 10 z^3 + 5 z^4 + z^5$, how does this computation "learn" that the sequence $1,5,10,10,5,1$ actually corresponds to the binomial coefficients and thus should be continued to the Beta function? If we perturb one of the coefficients, say change $5 z^4$ to $2 z^4$, the result would look like a "perturbed" Beta function (blue below vs the Beta function in orange).

Any explanation (or references in the literature) for what is going on here?
 A: While reordering the database of my library, I found a reference to paper [1], which I think is somewhat pertaining to the question. The Authors, in order to prove a result on the localization of singularities of lacunary power series on the boundary of their circle of convergence, establish a necessary and sufficient condition for the analytic continuation across boundary arcs, which involves exactly the possibility of "continuing" the function $a:\Bbb N\to\Bbb C$ (in our notation) to the whole complex plane. Precisely, they prove the following result
([1], theorem 2, §2.2, pp. 564-567) If $\sigma \in [0,\pi )$, then the open arc
$$
\gamma_\sigma\triangleq \{e^{i \theta}: \alpha < \theta < 2\pi - \alpha\}
$$ is an arc for regularity of the power series $f$ (i.e. it can be analytically continued across it) if and only if there exists an entire function of exponential type $\varphi$ satisfying the following conditions

*

*$\varphi (n) = a(n)$ for $n \in\mathbb N$,

*$h_\varphi (0) = 0$ and
$$
\limsup_{\theta \to 0} \frac{h_\varphi (\theta )}{| \theta |} \leq \alpha
$$
where $h_\varphi (\theta )$ is the indicator function of $\varphi$, i.e.
$$
h_\varphi (\theta )\triangleq\limsup_{r \to \infty}\big (r^{-1} \log | \varphi (r e^{i \theta})| \big)
$$
The proof is constructive: the entire function is constructed in the form
$$
\varphi(\zeta)=\frac{1}{2\pi i}\oint\limits_{\Gamma_{\gamma_\sigma}} f(z) z^{-\zeta-1}\mathrm{d} z,
$$
where $\Gamma_{\gamma_\sigma}$ is a properly constructed simple closed path (and this possibly accounts for the nice observation of Igor Khavine).
Final note
Theorem 2 of [1] above states that the possibility of representing the coefficients of $f$ as the evaluation over the integers of an exponential entire function belonging to a particular class is equivalent to its analytic continuablity across an open arc of the circumference of its convergence disk. However, if $f$ is non continuable across the boundaries of convergence disk (i.e. its circumference is a "cut" for $f$) the above theorem is not applicable: I'd like to know if, in such case, the Fourier integral approach described by prof. Tao in his answer could be used to say something on the structure of this kind of entire functions.
Reference
[1] Norair Arakelian and Wolfgang Luh and Jürgen Müller, "On the localization of singularities of lacunary power series", (English)
Complex Variables and Elliptic Equations 52, No. 7, 561-573 (2007), MR2340942, Zbl 1123.30001.
A: Edit (1/21/21): (Start)
For your first example, a classic method of interpolation is related to fractional calculus.
$$k!\;  a(k) = k! \; \oint_{|z|=r} \frac{e^z}{z^{k+1}} \; dz = e^{-1}k! \; \oint_{{|z|=r}} \frac{e^{z+1}}{z^{k+1}} \; dz $$
$$=  e^{-1}k! \; \oint_{|z-1|=1} \frac{e^{z}}{(z-1)^{k+1}} \; dz =e^{-1} D^k_{z=1} e^z.$$
Interpolating using a standard fractional integroderivative, of which there are several reps,
$$\lambda! \; a(\lambda) = \; e^{-1} D_{z=1}^{\lambda} \; e^z = e^{-1} \; \sum_{n \ge 0} \frac{z^{n-\lambda}}{(n-\lambda)!}\; |_{z=1}$$
$$ = e^{-z} \; \sum_{n \ge 0} \frac{z^{n-\lambda}}{(n-\lambda)!}\; |_{z=1} = e^{-z} z^{-\lambda}\; E_{1,-\lambda}(z) \; |_{z=1},$$
where $E_{\alpha,\beta}(z)$ is the Mittag-Leffler function (general definition in Wikipedia, MathWorld; some history and applications), encountered very early on by anyone exploring fractional calculus.
This method of interpolation gives the entire function
$$ a(\lambda) = e^{-1} \; E_{1,-\lambda}(1) \frac{1}{\lambda!}  =  e^{-1} \; \sum_{n \ge 0} \frac{1}{(n-\lambda)!} \; \frac{1}{\lambda!}, $$
which matches the OP's first graph and extends it to any real or complex argument, giving, of course, $a(k) = 1/k!$ for $k=0,1,2, ...$.
Similarly, the fractional calculus can be applied to usefully interpolate coefficients generated by $ D_z^n \; f(z)$ to those of $D_z^{\lambda} f(z)$, with care taken with defining branch cuts and contours of integration for the complex function for the point of evaluation. The fundamental example is Euler's integral rep of the beta function as depicted below. Another is the interpolation of the associated Laguerre polynomials to the confluent hypergeometric functions and, therefore, interpolations of the associated coefficients. (See links below.)
Other productive methods of interpolation--all can be related to fractional calculus--are sketched below.
(End)
Edit 1/23/21: (Start)
The interpolation you desire is obtained by applying Ramanujan's favorite Master Formula, a.k.a. Mellin transform interpolation, as discussed below and illustrated in four other MO-Qs--Q1, Q2, Q3, and Q4. This, naturally, amounts to replacing $k$ by $s$ only outside your Cauchy integral rather than within for the reason Terry Tao notes, which is equivalent to noting the Cauchy contour integration is not equal to the Mellin transform integration.
(End)
The Mellin transform pair allows for interpolation of the coefficients of generating functions, often directly connected to sinc and/or Newton interpolation. Basically, the following is a sketch of the analytics of Ramanujan's Master Formula/Theorem, which he so profoundly and elegantly wielded.
Here $(a.)^n := a_n$ in the Taylor series $e^{a.x}$ is interpolated as $a_{-s}$ via the Mellin transform of the Taylor series $f(x) = e^{-a.x}$.
First consider the normalized Mellin transform and its inverse
$$F(s) = MT[f(x)] = \int_{0}^{\infty} f(x) \; \frac{x^{s-1}}{(s-1)!} \; dx$$
$$f(x) = MT^{-1}[F(s)] = \frac{1}{2 \pi i} \int_{\sigma - i \infty}^{\sigma + i \infty} \frac{\pi}{\sin(\pi s)} F(s) \frac{x^{-s}}{(-s)!} \; ds .$$
Then the RMT holds for a class of functions such that the simple poles of the inverse sine in the inverse Mellin transform give
$$f(x) = e^{-a.x} = \sum_{n \geq 0} \frac{(-a.x)^n}{n!} = \sum_{n=0} a_n \frac{(-x)^n}{n!} =  
 \sum_{n=0} F(-n) \frac{(-x)^n}{n!} \;  ,$$
that is, such that we may close the complex contour to the left (e.g., in the sense of the limit of a semicircle with its radius expanding to infinity) for $0 < \sigma < 1$ and $0 < x < 1$ when $F(s)$ has no singularities/poles within the contour. This rep allows an extension of the RMT (and the Mellin transform) to cases in which poles are present in $F(s)$ and to other ranges of $x$.
Also note (see, e.g., Gelfand and Shilov's "Generalized Functions") the relation
$$D_x^{m+n+1} \; H(x) \frac{x^m}{m!} = H(x) \frac{x^{-n-1}}{(-n-1)!} = \delta^{(n)}(x),$$
reflected in the two (of several) reps of the fractional differintegro op equivalent under analytic continuation
$$\frac{x^{\alpha-\beta}}{(\alpha-\beta)!} = \frac{d^{\beta}}{dx^\beta}\frac{x^{\alpha}}{\alpha!}=\int_{0}^{x}\frac{z^{\alpha}}{\alpha!}\frac{(x-z)^{-\beta-1}}{(-\beta-1)!} dz  = \frac{1}{2\pi i} \oint_{|z-x|=|x|}\frac{z^{\alpha}}{\alpha!}\frac{\beta!}{(z-x)^{\beta+1}}dz ,$$
with $H(x)$ the Heaviside step function and $\delta(x)$, the Diraac delta. (This is equivalent through binomial expansion to a sinc function/cardinal series interpolation of the binomial coefficients $\binom{s}{\alpha} = \sum_{n \geq 0} \binom{s}{n} \frac{\sin(\pi(\alpha-n))}{\pi(\alpha-n)}=\sum_{n \geq 0} \binom{s}{n} \binom{0}{\alpha-n} $, an instance also of the Chu-Vandermonde identity. The contour integral for this beta function integral is easily transformed into a bandlimited Fourier transform. See, e.g., the relation to Newton interpolation in the MSE_Q "Why is the Euler gamma function the best extension of the factorial function to the reals?" For more detail and connections, see my post "Fractional calculus and interpolation of generalized binomial coefficients." The integral reps of the confluent (see this MO-Q) and non-degenerate hypergeometric functions can be expressed effectually as this or easily related differintegral ops acting on simple functions. Years ago I found the Danish masters Niels Nielsen and Niels Norlund to be most informative on contour integrals in interpolation. There is also a connection to Riemann surfaces via Pochhammer's contour for the beta function integral.)
So, under the conditions above,
$$F(-n) = \int_{0}^{\infty} f(x) \; \frac{x^{-n-1}}{(-n-1)!} \; dx = \int_{0}^{\infty} e^{-a. x} \; \delta^{(n)}(x) \; dx = a_n,$$
and this suggests the analytic continuation and relation to umbral calculus
$$F(s) = \int_{0}^{\infty} f(x) \; \frac{x^{s-1}}{(s-1)!} \; dx = \int_{0}^{\infty} e^{-a.x} \; \frac{x^{s-1}}{(s-1)!} \; dx = (a.)^{-s} = a_{-s}.$$
The iconic guiding example is the  Euler gamma function integral rep with $(a.)^n = a_n = c^n$
$$ (a.)^{-s} = a_{-s}  = c^{-s} = F(s) = MT[f(x)= e^{-c\; x}] = \int_{0}^{\infty} e^{-c \; x} \; \frac{x^{s-1}}{(s-1)!} \; dx = \frac{1}{c^{s}},$$
giving the interpolation of the coefficients of the Taylor series of $e^{cx}$, i.e., $a_n = c^n$, as $a_{-s}=c^{-s}=F(s)$.
Another useful example, which vividly illustrates the relation to the Appell Sheffer sequences of umbral calculus (of which the $x^n$ with e.g.f. $e^{x}$ is the basic example), is the integral rep for (what I call) the Bernoulli function, simply related to the Hurwitz zeta function and generalizing the Bernoulli polynomials,
$$ B_{-s}(z) = (B.(z))^{-s} = \int_{0}^{\infty} e^{-B.(z)t} \; \frac{t^{s-1}}{(s-1)!} \; dt $$
$$ = \int_{0}^{\infty} \frac{-t}{e^{-t}-1} \; e^{-zt} \frac{t^{s-1}}{(s-1)!} \; dt = s \; \zeta(s+1,z)$$
where the e.g.f. for the Bernoulli polynomials with $(b.)^n = b_n$ the Bernoulli numbers is
$$e^{B.(x)t} = e^{(b.+x)t} = e^{b.t} e^{xt} = \frac{t}{e^t-1} \; e^{xt}.$$
Note that
$$B_n(z) = -n \; \zeta(1-n,z),$$
$$B_n(1) = -n \; \zeta(1-n,1) =-n \; \zeta(1-n) (Riemann) =  (-1)^n B_n(0) = (-1)^n b_n.$$
Through this characterization, it is not too difficult to show that the Bernoulli function inherits all the elegant properties of a regular Appell sequence, such as $D_z \; B_{s}(z) =  s \; B_{s-1}(z)$.

Hankel contour deformations, Hadamard finite part regularization of the Mellin transform guided by the inverse Mellin transform, the Mellin-Barnes contour integral, and other methods of analytic continuation can be used to extend the range of $s$ and other parameters, as have been done for the integral reps of the Euler gamma and beta functions and Riemann and other zeta functions and their generalizations. In addition, interchanging summation and integration often gives rise to useful asymptotic expansions of functions a la Borel, Heaviside, Hardy, and Poincare.
Riemann knew all this stuff. Ramanujan intuited it. Hardy formalized it. (I stumbled across it on a journey starting from the ladder ops of QM and a brief comment by my old math prof Stallybrass about the sequence $D^{m+n} H(x) \frac{x^m}{m!}$ in his integral transforms class an eon ago.)
For application to defining fractional powers of operators, see my answer and comments therein to the MO-Q "What does the inverse Mellin transform really mean?" and several of my blog posts, such as "The Creation / Raising Operators for Appell Sequences."
A: The usual way to interpret the integral $$\frac{1}{2\pi i} \oint \exp(z) z^{-s} \, \mathrm{d}z$$ for nonintegral $s$ is to take a Hankel contour, and in this case the value is indeed $\frac{1}{\Gamma(s)}$ (see 5.9.2). Without knowing how you produced these plots, I'm not sure what else to say about your first example.
This is also related to Ramanujan's master theorem, which expresses this analytic continuation as a Mellin transform. Hardy gave some conditions under which this identity can be made rigorous.
A: It is easy to see that for every convergent power series
$$f(z)=\sum_{n=0}^\infty a_nz^n$$
there exists an entire function of exponential type $F$ which interpolates
the coefficients: $F(n)=a_n$. This function is not unique.
But existence of such a function with some additional properties is related to analytic continuation properties of the power series.
Interpolating coefficients of a power series by an entire function of $n$ is a common powerful method of the study of singularities of power series. There is a close relation between the growth of this entire function and analytic continuation of the power series. See, for example,
L. Bieberbach, Analytische Fortsetzung, Springer 1955,
or V. Bernstein, Lecons sur les progres recents de la theorie de series de Dirichlet, Paris, 1933.
Sorry, I do not know any books in English, but you can look at this paper
arXiv:0709.2360, Theorem C on p. 7 for a typical result.
Another nice result is the theorem of Leau:
A convergent power series
$$\sum_0^\infty a_nz^n$$
has an analytic continuation to $\overline{C}\backslash \{1\}$ if and only if
$a_n=f(n)$ where $f$ is an entire function of minimal exponential type.
And there are many theorems in between.
A: As pointed out in comments, your extension of $a(k)$ depends on the choice of branch used for $z^{-k-1}$.  All such extensions, though, are precisely the band-limited extensions of $a(k)$ given by the Nyquist-Shannon sampling theorem.  For instance, if one uses the standard branch cut on the negative real axis, then the extension
$$ a(x) = \frac{1}{2\pi} \int_{-\pi}^{\pi} f(e^{i\theta}) e^{-ix\theta}\ d\theta$$
is the unique extension of
$$ a(k) = \frac{1}{2\pi} \int_{-\pi}^{\pi} f(e^{i\theta}) e^{-ik\theta}\ d\theta$$
from the integers to the reals which is band-limited to the spectral interval $[-1/2,1/2]$ (or $[-\pi,\pi]$, depending on one's choices of convention for Fourier transform).  The extension can be given explicitly by the Whittaker-Shannon interpolation formula
$$ a(x) = \sum_{k \in {\bf Z}} a(k) \mathrm{sinc}(\pi(x-k))$$
where $\mathrm{sinc}(x) = \frac{\sin(x)}{x}$ is the sinc function (with convention $\mathrm{sinc}(0)=1$).  Other choices of branch cut lead to slightly different formulae of this type.
The reason why this procedure gives the "wrong" extension for the reciprocal Gamma function $1/\Gamma(k+1)$ is simply that this function is not band-limited.  From the Hankel contour formula mentioned in another answer one instead has
$$ \frac{1}{\Gamma(k+1)} = \frac{1}{2\pi i} \oint e^z z^{-k-1}\ dz + \frac{\sin(\pi k)}{\pi} \int_1^\infty e^{-t} t^{-k-1}\ dt$$
(I have not checked for sign errors etc.); the second term is not band-limited and is causing the discrepancy between your band-limited interpolation and the interpolation you were hoping for.
