The Euler sums are given by $$S_{p,q} = \sum_{n = 1}^{\infty} \frac{H_{n}^{(p)}}{n^q},$$ where $$H_{n}^{(p)} = \sum_{j = 1}^{n} \frac{1}{j^p}.$$ According to Wolfram, Eq. (19), the following special case holds: $$S_{1,3} \equiv \sum_{n = 1}^{\infty} \frac{H_{n}}{n^3} = \frac{5}{4}\zeta(4),$$ where $H_n = H_n^{(1)}$ is the $n$-th harmonic number. I am interested in the following "modified Euler sum": $$\sum_{n = 1}^{\infty} \frac{H_{2n}}{(2n)^3} \equiv \frac{1}{8} \sum_{n = 1}^{\infty} \frac{H_{2n}}{n^3}.$$ Can it as well be expressed in terms of $\zeta$-values? I have not been able to locate any relevant literature.
Update I: For what it may be worth (it is not even strictly related to the original question of this post, although it does feature sums containing $H_{2n}$), I have been able to derive the following identity: $$\frac{124}{5}\zeta(5) = 2\pi^{2}\ln2 + \frac{16}{\pi^{2}} \sum_{n=1}^{\infty}[\frac{6H_{n}}{(2n+1)^{4}} + \frac{H_{n}^{(2)}}{(2n+1)^{3}} - \frac{12H_{2n}}{(2n+1)^{4}} - \frac{4H_{2n}^{(2)}}{(2n+1)^{3}}],$$ where $H_{n}^{(2)} \equiv \sum_{k=1}^{n} 1/k^{2}$, as standardly. This has been done by considering a contour integral at infinity (as a limit, of course) of the function $$f(s) = \frac{\Psi(s)}{s^{3}\cos(\pi s)^{2}},$$ the corresponding sum of all residues having to vanish. However, that identity leaves me no less stuck.
Update II: Using the following family of functions: $$F_{k}(s) = \frac{\Psi(s)}{s^{k}\sin(\pi s)}[1 + \cos(\pi s)],$$ where $k$ is a natural number, I have been able to establish the following identity [note the condition on $k$]: $$\sum_{n=1}^{\infty}\frac{H_{2n}}{(2n)^{k}} = 2^{-k} \gamma \zeta(k) + 2^{-(k+2)}(k+1)\zeta(k+1) - \frac{\pi}{4} \rm{Res}(F_{k};0),\text{ for } k \text{ even},$$ where $\rm{Res}(F_{k};0)$ denotes, of course, the residue of $F(s)$ at $s = 0$. It implies the following specific cases [note the presence of $n^{k}$ in the denominators of the sums, not $(2n)^{k}$ as above] \begin{align} \sum_{n=1}^{\infty}\frac{H_{2n}}{n^{2}} &= \frac{11 \zeta(3)}{4},\\ \sum_{n=1}^{\infty}\frac{H_{2n}}{n^{4}} &= \frac{37 \zeta(5)}{4} - 4\zeta(2)\zeta(3),\\ \sum_{n=1}^{\infty}\frac{H_{2n}}{n^{6}} &= \frac{135 \zeta(7)}{4} - 16\zeta(2)\zeta(5) - 4\zeta(4)\zeta(3). \end{align} Unfortunately, the method fails for k odd, in the sense that the resulting condition (of the sum of residues having to vanish) does not contain any harmonic numbers.
An attempt to get an expression through the integrals.
$$\sum_{n \geq 1}\frac{H_{2n}}{n^3}=\sum_{n \geq 1}\frac{1}{n^3} \int_0^1 \frac{1-x^{2n}}{1-x}dx=$$
$$\int_0^1 \frac{dx}{1-x} \left(\zeta(3)-\sum_{n \geq 1}\frac{x^{2n}}{n^3} \right)=\int_0^1 \frac{dx}{1-x} \left(\zeta(3)-\text{Li}_3 (x^2)\right)=$$
Let's change the variable to $x=e^{-t}$, then we have:
$$=\int_0^\infty \frac{ dt}{e^t-1} \left(\zeta(3)-\text{Li}_3 (e^{-2t})\right)=$$
There are known integral expressions:
$$\zeta(3)=\frac{1}{2} \int_0^\infty \frac{u^2 du}{e^u-1}$$
$$\text{Li}_3 (e^{-2t})=\frac{1}{2} \int_0^\infty \frac{u^2 du}{e^{2t+u}-1}$$
So we can write:
$$=\frac{1}{2} \int_0^\infty \int_0^\infty\frac{ dt}{e^t-1}\left( \frac{1}{e^u-1}-\frac{1}{e^{2t+u}-1} \right) u^2 du=$$
$$=\frac{1}{2} \int_0^\infty \int_0^\infty\frac{ dt}{e^t-1} \frac{e^u (e^{2t}-1)u^2 du}{(e^u-1)(e^{2t+u}-1)} = \frac{1}{2} \int_0^\infty \int_0^\infty \frac{e^u (e^t+1)u^2 du~dt }{(e^u-1)(e^{2t+u}-1)}$$
So we have:
$$\sum_{n \geq 1}\frac{H_{2n}}{n^3}=\frac{1}{2} \int_0^\infty \int_0^\infty \frac{e^u (e^t+1)u^2 du~dt }{(e^u-1)(e^{2t+u}-1)}$$
We could set:
$$e^{-t}=x \\ e^{-u}=y$$
Then
$$\sum_{n \geq 1}\frac{H_{2n}}{n^3}=\frac{1}{2} \int_0^1 \int_0^1 \frac{ (1+x) \ln^2 y ~dx~dy }{(1-y)(1-x^2 y)}$$
Which could be separated into two distinct parts:
$$\sum_{n \geq 1}\frac{H_{2n}}{n^3}=\frac{1}{2} \int_0^1 \int_0^1 \frac{ \ln^2 y ~dx~dy }{(1-y)(1-x^2 y)}+\frac{1}{4} \int_0^1 \int_0^1 \frac{ \ln^2 y ~dx~dy }{(1-y)(1-x y)}=$$
$$=\frac{1}{2} \int_0^1 \int_0^1 \frac{ \ln^2 y ~dx~dy }{(1-y)(1-x^2 y)}-\frac{1}{4} \int_0^1 \frac{ \ln^2 y \ln(1-y) ~dy }{y(1-y)}$$
Mathematica gives for that last integral:
$$-\frac{1}{4} \int_0^1 \frac{ \ln^2 y \ln(1-y) ~dy }{y(1-y)}=\frac{\pi^4}{144}$$
So we can write (also taking the first integral w.r.t. $x$):
$$\sum_{n \geq 1}\frac{H_{2n}}{n^3}=\frac{1}{2} \int_0^1 \frac{ \tanh^{-1} \sqrt{y} \ln^2 y ~dy }{\sqrt{y}(1-y)}+\frac{\pi^4}{144}$$
The first integral can be simplified:
No idea what to do with the last integral, Mathematica can't take it. However, the integrand is nice, there's not even any singularities.
We could try integration by parts with $$u=\tanh^{-1} y \ln^2 y \\ dv=\frac{dy }{1-y^2}$$
$$I=\int_0^1 \frac{ \tanh^{-1} y \ln^2 y ~dy }{1-y^2}=(\tanh^{-1} y)^2 \ln^2 y \bigg|_0^1-I-2 \int_0^1 (\tanh^{-1} y)^2 \ln y ~\frac{dy}{y}$$
So:
$$I=-\int_0^1 (\tanh^{-1} y)^2 \ln y ~\frac{dy}{y}= \int_0^\infty t (\tanh^{-1} e^{-t})^2 dt$$
Not really a simplification, but also an interesting expression.
Moreover, it is very fitting for numerical computation of the integral, because for moderately large $t$ we have $$\tanh^{-1} e^{-t} \asymp e^{-t}$$ and so the integral becomes elementary.
Hence, we can write:
$$\sum_{n \geq 1}\frac{H_{2n}}{n^3} = \frac{\pi^4}{144}+4\int_0^\infty t (\tanh^{-1} e^{-t})^2 dt$$
$$\sum_{n \geq 1}\frac{H_{2n}}{n^3} \approx \frac{\pi^4}{144}+(1+2 \tau) e^{-2 \tau}+ 4\int_0^\tau t (\tanh^{-1} e^{-t})^2 dt$$
Where the last integral can be evaluated with any suitable numerical quadrature. $\tau=5$ already gives $6$ accurate digits for the series.