I'm interested in finding two sets of $N$ unitary $N \times N$ matrices $U_{1}, \ldots, U_{N}$, $V_{1}, \ldots, V_{N}$ such that:

$ \sup\limits_{X, Y}\sum\limits_{j,k = 1}^{N} |\mathrm{Tr}(YU_{j}XV_{k}^{\ast})|^2 \ll N$

where the sup is over all Hermitian traceless matrices with Frobenius norm $1$ nad $\ll$ means that the sum is asymptotically smaller than $N$.

The motivation comes from a problem in quantum information theory. An explicit construction of such sets of $U_{j}$ and $V_{k}$ would be the most desirable result, although pseudorandom constructions would also be very interesting. This statement is true if instead of taking unitaries we take matrices with independent Gaussian entries (rescaled properly so that on average each column has norm 1 etc.); however, I don't know how to show that random unitaries satisfy this property, so a hint for such a proof would also be welcome.

I would like to think of such unitaries as having some sort of "spread" property similar to objects showing up in the study of pseudorandomness like mutually unbiased bases, randomness extractors etc. However, I'm not sure if there is any definite connection to those notions.