Las continuaciones analíticas de Riemann y la derivación de las ecuaciones funcionales para $\zeta$ y $\xi$ parecen bastante naturales e intuitivos desde la perspectiva del análisis complejo básico.
Riemann en la segunda ecuación de su trabajo clásico Sobre el número de números primos menores que una cantidad dada (1859) escribe la transformada de Laplace (1749-1827)
$$\int_{0}^{+\infty}e^{-nx}x^{s-1}dx=\frac{(s-1)!}{n^s},$$ válido para $real(s)>0.$
Con $n=1$ esta es la representación integral icónica de Euler (1707-1783) de la función gamma, y observando que
$(s-1)!=\frac{\pi}{sin(\pi s)}\frac{1}{(-s)!}$ de la relación simétrica $\frac{sin(\pi s)}{\pi s}=\frac{1}{s!(-s)!},$
esto se puede reescribir como
$$\frac{sin(\pi s)}{\pi}\int_{0}^{\infty}e^{-x}x^{s-1}dx=\frac{1}{(-s)!},$$
sugiriendo con toda naturalidad a alguien tan familiarizado con la continuación analítica como Riemann que
$$\frac{-1}{2\pi i}\int_{+\infty}^{+\infty}e^{-x}(-x)^{s-1}dx=\frac{1}{(-s)!},$$
válido para todos $s$ , donde la integral de línea se amplía alrededor del eje real positivo en el plano complejo para intercalarla con un corte de rama para $x>0$ y hacer un bucle en el origen en sentido positivo desde el infinito positivo hasta el infinito positivo. Al desinflar el contorno hacia el eje real se introduce un $-exp(i\pi s)+exp(-i\pi s)=-2isin(\pi s)$ . (Este contorno especial se llama ahora el Contorno de Hankel después de Hermann Hankel (1839-1873), que se convirtió en alumno de Riemann en 1860 y publicó esta integral para el fct gamma recíproco. en su habilitación de 1863. Lo más probable es que Riemann le introdujera en esta maniobra).
Riemann en su tercera ecuación observa que
$$(s-1)!\zeta(s)=(s-1)!\sum_{n=1}^{\infty }\frac{1}{n^s}=\int_{0}^{+\infty}\sum_{n=1}^{\infty }e^{-nx}x^{s-1}dx=\int_{0}^{+\infty}\frac{1}{e^x-1}x^{s-1}dx$$
y luego escribe inmediatamente como su cuarta igualdad la continuación analítica
$$2sin(\pi s)(s-1)!\zeta(s)=i\int_{+\infty}^{+\infty}\frac{(-x)^{s-1}}{e^x-1}dx,$$
válido para todos $s$ que se puede reescribir como
$$\frac{\zeta(s)}{(-s)!}=\frac{-1}{2\pi i}\int_{+\infty}^{+\infty}\frac{(-x)^{s-1}}{e^x-1}dx$$
como sugiere naturalmente la continuación analítica del recíproco de gamma anterior.
Para $m=0,1,2, ...,$ esto da
$$\zeta(-m)=\frac{(-1)^{m}}{2\pi i}\oint_{|z|=1}\frac{m!}{z^{m+1}}\frac{1}{e^z-1}dz=\frac{(-1)^{m}}{m+1}\frac{1}{2\pi i}\oint_{|z|=1}\frac{(m+1)!}{z^{m+2}}\frac{z}{e^z-1}dz$$
a partir de la cual se puede ver, si se conoce la fct. generadora de exponentes (f.g.) para los números de Bernoulli, que la integral desaparece para los números pares $m$ . Euler publicó el e.g.f. en 1740 ( MSE-Q144436 ), y Riemann ciertamente estaba familiarizado con estos números y afirma que su cuarta igualdad implica la desaparición de $\zeta(s)$ para $m$ incluso (pero no da ninguna prueba explícita). Ciertamente también conocía la ecuación funcional heurística de Euler para valores enteros de la fct. zeta, y Edwards en Función Zeta de Riemann (pg. 12, ed. de Dover) incluso especula que ".. bien puede haber sido este problema de derivar (2) [la fórmula de Euler para $\zeta(2n)$ para que sea positivo $n$ ] de nuevo que llevó a Riemann al descubrimiento de la ecuación funcional ...."
Riemann procede entonces a derivar la eqn. funcional para zeta a partir de su igualdad utilizando las singularidades de $\frac{1}{e^z-1}$ para obtener básicamente
$$\zeta(s)=2(2\pi)^{s-1}\Gamma(1-s)\sin(\tfrac12\pi s)\zeta(1-s),$$
y dice tres líneas más adelante esencialmente que puede expresarse simétricamente alrededor de $s=1/2$ como
$$\xi(s) = \pi^{-s/2}\ \Gamma\left(\frac{s}{2}\right)\ \zeta(s)=\xi(1-s).$$
Riemann dice entonces: "Esta propiedad de la función [ $\xi(s)=\xi(1-s)$ me indujo a introducir, en lugar de $(s-1)!$ la integral $(s/2-1)!$ en el término general de la serie $\sum \frac{1}{n^s}$ con lo que se obtiene una expresión muy conveniente para la función $\zeta(s)$ ." Y luego procede a introducir lo que Edwards llama una segunda prueba de la eqn. funcional utilizando la función theta de Jacobi.
se pregunta Edwards:
"Dado que la segunda prueba hace totalmente innecesaria la primera, cabe preguntarse por qué Riemann incluyó la primera prueba. Tal vez la primera prueba muestra el argumento por el que descubrió originalmente la ecuación funcional o tal vez exhibe algunas propiedades que fueron importantes en su comprensión de la misma."
Me pregunto si, al evolucionar sus ideas antes de escribir el documento construyó por primera vez $\xi(s)$ observando que al multiplicar $\zeta(s)$ par $\Gamma(\frac{s}{2})$ introduce un polo simple en $s=0$ reflejando así el polo de $\zeta(s)$ en $s=1$ a través de la línea $s=1/2$ y que los otros polos simples de $\Gamma(\frac{s}{2})$ son eliminados por los ceros de la recta real de la función zeta. El $\pi^{-s/2}$ puede determinarse fácilmente como una normalización por una función entera $c^s$ donde $c$ es una constante, utilizando la simetría compleja conjugada de las fct. gamma y zeta alrededor del eje real. Riemann tenía una fina intuición física y habría pensado de forma holística en términos de los ceros de una función (véase la prueba de Euler del problema de Basilea) y sus polos, cuya importancia ciertamente destacó.
Extendamos el razonamiento anterior para la función theta de Jacobi
$$\vartheta (0;ix^2)=\sum_{n=-\infty,}^{\infty }exp(-\pi n^{2}x^2).$$
Considerando una transformada de Mellin modificada como una interpolación de los coeficientes de la serie de Taylor ( MO-Q79868 ), es fácil adivinar (nótese los ceros de los coeficientes) que
$$\int^{\infty}_{0}\exp(-x^2)\frac{x^{s-1}}{(s-1)!} dx = \cos(\pi\frac{ s}{2})\frac{(-s)!}{(-\frac{s}{2})!} = \frac{1}{2}\frac{(\frac{s}{2}-1)!}{(s-1)!},$$
y, por lo tanto,
$$\int^{\infty}_{0}\exp(-\pi (n x)^2)x^{s-1} dx = \frac{1}{2}\pi^{-s/2}(\frac{s}{2}-1)! \frac{1}{n^s}.$$
A estas alturas deberías ser capaz de completar la línea de razonamiento para obtener, por $real(s)>1,$
$$\xi(s)=\int_{0^+}^{\infty }[\vartheta (0;ix^2)-1)]x^{s-1}dx=\pi^{-s/2}\ \Gamma\left(\frac{s}{2}\right)\ \zeta(s).$$
Hacer una continuación analítica como se hace para la función gamma en MSE-Q13956 para obtener, por 0<real(s)<1,
$$\xi(s)=\int_{0^+}^{\infty }[\vartheta (0;ix^2)-(1+\frac{1}{x})]x^{s-1}dx.$$
Entonces utiliza las simetrías de la transformada de Mellin y el hecho de que $\xi(s)=\xi(1-s)$ (como se explica en MSE-Q28737 ) para obtener la ecuación funcional
$$\vartheta (0;ix^2)=\frac{1}{x}\vartheta (0;\frac{i}{x^{2}}).$$
Actualización (5 de enero de 2021):
Esto quizá responda a la pregunta de Edwards y proporcione un camino sencillo hacia la identidad funcional.
Euler adquirió inicialmente una reputación con la resolución del problema de Basilea para establecer el valor de $\zeta(2)$ y además las identidades
$$\frac{2}{(2\pi)^{2n}}\:(2n-1)!\:\zeta(2n)=(-1)^{n+1}\frac{B_{2n}}{2n}=(-1)^{n}\zeta(1-2n).$$
Como se ha señalado anteriormente, Riemann incorporó la f.g.e. de los números de Bernoulli en su transformada normalizada de Mellin/Laplace para la función zeta. No es un gran salto de fe creer que Riemann (de las superficies y derivadas epónimas) comprendió la propiedad de interpolación de la transformada de Mellin. Dada esa suposición, podría haber anotado fácilmente y quizás formulado inicialmente la representación integral de Mellin a partir de
$$b_n = D^m_{z=0} \; e^{b.z} = D^m_{z=0} \; \frac{z}{e^z-1}$$
$$ = (-1)^{m} \; \frac{1}{2\pi i}\oint_{|z|=1}\frac{m!}{z^{m+1}} \; \frac{z}{e^z-1} \; dz $$
$$= (-1)^{m} \; \frac{1}{2\pi i}\oint_{|z|=1}\frac{m!}{z^{m+1}} \; e^{b.z} \; dz$$
$$ =(-1)^{m} \frac{1}{2\pi i} \; \oint_{|z|=1}\frac{m!}{z^{m+1}} \; [1 - \frac{z}{2}+ \sum_{n \geq 2} \; \cos(\frac{\pi n}{2}) (-2) \; (2\pi)^{-n} \; n! \; \zeta(n) \; \frac{z^n}{n!}] \; dz $$
$$ =(-1)^{m} \frac{1}{2\pi i} \; \oint_{|z|=1}\frac{m!}{z^{m+1}} \; [1 - \frac{z}{2}+ \sum_{n \geq 2} \; (-n) \;\zeta(1-n) \; \frac{z^n}{n!}] \; dz . $$
Las maniobras de deformación de Hankel anteriores podrían invertirse para obtener las transformadas de Mellin equivalentes a estas integrales de contorno (véase este MO_A ), pero esto no cambiaría la equivalencia de los coeficientes en las dos f.e.g. para los números de Bernoulli, por lo que a partir de la propiedad interpoladora de la transformada de Mellin (y por tanto de las integrales de Cauchy), Riemann podría haber continuado simplemente de forma analítica $n$ a $1-s$ para conjeturar la identidad del objetivo
$$ \cos(\frac{\pi n}{2}) \; 2 \; (2\pi)^{-n} \; n! \; \zeta(n)] \; |_{n \to 1-s} = n \;\zeta(1-n) \; |_{n \to 1-s},$$
dando la identidad funcional de reflexión
$$\cos(\frac{\pi (1-s)}{2}) \; 2 \; (2\pi)^{s-1} \; (1-s)! \; \zeta(1-s) = (1-s) \;\zeta(s) , $$
o
$$2 \; (2\pi)^{s-1} \; \sin(\frac{\pi s}{2}) \; (-s)! \; \zeta(1-s) = \zeta(s) . $$
Edición 1/24/21: Una cuarta forma bastante sencilla de derivar la ecuación funcional que implica series de Fourier y la distribución de suma del núcleo de Poisson, pero evitando la función theta, se da en mi respuesta a este reciente MO-Q .