I managed to find$$\int\limits_0^\infty \frac{\ln^{2a}(x)\ln(1+x)}{\sqrt{x}(1+x)}\mathrm{d}x=-\pi\lim_{m\to \frac12 }\frac{d^{2a}}{d m^{2a}} \frac{\psi(1-m) + \gamma}{\sin(m\pi)}.$$
To show this relation, we follow the same approach here:
Reduce $n$ by $m$ in the beta function:
$$\int_0^\infty\frac{x^{m-1}}{(1+x)^{m+n}}\mathrm{d}x=\operatorname{B}(m,n)=\frac{\Gamma(m)\Gamma(n)}{\Gamma(m+n)},$$
we have
$$\int_0^\infty\frac{x^{m-1}}{(1+x)^{n}}\mathrm{d}x=\frac{\Gamma(m)\Gamma(n-m)}{\Gamma(n)}.$$
Differentiate both sides $2a$ times with respect to $m$ and once with respect to $n$,
\begin{gather*}
\frac{\partial^{2a}}{\partial m^{2a}} \frac{\partial}{\partial n} \frac{\Gamma(m)\Gamma(n-m)}{\Gamma(n)}=\frac{\partial^{2a}}{\partial m^{2a}} \frac{\partial}{\partial n}\int\limits_0^\infty \frac{x^{m-1}}{(1+x)^n}\mathrm{d}x\\
\{\text{use differentiation under the integral sign theorem}\}\\
=\int\limits_0^\infty \frac{\partial^{2a}}{\partial m^{2a}} \frac{\partial}{\partial n}\frac{x^{m-1}}{(1+x)^n}\mathrm{d}x\\
=-\int\limits_0^\infty \frac{\ln^{2a}(x)\ln(1+x)x^{m-1}}{(1+x)^n}\mathrm{d}x.
\end{gather*}
Now take the limit on both sides letting $m\to 1/2$ and $n\to1$,
\begin{gather*}
-\int\limits_0^\infty \frac{\ln^{2a}(x)\ln(1+x)}{\sqrt{x}(1+x)}\mathrm{d}x=\lim_{\substack{m\to 1/2 \\ n \to 1}}\frac{\partial^{2a}}{\partial m^{2a}} \frac{\partial}{\partial n} \frac{\Gamma(m)\Gamma(n-m)}{\Gamma(n)}\\
=\lim_{\substack{m\to 1/2 \\ n \to 1}}\frac{\partial^{2a}}{\partial m^{2a}}\Gamma(m)\left( \frac{\partial}{\partial n} \frac{\Gamma(n-m)}{\Gamma(n)}\right)\\
=\lim_{\substack{m\to 1/2 \\ n \to 1}}\frac{\partial^{2a}}{\partial m^{2a}} \Gamma(m)\left(\frac{\Gamma(n-m)[{\psi}(n-m) -\psi(n)]}{\Gamma(n)}\right)\\
\{\text{evaluate the limit when $n\to 1$ and use $\psi(1)=-\gamma$}\}\\
=\lim_{m\to 1/2 }\frac{\partial^{2a}}{\partial m^{2a}} \Gamma(m)\Gamma(1-m)[\psi(1-m) + \gamma]\\
\left\{\text{use $\Gamma(m)\Gamma(1-m)=\frac{\pi}{\sin(m\pi )}$}\right\}\\
\left\{\text{and write $\frac{\partial}{\partial m}$ as $\frac{d}{dm}$, since we have one variable left}\right\}\\
=\pi\lim_{m\to \frac12 }\frac{d^{2a}}{d m^{2a}} \frac{\psi(1-m) + \gamma}{\sin(m\pi)}.
\end{gather*}
The question here is can we write the limit in terms of finite summation?
Best Answer
This problem deals with the relation
$$\int\limits_0^\infty \frac{\ln^{2a}(x)\ln(1+x)}{\sqrt{x}(1+x)}\mathrm{d}x=-\pi\lim_{m\to \frac12 }\frac{\mathrm{d}^{2a}}{\mathrm{d} m^{2a}} \frac{\psi(1-m) + \gamma}{\sin(m\pi)}\tag{*}$$
The solution to the problem of the OP (write the r.h.s. as a finite sum) is found in the section "finite sum" below.
Optionally, we also try to verify the relation $(*)$. The main task, the transformation of the integral to a sum (plus explicit terms), is done in the first section.
Transformation of the integral
EDIT 17.05.21
I just discovered that Mathematica solves the generating integral (version 8.0 immediately, 10.1 via the antiderivative)
$$\int_0^{\infty } \frac{x^z \log (x+1)}{x+1} \, dx=\pi \csc (\pi z) (\psi ^{(0)}(-z)+\gamma )$$
The powers of $\log(x)$ unter the integral can be generated by differentiation.
End EDIT
We wish to evaluate the integral
$$i = \int_0^{\infty}\frac{ \log(x)^{2a}\log(1+x)}{\sqrt{x} (1+x)}\,dx\tag{1}$$
Where $2a$ is specified implicitly in the OP as a positive integer.
Splitting the integration region we can write
$$i=i_1 + i_2\tag{2}$$
with
$$i_1 = \int_0^{1}\frac{ \log(x)^{2a}\log(1+x)}{\sqrt{x} (1+x)}\,dx\tag{3a}$$
$$i_2 = \int_1^{\infty}\frac{ \log(x)^{2a}\log(1+x)}{\sqrt{x} (1+x)}\,dx\tag{3b}$$
Now letting $x\to \frac{1}{y}$ in $i_2$ we can write (skipping some steps)
$$i_2 = (-1)^{2a} i_1 + i_3\tag{4}$$
where
$$\begin{align}i_3 & = \int_0^1 \frac{(-\log (y))^{2 a+1}}{\sqrt{y} (y+1)} \, dy \\ & = 2^{-2 (a+1)} \left(\zeta \left(2 a+2,\frac{1}{4}\right)-\zeta \left(2 a+2,\frac{3}{4}\right)\right) \Gamma (2 a+2)\end{align}\tag{5}$$
Observing now that
$$\frac{\log(1+x)}{1+x} = -\sum_{n=1}^{\infty} H_n (-x)^n\tag{6}$$
We can write
$$i_1 =i_{1s}:=(-2)^{2 a+1}\Gamma (2 a+1) \sum_{n=1}^{\infty}(-1)^n \frac{H_n}{(2 n+1)^{2 a+1}}\tag{7a}$$
so that finally the integral $(1)$ can be written as
$$i=\left((-1)^{2 a}+1\right) i_{1s}+i_3$$
Hence the task of the OP reduces write the following sum $s$ as a sum of a finite number of terms
$$s = \sum_{n=1}^{\infty} (-1)^n\frac{ H_n}{ (1+2n)^{2a+1}}\tag{8}$$
Observation: this sum is evaluated in [1] identity $(4.91)$ .
Finite sum
Here we show that the r.h.s. of $(*)$ can be written an as finite sum.
Starting directly from the formula in the OP, letting $2a = k$ I get for the first few limits these expressions (format {k, limit})
$$ \begin{array}{c} \left\{0,-\pi \left(\psi ^{(0)}\left(\frac{1}{2}\right)+\gamma \right)\right\} \\ \left\{1,\frac{\pi ^3}{2}\right\} \\ \left\{2,-\pi \left(\pi ^2 \psi ^{(0)}\left(\frac{1}{2}\right)+\gamma \pi ^2+\psi ^{(2)}\left(\frac{1}{2}\right)\right)\right\} \\ \left\{3,\frac{5 \pi ^5}{2}\right\} \\ \left\{4,-\pi \left(5 \pi ^4 \psi ^{(0)}\left(\frac{1}{2}\right)+5 \gamma \pi ^4+6 \pi ^2 \psi ^{(2)}\left(\frac{1}{2}\right)+\psi ^{(4)}\left(\frac{1}{2}\right)\right)\right\} \\ \left\{5,\frac{61 \pi ^7}{2}\right\} \\ \end{array}\tag{9}$$
Now we could try to guess the rule for the polygamma functions and write the powers of $\pi$ in terms of $\zeta$-functions. But it is easier to proceed systematically.
We expand the $k$-th derivative of the two factor expression into a finite binomial sum
$$(\frac{d}{dm})^k (A B) = \sum_{j=0}^{k}\binom{k}{j} A^{(k-j)} B^{(j)}\tag{10}$$
where
$$A = \psi (1-m)+\gamma, B = \frac{1}{\sin{\pi m}}\tag{11}$$
Now using
$$\frac{\partial \psi ^{(p)}(1-m)}{\partial m}=-\psi ^{(p+1)}(1-m)\tag{12}$$
and the partial fraction decomposition
$$\frac{\pi}{\sin (\pi m)}=\sum _{n=-\infty }^{\infty } \frac{(-1)^n}{m-n}\tag{13}$$
we obtain for the j-th derivative the following expression
$$\begin{align}\frac{\partial ^j}{\partial m^j}\left(\frac{\pi }{\sin (\pi m)}\right)|_{m\to \frac{1}{2}}=-2^{-j-1}j! \left((-1)^k+1\right) \\ \times \left(-\zeta \left(j+1,\frac{1}{4}\right)+\zeta \left(j+1,-\frac{1}{4}\right)+4^{k+1}\right)\end{align}\tag{14}$$
With $(10)$, $(11)$, $(12)$, and $(14)$ we have all ingredients to expand the r.h.s. of $(*)$ into a finite sum.
I obtain for the finite sum requested in the OP (letting $2a \to k$) the following expression
$$r(k) := -\pi\lim_{m\to \frac12 }\frac{\mathrm{d}^{k}}{\mathrm{d} m^{k}} \frac{\psi(1-m) + \gamma}{\sin(m\pi)} \\=U(k) + \sum_{j=1}^{k} \binom{k}{j} U(k-j) V(j)\tag{15}$$
where
$$U(p) = -\pi \left(\gamma \delta _{0,p}+(-1)^p \psi ^{(p)}\left(\frac{1}{2}\right)\right) \tag{15a}$$
$$\begin{align}V(p) =-\frac{1}{\pi} \left(2^{-p-1} \left((-1)^p+1\right) p! \\ \times \left(-\zeta \left(p+1,\frac{1}{4}\right)+\zeta \left(p+1,-\frac{1}{4}\right)+4^{p+1}\right)\right) \end{align}\tag{15b}$$
Here $\delta _{0,p}$ is the Kronecker delta.
The first few $r(k)$ after simplification (FullSimplify) are in the format ${k,r(k)}$
for even $k$
$$ \begin{array}{c} \{0,\pi \log (4)\} \\ \left\{2,\pi ^3 \log (4)+14 \pi \zeta (3)\right\} \\ \left\{4,5 \pi ^5 \log (4)+84 \pi ^3 \zeta (3)+744 \pi \zeta (5)\right\} \\ \left\{6,61 \pi ^7 \log (4)+1050 \pi ^5 \zeta (3)+11160 \pi ^3 \zeta (5)+91440 \pi \zeta (7)\right\} \\ \left\{8,1385 \pi ^9 \log (4)+23912 \pi ^7 \zeta (3)+260400 \pi ^5 \zeta (5)+2560320 \pi ^3 \zeta (7)+20603520 \pi \zeta (9)\right\} \\ \end{array}\tag{16a} $$
for odd $k$
$$ \begin{array}{l} \left\{1,\frac{\pi ^3}{2}\right\} \\ \left\{3,\frac{5 \pi ^5}{2}\right\} \\ \left\{5,\frac{61 \pi ^7}{2}\right\} \\ \left\{7,\frac{1385 \pi ^9}{2}\right\} \\ \left\{9,\frac{50521 \pi ^{11}}{2}\right\} \\ \end{array}\tag{16b} $$
This finalizes the solution of the problem in the OP.
Discussion
The numbers $c_{n}$ appearing here for both even and odd numbers
$$\{1, 1, 5, 61, 1 385, 50 521, 2 702 765, 199 360 981, ... \}$$
can be found in the impressivly long entry https://oeis.org/A000364:
A000364 Euler (or secant or "Zig") numbers.
They are related to the Euler polynomials: $c_{n}=2^n E_{n}(\frac{1}{2})$
References
[1] Ali Shadhar Olaikhan, "An introduction to harmonic series and logarithmic integrals", April 2021, ISBN 978-1-7367360-0-5