  A set of $k$ orthonormal bases of $\mathbb C^d$ is called mutually unbiased
if $|\langle e,f\rangle |^2 = 1/d$ whenever $e$ and $f$ are basis vectors in
distinct bases. A natural question is for which pairs $(d,k)$ there exist~$k$
mutually unbiased bases in dimension $d$. The (well-known) upper bound $k \leq
d+1$ is attained when~$d$ is a power of a prime. For all other dimensions it is
an open problem whether the bound can be attained. Navascu\'es, Pironio, and
Ac\'in showed how to reformulate the existence question in terms of the
existence of a certain $C^*$-algebra. This naturally leads to a noncommutative
polynomial optimization problem and an associated hierarchy of semidefinite
programs. The problem has a symmetry coming from the wreath product of $S_d$
and $S_k$.
  We exploit this symmetry (analytically) to reduce the size of the
semidefinite programs making them (numerically) tractable. A key step is a
novel explicit decomposition of the $S_d \wr S_k$-module $\mathbb C^{([d]\times
[k])^t}$ into irreducible modules. We present numerical results for small $d,k$
and low levels of the hierarchy. In particular, we obtain sum-of-squares proofs
for the (well-known) fact that there do not exist $d+2$ mutually unbiased bases
in dimensions~$d=2,3,4,5,6,7,8$. Moreover, our numerical results indicate that
a sum-of-squares refutation, in the above-mentioned framework, of the existence
of more than $3$ MUBs in dimension $6$ requires polynomials of total degree at
least~$12$.
