I'm looking for $d^*$, the smallest value of $d$ such that the following holds: $$\frac{\left(\sum_{i=1}^d i^{-p}\right)^2}{\sum_{i=1}^d i^{-2p}}\ge \frac{1}{2}\frac{\zeta(p)^2}{\zeta(2p)}$$
with $\zeta$ referring to Riemann Zeta function.
Empirically this seems to scale as the function below, why?
$$d^*\approx c_1 \exp\left(\frac{c_2}{p-1}\right)$$
Edit
Mathematica user found a way to express $c_1$ and $c_2$ in terms of Euler gamma constant and $\sqrt{2}$
Motivation: $\frac{\zeta(p)^2}{\zeta(2p)}$ gives "effective rank" of infinite-dimensional linear regression problem with eigenvalues $1^{-p},2^{-p},\ldots$ (ie, how many observations are needed for prediction). This question gives $d$ such that finite-dimensional truncation after $d$ terms has similar effective rank.
ùToo long for comments.*
The problem is to find, for a given value of $p$, the zero of function $$f(d)=2 \zeta (2 p) \left(H_d^{(p)}\right){}^2-\zeta (p)^2 H_d^{(2 p)}$$
Making the problem continuous, let $d=e^x$ and the function is quite nice (except that it shows a minimum value for small values of $x$.
Consider instead $$g(x)=\log \left(2 \zeta (2 p) \left(H_{e^x}^{(p)}\right){}^2\right)-\log \left(\zeta (p)^2 H_{e^x}^{(2 p)}\right)$$ which is very smooth and Newton method converges quitae fast.
Trying for $p=1.1$ and the ridculous $x_0=1$, the iterates are $$\left( \begin{array}{cc} n & x_n \\ 0 & 1.00000 \\ 1 & 4.90872 \\ 2 & 8.71146 \\ 3 & 11.1061 \\ 4 & 11.6858 \\ 5 & 11.7114 \\ 6 & 11.7115 \\ \end{array} \right)$$
Repeating you calculations for $1.001 \leq p \leq 1.200$, the fit you proposed is almost perfect; what I obtained for the empirical model $$x=-a + \frac b {p-1}$$ is $$\begin{array}{l|lll} \text{} & \text{Estimate} & \text{Std Error} & \text{Confidence Interval} \\ \hline a & 0.56830 & 5.95\times 10^{-4} & \{0.56712,0.56947\} \\ b & 1.22797 & 6.57\times 10^{-6} & \{1.22796,1.22799\} \\ \end{array}$$
I have not been able to find any way to conclude using asymptotic analysis.