Suppose $H$ is a random positive definite matrix normalized to have norm 1. Can we say anything about behavior of $f(s)$?

$$f(s)=\operatorname{Tr}[H(I-H)^s]$$

Taking $H=A^T A$ with entries of $A$ sampled from standard normal or uniform(-1,1) produces essentially identical curves, which suggests universality. What is the shape of this curve?

If we represent $i$th eigenvalue of $H$ as $y=h(i)$, we can approximate the sum in terms of integral which reduces to Laplace transform

$$f(t)=\int_0^\infty h(i) \exp(-t h(i))=\mathcal{L}(yh^{-1}(y)')$$

However, plugging this into Mathematica for $h(i)$ representing semicircle law gives something unintelligible in terms of Bessel and Struve functions.

Notebook

Motivation: $f(2s)$ gives the loss observed after $s$ steps of gradient descent on a quadratic $H$ (math). I'm interested in understanding behavior of this decay for some families of random matrices so that I could try fitting them to existing loss trajectories on real problems (unknown $H$). Theoretically, all eigenvalues of $H$ can be reconstructed from enough values of $f(t)$ by using inverse Laplace transform, but it needs a lot of values and not numerically stable