A series involves harmonic number How do we get a closed form for
$$\sum_{n=1}^\infty\frac{H_n}{(2n+1)^2}$$
 A: $\newcommand{\+}{^{\dagger}}
 \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
 \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
 \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
 \newcommand{\dd}{{\rm d}}
 \newcommand{\down}{\downarrow}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
 \newcommand{\fermi}{\,{\rm f}}
 \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
 \newcommand{\half}{{1 \over 2}}
 \newcommand{\ic}{{\rm i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\isdiv}{\,\left.\right\vert\,}
 \newcommand{\ket}[1]{\left\vert #1\right\rangle}
 \newcommand{\ol}[1]{\overline{#1}}
 \newcommand{\pars}[1]{\left(\, #1 \,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\pp}{{\cal P}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
 \newcommand{\sech}{\,{\rm sech}}
 \newcommand{\sgn}{\,{\rm sgn}}
 \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}
 \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
 \newcommand{\wt}[1]{\widetilde{#1}}$
$\ds{\sum_{n = 1}^{\infty}{H_{n} \over \pars{2n + 1}^{2}}:\ {\large ?}}$.

Lets consider
  $\ds{\fermi\pars{x}\equiv
     \sum_{n = 1}^{\infty}{H_{n} \over \pars{2n + 1}^{2}}\,x^{2n + 1}.
     \qquad\fermi\pars{1}={\large ?}\,,\quad \fermi\pars{0} = 0}$.

\begin{align}
\fermi'\pars{x}&=\sum_{n = 1}^{\infty}{H_{n} \over 2n + 1}\,x^{2n}
\ \imp\
\bracks{x\fermi'\pars{x}}'=\sum_{n = 1}^{\infty}H_{n}\,x^{2n}
=-\,{\ln\pars{1 - x^{2}} \over 1 - x^{2}}\,,\qquad\fermi'\pars{0} = 0
\end{align}
where we used the
Harmonic Number Generating Function.

Then
  \begin{align}
&x\fermi'\pars{x}=-\int_{0}^{x}{\ln\pars{1 - t^{2}} \over 1 - t^{2}}\,\dd t
\\[3mm]&\imp
\fermi\pars{1}=-\int_{0}^{1}{\dd x \over x}\int_{0}^{x}
{\ln\pars{1 - t^{2}} \over 1 - t^{2}}\,\dd t
=-\int_{0}^{1}{\ln\pars{1 - t^{2}} \over 1 - t^{2}}\int_{t}^{1}{\dd x \over x}
\,\dd t
\end{align}

$$\begin{array}{|c|}\hline\\
\quad\sum_{n = 1}^{\infty}{H_{n} \over \pars{2n + 1}^{2}}
=\int_{0}^{1}{\ln\pars{t}\ln\pars{1 - t^{2}} \over 1 - t^{2}}\,\dd t\quad
\\ \\ \hline
\end{array}
$$

\begin{align}
&\color{#c00000}{\sum_{n = 1}^{\infty}{H_{n} \over \pars{2n + 1}^{2}}}
=\int_{0}^{1}{\ln\pars{t^{1/2}}\ln\pars{1 - t} \over 1 - t}\,\half\,t^{-1/2}
\,\dd t
={1 \over 4}\int_{0}^{1}{t^{-1/2}\ln\pars{t}\ln\pars{1 - t} \over 1 - t}\,\dd t
\\[3mm]&={1 \over 4}\lim_{\mu\ \to\ 0 \atop{\vphantom{\LARGE A}\nu\ \to\ 0}}
\partiald{}{\mu}\partiald{}{\nu}\int_{0}^{1}t^{\mu - 1/2}
\pars{1 - t}^{\nu - 1}\,\dd t
={1 \over 4}\lim_{\mu\ \to\ 0 \atop{\vphantom{\LARGE A}\nu\ \to\ 0}}
\partiald{}{\nu}\Gamma\pars{\nu}\partiald{}{\mu}
\bracks{\Gamma\pars{\mu + 1/2} \over \Gamma\pars{\mu + \nu + 1/2}}
\\[3mm]&={1 \over 4}\lim_{\nu\ \to\ 0}
\partiald{}{\nu}\braces{%
\Gamma\pars{\nu}\,{\Gamma\pars{1/2} \over \Gamma\pars{\nu + 1/2}}
\bracks{\Psi\pars{\half} - \Psi\pars{\nu + \half}}}
\\[3mm]&=-\,{1 \over 4}\,\Gamma\pars{\half}\lim_{\nu\ \to\ 0}
\partiald{}{\nu}\bracks{%
{\Gamma\pars{\nu + 1} \over \Gamma\pars{\nu + 1/2}}\,
{\Psi\pars{1/2 + \nu} - \Psi\pars{1/2} \over \nu}}
\\[3mm]&=-\,{1 \over 4}\,\Gamma\pars{\half}\lim_{\nu\ \to\ 0}
\partiald{}{\nu}\braces{%
{\Gamma\pars{\nu + 1} \over \Gamma\pars{\nu + 1/2}}\,
\bracks{\Psi'\pars{\half} + \half\,\Psi''\pars{\half}\nu}}
\\[3mm]&={\pi^{2}\gamma + \pi^{2}\Psi\pars{1/2} + 14\zeta\pars{3} \over 8}
\quad\mbox{where we used}\quad\Psi\pars{1} = -\gamma\,,\quad
\Psi''\pars{\half} = -14\zeta\pars{3}.
\end{align}

With $\ds{\Psi\pars{\half} = -2\ln\pars{2} - \gamma}$:
$$
\color{#66f}{\large\sum_{n = 1}^{\infty}{H_{n} \over \pars{2n + 1}^{2}}
={1 \over 4}\,\bracks{7\zeta\pars{3} - \pi^{2}\ln\pars{2}}}
\approx {\tt 0.3933}
$$
A: Here's another solution. I'll denote various versions of the sum
$$
\sum_{k=1}^\infty\sum_{j=1}^k\frac1j\frac1{k^2}
$$
by an $S$ with two subscripts indicating which parities are included, the first subscript referring to the parity of $j$ and the second to the parity of $k$, with '$\mathrm e$' denoting only the even terms, '$\mathrm o$' denoting only the odd terms, '$+$' denoting the sum of the even and odd terms, i.e. the regular sum, and '$-$' denoting the difference between the even and the odd terms, i.e. the alternating sum. Then
$$
\begin{align}
\sum_{n=1}^\infty\frac{H_n}{(2n+1)^2}
&=
2\sum_{n=1}^\infty\sum_{i=1}^n\frac1{2i}\frac1{(2n+1)^2}
\\
&=
2S_{\mathrm{eo}}
\\
&=
2(S_{++}-S_{\mathrm o+}-S_{\mathrm{ee}})
\\
&=
2\left(S_{++}-S_{\mathrm o+}-\frac18S_{++}\right)
\\
&=
2\left(\frac38S_{++}+\left(\frac12S_{++}-S_{\mathrm o+}\right)\right)
\\
&=
\frac34S_{++}+S_{-+}
\\
&=
\frac32\zeta(3)+\sum_{k=1}^\infty\sum_{j=1}^k\frac{(-1)^j}j\frac1{k^2}\;,
\end{align}
$$
where I used the result $\sum_nH_n/n^2=2\zeta(3)$ from the blog post Aeolian linked to and reduced the present problem to finding the analogue of that result with the sign alternating with $j$, which we can rewrite as
$$
\begin{align}
\sum_{k=1}^\infty\sum_{j=1}^k\frac{(-1)^j}j\frac1{k^2}
&=
\sum_{k=1}^\infty\sum_{j=1}^\infty\frac{(-1)^j}j\frac1{k^2}-\sum_{k=1}^\infty\sum_{j=k+1}^\infty\frac{(-1)^j}j\frac1{k^2}
\\
&=
-\zeta(2)\log2+\sum_{j=1}^\infty\frac{(-1)^j}{j+1}\sum_{k=1}^j\frac1{k^2}\;.
\end{align}
$$
This last double sum can be evaluated by the method applied in the blog post, making use of the fact that summing the coefficients of a power series in $x$ corresponds to dividing it by $1-x$:
$$
\begin{align}
\sum_{j=1}^\infty x^j\sum_{k=1}^j\frac1{k^2}=\def\Li{\operatorname{Li}}\frac{\Li_2(x)}{1-x}\;,
\end{align}
$$
where $\Li_2$ is the dilogarithm. Thus
$$
\begin{align}
\sum_{j=1}^\infty\frac{(-1)^j}{j+1}\sum_{k=1}^j\frac1{k^2}
&=
\int_0^1\sum_{j=1}^\infty (-x)^j\sum_{k=1}^j\frac1{k^2}\mathrm dx
\\
&=
\int_0^1\frac{\Li_2(-x)}{1+x}\mathrm dx
\\
&=
\left[\Li_2(-x)\log(1+x)\right]_0^1+\int_0^1\frac{\log^2(1+x)}x\mathrm dx
\\
&=-\frac{\zeta(2)}2\log2+\frac{\zeta(3)}4\;,
\end{align}
$$
where the boundary term is evaluated using $\Li_2(-1)=-\eta(2)=-\zeta(2)+2\zeta(2)/4=-\zeta(2)/2$ and the integral in the second term is evaluated in this separate question. Putting it all together, we have
$$
\begin{align}
\sum_{n=1}^\infty\frac{H_n}{(2n+1)^2}
&=
\frac74\zeta(3)-\frac32\zeta(2)\log2
\\
&=
\frac74\zeta(3)-\frac{\pi^2}4\log2\;.
\end{align}
$$
I believe all the rearrangements can be justified, despite the series being only conditionally convergent in $j$, by considering the partial sums with $j$ and $k$ both going up to $M$; then all the rearrangements can be carried out within that finite square of the grid, and the sums of the remaining terms go to zero with $M\to\infty$.
A: I gave an integral representation for a more general form. Here is an integral representation for your sum

$$\sum_{n=1}^\infty\frac{H_n}{(2n+1)^2}= \frac{1}{4}\,\int_{0}^{1}\!{\frac {\ln  \left( 1-z \right) \ln  \left( z\right) }{z\sqrt {1-z}}}{dz}= \frac{1}{4}(7\,\zeta  \left( 3 \right) -{\pi }^{2}\ln  \left( 2 \right))\sim 0.393327464. $$

The above integral can be evaluated through beta function. Here is the technique from previous problems. Basically, you need to consider the integral

$$ \int_{0}^{1} z^s (1-z)^{w-1/2} dz. $$

A: using the following identity proved by Random Variable here
$$S= \sum_{n=1}^{\infty} \frac{H_{n}}{ (n+a)^{2}}= \left(\gamma + \psi(a) \right) \psi_{1}(a) - \frac{\psi_{2}(a)}{2} \, , \quad a >0.$$
take $\ a=1/2$
$$S= \sum_{n=1}^{\infty} \frac{H_{n}}{ (2n+1)^{2}}=\frac74\zeta(3)-\frac32\ln2\zeta(2)$$
a similar identity was proved here by the mathematician Anthony Sofo, when he published some related work in 2011.
