Why does the Gamma-function complete the Riemann Zeta function? Defining   $$\xi(s) := \pi^{-s/2}\ \Gamma\left(\frac{s}{2}\right)\ \zeta(s)$$  yields      $\xi(s) = \xi(1 - s)$ (where $\zeta$ is the Riemann Zeta function).
Is there any conceptual explanation - or intuition, even if it cannot be made into a proof - for this? Why of all functions does one have to put the Gamma-function there?
Whoever did this first probably had some reason to try out the Gamma-function. What was it?
(Best case scenario) Is there some uniform way of producing a factor out of a norm on the rationals which yields the other factors for the p-adic norms and the Gamma factor for the absolute value?
 A: One way to get started is to look at the integral for the gamma function:
$$\Gamma(s) = \int_0^\infty t^{s-1} e^{-t}\,dt$$
Subsitute $t=nx$ in the integral to arrive at
$$\frac{\Gamma(s)}{n^s} = \int_0^\infty e^{-nx}x^{s-1}\,dx$$
which we then sum up to get
$$\Gamma(s)\zeta(s)=\int_0^\infty \frac{x^{s-1}}{e^x-1}\,dx$$
which already shows that there is some connection between the gamma and zeta functions, and it does in fact allow us to extend the definition of the zeta function into the critical strip.
What comes next is far less obvious, but the idea is to introduce a branch cut for $x^{s-1}$ along the positive real axis, and to replace the above integral by one running from $+\infty$ along the bottom of the positive real axis, around the origin, and back to $+\infty$ along the top of the real axis. This introduces an extra factor $1-e^{2\pi i s}$. Now start expanding the circle around the origin, taking account of the poles of the integrand along the imaginary axis as we go, and end up with
$$\Gamma(s)\zeta(s)=(2\pi)^{s-1}\Gamma(1-s)\sin(\tfrac12\pi s)\zeta(1-s).$$
From there, some cleanup still remains. As I said, this is not terribly intuitive, so it doesn't answer your question, but the first paragraph should at least give you a notion how the gamma and zeta functions are interrelated.
A: To the best of my understanding, the answer is yes, and this uniform way consists of some integration over the local field.  This is explained in John Tate's dissertation.  One starts with a certain smooth rapidly decreasing function, for which one takes the characteristic function of the p-adic integers in the nonarchimedean case and the function $e^{-|x|^2}$ for an archimedean field.  This is being multiplied with $|x|^s$ (approximately) and integrated over the Haar measure of the additive group of the field.  This produces the $\Gamma$-factor for an archimedean field and $(1-p^{-s})^{-1}$ for a p-adic field.
A: As has been explained above, the zeta function has a factor for each completion of $\mathbb{Q}$. The factor at $\mathbb{R}$ has to do with integrating $e^{- \pi x^2}$ and the factor at $\mathbb{Q}_p$ has to do with integrating the characteristic function of $\mathbb{Z}_p$.
Some people might wonder why these two functions were chosen. The answer is simple: they are both their own Fourier transforms.
Also, I don't think anyone has recommended Terry Tao's expository post on this material yet. It is quite good.
A: I'm not sure of the history of the gamma factor, though I would suggest that no one "tried it out", but rather it simply arose in trying to prove of the functional equation. Riemann was the first to prove the functional equation, and his proof essentially follows that in Harald Hanche-Olsen's answer, which makes my explanation plausible. Alternatively, the functional equation of the zeta function comes out of the functional equation of a theta series, and the Mellin transform of a theta series gives rise to a Gamma function. This latter explanation arises more naturally for modular forms: the L-function of a modular form is also completed by a gamma factor to obtain a functional equation; in this case, the completed L-function is simply the Mellin transform of the modular form itself.
Furthermore, as Leonid Positselski answers, it is indeed true that Tate's thesis provides a uniform way of obtain the gamma factors at infinity in the same manner as one obtains the local L-factors at finite places.
More generally, there is a recipe given an arbitrary motive for the expected gamma factors that should give a functional equation for the motivic L-functions. These are due to Deligne and Serre (I believe) and are determined by the Hodge structure of the motive (see Deligne's corvallis article "Valeurs de fonctions L..."). This shows that there's a uniform way of obtaining the gamma factors as one varies the L-function one is studying, an orthogonal question to the one Leonid Positselski answered.
A: Gamma function arises when we consecutively differentiate an Appell sequence. An example of Appell polynomials are Bernoulli polynomials. When we differentiate it, the factors combine with themselves:
$$B_n'(x)=nB_{n-1}(x)$$
$$B_n''(x)=n(n-1)B_{n-2}(x)$$
$$B_n'''(x)=n(n-1)(n-2)B_{n-3}(x)$$
They are just another name for Hurwitz Zeta function:
$$B_n(x) = -n \zeta(1-n,x)$$
Thus, for $f(s,q)=\zeta(s,-q)$
$$\frac\partial{\partial q}f(s,q)= s f(s+1,q)$$
$$\frac{\partial^2}{\partial q^2}f(s,q)= s(s+1) f(s+2,q)$$
$$\frac{\partial^3}{\partial q^3}f(s,q)= s(s+1)(s+2) f(s+3,q)$$
Since Reihmann zeta is Hurwitz zeta evaluated at $q=1$, the expression you give is apparently consecutive derivative of Hurwitz Zeta, with factor $\pi^{-s}$ appearing if we normalize Hurwitz Zeta by stretching it horizontally by factor of pi.
Consecutive derivatives of Hurwitz Zeta in turn are nothing more than just polygamma function.

For instance, here is the function $-1/x$:


If we add infinitely many similar functions with a shift of pi/2 each in both directions, we get $\tan x$. But if we do the same only in one direction, we get "incomplete tangent":

The yellow one is $\operatorname{pg}(x)=\frac 1\pi \psi (\frac x\pi)$, the blue one is $\operatorname{cpg}(x)=-\frac 1\pi \psi (1-\frac x\pi)$. They obey $\operatorname{cpg}(x)+\operatorname{pg}(x)=-\cot(x)$.
Now if we differentiate cpg(x) we get:
$$(\operatorname{cpg}(x))^{(s-1)}=\pi^{-s}\Gamma(s)\zeta(s,1-\frac x\pi)$$
Compare it with yours formula:
$$\xi(2s) = \pi^{-s}\Gamma\left(s\right)\zeta(2s)$$
A: There may also be some interest in the point of the "local functional equation", namely, that in fact the Gamma function (with the power of $\pi$) is just one (optimized) possibility, and somehow making a suboptimal choice doesn't really matter:
For a Schwartz function $f$ on $\mathbb R$, let $\Gamma(f,s)=\int_{\mathbb R^\times} |x|^s\,f(x)\;{dx\over |x|}$. The usual Gamma factor is obtained by taking a Gaussian. The local functional equation (proven by changing variables in the defining integrals, in the range $0<{\rm Re}(s)<1$, is
$$\Gamma(f,s)\cdot \Gamma(\hat{g},1-s)\;=\; \Gamma(\hat{f},1-s)\cdot \Gamma(g,s)$$
for any two Schwartz functions $f,g$. And Riemann's argument proves
$$
\Gamma(f,s)\cdot \zeta(s) \;=\; \Gamma(\hat{f},1-s)\cdot \zeta(1-s)
$$
for any Schwartz $f$.
A: Multiple answers and comments have already pointed out that the conceptual role of $\pi^{-s/2}\Gamma(s/2)$ comes from the viewpoint of Iwasawa and Tate, which for $\text{Re}(s) > 1$ creates this function as $\int_{\mathbf R^\times} e^{-\pi x^2}|x|^s\,dx/|x|$, an integral over the multiplicative group $\mathbf R^\times$ of the function $e^{-\pi x^2}$ that is self-dual for the Fourier transform on the additive group $\mathbf R$ relative to the self-duality $\langle x,y\rangle = e^{2\pi ixy}$ or $\langle x,y\rangle = e^{-2\pi ixy}$ on $\mathbf R$. (If we use another self-duality of $\mathbf R$ then $e^{-ax^2}$ would be self-dual for some $a \not= \pi$ instead.)
It's also been said elsewhere on this page that there are many self-dual Schwartz functions on $\mathbf R$, or more specifically many even self-dual Schwartz functions on $\mathbf R$: for Schwartz $f$ on $\mathbf R$ and $\text{Re}(s) > 0$, we have $\int_{\mathbf R^\times} f(x)|x|^s\,dx/|x| = \int_{0}^\infty (f(x) + f(-x))x^s\,dx/x$ and this is $0$ when $f$ is odd, so we may as well assume $f$ is even since $f(x) + f(-x)$ is even anyway and we want to avoid the silly equation $0=0$ even if it is a valid equation.
For arbitrary Schwartz $f$ on $\mathbf R$, set $\Gamma_f(s) = \int_{0}^\infty f(x)x^s\,dx/x$, which is a mild modification of the function $\Gamma(f,s)$ in Paul Garrett's answer (his $\Gamma(f,s)$ is my $\Gamma_{f(x)+f(-x)}(s)$ by a formula I wrote in the previous paragraph). This function converges absolutely and is analytic for $\text{Re}(s) > 0$, and it extends meromorphically to $\mathbf C$ by repeated integration by parts (the same way the $\Gamma$-function can be extended to $\mathbf C$ from its integral definition for $\text{Re}(s) > 0$), and Tate's thesis shows there is a general functional equation $\Gamma_f(s)\zeta(s) = \Gamma_{\hat{f}}(1-s)\zeta(1-s)$ where $\hat{f}$ is the Fourier transform of $f$ (for the self-duality on $\mathbf R$ given by $\langle x,y\rangle = e^{-2\pi ixy}$), so if $f$ is self-dual then we get $$\Gamma_f(s)\zeta(s) = \Gamma_{f}(1-s)\zeta(1-s),$$ a very nice functional equation indeed, especially if we use even $f$ to avoid $0 = 0$. 
All of what I wrote so far has appeared explicitly or implicitly in some of the other comments or answers. Since there are many self-dual even Schwartz functions $f$ on $\mathbf R$, what is it about the choice $f(x) = e^{-\pi x^2}$, leading to $\Gamma_f(s) = (1/2)\pi^{-s/2}\Gamma(s/2)$ (an extra $1/2$ on both sides of the functional equation can be cancelled) that is so nice? I have not seen the following property pointed out yet: with this choice of $f$ and familiarity with the $\Gamma$-function we know $\Gamma_f(s) \not= 0$ for $\text{Re}(s) > 1$ (in fact for $\text{Re}(s) > 0$), so therefore $\Gamma_f(s)\zeta(s) \not= 0$ for $\text{Re}(s) > 1$ from $\zeta(s)$ being nonvanishing there, and then by the functional equation $\Gamma_f(s)\zeta(s) \not= 0$ for $\text{Re}(s) < 0$, which means all zeros of $\Gamma_f(s)\zeta(s)$ have $0 \leq \text{Re}(s) \leq 1$.  If you want to use a totally random even Schwartz function for $f$ in order to define a factor $\Gamma_f(s)$ that completes the Riemann zeta-function, you will get the nice-looking nontrivial functional equation displayed above but how are you going to use $\Gamma_f(s)\zeta(s)$ to analyze the location of zeros of $\zeta(s)$ (including discovering its  trivial zeros, whether or not you consider those important) if you do not know where $\Gamma_f(s)$ has its zeros and poles? 
So although there are many even Schwartz functions $f$ on $\mathbf R$ besides $e^{-\pi x^2}$ that you could use to get a nice functional equation by multiplying $\zeta(s)$ by $\Gamma_f(s)$, the reason that the choice $f(x) = e^{-\pi x^2}$ is so convenient is that we actually know the zeros and poles of $\Gamma_f(s) = (1/2)\pi^{-s/2}\Gamma(s/2)$: it has no zeros in $\mathbf C$ and it has simple poles at $0, -2, -4, \ldots$. For even self-dual Schwartz $f$ on $\mathbf R$ that are not simple modifications of $e^{-\pi x^2}$, how feasible is it to determine whether or not $\Gamma_f(s) \not= 0$ for $\text{Re}(s) > 1$ (or $\text{Re}(s) > 0$)?  The method of meromorphically continuing $\Gamma_f(s)$ from the half-plane $\text{Re}(s) > 0$ where it is analytic to all of $\mathbf C$ shows that its only possible poles are at $0, -1, -2, -3, \ldots$ with orders at most $1$ and the residue at $s = -n$ is $(-1/n!)\int_0^\infty f^{(n+1)}(x)\,dx$, which by the Fundamental Theorem of Calculus is $(-1/n!)(f^{(n)}(\infty) - f^{(n)}(0)) = f^{(n)}(0)/n!$.
Therefore you could determine the poles of $\Gamma_f$ by seeing when $f^{(n)}(0)$ is 0 and not 0, but how are you going to determine where the zeros of $\Gamma_f$ are or that there are no zeros? (EDIT: for even $f$, its odd-order derivatives vanish at $0$, so the residue at $-n$ vanishes when $n$ is odd, which means the poles of $\Gamma_f(s)$ can only be at $n = 0, -2, -4, -6, \ldots$. Those are all simple poles of $\pi^{-s/2}\Gamma(s/2)$, which has no zeros, so $G(s) := \Gamma_f(s)/(\pi^{-s/2}\Gamma(s/2))$ is an entire function. Thus $\Gamma_f(s) = G(s)\pi^{-s/2}\Gamma(s/2)$ with $G$ entire, so $\pi^{-s/2}\Gamma(s/2)$ a "holomorphic gcd" of all $\Gamma_f(s)$ for even Schwartz functions $f$ on $\mathbf R$. The exponential factor $\pi^{-s/2}$ was kind of irrelevant to drag through the calculation since it has no zeros or poles, but it's traditionally seen alongside $\Gamma(s/2)$ so I used it. This addresses comments below by Will Sawin and Venkataramana.)
Example: the function $f(x) = 1/(e^{\pi x} + e^{-\pi x})$ is an even self-dual Schwartz function on $\mathbf R$. Can someone determine in a self-contained way (i.e., not using $\zeta(s)$) where $\Gamma_f(s)$ has its zeros on $\mathbf C$, or determine if it has no zeros? 
Edit: Ignoring the wacky example just above, in some comments below I work out an example with $f(x)$ being a 4th degree Hermite polynomial times a Gaussian and find that $\Gamma_f(s)$ has two zeros with positive real part, at $s = (1\pm \sqrt{-2})/2$.
A: "Why of all functions does one have to put the Gamma-function there?"
$\zeta(s)$ has trivial zeroes at $-2, -4, -6$, etc.
$\zeta(1-s)$ thus has trivial zeroes at $s=3, 5, 7$, etc - a completely different set of zeroes.
To make a reflection formula where $\zeta(s)$ is somehow equal $\zeta(1-s)$, you have to get rid of the two differing sets of trivial zeroes. Multiplying by the gamma is perfect for this since its poles will cancel out those zeroes. For example, $\Gamma(s/2)$ has poles at $0, 2, 4, 6$, etc. and should go with $\zeta(s)$. $\Gamma((1-s)/2)$ has poles at $s=1, 3, 5$, etc. and should go with $\zeta(1-s)$.
It's possible to prove that gamma is the right choice, but Euler no doubt discovered that gamma is the right function through numerical experimentation - when he discovered the zeta reflection formula like 250 years ago.
A: 
Whoever did this first probably had
  some reason to try out the
  Gamma-function. What was it?

The first one to do this was, precisely, Riemann in his famous (and 150 years old) paper: Ueber die Anzahl der Primzahlen unter einer gegebenen Grösse.  There he proved the functional equation as well, with the method that Harald explained above.
A: Although there is already an answer of mine, I want to add another answer.
This is TL;DR.
Short answer. This is because logarithmic function lacks factorial in its Taylor expansion.
Medium answer. Riemann's functional equation links exponential and trigonometric functions with logarithms and inverse trigonometric. It contains everything what you need to make an exponent from a logarithm.
Long answer.
This is Taylor series for logarithm:
$$\ln(z+1)=z-\frac{z^2}{2}+\frac{z^3}{3}-\frac{z^4}{4}+\frac{z^5}{5}-\frac{z^6}{6}+\frac{z^7}{7}-\frac{z^8}{8}+\frac{z^9}{9}-\frac{z^{10}}{10}+O\left(z^{11}\right)$$
This is Taylor series for exponent:
$$\exp (z)-1=z+\frac{z^2}{2!}+\frac{z^3}{3!}+\frac{z^4}{4!}+\frac{z^5}{5!}+\frac{z^6}{6!}+\frac{z^7}{7!}+\frac{z^8}{8!}+\frac{z^9}{9!}+\frac{z^{10}}{10!}+O\left(z^{11}\right)$$
What should we add to the former to get the later? Well, we have to add the factorial and remove the counter from the denominator.
Consider such algebraic element $\omega_+$ (not a real number) on which a function "standard part" is implemented in such a way, that $\operatorname{st} \omega_+^n=B_n^*$ where $B_n^*$ are Bernoulli numbers (with $B_1^*=1/2$), or more generally, $\operatorname{st}\omega_+^x=-x\zeta(1-x)$.
Now consider the function
$$\frac{z}{2\pi} \log \left(\frac{\omega _+-\frac{z}{2 \pi }}{\omega _++\frac{z}{2 \pi }}\right)$$
Its Taylor series is 
$$-\frac{z^2}{2 \left(\pi ^2 \omega _+\right)}-\frac{z^4}{24 \left(\pi ^4 \omega _+^3\right)}-\frac{z^6}{160 \left(\pi ^6 \omega _+^5\right)}-\frac{z^8}{896 \left(\pi ^8 \omega _+^7\right)}-\frac{z^{10}}{4608 \left(\pi ^{10} \omega _+^9\right)}+O\left(z^{11}\right)$$
Following Riemann's functional equation and our definition, we have:
$$\operatorname{st}\omega_+^{-x}=\operatorname{st}\frac{-\omega_+^{x+1} 2^x\pi^{x+1}}{\sin(\pi x/2)\Gamma(x)(x+1)}$$
So we can substitute the negative powers of $\omega_+$ with positive powers without changing the standard part of the whole expression.
The non-zero terms are 
$$\frac{2 \left(-\frac{1}{2 \pi }\right)^n \left(-\omega _+\right){}^{1-n}}{n-1}$$
and after substitution we have
$$\frac{\omega _+^n \sec \left(\frac{\pi  n}{2}\right)}{\Gamma (n+1)}$$
The resulting series is 
$$\frac{1}{2} \omega _+^2 z^2+\frac{1}{24} \omega _+^4 z^4-\frac{1}{720} \omega _+^6 z^6+\frac{\omega _+^8 z^8}{40320}-\frac{\omega _+^{10} z^{10}}{3628800}+O\left(z^{11}\right)$$
oh, wait... is not it similar to
$$\cos \left(\omega _+ z\right)=1-\frac{1}{2} \omega _+^2 z^2+\frac{1}{24} \omega _+^4 z^4-\frac{1}{720} \omega _+^6 z^6+\frac{\omega _+^8 z^8}{40320}-\frac{\omega _+^{10} z^{10}}{3628800}+O\left(z^{11}\right)$$
Well, we got:
$$\operatorname{st}\frac{z}{2 \pi } \log \left(\frac{\omega _+-\frac{z}{2 \pi }}{\omega _++\frac{z}{2 \pi }}\right)=\operatorname{st}(\cos \left(\omega _+ z\right)-1)$$
In a similar way one can establish other impressive relations:
$$\operatorname{st}(\exp \left(\omega _+ z\right)-\omega _+ z-1)=\operatorname{st}\frac{i z}{2 \pi } \log \left(\frac{\omega _+-\frac{i z}{2 \pi }}{\omega _++\frac{i z}{2 \pi }}\right)$$
$$\operatorname{st}\cos \left(\omega _+ z\right)=\operatorname{st}\frac{ z}{2 \pi } \log \left(\frac{\omega _+-\frac{ z}{2 \pi }}{\omega _-+\frac{ z}{2 \pi }}\right)$$
$$\operatorname{st}\cosh \left(\omega _+ z\right)=\operatorname{st}\frac{i z}{2 \pi } \log \left(\frac{\omega _+-\frac{i z}{2 \pi }}{\omega _-+\frac{i z}{2 \pi }}\right)$$
(where $\omega_-=\omega_+-1$).
In other words, Riemann's functional equation is a direct bridge that connects exponential function to logarithm, trigonometric functions to inverse trigonometric, transforming each term of the series separately.
