Deriving the functional equation for $\zeta(s)$ from summing the powers of the zeros required to count the integers When counting the number of integers $n(x)$ below a certain non-integer number $x$, the following series could be used:
$$n(x) = x-\frac12 + \sum_{n=1}^{\infty} \left(\frac{e^{x \mu_n}} {\mu_n}+\frac{e^{x \overline{\mu_n}}} {\overline{\mu_n}}\right)$$
where $\mu_n = 2\pi n i$ which are the zeros of the function $\xi_i(s) = \frac{2}{s}\sinh\left(\frac{s}{2}\right)$ that has the simple Hadamard product:
$$\displaystyle  \xi_i(s) = \prod_{n=1}^\infty \left(1- \frac{s}{2 \pi ni} \right) \left(1- \frac{s}{{-2 \pi ni}} \right)$$
Note that $\xi_i(0)=1$ just like $\xi(0)=1$ in the Hadamard product of the non-trivial zeros of Riemann $\xi$-function when ignoring its probably superfluous factor $\frac12$.
Summing the powers of these paired zeros as follows yields ($B_r$ = Bernoulli number):
$$\hat{\sigma}_r = \sum_{n=1}^\infty \left(\frac{1}{(2\pi ni)^r}+ \frac{1}{(-2\pi ni)^r}\right) = -\frac{B_{r}}{r\,\Gamma(r)} \qquad r \in \mathbb{N}, r \gt 1\tag{1}$$
The domain of the series could be expanded as follows:
$$\hat{\sigma}_s = \frac{1}{(2\pi i)^s}\,\left(1+e^{\pi s i}\right)\sum_{n=1}^\infty \frac{1}{n^s}\qquad s \in \mathbb{C}, \Re(s) \gt 1 \tag{2}$$
$$\hat{\sigma}_s = 2^{1-s}\,\pi^{-s}\cos\left(\frac{\pi s}{2}\right)\sum_{n=1}^\infty\frac{1}{n^s}\qquad s \in \mathbb{C}, \Re(s) \gt 1 \tag{3}$$
Transferring the $\Gamma(r)$ from the RHS of (1) and $r \mapsto s$ gives:
$$2^{1-s}\,\pi^{-s}\cos\left(\frac{\pi s}{2}\right)\,\Gamma(s)\,\zeta(s) = \,\,? \tag{4}$$
which is 5/6-th of the famous functional equation. We know through various proofs (e.g. 7 different ones are listed in Titchmarsh book on the Zeta functions) that the ? $= \zeta(1-s)$ and that this provides the full analytic continuation of $\zeta(s)$ towards $s \in \mathbb{C}\,\, /\,\, {1}$.
Question: (I hope not too trivial...)
I know the Euler product reflects the multiplicative structure of the integers, whereas the functional equation reflects the additive structure, but is there an intuitive explanation on why the functional equation should emerge from summing the powers of the zeros required for the oscillating term to count the integers?
P.S.:
I read this interesting discussion, but could not derive the answer from it.
 A: The intermediary seems to be the Bernoulli number sequence which was originally birthed in summing up powers of the integers and in turn eventually gave birth, via the midwife the Mellin transform, to the Riemann and Hurwitz zeta functions. The MO-Q to which you link on motivating derivations of the functional equation for the Riemann zeta has an analytic continuation of the coefficients of the e.g.f. for the Bernoullis (the AC in fact gives the Riemann zeta function) with the numbers expressed two different ways, from which the FE of the Riemann zeta falls out. Your Eqn. 1 could be used to replace one of those reps for the Bernoullis--the one containing $\cos(\frac{\pi n}{2})$--giving the same end result, the FE. (Another perspective on the AC of the Bernoulli numbers to the Hurwitz and Riemann zeta functions is presented in this MO-Q.)
If you take the derivative of your initial equation, you obtain the Dirac delta function/operator comb on the left and a sum of cosines of the right, giving the core Poisson summation identity. The Mellin transform of the Dirac comb will give you the Riemann zeta function. For more on this, see "The Correspondence Principle" by Hughes and Ninham.
Edit 1/23-4/21:
Let me be elaborate on the last paragraph.
As you depict in your associated MSE-Q, a doubly infinite staircase function is obtained by adding $x$ to the Fourier series rep of the sawtooth wave. For $x > 0$, you can write the piecewise continuous semi-infinite staircase function as
$$H(x) \; n(x) = \sum_{n \geq 1} H(x-n) = H(x) [ \; x - \frac{1}{2} + 2 \sum_{n \geq 1} \frac{\sin(2 \pi n x)}{2 \pi n} \; ],$$
where $H(x)$ is the Heaviside step function (Heaviside knew all this).
Taking the derivative of both sides gives, for $x > 0 $, half of the core of the Poisson-summation distribution formula
$$ \sum_{n \geq 1} \delta(x-n) = H(x) [\;1 + 2 \sum_{n \ge 1} \cos(2 \pi n x) \;],$$
and, since
$$
\int_{0^{+}}^{\infty} x^{s-1} \delta(x-n) \; dx 
= n^{s-1}$$
and
$$
2 \;\int_{0^{+}}^{\infty} x^{s-1} \cos(2\pi n x) dx 
= 2 \; (2\pi n)^{-s} \int_{0}^{\infty} x^{s-1} \cos(x) \; dx$$
$$= 2\; (2\pi n)^{-s} \; (s-1)!\; \cos(\frac{\pi}{2}s)$$
for $0 < Re(s) < 1$, taking the RHS as the analytic continuation for all $s$, we have a rudimentary form of the zeta FE crystallizing.
The Mellin transform term by term of the Dirac comb gives the Riemann zeta function series rep
$$ \zeta(1-s) = \sum_{n \ge 1} \frac{1}{n^{1-s}}$$
for $Re(s) < 0$. However, the $n =0$ term, i.e., the constant term, in the cosine series poses a problem in the term by term Mellin transformation of the series. Tossing it out--regularizing through the Hadamard finite-part scheme, justified by an inverse Mellin transform rep just as for the AC of the integral for the Euler gamma function--and equating the analytically continued Mellin transforms of the two reps gives the Riemann zeta functional symmetry equation
$$\zeta(1-s) = 2 \; (2\pi)^{-s} \; (s-1)! \; \cos(\frac{\pi}{2}s) \; \zeta(s).$$

Note how Mellin interpolation (MI) of the coefficients of an e.g.f. (a.k.a. Ramanujan's favorite Master Formula) underlies these transforms:
$$ \cos(2\pi n x) = \sum_{k \ge 0} \cos(\pi \frac{k}{2}) (2\pi n)^k \frac{x^k}{k!} = \sum_{k \ge 0} c_k \frac{x^k}{k!} = e^{c. x} ,$$
so to MI the coefficients, apply the normalized Mellin transform to the e.g.f. with the argument negated (in this case negation returns the same function)
$$\int_{0}^{\infty} e^{-c.x} \; \frac{x^{s-1}}{(s-1)!} \;  dx = (c.)^{-s} = c_{-s} $$
$$ = \int_{0}^{\infty} \frac{x^{s-1}}{(s-1)!} \; \cos(-2\pi n x) \; dx = \cos(\pi \frac{k}{2}) (2\pi n)^k  \; |_{k \to -s}. $$
For completeness, playing fast and loose with Dirac delta function/op reps, we can again apply MI via
$$
\int_{0}^{\infty} \frac{x^{s-1}}{(s-1)!} \; \delta(x-n) \; dx =\int_{0}^{\infty} \frac{x^{s-1}}{(s-1)!} \; \frac{1}{n} \delta(1-\frac{x}{n}) \; dx $$
$$
=\int_{0}^{\infty} \frac{x^{s-1}}{(s-1)!} \; \frac{1}{n} \frac{(1-\frac{x}{n})^{-1}}{(-1)!}  \; dx  = \int_{0}^{\infty} \frac{x^{s-1}}{(s-1)!} \;  \sum_{k \geq 0}(-1)^k \frac{1}{n^{k+1}} \; \frac{1}{(-k-1)!} \; \frac{x^k}{k!}  \; dx$$
$$ =\frac{1}{n^{k+1}} \; \frac{1}{(-k-1)!} \; |_{k \to -s} = \frac{1}{(s-1)!} \; n^{s-1} .$$
This is consistent with the limiting case of $H(1-x) \; \frac{(1-x)^{\omega}}{\omega!}$ as $\omega$ tends to $-1$ for the analytically continued integral rep of the Euler beta function, with $H(x)$ the Heaviside step function, and, therefore, fractional calculus. Being cautiously semi-conservative, one could look at the inverse Mellin transform rep of $\delta(x-n)$.
