So it might help to just start with one state. We diagonalize $\rho$ and we find out that it is
$$\rho =\rho_1 |\psi_1\rangle\langle\psi_1| +\dots+\rho_d|\psi_d\rangle\langle\psi_d| $$.
Define $$|e_0\rangle = \sqrt{\frac1d} \big(|\psi_1\rangle+ |\psi_2\rangle + \dots + |\psi_d\rangle\big).$$
Now observe, the occupation in this state is $$\operatorname{Tr}\big(|e_0\rangle\langle e_0|~\rho\big)=\langle e_0|\rho|e_0\rangle = \frac1d \sum_{k=1}^d \rho_d.$$
That sum is conveniently 1, so we have our first vector.
What about our next vector? Well suppose $d$ is even, consider
$$|e_1\rangle = \sqrt{\frac1d} \big(|\psi_1\rangle -|\psi_2\rangle + \dots - |\psi_d\rangle\big).$$
This is clearly orthogonal to $|e_1\rangle$, but it has the exact same sum when we form $
\operatorname{Tr}\big(|e_1\rangle\langle e_1|~\rho\big)
$ so we have our second vector.
Can also try going $1,i,-1,-i$ and if $d$ is divisible by 4, this will be orthogonal to both of the above and will generate the same $1/d$.
In fact, any phase factors leading the terms will allow us to get the same $1/d$ and so the only thing we have to focus on is orthogonality. So the basic idea there is just that this discrete Fourier transform arrangement does the job well. Indeed both of the above are special cases of this. If you have a sum of nth roots of unity, the finite geometric series formula automatically forces the sum to be zero, which you can kind of intuit because that sum calculates the center of mass of a bunch of masses evenly spread out on a ring of radius one around an origin. This basic rule gives you your Kronecker delta for your overlap integrals $\langle e_i | e_k\rangle =\delta_{ik}$, so you have orthogonality and you have the trivial probability $1/d$ and everything is good.
