Question
I would first like to show that the equality below holds: $$ -1 + \frac{1}{2}\left( 2k+1 \right)\ln\left( \frac{1 + \frac{1}{2k+1}}{1 - \frac{1}{2k+1}} \right) = \sum^\infty_{l=1} \frac{1}{\left( 2l+1\right)\left( 2k+1 \right)^{2l}}. $$
Then the following inequality for the lower bound (the upper bound is easily shown by doing the infinite sum, but I am having trouble checking the first inequality for the lower bound): $$ \frac{1}{12\left(n+\frac{1}{2}\right)} < \sum^\infty_{k=n} \frac{1}{3}\frac{1}{\left(2k+1\right)^2}<\epsilon_n < \sum^\infty_{k=n}\sum^\infty_{l=1}\frac{1}{3}\frac{1}{\left(2k+1\right)^{2l}} = \frac{1}{12n}. $$
Details of the derivation
The logarithm of the factorial can be written as $$ \ln n! = \left( n + \frac{1}{2} \right) \ln n - n + \left( 1 + \int^n_1\frac{x - \lfloor x \rfloor - \frac{1}{2}}{x}\mathrm{d}x \right). $$
Let
$$C = \exp\left(1 + \int^\infty_1 \frac{x - \lfloor x \rfloor - \frac{1}{2}}{x}\mathrm{d}x \right) $$
and
$$\epsilon_n = -\int^\infty_n \frac{x - \lfloor x \rfloor - \frac{1}{2}}{x}\mathrm{d}x. $$
It can be shown that as $n$ approaches infinity,
$$ C \sim \sqrt{2\pi}. $$
The original factorial expression is then $$ n! = Cn^{n+\frac{1}{2}}e^{-n+\epsilon_n} \sim \sqrt{2\pi}n^{n+\frac{1}{2}}e^{-n+\epsilon_n}. $$
Knowing this and the series form from (1), and using the inequality from (2), we can bound the approximation as: $$ \sqrt{2\pi}n^{n+\frac{1}{2}}\exp\left[-n+\frac{1}{12\left(n+\frac{1}{2}\right)}\right] < n! < \sqrt{2\pi}n^{n+\frac{1}{2}}\exp\left[-n+\frac{1}{12n}\right], $$ and $$ n! \sim \sqrt{2\pi}n^{n+\frac{1}{2}} e^{-n+\frac{1}{12n}}. $$
For the first question, working the rhs, let $x=\frac 1{2k+1}$
$$\sum^\infty_{l=1} \frac{1}{\left( 2l+1\right)\left( 2k+1 \right)^{2l}}=\sum^\infty_{l=1} \frac{x^{2l}}{\left( 2l+1\right)}=\frac 1 x\sum^\infty_{l=1} \frac{x^{2l+1}}{\left( 2l+1\right)}$$ $$\sum^\infty_{l=1} \frac{x^{2l+1}}{\left( 2l+1\right)}=-x+\sum^\infty_{l=0} \frac{x^{2l+1}}{\left( 2l+1\right)}=\tanh ^{-1}(x)-x$$
Use the logarithmic representation of the hyperbolic arctangent function, replace $x$ by its definition to obtain the lhs.
For the second one, you can improve the bounds (assuming that you do not want to use Stirling series).
Consider $$F=\log(n!)-\log\left(\sqrt{2\pi}n^{n+\frac{1}{2}}\exp\left[-n+\frac{1}{12\left(n+a\right)}\right] \right)$$ and expand $$F=\frac{a}{12 n^2}-\frac{30 a^2+1}{360 n^3}+\frac{a^3}{12 n^4}-\frac{105 a^4-1}{1260 n^5}+O\left(\frac{1}{n^6}\right)$$
Let $a_2$ be the solution of $$\frac{a}{12 n^2}-\frac{30 a^2+1}{360 n^3}=0$$ that is to say $$a_2=\frac{n}{2}-\frac{1}{2} \sqrt{n^2-\frac{2}{15}}$$ and then $$F_2=\frac{1}{1260 n^5}+O\left(\frac{1}{n^7}\right)$$ Let $a_3$ be the solution of $$\frac{a}{12 n^2}-\frac{30 a^2+1}{360 n^3}+\frac{a^3}{12 n^4}=0$$ that is to say $$a_3=\frac{n}{3} \left(2 \sqrt{2} \sinh \left(\frac{1}{3} \text{csch}^{-1}\left(-\frac{40 \sqrt{2} n^2}{70 n^2-9}\right)\right)+1\right)$$ which also leads to $$F_3=\frac{1}{1260 n^5}+O\left(\frac{1}{n^7}\right)$$
So, stay with $a_2$ which is (as one could expect)
$$a_2=\frac{n}{2}-\frac{1}{2} \sqrt{n^2-\frac{2}{15}}=\frac{1}{30 n}+\frac{1}{900 n^3}+O\left(\frac{1}{n^5}\right)$$
and $e^{F_2}$ is a very tight lower bound of $n!$