Norms of Wigner matrices under power law decay

Question

Suppose $\Sigma=\operatorname{diag}(h)$ where $h=(1^{-p},2^{-p},3^{-p},\ldots,d^{-p})$ and $p> 1$

$X$ is a matrix with $b$ rows sampled independently from $\operatorname{Normal}(0,\Sigma)$

Suppose $d$ is large. How does $\left\|\frac{1}{b}XX^T\right\|_{\operatorname{op}}$ depend on $b$?

For $p=1,d=30$, the formula below gives a nice fit to the average value observed in simulation, can it be improved/generalized?

$$\left\|\frac{1}{b}XX^T\right\|_{\operatorname{op}}\approx\frac{1}{b}\left(\|h\|_1+(1+b)\|h\|_\infty\right)$$

Notebook

The formula comes from analogy with Frobenius norm where Wick's theorem can be applied.