# Black hole encircled by a thin disk: fully relativistic solution111Dedicated
to our teacher, colleague and friend, Professor Oldřich Semerák, on the
occasion of his 60th orbit around the Sun. Oldřich’s well-known lectures –
famous not only among us (former) students of ‘math-phys’ in Prague – were the
first impulse that led us to devote our studies to general relativity.

Petr Kotlařík Institute of Theoretical Physics, Faculty of Mathematics and
Physics, Charles University, Prague, Czech Republic  David Kofroň Institute of
Theoretical Physics, Faculty of Mathematics and Physics, Charles University,
Prague, Czech Republic

###### Abstract

We give a full metric describing the gravitational field of a static and
axisymmetric thin disk without radial pressure encircling a Schwarzschild
black hole. The disk density profiles are astrophysically realistic,
stretching from the horizon to radial infinity, yet falling off quickly at
both these locations. The metric functions are expressed as finite series of
Legendre polynomials. Main advantages of the solution are that (i) the disks
have no edges, so their fields are everywhere regular (outside the horizon),
and that (ii) all non-trivial metric functions are provided analytically and
in closed forms. We examine and illustrate basic properties of the black-
hole—disk space-times.

Gravitation – Relativity – black-hole physics – accretion disks

††software: Wolfram Mathematica 13

##  1 Introduction

Black holes are once again in the spotlight of astrophysical interest thanks
to the unprecedented improvement of imaging techniques in the last couple of
years. The key aspect in imaging a black hole shadow is the presence of
accreting matter. Having nonvanishing angular momentum, such matter usually
forms an accretion disk inspiraling into the black hole. The usual
simplification of this model is neglecting the disk self-gravitation because
the real accretion disks are supposed to be much lighter than the central
black hole. However, in many situations, this is not necessarily true, and
many disk properties such as stability and positions of significant orbits are
sensible to the gravitational field. Thus, we should include also the disk’s
own gravity in our considerations.

Nevertheless, the complexity of Einstein’s equations makes the treatment in
the most general case almost impossible using just analytical methods. We have
to simplify the problem at some point. A great source of difficulties is the
presence of overall net rotation and the resulting dragging effects. But if
the rotation can be neglected or if it is compensated by e.g. two equal
counter-orbiting streams, the situation is much simpler. Namely, any static
and axially symmetric vacuum spacetime (or any non-empty spacetime where the
condition $T^{\rho}_{\rho}+T^{z}_{z}=0$ for the stress-energy tensor holds)
can always be described by a Weyl-type metric

$\operatorname{d\\!}s^{2}=-e^{2\nu}\operatorname{d\\!}t^{2}+\rho^{2}e^{-2\nu}\operatorname{d\\!}\phi^{2}+e^{2\lambda-2\nu}(\operatorname{d\\!}\rho^{2}+\operatorname{d\\!}z^{2})\;,$
(1)

where $t,\rho,\phi,z$ are Weyl coordinates and functions $\nu,\lambda$ depend
solely on $\rho,z$. Outside sources, the function $\nu$ satisfies 3
dimensional Laplace equation (therefore $\nu$ is the counterpart of Newtonian
gravitational potential), and for $\lambda$ we have

$\lambda_{,\rho}=\rho(\nu_{,\rho}^{2}-\nu_{,z}^{2})\;,\qquad\lambda_{,z}=2\rho\nu_{,\rho}\nu_{,z}\;.$
(2)

Just like in electrostatics and Newtonian gravity, in order to find the
external gravitational field in GR, we have to solve the Laplace (or Poisson)
equation. But that is only the (Newtonian) part of the whole story. The field
is not fully determined only by the potential – there is also the second
metric function $\lambda$ which influences the geometry in the meridional
$(\rho,z)$ plane. This feature can significantly deviate from the Newtonian
picture. Moreover, while the linearity of the Laplace equation means a
considerate simplification in obtaining the field generated by more than one
individual source (the potentials simply add up), the non-linear nature of
Einstein’s equations exhibits itself in the second metric function $\lambda$.

The first static superposition of a black hole and a thin disk (inverted
solution of Morgan & Morgan (1969)) was considered by Lemos & Letelier (1993)
and further studied in e.g. Semerák (2003). Later, Semerák (2004) found a
potential of a disk with a general power-law density profile in terms of
infinite series and studied its properties when superposed with a black hole.
In a recent revision of this topic in Kotlařík et al. (2022), we provided this
result in closed-form. All these disks are thin, infinite, and have an inner
rim, where higher derivatives of curvature are singular. Also, only the
Newtonian part (i.e. gravitational potential) of the superposition is given
analytically; the second metric function $\lambda$, if needed, has to be
solved numerically using (2). Here, we provide exact full relativistic
solutions describing superposition of a black hole with a disk which extends
from the horizon to infinity. The fields are everywhere regular, and both
metric functions of the superposition are obtained in closed-forms. It may
seem inappropriate to consider disks which reach down to the horizon, since
the stationary horizons cannot host any matter, and since quasi-stationary
accretion disks are assumed to end somewhere about the innermost stable
circular orbit. However, (i) our disk densities drop to zero towards the black
hole, so there is really no matter at the very horizon, and (ii) the accretion
disks around astrophysical black holes indeed continue towards the horizon (at
least partially), though their matter is infalling there rather than orbiting
in a quasi-circular regime. (Below the last stable orbit, the disk should not
be interpreted in terms of stationary circular motion.) Hence, it is in fact
more realistic to model the gravitational field of an accreting black hole
employing a disk that does have a certain modest density going down to the
horizon.

Before we proceed further, let us also mention other exact static black-
hole–disk superposition, where both non-trivial metric functions have been
obtained. González et al. (2009) found a family of infinite disks with an
inner edge suitable for the superposition with a black hole. However, Gleiser
(2012) showed that the disks have problematic physical interpretation. More
realistic result was obtained by Vieira (2020) using ‘displace, cut and
reflect’ method applied to solutions of $N$ black holes (singular rods in Weyl
coordinates) arranged in a linear chain on the symmetry axis.

The paper is organized as follows. First, in the Section 2, we review a family
of thin infinite (in extent) – but of finite mass – disks due to Kuzmin (1956)
and Toomre (1963), which were later studied in the context of GR by Bičák et
al. (1993). Then, by inversion (Kelvin transformation) of the potential with
respect to some Weyl radius, we find the potential describing, again, thin
infinite disks of finite mass, but in this case, their surface densities are
zero at $\rho=0$. The same potential has been found before, as a special case,
by Vogt & Letelier (2009) via superposition within the Kuzmin-Toomre family.
Here, we contribute to this result with the second metric function $\lambda$
and thus obtain full relativistic solutions describing inverted Kuzmin-Toomre
disks. Moreover, in Section 3, we also provide exact solutions of $\lambda$ to
all ‘single rings’ found in Vogt & Letelier (2009) – disks similar to the
inverted Kuzmin-Toomre ones, but their surface density falls faster. All these
disks are suitable for the superposition with a central black hole, thanks to
their ‘annular’ character. This is done in Sec. 4 where, again, both metric
functions of the black-hole–disk superposition are derived. In the subsequent
Sections 5 and 6 we study some properties of the resulting gravitational
field, namely the disk’s influence on the black-hole horizon, radial profiles
of densities, pressures and circular speeds in the equatorial plane. Finally,
in Section 7 we make several concluding remarks.

##  2 Kuzmin-Toomre disks and their inversion

The $n$-order Kuzmin-Toomre disk, $n\geq 0$, is described by the potential
(Kuzmin, 1956; Toomre, 1963; Evans & de Zeeuw, 1992; Bičák et al., 1993)

$\nu_{n}=-\frac{\mathcal{M}}{(2n-1)!!}\sum_{k=0}^{n}\frac{(2n-k)!}{2^{n-k}(n-k)!}\frac{b^{k}}{r_{b}^{k+1}}P_{k}\left(|\cos\theta_{b}|\right)\;,$
(3)

where we have denoted

$r_{b}^{2}\mathrel{\mathop{\mathchar
58\relax}}=\rho^{2}+(|z|+b)^{2}\;,\qquad|\cos\theta_{b}|\mathrel{\mathop{\mathchar
58\relax}}=\frac{|z|+b}{r_{b}}\;,$ (4)

$\mathcal{M}$ stands for the total mass of the disk, and $P_{k}$ are Legendre
polynomials. The corresponding Newtonian density222$\nu$ is the solution of
$\Delta\nu(\rho,z)=4\pi w(\rho)\delta(z)$, where $\Delta$ is standard Laplace
operator in cylindrical coordinates, $\delta$ is delta distribution. reads

$w_{n}=\frac{(2n+1)b^{2n+1}}{2\pi}\frac{\mathcal{M}}{(\rho^{2}+b^{2})^{n+3/2}}\;.$
(5)

Such a potential can be constructed by considering some mass distribution
along the negative half of the axis $z<0$, then cutting the solution along the
equatorial plane $z=0$ and reflecting its upper part $z>0$ to the negative
part $z<0$ of the axis – the so-called ‘displace, cut and reflect’ method
(Kuzmin, 1956; Vieira, 2020). The resulting field is then symmetric with
respect to the equatorial plane $z=0$ and given by

$\nu=-\int_{-\infty}^{\infty}\frac{W(b^{\prime})\operatorname{d\\!}b^{\prime}}{\sqrt{\rho^{2}+(|z|+b^{\prime})^{2})}}\;.$
(6)

For Kuzmin-Toomre solution, the weight function $W$ is a distribution given by
some combination of the derivatives of delta distribution (Bičák et al.,
1993). Note that in paper Bičák et al. (1993), the weight function given in
(2.23) does not in fact correspond to the Kuzmin-Toomre disks, but rather to
the solution obtained by inversion with respect to $\rho=b$, only with
different normalization – see the explicit treatment below.

Figure 1: Radial (Newtonian) surface density profiles of the $n$th order
Kuzmin-Toomre disks (left plot) and their inversions (right plot)
corresponding to $b=10\mathcal{M}$ in both plots. The horizontal axes are in
units of the disk masses $\mathcal{M}$, the vertical axes are in units of
$\mathcal{M}^{-2}$. In the left plot the maxima lie at
$\rho_{\text{max}}=\sqrt{2n/3}b$.

Inverted Kuzmin-Toomre disks are obtained by performing an inversion with
respect to333Performing the inversion around an arbitrary radius $a$ leads to
the same density profile just with different Weyl distance $b\longrightarrow
a^{2}/b$. $b$ (Kelvin transformation), i.e. the coordinate transformation

$\rho\longrightarrow\frac{b^{2}\rho}{\rho^{2}+z^{2}}\;,\qquad
z\longrightarrow\frac{b^{2}z}{\rho^{2}+z^{2}}\;,$ (7)

and the respective potential and density transformation

$\nu_{n}(\rho,z)\longrightarrow\frac{b}{\sqrt{\rho^{2}+z^{2}}}\nu_{n}\left(\frac{b^{2}\rho}{\rho^{2}+z^{2}},\frac{b^{2}z}{\rho^{2}+z^{2}}\right)\quad\Longrightarrow\quad
w_{n}(\rho)\longrightarrow\frac{b^{3}}{\rho^{3}}w_{n}(b^{2}/\rho)\;.$ (8)

Figure 2: Coordinate plots of the inverted $n$th order Kuzmin-Toomre disk
potentials in the equatorial plane (left plot) and on the axis (right plot)
corresponding to $b=3\mathcal{M}$. The horizontal axes are in units of the
disk mass $\mathcal{M}$ while the vertical axes are dimensionless.

The inversion of (3), after a suitable rearrangement of the sum, is again
separable in $r_{b}$ and $\cos\theta_{b}$ and can be easily derived from (3)
by (i) replacing $b\rightarrow-b$, (ii) multiplying each term by the
hypergeometric function
$\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}(1+k,k-n;k-2n;2)$,
(iii) fixing the normalization, so the total mass of the disk is still
$\mathcal{M}$. The inverted potential reads444Note that
$\frac{(2n-k)!}{2^{n-k}(n-k)!}\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}(1+k,k-n;k-2n;2)=\sum_{j=k}^{n}\binom{j}{k}\frac{(2n-j)!}{2^{n-j}(n-j)!}\;.$

$\nu_{n}^{\text{i}}=-\binom{n+1/2}{n}\frac{\mathcal{M}}{(1+2n)!!}\sum_{k=0}^{n}\frac{(2n-k)!}{2^{n-k}(n-k)!}\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}(1+k,k-n;k-2n;2)\frac{(-b)^{k}}{r_{b}^{k+1}}P_{k}\left(|\cos\theta_{b}|\right)\;,$
(9)

and the corresponding Newtonian density profile

$w_{n}^{\text{i}}=\binom{n+1/2}{n}\frac{\mathcal{M}b}{2\pi}\frac{\rho^{2n}}{(\rho^{2}+b^{2})^{n+3/2}}\;.$
(10)

Similarly as the original Kuzmin-Toomre disks (3), the $n$th-order potential
follows also from the recurrence relation

$\nu_{n+1}^{\text{i}}=\nu_{n}^{\text{i}}+\frac{b}{2(n+1)}\mathinner{\dfrac{\partial{}}{\partial{b}}}\nu_{n}^{\text{i}}\;,\qquad\nu^{\text{i}}_{0}=-\frac{\mathcal{M}}{r_{b}}\;.$
(11)

The shape of the disk potential along the equatorial plane and along the axis
is illustrated in Fig. 2. The Newtonian surface densities are depicted in Fig.
1.

The second metric function $\lambda_{n}$ can be also obtained following the
same steps as Bičák et al. (1993). They consist in rewriting the derivatives
in Eq. (2) in terms of $(r_{b},\theta_{b})$ and integrating from $r_{b}$ to
the infinity requiring $\lambda\rightarrow 0$ when $r_{b}\rightarrow\infty$.
Then

$\lambda_{n}^{\text{i}}=-\binom{n+1/2}{n}^{2}\frac{\mathcal{M}^{2}\sin^{2}\theta_{b}}{\left[(1+2n)!!\right]^{2}}\sum_{k,l=0}^{n}\mathcal{B}_{k,l}\frac{(-b)^{k+l}}{r_{b}^{k+l+2}}\mathcal{P}_{k,l}(\theta_{b})\;,$
(12)

where

$\displaystyle\mathcal{B}_{k,l}$
$\displaystyle\equiv\frac{(2n-k)!(2n-l)!}{2^{2n-k-l}(n-k)!(n-l)!(k+l+2)}\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}(1+k,k-n;k-2n;2)\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}(1+l,l-n;l-2n;2)\;,$
(13) $\displaystyle\mathcal{P}_{k,l}$
$\displaystyle\equiv(k+1)(l+1)P_{k}P_{l}+2(k+1)|\cos\theta_{b}|P_{k}P^{\prime}_{l}-\sin^{2}\theta_{b}P_{k}^{\prime}P_{l}^{\prime}\;,$
(14) $\displaystyle P^{\prime}_{k}$
$\displaystyle\equiv\mathinner{\dfrac{\operatorname{d\\!}{}}{\operatorname{d\\!}{|\cos\theta_{b}|}}}P_{k}(|\cos\theta_{b}|)\;.$
(15)

##  3 Vogt-Letelier disks

Figure 3: Newtonian surface densities of Vogt-Letelier disks for $n=1$ and
$b=10\mathcal{M}$. Maxima lie at $\rho_{\text{max}}=2b\sqrt{\frac{n}{m+3/2}}$.
The horizontal axis is in units of $\mathcal{M}$ while the vertical axis is in
units of $\mathcal{M}^{-2}$. Figure 4: Coordinate plots of the Vogt-Letelier
disks’ potential ($n=1$) in the equatorial plane (left plot) and on the axis
(right plot). In both plots, $b=10\mathcal{M}$, the horizontal axes are in
units of $\mathcal{M}$, while the vertical axes are dimensionless.

Vogt & Letelier (2009) obtained a broader family of disks by a superposition
within the original Kuzmin-Toomre disks. They took

$\nu^{(m,n)}=W^{(m,n)}\sum_{k=0}^{n}(-1)^{k}\binom{n}{k}\frac{\nu_{m+k}}{2m+2k+1}\;,$
(16)

where $W^{(m,n)}$ is a normalization factor dependent only on $m,n$. Resulting
disks are again of infinite extent with the surface density

$w^{(m,n)}=W^{(m,n)}\frac{\mathcal{M}b^{2m+1}}{2\pi}\frac{\rho^{2n}}{(\rho^{2}+b^{2})^{m+n+3/2}}\;,$
(17)

which behaves as $\mathcal{O}(\rho^{-2m-3})$ at infinity, thus the total mass
of each disk is again finite. We fix $W^{(m,n)}$, so the total mass of the
disk is still $\mathcal{M}$, more precisely

$W^{(m,n)}=(2m+1)\binom{m+n+1/2}{n}\Rightarrow 2\pi\int_{0}^{\infty}\rho
w^{(m,n)}\operatorname{d\\!}\rho=W^{(m,n)}\frac{\mathcal{M}b^{2m+1}}{2\pi}\int_{0}^{\infty}\frac{\rho^{2n}}{(\rho^{2}+b^{2})^{m+n+3/2}}\operatorname{d\\!}\rho=\mathcal{M}\;.$
(18)

For $m=0$, we get the inverted Kuzmin-Toomre disks, which we diskussed in the
previous section. disks of higher $m$ can be obtained using a recurrence
relation

$\frac{(2m+1)(2n+3)}{2m+2n+3}\nu^{(m+1,n)}=\nu^{(m,n)}+\frac{4m(n+1)}{2m+2n+3}\nu^{(m,n+1)}-b\mathinner{\dfrac{\partial{}}{\partial{b}}}\nu^{(m,n)}\;.$
(19)

After substitution of the Kuzmin-Toomre potential into (16) the double sum can
be rearranged and performed over the $k$ index. It leads to

$\nu^{(m,n)}=-(2m+1)\mathcal{M}\binom{m+n+1/2}{n}\sum_{j=0}^{m+n}\mathcal{Q}_{j}^{(m,n)}\frac{b^{j}}{r_{b}^{j+1}}P_{j}(\cos\theta_{b})\;,$
(20)

where we have already substituted $W^{(m,n)}$ from (18) and

$\mathcal{Q}_{j}^{(m,n)}=\begin{cases}\sum_{k=0}^{n}(-1)^{k}\binom{n}{k}\frac{2^{j-k-m}(2m+2k-j)!}{(m+k-j)!(2m+2k+1)!!}=\frac{2^{j-m}(2m-j)!}{(2m+1)!!(m-j)!}\operatorname{{\vphantom{\mathnormal{F}}}_{3}\kern-0.5pt\mathnormal{F}_{2}}\left(\frac{2m-j+1}{2},\frac{2m-j+2}{2},-n;\frac{2m+3}{2},m-j+1;1\right)&\text{
if }j\leq m\\\
\sum_{k=j}^{m+n}(-1)^{k-m}\binom{n}{k-m}\frac{2^{j-k}(2k-j)!}{(k-j)!(2k+1)!!}=\frac{(-1)^{j-m}j!}{(2j+1)!!}\binom{n}{j-m}\operatorname{{\vphantom{\mathnormal{F}}}_{3}\kern-0.5pt\mathnormal{F}_{2}}\left(\frac{j+1}{2},\frac{j+2}{2},j-m-n;\frac{2j+3}{2},j-m+1;1\right)&\text{
if }j>m\;,\end{cases}$

where
$\operatorname{{\vphantom{\mathnormal{F}}}_{3}\kern-0.5pt\mathnormal{F}_{2}}$
is the generalized hypergeometric function.

The second metric function $\lambda$ is again of the same structure (12) as
the original or inverted Kuzmin-Toomre disks, but now

$\displaystyle\lambda^{(m,n)}=-(2m+1)^{2}\binom{m+n+1/2}{n}^{2}\mathcal{M}^{2}\sin^{2}\theta_{b}\sum_{k,l=0}^{m+n}\mathcal{B}^{(m,n)}_{k,l}\frac{b^{k+l}}{r_{b}^{k+l+2}}\mathcal{P}_{k,l}(\theta_{b})\;,$
(21)
$\displaystyle\mathcal{B}^{(m,n)}_{k.l}=\frac{\mathcal{Q}_{l}^{(m,n)}\mathcal{Q}_{k}^{(m,n)}}{k+l+2}\;,$
(22)

while the polynomials $\mathcal{P}_{k,l}$ are the same as in (12). A shorter
expression can be found by the off-diagonal re-summation

$\lambda^{(m,n)}=-(2m+1)^{2}\binom{m+n+1/2}{n}^{2}\mathcal{M}^{2}\sum_{v=0}^{2(m+n)}\sum_{u=0}^{v}\frac{(1+u)(1+v-u)}{v+2}\frac{b^{v}}{r_{b}^{-2-v}}\mathcal{Q}_{u}^{(m,n)}\mathcal{Q}_{v-u}^{(m,n)}(P_{u}P_{v-u}-P_{1+u}P_{1+v-u})\;,$
(23)

where $P_{u}\equiv P_{u}(|\cos\theta_{b}|)$ are the Legendre polynomials and
we set $\mathcal{Q}_{u}^{(m,n)}=0$ for $u>m+n$. For the inverted Kuzmin-Toomre
disks (i.e. $m=0$) the coefficients are reduced
$\mathcal{Q}_{u}^{(0,n)}=\frac{(-1)^{u}}{(2n+1)!!}\frac{(2n-u)!}{2^{n-u}(n-u)!}\operatorname{{\vphantom{\mathnormal{F}}}_{2}\kern-0.5pt\mathnormal{F}_{1}}(1+u,u-n;u-2n;2)$.

Newtonian density profiles of several Vogt-Letelier disks are plotted in Fig.
3 and the behavior of disks’ potentials in the equatorial plane and on the
axis is illustrated in Fig. 4.

##  4 Superposition with a Schwarzschild black hole

All Vogt-Letelier disks have a clear ‘annular’ character – their densities are
exactly zero in the centre $(\rho=0,z=0)$. Hence, it is physically reasonable
to make a static superposition of such a disk with a black hole. In Weyl
coordinates, the Schwarzschild black hole of the mass $M$ is a singular rod
placed on the axis $(\rho=0,|z|\leq M)$ producing gravitational field
described by

$\displaystyle\nu_{\text{Schw}}=\frac{1}{2}\ln\frac{d_{1}+d_{2}-2M}{d_{1}+d_{2}+2M}\;,\qquad\lambda_{\text{Schw}}=\frac{1}{2}\ln\frac{(d_{1}+d_{2})^{2}-4M^{2}}{4d_{1}d_{2}}\;,$
(24)

where $d_{1,2}\equiv\sqrt{\rho^{2}+(|z|\mp M)^{2}}$. Then, the disk starts at
the black-hole horizon, but the density drops to zero sufficiently quickly
there. Due to the linearity of Laplace equation, the superposition of
gravitational potentials is a simple sum of the individual sources, i.e.
$\nu=\nu_{\text{Schw}}+\nu_{\text{disk}}$. The non-linearity of the Einstein
equations manifests itself in the second metric function $\lambda$ which does
not superpose that simply. In fact, for two individual sources, a black hole
and a disk, we can write
$\lambda=\lambda_{\text{Schw}}+\lambda_{\text{disk}}+\lambda_{\text{int}}$,
where $\lambda_{\text{Schw}}$ and $\lambda_{\text{disk}}$ are contributions
from the black hole and the disk alone (thus satisfying (2) with just
$\nu_{\text{Schw}}$ and $\nu_{\text{disk}}$ respectively) and
$\lambda_{\text{int}}$ is an ‘interaction’ term for which

$\displaystyle\lambda_{\text{int},\rho}$
$\displaystyle=2\rho(\nu_{\text{Schw},\rho}\nu_{\text{disk},\rho}-\nu_{\text{Schw},z}\nu_{\text{disk},z})\;,$
(25) $\displaystyle\lambda_{\text{int},z}$
$\displaystyle=2\rho(\nu_{\text{Schw},\rho}\nu_{\text{disk},z}+\nu_{\text{Schw},z}\nu_{\text{disk},\rho})\;.$
(26)

Notice that these conditions are already linear in $\nu_{\text{disk}}$, i.e.
if the potential of the disk is a sum of two components then also $\lambda$
must add in the same way. Hence, $\lambda_{\text{int}}$ satisfies the same
recurrence relations (11) for the inverted Kuzmin-Toomre disks and (19) for
Vogt-Letelier disks of higher $m$, namely

$\displaystyle\lambda^{(0,n+1)}_{\text{int}}=\lambda^{(0,n)}_{\text{int}}+\frac{b}{2(n+1)}\mathinner{\dfrac{\partial{}}{\partial{b}}}\lambda^{(0,n)}_{\text{int}}\;,\qquad\lambda^{(0,0)}_{\text{int}}=-\frac{\mathcal{M}}{r_{b}}\left(\frac{d_{1}}{b+M}-\frac{d_{2}}{b-M}\right)-\frac{2\mathcal{M}M}{b^{2}-M^{2}}\;,$
(27)
$\displaystyle\frac{(2m+1)(2n+3)}{2m+2n+3}\lambda_{\text{int}}^{(m+1,n)}=\lambda_{\text{int}}^{(m,n)}+\frac{4m(n+1)}{2m+2n+3}\lambda_{\text{int}}^{(m,n+1)}-b\mathinner{\dfrac{\partial{}}{\partial{b}}}\lambda_{\text{int}}^{(m,n)}\;.$
(28)

Superposition with a black hole is best represented in Schwarzschild
coordinates555Note that Schwarzschild coordinates $(r,\theta)$ are different
from those coordinates with subscript $b$ used in disks’ potentials.
$(r,\theta)$ introduced by

$\rho=\sqrt{r(r-2M)}\sin\theta\;,\qquad z=(r-M)\cos\theta\;.$ (29)

In these coordinates the black-hole horizon is rendered really spherical (at
$r=2M$), the Schwarszchild contribution becomes

$\nu_{\text{Schw}}=\frac{1}{2}\ln\left(1-\frac{2M}{r}\right)\;,\qquad\lambda_{\text{Schw}}=\frac{1}{2}\ln\frac{r^{2}-2Mr}{r^{2}-2Mr+M^{2}\sin^{2}\theta}\;,$
(30)

and the metric of the superposition reads

$\operatorname{d\\!}s^{2}=-\left(1-\frac{2M}{r}\right)e^{2\nu_{\text{disk}}}\operatorname{d\\!}t^{2}+e^{2\lambda_{\text{ext}}-2\nu_{\text{disk}}}\frac{\operatorname{d\\!}r^{2}}{1-\frac{2M}{r}}+r^{2}e^{-2\nu_{\text{disk}}}(e^{2\lambda_{\text{ext}}}\operatorname{d\\!}\theta^{2}+\sin^{2}\theta\operatorname{d\\!}\phi^{2})\;,$
(31)

where
$\lambda_{\text{ext}}\equiv\lambda-\lambda_{\text{Schw}}=\lambda_{\text{disk}}+\lambda_{\text{int}}$.
The resulting field is illustrated in Fig. 6 on several contour plots. In Fig.
5, we show radial plots of $\lambda_{\text{ext}}$ in the equatorial plane.

Figure 5: Equatorial plots (in Weyl radius) of
$\lambda_{\text{ext}}=\lambda-\lambda_{\text{Schw}}$ of the superposition of a
black hole and a disk (inverted Kuzmin-Toomre disks in the left plot; Vogt-
Letelier disks in the right plot). Parameters of the disks are
$\mathcal{M}=M$, $b=3M$ in the left plot, and $\mathcal{M}=3.8M$, $b=10M$ in
the right plot. The horizontal axes are in units of $M$, while the vertical
axes are dimensionless. In the right plot, the dashed black lines correspond
to $m=1,2,3,4$, while the solid red lines correspond to $m=5,6,7,8$. The
(color) line differentiation is chosen only for the sake of clarity. Figure 6:
Meridional plane ($\phi=\text{const}$) contour plots of the superposition of
the black hole and the inverted 3rd order Kuzmin-Toomre disk. Namely, the
gravitational potential $\nu=\nu_{\text{Schw}}+\nu_{\text{disk}}$ (top plot),
the second metric function
$\lambda_{\text{ext}}=\lambda-\lambda_{\text{Schw}}$ (middle plot) and the
difference $\lambda_{\text{ext}}-\nu_{\text{disk}}$ (bottom plot), are shown
in Schwarzschild coordinates. The disk lies in the equatorial plane $z=0$ and
is highlighted by a thick black line. The thick blue line shows the Newtonian
density profile of the disk; the white circle in the origin shows the black
hole. To illustrate how the field of the black hole is distorted by the
presence of the disk, rather extreme parameters are chosen in all plots – the
mass of the disk $\mathcal{M}=10M$ and $b=8M$. Both axes are in units of the
black hole mass $M$.

##  5 Isometric embedding of the horizon

Even though the horizon of the black hole always lies at the Schwarzschild
radius $r=2M$ and remains ‘spherical’ in the Schwarszchild coordinates, its
intrinsic geometry changes due to the presence of the surrounding matter
distribution. At any given (coordinate) time, the horizon is a 2D-surface
$(t=\text{konst},r=2M)$ with an induced metric

$\operatorname{d\\!}s^{2}_{\text{H}}=4M^{2}e^{-2\nu^{\text{H}}_{\text{disk}}(\theta)}\left[e^{4\nu^{\text{H}}_{\text{disk}}(\theta)-4\nu^{\text{H}}_{\text{disk}}(\theta=0)}\operatorname{d\\!}\theta^{2}+\sin^{2}\theta\operatorname{d\\!}\phi^{2}\right]\;,$
(32)

where $\nu^{\text{H}}_{\text{disk}}$ is the limit of the disk potential at the
horizon ($r\rightarrow 2M$), and where we have already used solution
$\lambda^{\text{H}}(\theta)=2\nu^{\text{H}}(\theta)-2\nu^{\text{H}}(\theta=0)$
valid on the horizon (for any static or even stationary axially symmetric
spacetime). For isometric embedding to Euclidean 3-space, we use the method by
Smarr (1973). It consists in rewriting the 2-metric in terms of the coordinate
$\mu=\cos\theta$, i.e.

$\operatorname{d\\!}s^{2}=\eta^{2}\left[f^{-1}(\mu)\operatorname{d\\!}\mu^{2}+f(\mu)\operatorname{d\\!}\phi^{2}\right]\;,$
(33)

where

$\eta=2Me^{-\nu_{\text{disk}}(\mu=1)}\;,\qquad
f=e^{2\nu_{\text{disk}}(\mu=1)-2\nu_{\text{disk}}(\mu)}(1-\mu^{2})\;.$ (34)

Then, the isometric embedding of the horizon 2-surface in three-dimensional
Euclidean space $(x,y,z)$ is given by

$x=\eta\sqrt{f}\cos\phi\;,\qquad y=\eta\sqrt{f}\sin\phi\;,\qquad
z=\eta\int_{0}^{\mu}\sqrt{\frac{1}{f}\left(1-\frac{1}{4}f_{,\mu}^{2}\right)}\operatorname{d\\!}\mu\;.$
(35)

See Fig. 7 for the numerical results. Another useful quantity for any 2D
surface is its Gauss curvature $\mathcal{K}\equiv\frac{{}^{(2)}R}{2}$, where
${}^{(2)}R$ is the corresponding 2D scalar curvature. For the metric (32), we
have

$\mathcal{K}(\theta)=\mathinner{\frac{1+3\nu_{\text{disk},\theta}\cot\theta+\nu_{\text{disk},\theta\theta}-2\nu_{\text{disk},\theta}^{2}}{4M^{2}e^{2\nu_{\text{disk}}(\theta)-4\nu_{\text{disk}}(0)}}\Biggr{\rvert}}_{r=2M}\;.$
(36)

The Gauss curvature is plotted against $\mathcal{M}$ and $b$ in Fig. 8. When
the Gauss curvature turns negative, the horizon is not globally embeddable
into Euclidean 3-space. Both, the embedding and the Gauss curvature, shows
flattened horizon in the direction of the disk.

Figure 7: Meridional section ($\phi=$ const) of the isometric embedding of the
black-hole horizon to Euclidean 3-space. The black hole is distorted due to
the inverted 3th order Kuzmin-Toomre disk with $b=2M$. The masses of the disk
range as $\mathcal{M}=5M,7.5M,10M,\dots 25M$. The horizon becomes more and
more flattened with increasing $\mathcal{M}$. For the last two cases, the
horizon is not globally embeddable into Euclidean 3-space. Both axes are in
units of $M$. Figure 8: Gauss curvature at the horizon on the axis is plotted
against the disk mass $\mathcal{M}$ (left plot) and $b$ (right plot). The
disks are inverted Kuzmin-Toomre of the 3rd order. In the left plot $b=0.75M$
is chosen, while in the right plot $\mathcal{M}=3M$ is chosen. Gauss curvature
turns negative if the disk is too massive or the density maximum is
sufficiently close to the black hole. Both horizontal axes are in units of $M$
while the vertical axes are in units of $M^{-2}$.

##  6 Physical properties of disks

Two rather simple physical interpretations of any static thin disks are (i) a
single component ideal fluid with a certain surface density $\sigma$ and an
azimuthal pressure $P$, or (ii) two identical counter-orbiting dust streams
with proper surface densities $\sigma_{+}=\sigma_{-}\equiv\frac{\sigma}{2}$
following circular geodesics. Both characteristics follow from the
diskontinuities of normal derivatives of the field over the equatorial plane

$\sigma+P=\frac{\nu_{,z}(z=0^{+})}{2\pi}=w(\rho)\;,\qquad
P=\frac{\lambda_{,z}(z=0^{+})}{4\pi}=\frac{\nu_{,z}(z=0^{+})}{2\pi}\rho\nu_{,\rho}=w(\rho)\rho\nu_{,\rho}\;,$
(37)

where $w(\rho)$ is the Newtonian density (10), or (17). The second
interpretation is possible only when $\sigma,P\geq 0$.

###  6.1 Radial profiles of densities and azimuthal pressures

Figure 9: Profiles of disk densities (top row), azimuthal pressures (middle
row) and circular velocities (bottom row) in terms of the radial coordinate
$\rho$ (left column), circumferential radius $r_{\text{cf}}$ (middle column)
and proper distance from the horizon $r_{\text{prop}}$ (right column). Eight
members of the inverted Kuzmin-Toomre family ($n=1,2,\dots,8$) are depicted in
each plot. We chose the mass of the disks $\mathcal{M}=M$ and $b=10M$, where
$M$ is the mass of the central black hole. The horizontal axes are in units of
$M$ while the densities and pressures are in units $M^{-1}$ and the velocities
in fractions of the speed of light.

Density and pressure (37) are not, yet, quantities measured by any physical
observer. Actually, an observer at rest with respect to the disk would measure
the density $e^{\nu-\lambda}\sigma$ and the pressure $e^{\nu-\lambda}P$ (Bičák
et al., 1993). We show radial profiles of $e^{\nu-\lambda}\sigma$ and
$e^{\nu-\lambda}P$ of the inverted Kuzmin-Toomre disks in Figs. 9, 10, and of
Vogt-Letelier disks ($n=1$) in Fig. 11, 12. The disks parameters in Fig. 12
are rather extreme (especially for higher orders in $m$, because the mass of
the disk is kept constant), but it shows the role of the multiplication factor
$e^{\nu-\lambda}$, where also the second metric function $\lambda$ is present.
However, the circular velocities in the disk plane (see the section below)
would be superluminal in some regions, so the double-stream interpretation of
such disks would not be possible.

###  6.2 Circular-velocity profiles

Another useful quantity is a physical speed $v$ of circular geodesics in the
equatorial plane measured locally by a static observer. Such a speed is given
by

$v^{2}\equiv\frac{P}{\sigma}=\frac{\rho\nu_{,\rho}}{1-\rho\nu_{,\rho}}\;,$
(38)

and its time-like condition $0\leq v^{2}<1$ covers both physical requirements
for the disk – the energy conditions and non-negativity of azimuthal pressure.
In Figs. 9, 10, 11, 13 we show velocity profiles for the inverted Kuzmin-
Toomre disks and Vogt-Letelier disks ($n=1$).

###  6.3 Coordinate and geometrical measures

Most of the statements about the spacetime is given in coordinate terms. Such
statements have to be taken with some caution, although the Weyl (or
Schwarzschild) coordinates represents some spacetime features adequately. In
particular, we should also check the physical properties of disks by employing
invariant measures like circumferential radius or proper radial distance. In
our case, the proper circumference corresponding to a certain $\rho$ (computed
along constant $t,\rho$ and $z$) reads

$\int_{0}^{2\pi}\sqrt{g_{\phi\phi}}\operatorname{d\\!}\phi=2\pi\sqrt{g_{\phi\phi}}=2\pi\rho
e^{-\nu}\qquad\Longrightarrow\qquad r_{\text{cf}}\mathrel{\mathop{\mathchar
58\relax}}=\rho e^{-\nu}\;,$ (39)

where we have denoted the circumferential radius $r_{\text{cf}}$ in such way
that the corresponding circumference is given as $2\pi r_{\text{cf}}$. The
proper radial distance from the black-hole horizon to a certain $\rho$,
calculated in the equatorial plane $z=0$ along the remaining coordinates
$(t,\phi)$ constant, is given by

$r_{\text{prop}}\mathrel{\mathop{\mathchar
58\relax}}=\int_{0}^{\rho}\sqrt{g_{\rho\rho}}\operatorname{d\\!}\rho=\int_{0}^{\rho}e^{\lambda-\nu}\operatorname{d\\!}\rho\;.$
(40)

The latter integral has to be solved numerically. All quantities in Figs. 9–13
are depicted against the coordinate $\rho$ as well as both geometrical
measures defined above. Note that the proper circumference is not zero at
$(\rho=0,z=0)$, instead it is the proper equatorial circumference of the
black-hole horizon. As illustrated in Sec. 5, this circumference must grow due
to the presence of the disk. Namely, the corresponding circumferential radius
(for any black-hole–Vogt-Letelier–disk spacetime) there reads

$r_{\text{cf}}(\rho=0,z=0)=\lim_{\rho\rightarrow 0}\rho
e^{-\nu(\rho,z=0)}=2M\exp\left(\frac{(2m+1)\mathcal{M}}{b}\binom{m+n+1/2}{n}\sum_{j=0}^{m+n}\mathcal{Q}_{j}^{(m,n)}\right)\;.$
(41)

Figure 10: Same profiles as in Fig. 9. The disks belong to the inverted
Kuzmin-Toomre family; the disk mass is chosen to be $\mathcal{M}=3M$, and
$b=10M$. Figure 11: Same profiles as in Fig. 9, but in this figure, eight
members of Vogt-Letelier disks $(m=1,2,\dots,8,n=1)$ are depicted in each
plot. The disks are of the same mass and $b$ as in Fig. 9, namely
$\mathcal{M}=M$ and $b=10M$. Figure 12: Profiles of disk densities (top row)
and azimuthal pressures (bottom row) of the Vogt-Letelier disks
$(m=1,2,\dots,8,n=1)$ in terms of radial distances, as in Fig. 9. The
rotational velocities are in separate figure, because rather extreme
parameters have been chosen, namely $\mathcal{M}=3.8M$ and $b=10M$. Meanwhile,
the densities and pressures are well behaving, the rotational velocities would
be superluminal in some regions of the disks, therefore the double-stream
interpretation would not be possible. Figure 13: Profiles of circular
velocities of the same disks as in Fig. 12, just with different disks mass
$\mathcal{M}=1.5M$.

##  7 Concluding remarks

We have found full metric describing static and axially symmetric
superposition of a black hole and a thin disk. The disks we considered have
finite mass, have no radial pressure and extend from the black-hole horizon to
infinity. They also have no edges, so the field is regular everywhere outside
the horizon, yet the disks surface densities fall off sufficiently quickly at
both extremes $\rho=0$ and $\rho\rightarrow\infty$. In particular, we have
used disks resulting from the inversion (Kelvin transformation) of Kuzmin-
Toomre solutions, as well as ‘simple’ disks due to Vogt & Letelier (2009) (of
which the inverted Kuzmin-Toomre disks are a distinct subclass). Both metric
functions are given analytically and in closed-forms. The gravitational field
is illustrated in plots of the metric functions and profiles of densities,
azimuthal pressures and circular velocities. The black-hole horizon was found
distorted due to the presence of the disk (flattened in the direction of the
disk) which is in agreement with many previous results, e.g. Semerák (2004,
2003); Semerák et al. (2001), or, in the stationary case, Kotlařík et al.
(2018).

We thank for support from GACR 21-11268S grant of the Czech Science Foundation
(D.K.) and from Grant schemes at Charles University, reg. n.
CZ.02.2.69/0.0/0.0/19_073/0016935 (P.K). We are also very grateful to Oldřich
Semerák for the initial idea of inverting Kuzmin-Toomre disks and for helpful
comments which improved the manuscript.

## References

  * Bičák et al. (1993) Bičák, J., Lynden-Bell, D., & Pichon, C. 1993, MNRAS, 265, 126, doi: 10.1093/mnras/265.1.126
  * Evans & de Zeeuw (1992) Evans, N. W., & de Zeeuw, P. T. 1992, MNRAS, 257, 152, doi: 10.1093/mnras/257.1.152
  * Gleiser (2012) Gleiser, R. J. 2012, Phys. Rev. D, 85, 028501, doi: 10.1103/PhysRevD.85.028501
  * González et al. (2009) González, G. A., Gutiérrez-Piñeres, A. C., & Viña-Cervantes, V. M. 2009, Phys. Rev. D, 79, 124048, doi: 10.1103/PhysRevD.79.124048
  * Kotlařík et al. (2022) Kotlařík, P., Kofroň, D., & Semerák, O. 2022, ApJ, 931, 161, doi: 10.3847/1538-4357/ac6027
  * Kotlařík et al. (2018) Kotlařík, P., Semerák, O., & Čížek, P. 2018, Phys. Rev. D, 97, doi: 10.1103/PhysRevD.97.084006
  * Kuzmin (1956) Kuzmin, G. G. 1956, Astr. Zh., 33, 27
  * Lemos & Letelier (1993) Lemos, J. P. S., & Letelier, P. S. 1993, CQGra, 10, L75, doi: 10.1088/0264-9381/10/6/003
  * Morgan & Morgan (1969) Morgan, T., & Morgan, L. 1969, Phys. Rev., 183, 1097, doi: 10.1103/PhysRev.183.1097
  * Semerák (2003) Semerák, O. 2003, CQGra, 20, 1613, doi: 10.1088/0264-9381/20/9/302
  * Semerák (2004) Semerák, O. 2004, CQGra, 21, 2203, doi: 10.1088/0264-9381/21/8/019
  * Semerák et al. (2001) Semerák, O., Zellerin, T., & Žáček, M. 2001, MNRAS, 322, 207, doi: 10.1046/j.1365-8711.2001.03814.x
  * Smarr (1973) Smarr, L. 1973, Phys. Rev. D, 7, 289, doi: 10.1103/PhysRevD.7.289
  * Toomre (1963) Toomre, A. 1963, ApJ, 138, 385, doi: 10.1086/147653
  * Vieira (2020) Vieira, R. S. 2020, CQGra, 37, doi: 10.1088/1361-6382/aba99b
  * Vogt & Letelier (2009) Vogt, D., & Letelier, P. S. 2009, MNRAS, 396, 1487, doi: 10.1111/j.1365-2966.2009.14803.x


# Mutually unbiased bases as a commuting polynomial optimisation problem

Luke Mortimer  ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute
of Science and Technology, 08860 Castelldefels, Spain

(August 28, 2024)

###### Abstract

We consider the problem of mutually unbiased bases as a polynomial
optimization problem over the reals. We heavily reduce it using known
symmetries before exploring it using two methods, combining a number of
optimization techniques. The first of these is a search for bases using
Lagrange-multipliers that converges rapidly in case of MUB existence, whilst
the second combines a hierarchy of semidefinite programs with branch-and-bound
techniques to perform a global search. We demonstrate that such an algorithm
would eventually solve the open question regarding dimension 6 with finite
memory, although it still remains intractable. We explore the idea that to
show the inexistence of bases, it suffices to search for orthonormal vector
sets of certain smaller sizes, rather than full bases. We use our two methods
to conjecture the minimum set sizes required to show infeasibility, proving it
for dimensions 3. The fact that such sub-problems seem to also be infeasible
heavily reduces the number of variables, by 66% in the case of the open
problem, potentially providing an large speedup for other algorithms and
bringing them into the realm of tractability.

##  I Introduction to mutually unbiased bases

A set of vectors form a complete orthonormal basis if each has a magnitude of
one, have pairwise inner products equal to zero, and completely span the space
in which they are defined. For example, the following is a complete
orthonormal basis in $\mathbb{C}^{2}$, generally referred to as the
computational basis:

$\displaystyle v_{11}=(1\leavevmode\nobreak\ \leavevmode\nobreak\ 0)\qquad
v_{12}=(0\leavevmode\nobreak\ \leavevmode\nobreak\ 1)$ (1)

Two orthonormal bases are considered mutually unbiased if the inner products
between their vectors each has a magnitude of $\frac{1}{\sqrt{d}}$. This
extends to sets of more than just two, such that each basis is mutually
unbiased with respect to every other basis, for example the following set of 3
mutually unbiased bases (MUBs) in dimension 2:

$\displaystyle v_{11}$ $\displaystyle=(1\leavevmode\nobreak\
\leavevmode\nobreak\ 0)\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
v_{12}=(0\leavevmode\nobreak\ \leavevmode\nobreak\ 1)$ (2) $\displaystyle
v_{21}$ $\displaystyle=\frac{1}{\sqrt{2}}(1\leavevmode\nobreak\
\leavevmode\nobreak\ 1)\qquad v_{22}=\frac{1}{\sqrt{2}}(1\leavevmode\nobreak\
\leavevmode\nobreak\ -1)$ (3) $\displaystyle v_{31}$
$\displaystyle=\frac{1}{\sqrt{2}}(1\leavevmode\nobreak\ \leavevmode\nobreak\
i)\qquad v_{32}=\frac{1}{\sqrt{2}}(1\leavevmode\nobreak\ \leavevmode\nobreak\
-i)$ (4)

Such bases have significant use in the field of quantum information theory,
providing measurements maximising Bell inequalities tabia2022bell ;
tavakoli2021mutually , their use in quantum state tomography rebon2013phase ;
yuan2016quantum , or their use in secure quantum communication protocols
mafu2013higher . This is due to the property that preparing in one basis and
measuring in a mutually unbiased basis results in a uniform probability of any
outcome, making it difficult for an attacker to gain information.

In any dimension $d$ it is known that there are at most $d+1$ MUBs. More
specifically, in a given dimension $d$ that is the power of a prime it is
known that one can always find the maximal set of $d+1$ MUBs
wootters1987wigner . It remains a decades-old open question exactly how many
MUBs exist in dimensions that are not the power of a prime, the lowest being
dimension 6. There exists a conjecture, known as Zauner’s conjecture, stating
that one can find at most 3 MUBs in dimension 6 zauner1999grundzuge . This
claim is supported by numerical searches which have failed to find 4 MUBs in
dimension 6 despite extensive effort colomer2022three ; raynal2011mutually .

##  II The MUB problem as polynomial optimisation

One of the most well-studied forms of problem is that of polynomial
optimisation, which in the general form consists of trying to minimize a real
polynomial $f(\vec{x})$ over a vector $\vec{x}$, subject to equality
constraints $g(\vec{x})=0$ and positivity constraints $h(\vec{x})>0$. These
polynomials can (and in this work will be) non-convex and non-linear, thus
allowing all problems in NP to be expressed jacquet2020reformulation . If the
polynomials were all convex, this type of problem could be solved easily using
convex optimisation techniques, many of which are known to converge to finite
precision with polynomial time complexity karmarkar1984new ;
laurent2005semidefinite . The idea of converting the MUB existence problem
into a polynomial optimisation problem has been studied before in the non-
commuting complex case gribling2021mutually  or in the case of studying MUB
constellations brierley2010mutually .

Here an objective function $f(\vec{x})$ is unnecessary, since we are only
trying to find any valid set of MUBs, rather the set of MUBs which minimizes
some function. The same is true for the positivity constraints, since in the
MUB problem we only need to constrain that various magnitudes of inner
products equal either 0, 1, or $\frac{1}{\sqrt{d}}$. So, for the problem of
finding $n$ MUBs in dimension $d$, we begin with the problem of finding
complex vectors $\vec{v}_{ij}$, where here the first index indicates the basis
(from $1\to n$) and the second indicates the vector within the basis (from
$1\to d$):

$\displaystyle\vec{v}_{ii}^{\dagger}\vec{v}_{ii}$ $\displaystyle=1$ (5)
$\displaystyle\vec{v}_{ik}^{\dagger}\vec{v}_{il}$ $\displaystyle=0$ (6)
$\displaystyle|\vec{v}_{ik}^{\dagger}\vec{v}_{jl}|$
$\displaystyle=\frac{1}{\sqrt{d}}$ (7) $\displaystyle i\neq j$
$\displaystyle\qquad k\neq l$ (8)

In order to simplify Equation 7 we introduce two new real variables
$c_{t},d_{t}$ for each of the pairs, constrained to the circle of magnitude
$\frac{1}{\sqrt{d}}$. This allows us to keep the equations second-order,
versus removing the norm by squaring the constraint which would allow for less
variables with the trade-off of a fourth-order equation set. Testing both
version we found it is significantly better regarding both the memory and the
time required for both of our algorithms if the order is kept as low as
possible. We also replace each element of each vector with a combination of
the real and complex parts, such that element $m$ of $\vec{v}_{ij}$ becomes
$v_{iim}=a_{iim}+ib_{iim}$, giving the following polynomial equation set over
the reals:

$\displaystyle\sum_{m=1}^{d}(a_{iim}-ib_{iim})(a_{iim}+ib_{iim})$
$\displaystyle=1$ (9)
$\displaystyle\sum_{m=1}^{d}(a_{ikm}-ib_{ikm})(a_{ilm}+ib_{ilm})$
$\displaystyle=0$ (10)
$\displaystyle\sum_{m=1}^{d}(a_{ikm}-ib_{ikm})(a_{jlm}+ib_{jlm})$
$\displaystyle=c_{i+dj}+id_{i+dj}$ (11) $\displaystyle
c_{i+dj}^{2}+d_{i+dj}^{2}$ $\displaystyle=\frac{1}{d}$ (12) $\displaystyle
i\neq j\qquad$ $\displaystyle k\neq l$ (13)

One can then expand the products, equating the real and complex sides of each
equation, as well as combining all variables into a single vector for
computational simplicity. Now we have a set of real equations in which every
feasible point is a set of MUBs, thus proving infeasibility of the equation
set for dimension $d$ and number of bases $n$ proves the inability to generate
MUBs for that $d$ and $n$. The number of variables quickly becomes large, as
shown in Table 1, which becomes a problem when discussing the global
optimisation required to prove infeasibility, .

n/d | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---
2 | 12 | 27 | 48 | 75 | 108
3 | 24 | 54 | 96 | 150 | 216
4 | 40 | 90 | 160 | 250 | 360
5 | - | 135 | 240 | 375 | -
6 | - | - | 336 | 525 | -
7 | - | - | - | 700 | -
Table 1: Number of variables when converting the mutually unbiased basis
problem into a real polynomial optimisation problem, without any reductions.
The formula for this is $nd^{2}+\frac{1}{2}d^{2}n(n-1)$. Note that due to the
non-convexity of the problem the worst-case complexity of such general
polynomial optimisation is exponential in the number of variables, thus here
the open question of $d=6$ and $n=4$ remains highly intractable.

##  III Reducing the number of variables using symmetries

Whilst one could use the algorithms presented in this work to optimise the
polynomial problem in its current state, we were unable to perform global
optimisation for even small sets in dimensions 2 and 3. Luckily for us, MUBs
have a large number of rotational and reflectional symmetries that we can
exploit to significantly reduce the number of variables and equations.

The first of these reductions can be derived from the fact that a unitary
rotation applied globally to a set of MUBs preserves the mutually unbiasedness
romero2005structure . As such, one can always choose the first basis to be any
valid orthonormal basis without loss of generality. For the sake of
simplicity, we choose the computational basis mentioned earlier, which as well
as removing a large number of variables from the optimisation also allows for
a further simplification: given that all inner products must have a fixed
inner product and given that each vector in the computational basis only
affects one element of each vector, one can now explicitly constraint the
magnitude of every element of every other vector to be $\frac{1}{\sqrt{d}}$,
for example in dimension 2:

$\displaystyle\text{if }\quad\vec{v}_{11}$
$\displaystyle=(1\leavevmode\nobreak\ \leavevmode\nobreak\
0)\quad\text{and}\quad\vec{v}_{12}=(0\leavevmode\nobreak\ \leavevmode\nobreak\
1)$ (14) $\displaystyle\implies\vec{v}_{21}$
$\displaystyle=\frac{1}{\sqrt{2}}(e^{i\theta_{1}}\leavevmode\nobreak\
\leavevmode\nobreak\ e^{i\theta_{2}})$ (15)

One can also choose to fix some of the rotational symmetries by fixing the
first vector of the second basis, for instance to the uniform vector of
$\vec{v}_{21}=\frac{1}{\sqrt{d}}(1\leavevmode\nobreak\ \leavevmode\nobreak\
1\leavevmode\nobreak\ \leavevmode\nobreak\ ...\leavevmode\nobreak\
\leavevmode\nobreak\ 1)$. A similar argument can be made for fixing the first
element of every vector to be $\frac{1}{\sqrt{d}}$, further reducing the
number of variables we need to optimise. These can both be done without loss
of generality due to the inner product being unaffected by the phases of the
bases. If one could find 4 MUBs in dimension 6 without such restrictions, one
could then simply add a phase to each basis/vector until we have MUBs that
follow the aforementioned restrictions.

Continuing with the idea of symmetry reductions, one can consider a number of
other symmetries which preserve the various inner products of the problem,
specifically certain permutations of variable. For instance, given a set of
MUBs, swapping the first and second bases is guaranteed to preserve the
mutually unbiasedness. The same can be said about swapping two vectors within
a basis, swapping the first and second elements of all vectors, as well as
taking the conjugate of everything. All of these are transformations which
preserve the inner product. In order to reduce the size of the search space we
define positivity constraints for each of the aforementioned symmetries, such
that one can always swap the sign of the constraint function by performing the
corresponding transformation on the equation set, whilst retaining mutual
unbiasedness.

For example, since it is known that we can swap vectors 1 and 2 in the 3rd
basis without loss of generality, once can define some function of the form
$h_{1}(\vec{v}_{31},\vec{v}_{32})=\sum_{m=1}^{d}(a_{31m}+b_{31m})-\sum_{m=1}^{d}(a_{32m}+b_{32m})$,
the difference between the sums of the components of each vector. This can
always be made positive by MUBs, since if a set of MUBs was to make this
negative one could simply swap the two bases, resulting in the negative of the
function. A similar function can be defined for the conjugation symmetry:
$h_{2}(\vec{b})=\sum_{m}b_{m}$. Given MUBs making this function negative one
can take the conjugate of the MUBs resulting in the negative of all of the
imaginary components $b_{m}$, thus flipping the sign of this function whilst
retaining the mutual unbiasedness. Adding each of these functions as
positivity constraints to the polynomial problem reduces the size of the
feasible set without loss of generality in the infeasibility case.

##  IV Reducing the number of variables using sub-bases

Whilst we have reduced the number of variables significantly with the
reductions so far, we are still considering the full equation set. If we now
shift the focus to not finding full sets of MUBs, but instead showing that
MUBs are infeasible, the question can be asked, is any subset of the equations
infeasible? If this were the case the infeasibility would extend to the full
equation set, whilst being much simpler to prove due to the reduction in
number of equations. Here we focus on the idea of optimising sub-bases, which
we define as a series of orthonormal vectors which does not contain enough
vectors to fully span the space, thus making it a subset of some basis. For
example, in dimension 3 one could consider the sub-basis of the computational
basis:

$\displaystyle v_{11}=(1\leavevmode\nobreak\ \leavevmode\nobreak\
0\leavevmode\nobreak\ \leavevmode\nobreak\ 0)\qquad
v_{12}=(0\leavevmode\nobreak\ \leavevmode\nobreak\ 1\leavevmode\nobreak\
\leavevmode\nobreak\ 0)$ (16)

Since finding a valid sub-basis is a requirement for finding a valid basis, if
a given problem does not permit sub-bases, it certainly doesn’t permit full
bases. As such we extend the notation with the mutually unbiased property,
turning our search from mutually unbiased bases (MUBs) to the sub-problem of
searching for mutually unbiased sub-bases (MUSBs), which we denote by listing
the number of vectors per set. For example in dimension 3, the MUSB problem
with sizes $\\{3,2,1\\}$ would imply the search for a full basis (of 3
vectors), a sub-basis of 2 vectors, plus a sub-basis of 1 vector, each being
mutually unbiased with each other. The choice of such sub-basis sizes, in
order to minimize the size of the problem whilst still remaining infeasible in
the cases that should be infeasible, is something that was determined by
starting with the smallest MUSB problem for a given $d$ and $n$, then
increasing the size until the problem appears to become infeasible as decided
by the upcoming methods. Combining all reductions results in a significant
drop in the number of variables, as given by table 2. For the open question of
4 MUBs in dimension 6 this results in a 66% reduction ($360\to 122$).

n/d | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---
2 | 0 | 0 | 6 | 8 | 50
3 | 2 | 4 | 14 | 28 | 119
4 | 6 | 10 | 24 | 42 | 122
5 | - | 18 | 36 | 58 | -
6 | - | - | 50 | 76 | -
7 | - | - | - | 96 | -
Table 2: Number of variables when converting the mutually unbiased basis
problem into a real polynomial optimisation problem, with all reductions. Here
the set sizes used are those with the lowest number of variables whilst still
retaining infeasibility in the presumed infeasible cases.

##  V Method 1: Search Space Without Local Minima

Given this reduced set of equations, we first considered methods to ensure
that this correctly conveys the problem, finding MUBs in the cases where they
should exist and failing to find MUBs in those where they are not. To do this
we use the idea of combining our many second-order equality constraints into
one giant fourth-order objective function, which could then be minimized
without constraints. Consider first some general equation set in which we want
to find a common root:

$\displaystyle h_{i}(\vec{x})=0\qquad\forall i$ (17)

which we can then combine by squaring each equation and summing:

$\displaystyle f(\vec{x})=\sum_{i}h_{i}^{2}(\vec{x})=0$ (18)

Since each term has to be non-negative, the only way such a function can reach
zero is if all of the squared equations also reach zero, thus transforming the
problem into an unconstrained optimisation. Here we just use the equality
constraints of the reduced equation set and ignore the positivity constraints
used for enforce some symmetries, since these only provide a slight
improvement to this method but add quite a bit of complexity. Our aim with
this optimisation is to converge rapidly to MUSBs if they exist, with a lack
of convergence serving as a suggestion that such a set of MUSBs (and thus a
corresponding set of MUBs) is unlikely to exist.

Based on the method of Lagrange multipliers, we begin by integrating our
polynomial $f(\vec{x})$ by one of our variables, $x_{m}$, specifically one
that is known to be able to reach zero without loss of generality. In the case
of the MUB problem this is fine, since we have a number of variables already
assumed to be zero due to rotational symmetries, however equally one could
create a new variable and integrate by it to create something that resembles a
traditional Lagrange multiplier. Regardless of the choice of variable, the
result is then a new fifth-order polynomial $f_{m}(\vec{x})$ with a number of
interesting properties. Considering the partial derivatives of this function
with respect to every variable, we find that one of these derivatives is our
original function, whilst every other derivative has our integrating variable
$x_{m}$ as a factor. As an arbitrary example, if our original function was:

$\displaystyle f(\vec{x})=x_{0}^{2}+x_{0}x_{1}+x_{2}^{4}$ (19)

and we choose to integrate by $x_{1}$:

$\displaystyle f_{1}(\vec{x})=\int f(\vec{x})\leavevmode\nobreak\
dx_{1}=x_{1}x_{0}^{2}+\frac{1}{2}x_{0}x_{1}^{2}+x_{1}x_{2}^{4}$ (20)

which has partial derivatives:

$\displaystyle\frac{\partial}{\partial x_{0}}f_{1}(\vec{x})$
$\displaystyle=2x_{1}x_{0}+\frac{1}{2}x_{1}^{2}$ (21)
$\displaystyle\frac{\partial}{\partial x_{1}}f_{1}(\vec{x})$
$\displaystyle=f(\vec{x})$ (22) $\displaystyle\frac{\partial}{\partial
x_{2}}f_{1}(\vec{x})$ $\displaystyle=4x_{1}x_{2}^{3}$ (23)

From this it can be seen that the set of partial derivatives can only be all
zero, a stationary point, if our original function is zero. The other
derivatives being zero can be guaranteed to exist if the variable $x_{1}$ can
also go to zero without loss of generality. As such, we have created a search
space in which every single stationary point is a solution to our problem. In
the case of MUBs this means a search space in which all stationary points are
MUBs, with no stationary points existing in the cases in which MUBs do not
exist.

Now the task is to find a stationary point, which we do using a Newton-style
iterative approach. This is done by calculating the matrix of second-
derivatives $H$ as well as the vector of partial derivatives $\vec{g}$ of our
function $f_{m}(\vec{x})$, then solving the system of linear equations
$H\vec{p}=\vec{g}$ to get the update direction $\vec{p}$. We then follow this
direction, scaled by some factor $0<\alpha\leq 1$. The result in an algorithm
which will eventually converge (given a small enough step size) if there are
MUBs, whilst endlessly traversing the search space if MUBs do not exist. As
such it does not constitute a formal proof of non-existence, simply serving as
suggestion as to whether the given problem is feasible or not.

For the numerical computation of this method we use the library Eigen eigen
for the linear system solving required for the Newton step. We tested a number
of solving methods but found that a partial-pivot LU decomposition works the
best in our case, offering an effective blend of speed and accuracy, as well
as offering a parallel implementation. Whilst we did experiment with a line-
search for choosing the parameter $\alpha$, this seemed to require more time
than it saved. We also tested a quasi-Newton update method in which we did not
explicitly calculate the full Hessian, but the loss of accuracy resulted in
very slow convergence. We also add a small ($10^{-10}$) term to the diagonals
of the Hessian to improve the numerical stability.

##  VI Method 2: Branch and bound combined with a hierarchy of SDPs

For our second method, we use another optimization technique: a hierarchy of
semidefinite programming relaxations similar to the Lasserre
lasserre2001global  or NPA navascues2008convergent  hierarchies for commuting
and non-commuting polynomial optimization, respectively. The idea behind such
techniques is to relax the problem and then add constraints of increasing
complexity (referred to as levels in some hierarchy) to tighten the
relaxation, often with guarantees of finite convergence. We begin by taking
our non-linear equations and creating new variables to linearize them, turning
each term $x_{i}x_{j}$ into a new variable $x_{ij}$, for example:

$\displaystyle x_{1}x_{2}+x_{2}^{2}+x_{1}\to x_{12}+x_{22}+x_{1}$ (25)

We start by assuming that the new variables have the same bounds as the
originals, meaning at this point the relaxation is quite bad, for instance
$x_{ij}$ could be negative despite $x_{i}$ and $x_{j}$ being positive since
for now there is no constraint linking them. It is, however, fast to compute,
since we only have linear equations and thus it has become a convex
optimisation problem. To then begin to constrain $x_{ij}$ to be closer to
$x_{i}x_{j}$ we add a hierarchy of positive-semidefinite (PSD) constraints,
which are also convex and thus the relaxation remains convex. The first level
of such a hierarchy considers the matrix whose first row contains all first-
order monomials:

$\displaystyle\begin{pmatrix}1&x_{1}&x_{2}\\\ x_{1}&x_{11}&x_{12}\\\
x_{2}&x_{12}&x_{22}\end{pmatrix}\succcurlyeq 0$ (26)

Going to the second level in our two variable example would look as follows,
with second order terms on the top row and then forming the corresponding
products:

$\displaystyle\begin{pmatrix}1&x_{1}&x_{2}&x_{11}&x_{12}&x_{22}\\\
x_{1}&x_{11}&x_{12}&x_{111}&x_{112}&x_{122}\\\
x_{2}&x_{12}&x_{22}&x_{112}&x_{122}&x_{222}\\\
x_{11}&x_{111}&x_{112}&x_{1111}&x_{1112}&x_{1122}\\\
x_{12}&x_{112}&x_{122}&x_{1112}&x_{1122}&x_{1222}\\\
x_{22}&x_{122}&x_{222}&x_{1122}&x_{1222}&x_{2222}\\\ \end{pmatrix}\succcurlyeq
0$ (27)

What we now have is an optimisation over a much larger set of variables with
linear equality constraints corresponding to the MUB problem and linear
inequality constraints to break symmetries, all to some chosen level of
semidefinite hierarchy with the hope that at some level this convex relaxation
of the problem will become infeasible in the cases of MUB non-existence.
However, we find that even high levels (>5) of this hierarchy alone do not
become infeasible, even for small cases. Some works dealing with the NPA
hierarchy suggest that one may need at least level 12 of that hierarchy
gribling2021mutually . If the level needed for the commuting case is similar,
it would require solving a semidefinite program of $110^{12}$ variables to
solve the $d=6$ $n=4$ case, which is at least six orders of magnitude larger
than what modern supercomputers can solve within a reasonable time. As such,
we ask the question: maybe we don’t need to solve it in a single convex
optimisation? Perhaps splitting it into a larger number of smaller convex
regions would be sufficient?

Consider now the first level of the first-order relaxation of a single
variable. Since it’s a 2x2 matrix, there are only 2 eigenvalues. For both to
be positive, the product of eigenvalues needs to be positive and thus the
determinant needs to be positive:

$\displaystyle\begin{pmatrix}1&x_{1}\\\ x_{1}&x_{11}\end{pmatrix}\succcurlyeq
0\qquad\implies\qquad x_{11}-x_{1}^{2}\geq 0$ (28)

Assuming we also use box constraints of $-1<x_{1}<1$, which we can assume in
the MUB case since all of the magnitudes are bounded, the search space looks
as given in Figure 1.

Figure 1: The search space created from the constraint $x_{11}-x_{1}^{2}\geq
0$ as well as the box constraint $-1\leq x_{1}\leq 1$

As one can see, this constraint still leaves many points such as $(0,1)$ well
within the search space whilst being very far from the $x_{11}=x_{1}^{2}$ line
to which we want to constrain. This shows the issue with trying to convexify
such an equality constraint. Our approach here is to stay in the first level
of the hierarchy whilst splitting such a convex region into two. Whilst for
many styles of constraint splitting a region into two would only half the size
of the region, in our case we can constrain further since we know the edge we
are trying to converge towards, as shown by Figure 2.

Figure 2: The two convex search spaces created from the constraint
$x_{11}-x_{1}^{2}>0$ with the box constraints $-1\leq x_{1}\leq 0$ and $0\leq
x_{1}\leq 1$

By first performing the first-level relaxation without any branching we obtain
some point in the search space. Checking the error of each variable versus its
square $\text{err}(x_{i})=(x_{ii}-x_{i}^{2})^{2}$ we get a heuristic for which
variables need constraining tighter. We found that a good branching method is
to choose the variable with the highest error, splitting it in two (as in
Figure 2) at the feasible point that was found by the previous iteration. As
such we have divided our whole search space into two, whilst also constraining
tighter towards the set of valid first and second-order monomials. One can
also choose to add an objective function to ensure that we find the "most-
feasible" point, done by adding a new variable $\lambda$ and asking that the
moment matrix plus $\lambda$ times the identity is positive-semidefinite, then
minimizing $\lambda$.

By repeating this process one can see that we will eventually reach regions
which are sufficiently small that the maximum error for any variable becomes
less than some $\epsilon$. In the MUB non-existence cases, under the
assumption that there exists some $\epsilon$ such that one cannot find MUBs
even allowing such error in the relations (e.g. that the various inner
products are now constrained within $\frac{1}{\sqrt{d}}\pm\epsilon$), such
small regions will result in convex infeasibility, which can be efficiently
detected. By checking that a large number of these convex regions are all
infeasible one can show that the original equation set is infeasible, thus
proving the MUB inexistence case.

Whilst in the above case we consider the first-level relaxation, moving to the
second level must be at least as constraining as the first, since the second-
level PSD matrix contains the first as a submatrix. The procedure works the
same as we move to higher-level relaxations, except now we have higher-order
terms in the first row of our PSD matrix, making it much larger. For example,
given 18 variables (as in the case of our reduction of the 5 MUBs in dimension
3 case), the first-level relaxation has a PSD matrix of 19x19, the second a
matrix of 343x343, the third 6175x6175 and so on. Given that SDPs in the size
of 300x300 are already difficult for desktop computers, doing a large number
of 6000x6000 SDPs is near intractable. The hope, however, would be that by
doing a tighter relaxation one does not need to perform so many branches to
prove full infeasibility, giving a trade-off between the number of iterations
needed and the level of the relaxation.

For the numerical computation of this method we have tested two SDP solvers:
MOSEK mosek  and Splitting Cone Solver (SCS) scs . For the small cases and the
first-level SDPs MOSEK works well, but has an incredibly high memory profile
due to its need to store various second-order derivatives, however this allows
it high accuracy and fast (logarithmic) convergence versus other methods. The
alternative is SCS, which as a first order method has a very minimal memory
usage, often using 100x less memory than MOSEK, at the cost of significantly
longer convergence times. For all second-level optimisation we use SCS, since
otherwise memory quickly becomes an issue, with MOSEK using around 80 GB for
proving dimension 3 non-existence whilst SCS uses only 100 MB. With both
libraries we exploit their sparsity and model parameterization features.

##  VII Results

We first begin by showing the amount of MUBs that seem to be feasible and
infeasible using our first method, shown in Table 3, given as the number of
iterations it required before we either found a satisfying vector or we
stopped the optimization. With this we show how easily the algorithm finds
solutions to contrast with the amount of time spend on the presumed non-
existence cases.

$n$ / $d$ | 2 | 3 | 4 | 5 | 6 | 10
---|---|---|---|---|---|---
2 | 1 | 1 | 55 | 200 | 84 | 175
3 | 59 | 55 | 61 | 63 | 73 | 226
4 | 100000 | 113 | 75 | 104 | 100000 | 100000
5 | 100000 | 100000 | 87 | 96 | 100000 | 100000
6 | 100000 | 100000 | 100000 | 235 | 100000 | 100000
7 | 100000 | 100000 | 100000 | 100000 | 100000 | 100000
Table 3: The number of iterations required to either find the MUSB or until
reaching the chosen iteration limit. These results, if we assume reaching the
limit of 100000 iterations implies infeasibility, agree with known results and
Zauner’s conjecture. Here the criteria for convergence was reaching a value
below $10^{-13}$ in the combined equation, roughly corresponding to a maximum
error of $10^{-7}$ in any of the original equality constraints.

The algorithm in general seems to find MUSBs with set sizes with less
deviation (i.e. $\\{4,4,3,3\\}$) faster than those that are heavily deviated
(i.e. $\\{6,5,2,1\\}$), presumably because there are more MUSBs to find due to
an increased number of symmetries. For highly uneven set sizes like
$\\{9,7,1,1,1,1,1,1,1,1\\}$ in dimension 9, the algorithm struggles, sometimes
taken several thousand iterations. For high dimensions (>15) there is some
instability that develops from the large linear systems, to combat this we
lower the step size $\alpha$ and in some cases add a larger diagonal term to
the Hessian to improve the condition number and thus numerical stability of
the solver.

Now applying the branch-and-bound algorithm, we correctly convergence to the
expected results in all cases in dimension 2 and 3. Whilst it was expected
that we would not find MUSBs of sizes {2,1,1,1} in dimension 2, as far as the
authors are aware it has not been shown that one cannot find MUSBs of sizes
{2,1,1,1,1} in dimension 3. We show in Table 4 the number of SDPs that were
required when branching with the first level of the hierarchy, whilst Table 5
shows the number using the second level. As expected, significantly less
branches were needed as higher levels of the hierarchy were used, however
using the third level or higher is intractable for anything higher than
dimension 2.

$n$ / $d$ | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---
2 | 0 | 0 | 11 | 43 | ?
3 | 2 | 40 | 81 | 19279 | ?
4 | 11 | 337 | 18004 | ? | ?
5 | - | 17153 | ? | ? | -
6 | - | - | ? | ? | -
7 | - | - | - | ? | -
Table 4: Number of SDPs required to reach either convergence to an MUSB or to prove infeasibility when using the branch-and-bound algorithm with the first level of the hierarchy. Convergence to an MUSB is defined as an SDP being feasible with a largest monomial error of at most $10^{-8}$. Bold text indicates the relevant problems for proving infeasibility. Here the 17153 iterations for dimension 3 takes around 3-4 minutes on a standard desktop. $n$ / $d$ | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---
2 | 0 | 0 | 13 | 13 | ?
3 | 1 | 18 | 35 | ? | ?
4 | 1 | 22 | ? | ? | ?
5 | - | 7 | ? | ? | -
6 | - | - | ? | ? | -
7 | - | - | - | ? | -
Table 5: Number of SDPs required to reach either convergence to MUSB or
proving infeasibility when using the branch-and-bound algorithm with the
second level of the hierarchy. Bold text indicates the relevant problems for
proving infeasibility. Here the 7 iterations for dimension 3 takes around a
minute on a standard desktop with SCS.

As this algorithm runs, it regularly declares regions as infeasible, each of
which has an associated area. Since in our problem the total search space is
bounded (each variable between $\pm\frac{1}{\sqrt{d}}$) we can obtain the
fraction of the search space that has been declared infeasible. If the entire
problem is infeasible then this total value increases until reaching 100%,
otherwise a solution is found before then. By considering the time taken to
reach the current percentage of infeasibility and extrapolating to 100%, the
algorithm progressively outputs an estimate for the time it will take to prove
infeasibility. For the $d=3$ and $n=5$ case (18 variables) this states 5
minutes at the start, the algorithm then finishes in around 4. When running
for $d=4$ and $n=6$ (50 variables) it estimates several months. For $d=6$ and
$n=4$ (118 variables), the open problem, it estimates something on the order
of hundreds of millions of years. Even assuming perfect parallelisation and
running the algorithm on the worlds best supercomputer (at the time of writing
the American “Frontier” rajaraman2023frontier ) it would still take hundreds
of years to solve this open problem with this algorithm. As such, further
reductions are needed, otherwise we may simply have to wait until classical
computational power reaches a level such that this becomes tractable.

Considering the results from both methods, we notice there appears to exist
some sort of pattern in the minimum sizes that are infeasible, for instance we
find that for dimension 6 $\\{6,3,3,3\\}$ does not seem to be feasible, whilst
$\\{6,3,3,2\\}$ is feasible and the algorithm finds a solution in very few
iterations. Although our first thoughts were that perhaps it relates to the
total number of vectors, it cannot be so simple due to the lack of freedom of
the final vector in an orthonormal set, for instance that in dimension 3,
$\\{2,2,2,2,2\\}$ and $\\{3,3,3,3,3\\}$ will both be equally infeasible
despite having very different total vector counts. We list a selection of
feasible and infeasible sets for dimensions 3, 4, 5 and 6 in Tables 6, 7, 8
and 9 respectively.

One possible conjecture is that the MUSB problem is infeasible if the number
of independent constraints is at least the number of variables minus one.
Whilst there are known theorems for the number of solutions of linear systems
that are overdetermined and underdetermined williams1990overdetermined , the
results do not extend trivially to non-linear systems. Whilst there are some
MUSBs such as $\\{6,6,1,1\\}$ which have more constraints (95) than variables
(92) despite being feasible, one can perform an analysis of the equations to
show that some are dependant on others. We have found a set of 10 constraints
for $\\{6,6,1,1\\}$ that even when removed, the optimisation seems to always
find a solution which satisfies them, suggesting that the other 85 equations
were sufficiently constraining. By extension this would also mean these
constraints are redundant for $\\{6,6,3,1\\}$, thus taking it from 149
equations to 139 equations vs 136 variables, leaving it still infeasible under
this conjecture and agreeing with the optimisation results.

set sizes | vecs | vars | eqns | result
---|---|---|---|---
{3,3,3,3} | 12 | 74 | 101 | proven feasible
{1,1,1,1,1} | 5 | 18 | 15 | proven feasible
{2,1,1,1,1} | 6 | 18 | 18 | proven infeasible
{3,1,1,1,1} | 7 | 18 | 18 | proven infeasible
Table 6: A list of some set sizes for dimension 3 as well as the conclusion made using our optimisation techniques. Here “proven feasible” means a valid MUSB was found to a tolerance of $10^{-13}$ using method 1 whilst “proven infeasible” means that we have branched the entire search space to infeasibility using method 2. set sizes | vecs | vars | eqns | result
---|---|---|---|---
{4,4,4,4,4} | 258 | 20 | 354 | proven feasible
{1,1,1,1,1,1} | 6 | 36 | 26 | proven feasible
{2,1,1,1,1,1} | 7 | 36 | 30 | proven feasible
{2,2,1,1,1,1} | 8 | 50 | 45 | proven feasible
{2,2,2,1,1,1} | 9 | 64 | 62 | proven feasible
{2,2,2,2,1,1} | 10 | 80 | 82 | seems infeasible
{4,1,1,1,1,1} | 9 | 36 | 34 | proven feasible
{4,2,1,1,1,1} | 10 | 50 | 50 | seems infeasible
Table 7: A list of some set sizes for dimension 4 as well as the conclusion made using our optimisation techniques. Here “proven feasible” implies a valid MUSB was found to a tolerance of $10^{-13}$ using method 1, whilst “seems infeasible” implies our search for a feasible point using method 1 was inconclusive even after 100000 iterations. set sizes | vecs | vars | eqns | result
---|---|---|---|---
{5,5,5,5,5,5} | 30 | 652 | 902 | proven feasible
{1,1,1,1,1,1,1} | 7 | 60 | 40 | proven feasible
{2,2,2,2,2,1,1} | 12 | 138 | 132 | proven feasible
{2,2,2,2,2,2,1} | 13 | 162 | 161 | seems infeasible
{3,2,2,1,1,1,1} | 11 | 96 | 90 | proven feasible
{3,3,2,1,1,1,1} | 12 | 116 | 114 | proven feasible
{3,3,3,1,1,1,1} | 13 | 136 | 140 | seems infeasible
{5,1,1,1,1,1,1} | 11 | 60 | 55 | proven feasible
{5,2,1,1,1,1,1} | 12 | 78 | 75 | proven feasible
{5,3,1,1,1,1,1} | 13 | 96 | 97 | seems infeasible
Table 8: A list of some set sizes for dimension 5 as well as the conclusion made using our optimisation techniques. Here “proven feasible” implies a valid MUSB was found to a tolerance of $10^{-13}$ using method 1, whilst “seems infeasible” implies our search for a feasible point using method 1 was inconclusive even after 100000 iterations. set sizes | vecs | vars | eqns | result
---|---|---|---|---
{6,6,6} | 18 | 170 | 206 | proven feasible
{1,1,1,1} | 4 | 22 | 7 | proven feasible
{2,1,1,1} | 5 | 22 | 9 | proven feasible
{2,2,1,1} | 6 | 36 | 18 | proven feasible
{2,2,2,1} | 7 | 50 | 29 | proven feasible
{2,2,2,2} | 8 | 66 | 43 | proven feasible
{3,2,2,1} | 8 | 50 | 33 | proven feasible
{3,2,2,2} | 9 | 66 | 48 | proven feasible
{3,3,3,2} | 11 | 102 | 86 | proven feasible
{3,3,3,3} | 12 | 122 | 109 | proven feasible
{4,2,1,1} | 8 | 36 | 24 | proven feasible
{4,3,3,3} | 13 | 122 | 117 | proven feasible
{4,4,3,3} | 14 | 144 | 144 | seems infeasible
{6,1,1,1} | 9 | 22 | 15 | proven feasible
{6,6,1,1} | 14 | 92 | 95 | proven feasible
{6,6,2,1} | 15 | 114 | 121 | proven feasible
{6,6,3,1} | 16 | 136 | 149 | seems infeasible
{6,3,3,2} | 14 | 102 | 100 | proven feasible
{6,3,3,3} | 15 | 122 | 125 | seems infeasible
{6,3,2,2,1} | 14 | 106 | 102 | proven feasible
{6,3,2,2,2} | 15 | 128 | 128 | seems infeasible
{6,2,2,1,1,1,1,1} | 15 | 134 | 131 | proven feasible
{6,2,2,2,1,1,1,1} | 16 | 158 | 160 | seems infeasible
Table 9: A list of some set sizes for dimension 6 as well as the conclusion
made using our optimisation techniques. Here “proven feasible” implies a valid
MUSB was found to a tolerance of $10^{-13}$ using method 1, whilst “seems
infeasible” implies our search for a feasible point using method 1 was
inconclusive even after 100000 iterations.

The algorithms used here are general and can be applied to any commuting
polynomial optimisation problem using our library, written as a single C++
header file. The library represents polynomials using hash tables, allowing
for high performance even with large equations. We release this code as an
open-source project mygithub , although it still needs work regarding user-
friendliness and documentation, for now only serving as a reference for the
methods used as well as containing the data files for each optimisation
referenced. The example relevant to this work can be compiled using “make mub”
and then run the “mub” executable as “./mub -h” to display the help text.

##  VIII Conclusion

We explored the polynomial representation of the open mutually unbiased bases
existence problem, showing a number of reductions in order to take full
advantage of the available symmetries. The result is then a second-order real
non-linear set of equations. We then show two methods of approaching such a
problem, a Lagrange-multiplier Newtonian-descent algorithm which we use for
quickly checking which of our relaxations are feasible, as well as a global
branch-and-bound algorithm. We introduce the notation of mutually unbiased
sub-bases and show that proving infeasibility can be done by proving
infeasibility of sub-bases, reducing the number of variables drastically.

Our reductions regarding the MUB problem may simplify future proofs and
generate new questions about why sub-bases cannot be found in certain
dimensions and why the minimum sub-basis size seems to follow a particular
pattern. Unfortunately, solving the open question in this framework by
searching the whole space still remains intractable, with our estimates
suggesting it is still beyond the capacity of even the most advanced
supercomputers at the time of writing. Nevertheless, given humanity’s ongoing
quest for greater computational power, we are confident that the technology of
the near future will bring the ability to solve this problem, accelerated by
any further reductions and simplifications we find.

##  IX Acknowledgements

Thanks to Máté Farkas, Leonardo Zambrano and Antonio Acín for the useful
discussions.

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 847517.

## References

  * [1] Gelo Noel M Tabia, Varun Satya Raj Bavana, Shih-Xian Yang, and Yeong-Cherng Liang.  Bell inequality violations with random mutually unbiased bases.  Physical Review A, 106(1):012209, 2022.

  * [2] Armin Tavakoli, Máté Farkas, Denis Rosset, Jean-Daniel Bancal, and Jedrzej Kaniewski.  Mutually unbiased bases and symmetric informationally complete measurements in bell experiments.  Science advances, 7(7):eabc3847, 2021.
  * [3] Lorena Rebón, C Iemmi, and S Ledesma.  Phase-measurement interferometry as a simulation of optimal quantum-state tomography.  Optik, 124(22):5548–5552, 2013.
  * [4] Hao Yuan, Zheng-Wei Zhou, and Guang-Can Guo.  Quantum state tomography via mutually unbiased measurements in driven cavity qed systems.  New Journal of Physics, 18(4):043013, 2016.
  * [5] Mhlambululi Mafu, Angela Dudley, Sandeep Goyal, Daniel Giovannini, Melanie McLaren, Miles J Padgett, Thomas Konrad, Francesco Petruccione, Norbert Lütkenhaus, and Andrew Forbes.  Higher-dimensional orbital-angular-momentum-based quantum key distribution with mutually unbiased bases.  Physical Review A, 88(3):032305, 2013.
  * [6] William K Wootters.  A wigner-function formulation of finite-state quantum mechanics.  Annals of Physics, 176(1):1–21, 1987.
  * [7] Gerhard Zauner.  Grundzüge einer nichtkommutativen designtheorie.  Ph. D. dissertation, PhD thesis, 1999.
  * [8] Maria Prat Colomer, Luke Mortimer, Irénée Frérot, Máté Farkas, and Antonio Acín.  Three numerical approaches to find mutually unbiased bases using bell inequalities.  Quantum, 6:778, 2022.
  * [9] Philippe Raynal, Xin Lü, and Berthold-Georg Englert.  Mutually unbiased bases in six dimensions: The four most distant bases.  Physical Review A, 83(6):062303, 2011.
  * [10] Stéphane Jacquet and Sylvain Hallé.  Reformulation of sat into a polynomial box-constrained optimization problem.  In Integrated Formal Methods: 16th International Conference, IFM 2020, Lugano, Switzerland, November 16–20, 2020, Proceedings 16, pages 387–394. Springer, 2020.
  * [11] Narendra Karmarkar.  A new polynomial-time algorithm for linear programming.  In Proceedings of the sixteenth annual ACM symposium on Theory of computing, pages 302–311, 1984.
  * [12] Monique Laurent and Franz Rendl.  Semidefinite programming and integer programming.  Handbooks in Operations Research and Management Science, 12:393–514, 2005.
  * [13] Sander Gribling and Sven Polak.  Mutually unbiased bases: polynomial optimization and symmetry.  arXiv preprint arXiv:2111.05698, 2021.
  * [14] Stephen Brierley and Stefan Weigert.  Mutually unbiased bases and semi-definite programming.  In Journal of Physics: Conference Series, volume 254, page 012008\. IOP Publishing, 2010.
  * [15] JL Romero, Gunnar Björk, AB Klimov, and LL Sánchez-Soto.  Structure of the sets of mutually unbiased bases for n qubits.  Physical Review A, 72(6):062310, 2005.
  * [16] Gaël Guennebaud, Benoît Jacob, et al.  Eigen v3.  http://eigen.tuxfamily.org, 2010.
  * [17] Jean B Lasserre.  Global optimization with polynomials and the problem of moments.  SIAM Journal on optimization, 11(3):796–817, 2001.
  * [18] Miguel Navascués, Stefano Pironio, and Antonio Acín.  A convergent hierarchy of semidefinite programs characterizing the set of quantum correlations.  New Journal of Physics, 10(7):073013, 2008.
  * [19] MOSEK ApS.  The MOSEK optimization toolbox for C++ manual. Version 10.0., 2022\.
  * [20] Brendan O’Donoghue, Eric Chu, Neal Parikh, and Stephen Boyd.  Conic optimization via operator splitting and homogeneous self-dual embedding.  Journal of Optimization Theory and Applications, 169(3):1042–1068, June 2016.
  * [21] V Rajaraman.  Frontier—world’s first exaflops supercomputer.  Resonance, 28(4):567–576, 2023.
  * [22] Gareth Williams.  Overdetermined systems of linear equations.  The American Mathematical Monthly, 97(6):511–513, 1990.
  * [23] Luke Mortimer.  poly.  https://github.com/Lumorti/poly, 2023.