Edit 5.2.2014:

Using this as a starting point:

$$\Lambda(n)=\lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n} \frac{\mu(d)}{d^{(s-1)}}.$$

Where $\Lambda(n)$ is the von Mangoldt function.

Edit 25.3.2018 All of the pictures I posted above are either slightly or almost completely wrong. I deleted the last two. The correct plot starts from the following symbolic relationships:

Let $\mu(n)$ be the Möbius function, then:

$$a(n) = \sum\limits_{d|n} d \cdot \mu(d)$$

$$T(n,k)=a(GCD(n,k))$$

$$T = \left( \begin{array}{ccccccc} +1&+1&+1&+1&+1&+1&+1&\cdots \\ +1&-1&+1&-1&+1&-1&+1 \\ +1&+1&-2&+1&+1&-2&+1 \\ +1&-1&+1&-1&+1&-1&+1 \\ +1&+1&+1&+1&-4&+1&+1 \\ +1&-1&-2&-1&+1&+2&+1 \\ +1&+1&+1&+1&+1&+1&-6 \\ \vdots&&&&&&&\ddots \end{array} \right)$$

$$\sum\limits_{k=1}^{\infty}\sum\limits_{n=1}^{\infty} \frac{T(n,k)}{n^c \cdot k^s} = \sum\limits_{n=1}^{\infty} \frac{\lim\limits_{z \rightarrow s} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c} \;\;\;\;\;\;\;\;(1)$$
$$\sum\limits_{n=1}^{\infty} \frac{\lim\limits_{z \rightarrow s} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c} = \frac{\zeta(s) \cdot \zeta(c)}{\zeta(c + s - 1)} \;\;\;\;\;\;\;\;\;(2)$$

which is part of the limit:

$$\frac{\zeta '(s)}{\zeta (s)}=\lim_{c\to 1} \, \left( \color{Red}{\zeta (c)}-\frac{\zeta (c) \zeta (s)}{\zeta (c+s-1)}\right) \;\;\;\;\;\;\;(3)$$
$$\int \frac{\zeta '(s)}{\zeta (s)} ds = \log(\zeta(s)) \;\;\;\;\;\;\;\;\;\;(4)$$

And the zeta zero counting function is roughly:

$$f(t)=\frac{\vartheta (t)+\Im\left(\log \left(\zeta \left(i t+\frac{1}{2}\right)\right)\right)}{\pi }\;\;\;\;\;\;\;\;\;\;(5)$$

where $\vartheta (t)$ is the Riemann Siegel Theta function.

Since we are looking for the Fourier series like plot we instead write right hand side of $(3)$ as:

$$\frac{\zeta '(s)}{\zeta (s)} \approx \sum\limits_{n=1}^{\infty} \left( \color{Red}{\frac{1}{n^c}}-\frac{\lim\limits_{z \rightarrow s} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c}\right)\;\;\;\;\;\;\;\;\;\;(6)$$
Integrating $(6)$ we achieve:

$$\log(\zeta(s))\approx \int\sum\limits_{n=1}^{\infty} \left( \color{Red}{\frac{1}{n^c}}-\frac{\lim\limits_{z \rightarrow s} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c}\right)ds\;\;\;\;\;\;\;\;\;\;(7)$$
The colored $\color{Red}{\zeta(c)}$ in $(3)$ and $\color{Red}{\frac{1}{n^c}}$ in $(6)$ is equivalent to the first column in matrix $T$.

This integral can be solved symbolically with the Euler-Maclaurin formula for the Riemann zeta function, and combined with Riemann Siegel Theta in $(5)$ it is something like:

$f(t)=\small \frac{\vartheta (t)+\Re\left(\sum _{\text{nn}=1}^{\text{nnn}} \left(\frac{1}{\text{nn}^c}-\frac{\sum _{d=1}^{\text{nn}} \text{If}\left[\text{nn} \bmod d=0,\sum _{r=1}^{q-1} \left(\sum _{i=0}^{r-2} -\frac{B_r d^{1-\left(\frac{1}{2}+i t\right)} k^{-r-\left(\frac{1}{2}+i t\right)} \left|S_r^{(i)}\right| \text{Expand}\left[\sum _{m=0}^i \left(\frac{1}{2}+i t\right)^m (\log (d)+\log (k))^m N\left[\frac{i!}{m!}\right]\right]}{r! (\log (d)+\log (k))^{i+1}}\right)+\mu (d) \sum _{n=1}^k \text{If}\left[n=1\lor (n=1\land d=1),\frac{1}{2}+i t,-\frac{d^{1-\left(\frac{1}{2}+i t\right)}}{n^{\frac{1}{2}+i t} (\log (d)+\log (n))}\right]+\text{Ei}\left(-\left(\left(i t+\frac{1}{2}\right)-1\right) \log (d)-\left(\left(i t+\frac{1}{2}\right)-1\right) \log (k)\right),0\right]}{\text{nn}^c}\right)\right)}{\pi } \;(8)$
where $s$ has been replaced with $1/2+it$. $f(t)$ in $(8)$ is not computationally feasible to evaluate so instead we resort to looking at the numerical integral.

If we instead look at the plot of the **numerical integral** in the right hand side in $(9)$ below, then for $c=1$;
$$1+\frac{\vartheta (t)+\Im\left(\log \left(\zeta \left(i t+\frac{1}{2}\right)\right)\right)}{\pi } \\
\approx \frac{\vartheta (t)+\Re\left(\int\sum\limits_{n=1}^{\infty} \left( \frac{1}{n^c}-\frac{\lim\limits_{z \rightarrow \frac{1}{2} + it} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c}\right)ds\right)}{\pi}\;\;\;\;\;(9)$$

something much more Fourier series like appears:

where there is a jump of size one at every Riemann zeta zero.

The Mathematica program for the plot of the numerical integral in $(9)$ is:

```
(*start*)
Clear[n, s, c, t, nn, delta];
c = 1;
delta = 40;
nn = 60;
a = Sum[1/N[n]^c -
Zeta[1/2 + I*t]*
Total[MoebiusMu[Divisors[n]]/Divisors[n]^(1/2 + I*t - 1)]/
n^c, {n, 1, 600}];
ListLinePlot[(Table[
delta*RiemannSiegelTheta[t], {t, 0, nn, 1/delta}] +
Accumulate[Table[Re[a], {t, 0, nn, 1/delta}]])/delta/Pi,
PlotStyle -> Thickness[0.003], DataRange -> {0, nn}]
(*end*)
```

Edit 15.8.2018:

The correct way to integrate the Euler-Maclaurin formula is:

$$\boxed{N(t)=\frac{1}{\pi}\left(\vartheta (t)-\Re\left(\sum _{n=1}^{\text{nn}} \frac{1}{n^c} \left(\underset{d \mid n}{\sum\limits_{d=1}^{n}} \left(f(\frac{1}{2}+it, d)-f(\frac{1}{2}+i0, d) \right)\right)\right)\right)}$$

where $\vartheta (t)$ is the Riemann-Siegel theta function,

and where:

$$f(s,d)=-i \mu (d) \left(\sum _{n=0}^{q-1} \left(\sum _{i=0}^{2 n+1} -\frac{d B_{2 (n+1)} k^{-2 n-1} \left|S_{2 n+1}^{(i)}\right| \log ^{-i-1}(d k) \Gamma (i+1,s \log (k d))}{(2 (n+1))!}\right)+\sum _{n=2}^k -\frac{d^{1-s} n^{-s}}{\log (d)+\log (n)}+\text{If}\left[d=1 \text{ then } 0\text{ else }-\frac{1^{-s} d^{1-s}}{\log (d)+\log (1)}\right]+\frac{d^{1-s} k^{-s}}{2 (\log (d)+\log (k))}+\text{Ei}(-(s-1) \log (d)-(s-1) \log (k))\right)$$

and where $\mu (d)$ is the Möbius function.

As a Mathematica program the **symbolic** integral is:

```
(* Integral *)
(*start*)
Clear[k, n, s, i, d, r, q, f, a, b, integral, c, nn, n, k];
nn = 30;
delta = 1/5;
k = 30;
q = 10;
c = 1;
rr = 60;
f[s_, d_] = -I*
MoebiusMu[
d]*(If[d == 1, 0, -((d^(1 - s) 1^-s)/(Log[d] + Log[1]))] +
Sum[-((d^(1 - s) n^-s)/(Log[d] + Log[n])), {n, 2, k}] +
ExpIntegralEi[-(-1 + s) Log[d] - (-1 + s) Log[k]] + (
d^(1 - s) k^-s)/(2 (Log[d] + Log[k])) +
Sum[Sum[-(BernoulliB[2*(n + 1)]/((2*(n + 1))!))*(Abs[
StirlingS1[2*n + 1, i]]) d k^(-2*n - 1)
Gamma[1 + i, s Log[k*d]] Log[k*d]^(-1 - i), {i, 0,
2*n + 1}], {n, 0, q - 1}]);
integral[a_, b_, d_] = f[b, d] - f[a, d];
Print["Counting to ", rr];
Monitor[ListLinePlot[
Table[(RiemannSiegelTheta[t] -
Re[Sum[Sum[
If[Mod[n, d] == 0, integral[1/2 + I*0, 1/2 + I*t, d]/n^c,
0], {d, 1, n}], {n, 1, nn}]])/Pi, {t, 0, rr, delta}],
PlotStyle -> Thickness[0.003], DataRange -> {0, rr},
ImageSize -> Large], Floor[t]]
(* end *)
```