These authors contributed equally to this work.

These authors contributed equally to this work.

# Three numerical approaches to find mutually unbiased bases using Bell
inequalities

Maria Prat Colomer  ICFO-Institut de Ciencies Fotoniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels, Spain  CFIS-Centre
de Formació Interdisciplinària Superior, UPC-Universitat Politècnica de
Catalunya, 08028 Barcelona, Spain  Luke Mortimer  ICFO-Institut de Ciencies
Fotoniques, The Barcelona Institute of Science and Technology, 08860
Castelldefels, Spain  Irénée Frérot  ICFO-Institut de Ciencies Fotoniques, The
Barcelona Institute of Science and Technology, 08860 Castelldefels, Spain
Univ Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 38000 Grenoble, France
Máté Farkas<EMAIL_ADDRESS>ICFO-Institut de Ciencies Fotoniques, The
Barcelona Institute of Science and Technology, 08860 Castelldefels, Spain
Antonio Acín  ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of
Science and Technology, 08860 Castelldefels, Spain  ICREA-Institucio Catalana
de Recerca i Estudis Avançats, Lluis Companys 23, 08010 Barcelona, Spain

###### Abstract

Mutually unbiased bases correspond to highly useful pairs of measurements in
quantum information theory. In the smallest composite dimension, six, it is
known that between three and seven mutually unbiased bases exist, with a
decades-old conjecture, known as Zauner’s conjecture, stating that there exist
at most three. Here we tackle Zauner’s conjecture numerically through the
construction of Bell inequalities for every pair of integers $n,d\geq 2$ that
can be maximally violated in dimension $d$ if and only if $n$ MUBs exist in
that dimension. Hence we turn Zauner’s conjecture into an optimisation
problem, which we address by means of three numerical methods: see-saw
optimisation, non-linear semidefinite programming and Monte Carlo techniques.
All three methods correctly identify the known cases in low dimensions and all
suggest that there do not exist four mutually unbiased bases in dimension six,
with all finding the same bases that numerically optimise the corresponding
Bell inequality. Moreover, these numerical optimisers appear to coincide with
the “four most distant bases” in dimension six, found through numerically
optimising a distance measure in [P. Raynal, X. Lü, B.-G. Englert, Phys. Rev.
A, 83 062303 (2011)]. Finally, the Monte Carlo results suggest that at most
three MUBs exist in dimension ten.

##  1 Introduction

Mutually unbiased bases (MUBs), on the one hand, are highly symmetric bases in
complex Hilbert spaces and, on the other hand, correspond to pairs of quantum
measurements. The defining property of a pair of MUBs is that the overlaps
between any two vectors from the two different bases is uniform. This property
translates to the corresponding measurements as follows: If a measurement
yields a definite outcome when measured on a quantum state, then a measurement
unbiased to it will yield a uniformly random outcome on the same state. This
feature makes MUBs widely useful in quantum information processing. MUBs were
originally introduced in the context of optimal state determination [1], but
since have been found to be useful in a variety of quantum information
processing tasks, such as quantum cryptography [2, 3, 4], quantum
communication tasks [5, 6], Bell inequalities [7, 8, 9] and so on (for a
review, see Ref. [10]).

While MUBs have been extensively studied both in the quantum information and
the mathematics community for decades, there are still open questions
regarding their structure. Most notably, the maximal number of bases that are
pairwise mutually unbiased is unknown for general Hilbert space dimension. A
general upper bound was shown by Wootters and Fields, stating that in
dimension $d$ there exist no more than $d+1$ MUBs [11]. In the same work, they
showed that this upper bound is saturated in prime power dimensions, by
providing an explicit construction. However, in composite dimensions the only
known generic lower bound on the number of MUBs is $p^{r}+1$, where $p^{r}$ is
the smallest prime power in the prime decomposition of the dimension (this
lower bound is shown using tensor products of the Wootters–Fields
construction). While in certain dimensions this lower bound has been improved
[12], there exists no composite dimension in which the exact number of MUBs is
known. For the smallest composite dimension, six, the number of MUBs is known
to be no more than seven and no less than three, from the general bounds.
However, which of the numbers in between is the exact number of MUBs in
dimension six is unknown (apart from the fact that it cannot be six, following
from a general result by Weiner [13]). It was first conjectured by Zauner in
1999 that there are no more than three MUBs in dimension six [14], and this
conjecture has not been resolved to date, despite substantial efforts.

There are numerous works trying to prove (or disprove) Zauner’s conjecture,
both analytically and numerically. While not providing an exhaustive list of
references here, let us note that on the analytic side, it has been shown that
Zauner’s conjecture is equivalent to a conjecture on orthogonal decompositions
of Lie algebras [15]. Furthermore, there exist various analytic constructions
of MUB triplet families (see Refs. [16, 17] and references therein), but thus
far there has not been found even a single vector that is unbiased to all the
vectors in any of these triplets. For notable recent developments on Zauner’s
conjecture see Refs. [18, 19].

On the numerical side, Bengtsson et al. introduced a distance measure of two
bases that is maximised if and only if the bases are mutually unbiased [20].
This construction turns the problem of finding a set of MUBs into an
optimisation problem, maximising all the pairwise distances within a set of
bases. Using this approach, Raynal, Lü and Englert later constructed a two-
parameter family of four bases in dimension six, such that for certain values
of the parameters, the bases coincide with the numerical maximiser of the
distance function [21]. These four bases are not MUBs, but based on the
numerical evidence the authors refer to them as “the four most distant bases
in dimension six”.

Since MUBs optimise various quantum information processing tasks, it is
natural to measure the closeness of a set of bases to MUBs in terms of some
quantum information processing protocol. This was studied formally by Aguilar
et al. [22], using the fact that MUBs optimise the success probability of a
communication task called quantum random access codes (QRACs). A slight
generalisation of the QRAC task is then optimised by a set of $n$ MUBs, and
finding $n$ MUBs in dimension $d$ corresponds to optimising the associated
success probability. With this method, Aguilar et al. managed to re-prove the
non-existence of $d+2$ MUBs in certain low dimensions using quantum
information theoretic tools. However, the case of dimension six remains open.

In this work, we employ similar ideas to tackle Zauner’s conjecture. Namely,
we study a recently introduced family of Bell inequalities, known to be
maximally violated by a pair of MUBs in dimension $d$ [9]. We then extend
these inequalities to new ones, maximally violated by a set of $n$ MUBs in
dimension $d$. Then, we apply three numerical methods for finding the maximal
value of these Bell inequalities in a fixed dimension. Namely, we apply see-
saw semidefinite programming (SDP), non-linear SDP, and Monte Carlo
techniques. While these methods are heuristic—in the sense that there is no
guarantee for finding a global maximum—they find the maximum in all the cases
where the maximum is known (i.e., it is known that $n$ MUBs exist in the given
dimension $d$). Furthermore, when applying these techniques to dimension six
and four bases, all the different numerical tools converge to the same bases,
and these four bases are—numerically—very close to the “four most distant
bases” of Ref. [21] (one should not expect exact equality, since the measure
optimised in Ref. [21] is different from the Bell inequalities we optimise).
Hence, our results support Zauner’s conjecture, even though no rigorous claim
can be made due to the heuristic nature of our methods. Finally, we were able
to implement the Monte Carlo algorithm for $d=10$, where—similarly to $d=6$—we
do not find more than three MUBs.

##  2 Preliminaries

In this section, we introduce the mathematical background and concepts
necessary for turning the MUB problem into an optimisation problem. Namely, we
formally introduce MUBs, Bell inequalities and the specific family of Bell
inequalities tailored for MUBs.

###  2.1 Mutually unbiased bases

Let us take a $d$-dimensional Hilbert space $\mathcal{H}\cong\mathbb{C}^{d}$,
and two orthonormal bases on it, $\\{|\hskip
0.7ptb^{1}_{j}\rangle\\}_{j=1}^{d}$ and $\\{|\hskip
0.7ptb^{2}_{k}\rangle\\}_{k=1}^{d}$. We say that these two bases are mutually
unbiased if

$|\langle b^{1}_{j}\hskip 0.7pt|\hskip
0.7ptb^{2}_{k}\rangle|^{2}=\frac{1}{d}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall
j,k\in[d],$ (1)

where $[d]\equiv\\{1,2,\ldots,d\\}$. A simple example in dimension two is the
computational basis $\\{|\hskip 0.7pt0\rangle,|\hskip 0.7pt1\rangle\\}$ and
the Hadamard basis $\\{\frac{1}{\sqrt{2}}(|\hskip 0.7pt0\rangle+|\hskip
0.7pt1\rangle),\frac{1}{\sqrt{2}}(|\hskip 0.7pt0\rangle-|\hskip
0.7pt1\rangle)\\}$.

One may also associate an orthonormal basis with a quantum measurement. In
general, a quantum measurement is described by a positive operator-valued
measure (POVM), which, in the $d$-dimensional, $d$-outcome case corresponds to
a set of $d$ positive semidefinite operators $B_{j}\geq 0$ on
$\mathbb{C}^{d}$, adding up to the identity operator $\mathbb{1}$. Given an
orthonormal basis $\\{|\hskip 0.7ptb_{j}\rangle\\}_{j=1}^{d}$ on
$\mathbb{C}^{d}$, one can define the corresponding POVM $\\{B_{j}=|\hskip
0.7ptb_{j}\rangle\langle b_{j}\hskip 0.7pt|\\}_{j=1}^{d}$, consisting of
rank-1 projections onto the basis elements. We say that two measurements are
MUBs if they correspond to a pair of orthonormal bases that are MUBs.

###  2.2 Bell inequalities

We look at MUBs in the context of Bell scenarios (see Ref. [23] for a review).
Bell scenarios describe physical experiments performed by two distant parties,
usually referred to as Alice and Bob. These parties share many copies of a
bipartite (quantum) state, and perform local measurements on these copies. The
experiment is described by the correlation, $p$, with elements $p(a,b|x,y)$
specifying the probability of Alice (Bob) observing locally the outcome $a$
($b$) upon choosing the measurement setting $x$ ($y$). In quantum theory, the
shared state is described by a density operator $\rho\geq 0$ with unit trace
($\operatorname{tr}\rho=1$) on a tensor product Hilbert space
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. The measurements are described by
local POVMs $\\{A^{x}_{a}\\}$ and $\\{B^{y}_{b}\\}$ on the Hilbert spaces
$\mathcal{H}_{A}$ and $\mathcal{H}_{B}$, respectively. For a fixed state and
measurements, the correlation is given by the Born rule,

$p(a,b|x,y)=\operatorname{tr}[\rho(A^{x}_{a}\otimes B^{y}_{b})].$ (2)

Note that for a pure state, $\rho=|\hskip 0.7pt\psi\rangle\langle\psi\hskip
0.7pt|$, with $|\hskip
0.7pt\psi\rangle\in\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ and
$\langle\psi\hskip 0.7pt|\hskip 0.7pt\psi\rangle=1$, the Born rule reduces to

$p(a,b|x,y)=\langle\psi\hskip 0.7pt|A^{x}_{a}\otimes B^{y}_{b}|\hskip
0.7pt\psi\rangle.$ (3)

Bell functionals are linear functionals of correlations, i.e., functionals of
the form

$W(p)=\sum_{a,b,x,y}c_{abxy}p(a,b|x,y),$ (4)

where $c_{abxy}$ are real coefficients. Non-trivial Bell inequalities are Bell
functionals for which $W(p)\leq\beta_{L}$ holds for every correlation of the
form

$p(a,b|x,y)=\int_{\Lambda}\mathrm{d}\mu(\lambda)p_{A}(a|x,\lambda)p_{B}(b|y,\lambda)$
(5)

(also called a local correlation, considered as the notion of classicality in
Bell scenarios), but for which there exists a quantum correlation $p$ of the
form (2) such that $W(p)>\beta_{L}$. In Eq. (5), $\Lambda$ is a measurable set
with a probability measure $\mu$, $\lambda\in\Lambda$, and $p_{A}$ and $p_{B}$
are conditional probability distributions.

While the original interest in Bell inequalities was precisely this separation
of local and quantum correlations, we will be interested in their quantum
maximum (or maximal quantum violation), i.e., the tight upper bound
$W(p)\leq\beta_{Q}$ satisfied by all correlations of the form (2).

###  2.3 Bell inequalities for mutually unbiased bases

In this work, we are interested in a family of Bell inequalities that was
introduced in Ref. [9], and is parametrised by an integer $d\geq 2$. For a
fixed $d$, Alice has $d^{2}$ measurement settings labelled as $x=x_{1}x_{2}$
with $x_{1},x_{2}\in[d]$. Each of these measurements has three outcomes,
$a\in\\{1,2,\perp\\}$. Bob, on the other hand, has two measurement settings,
$y\in\\{1,2\\}$, with $d$ outcomes each, $b\in[d]$. The Bell inequality then
reads

$\begin{split}W_{d}(p)&\left.=\sum_{x_{1},x_{2},y}\big{[}p(y,x_{y}|x_{1}x_{2},y)-p(\bar{y},x_{y}|x_{1}x_{2},y)\big{]}\right.\\\
&\left.-\frac{1}{2}\sqrt{\frac{d-1}{d}}\sum_{x_{1},x_{2}}\big{[}p_{A}(1|x_{1}x_{2})+p_{A}(2|x_{1}x_{2})\big{]},\right.\end{split}$
(6)

where $\bar{y}$ flips the value of $y\in\\{1,2\\}$, and
$p_{A}(a|x)=\sum_{b}p(a,b|x,y)$ is the marginal probability distribution of
Alice (which is independent of $y$).

This is a non-trivial Bell inequality with maximal quantum violation
$\beta_{Q}=\sqrt{d(d-1)}$ [9]. The maximal violation can be achieved with the
maximally entangled $d$-dimensional state $|\hskip
0.7pt\phi^{+}_{d}\rangle\equiv\frac{1}{\sqrt{d}}\sum_{j=1}^{d}|\hskip
0.7ptj\rangle\otimes|\hskip 0.7ptj\rangle$, and any pair of MUB measurements
on Bob’s side. Moreover, if the dimension is fixed to be $d$, this is the only
way in which the maximal quantum violation can be achieved, up to local
unitary freedom [9]. This property of the Bell inequality (6) forms the core
of our numerical approaches to construct MUBs.

The above Bell inequality can be straightforwardly extended to a set of $n$
measurements on Bob’s side. The Bell inequality for $n$ measurements is a sum
of Bell inequalities of the form (6). For each pair $y,z\in[n]$ such that
$y<z$ (denoted in the following as $(y,z)\in\text{Pairs}[n]$), we introduce
$d^{2}$ settings for Alice, labelled as $x=(y,z)x_{y}x_{z}$ with
$x_{y},x_{z}\in[d]$, and take a copy of the Bell inequality in Eq. (6),
defined as

$\begin{split}W^{(y,z)}_{d}(p)&\left.=\sum_{x_{y},x_{z},w}\big{[}p(a_{w},x_{w}|(y,z)x_{y}x_{z},w)\right.\\\
&\left.-p(\bar{a}_{w},x_{w}|(y,z)x_{y}x_{z},w)\big{]}\right.\\\
&\left.-\frac{1}{2}\sqrt{\frac{d-1}{d}}\sum_{x_{y},x_{z}}\big{[}p_{A}(1|(y,z)x_{y}x_{z})\right.\\\
&\left.+p_{A}(2|(y,z)x_{y}x_{z})\big{]},\right.\end{split}$ (7)

where $w\in\\{y,z\\}$, $a_{y}=1$, $a_{z}=2$, and $\bar{a}_{w}$ flips the value
of $a_{w}$. The final Bell inequality then reads

$W_{d,n}(p)=\sum_{(y,z)\in\text{Pairs}[n]}W^{(y,z)}_{d}(p).$ (8)

It is clear that $W_{d,n}(p)\leq\binom{n}{2}\sqrt{d(d-1)}$, by applying the
known bound to each individual term in the above sum. Moreover, if the
dimension is $d$, the only way to reach this maximum (up to local unitary
freedom) is by using the maximally entangled state, and if the $n$
measurements on Bob’s side correspond to MUBs. Hence, we can reformulate the
MUB problem in terms of these Bell inequalities:

######  Proposition 1.

$W_{d,n}(p)=\binom{n}{2}\sqrt{d(d-1)}$ can be achieved in dimension $d$ if and
only if $n$ MUBs exist in dimension $d$.

###  2.4 The optimisation problem

According to the above proposition, finding $n$ MUBs in dimension $d$ can be
cast as an optimisation problem, maximising $W_{d,n}(p)$ over $d$-dimensional
quantum states and measurements. To see how to do this explicitly, let us
first write out the Bell inequality $W_{d}$ in Eq. (6) in terms of a quantum
state $|\hskip 0.7pt\psi\rangle$ and measurements $\\{A^{x}_{a}\\}$,
$\\{B^{y}_{b}\\}$, using Born’s rule in Eq. (3):

$\begin{split}W_{d}(|\hskip
0.7pt\psi\rangle,&\left.\\{A^{x}_{a}\\},\\{B^{y}_{b}\\})=\sum_{x_{1},x_{2},y}\big{(}\langle\psi\hskip
0.7pt|A^{x_{1}x_{2}}_{y}\otimes B^{y}_{x_{y}}|\hskip
0.7pt\psi\rangle\right.\\\ &\left.-\langle\psi\hskip
0.7pt|A^{x_{1}x_{2}}_{\bar{y}}\otimes B^{y}_{x_{y}}|\hskip
0.7pt\psi\rangle\big{)}\right.\\\
&\left.-\frac{1}{2}\sqrt{\frac{d-1}{d}}\sum_{x_{1},x_{2}}\big{(}\langle\psi\hskip
0.7pt|A^{x_{1}x_{2}}_{1}\otimes\mathbb{1}|\hskip 0.7pt\psi\rangle\right.\\\
&\left.+\langle\psi\hskip 0.7pt|A^{x_{1}x_{2}}_{2}\otimes\mathbb{1}|\hskip
0.7pt\psi\rangle\big{)}.\right.\\\ &\left.=\sum_{j,k}\langle\psi\hskip
0.7pt|\big{[}(A^{jk}_{1}-A^{jk}_{2})\otimes(B^{1}_{j}-B^{2}_{k})\right.\\\
&\left.-\frac{1}{2}\sqrt{\frac{d-1}{d}}(A^{jk}_{1}+A^{jk}_{2})\otimes\mathbb{1}\big{]}|\hskip
0.7pt\psi\rangle,\right.\end{split}$ (9)

where in the second equality we wrote out the summation over $y$ and switched
to the notation $x_{1}x_{2}\to jk$ with $j,k\in[d]$. Maximising the inequality
(9) in terms of the state and the measurements can then be written as the
optimisation problem

$\begin{split}\max_{|\hskip
0.7pt\psi\rangle,\\{A^{x}_{a}\\},\\{B^{y}_{b}\\}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ &\left.W_{d}(|\hskip
0.7pt\psi\rangle,\\{A^{x}_{a}\\},\\{B^{y}_{b}\\})\right.\\\
\text{s.t.}\quad\quad\quad&\left.|\hskip
0.7pt\psi\rangle\in\mathbb{C}^{d}\otimes\mathbb{C}^{d},\leavevmode\nobreak\
\leavevmode\nobreak\ \langle\psi\hskip 0.7pt|\hskip
0.7pt\psi\rangle=1\right.\\\
&\left.A^{x}_{a},B^{y}_{b}\in\mathcal{L}_{\text{sa}}(\mathbb{C}^{d})\leavevmode\nobreak\
\leavevmode\nobreak\ \forall a,b,x,y\right.\\\ &\left.A^{x}_{a}\geq
0\leavevmode\nobreak\ \leavevmode\nobreak\ \forall a,x,\leavevmode\nobreak\
\leavevmode\nobreak\ B^{y}_{b}\geq 0\leavevmode\nobreak\ \leavevmode\nobreak\
\forall b,y\right.\\\ &\left.\sum_{a}A^{x}_{a}=\mathbb{1}\leavevmode\nobreak\
\leavevmode\nobreak\ \forall x,\leavevmode\nobreak\ \leavevmode\nobreak\
\sum_{b}B^{y}_{b}=\mathbb{1}\leavevmode\nobreak\ \leavevmode\nobreak\ \forall
y,\right.\end{split}$ (10)

where $\mathcal{L}_{\text{sa}}(\mathbb{C}^{d})$ is the set of self-adjoint
linear operators on $\mathbb{C}^{d}$. Optimising $W_{d,n}$ can be written in a
similar fashion, with:

$\begin{split}W_{d,n}&\left.(|\hskip
0.7pt\psi\rangle,\\{A^{x}_{a}\\},\\{B^{y}_{b}\\})=\sum_{(y,z)\in\text{Pairs}[n]}\Bigg{\\{}\right.\\\
&\left.\sum_{j,k}\langle\psi\hskip
0.7pt|\big{[}(A^{(y,z)jk}_{1}-A^{(y,z)jk}_{2})\otimes(B^{y}_{j}-B^{z}_{k})\right.\\\
&\left.-\frac{1}{2}\sqrt{\frac{d-1}{d}}(A^{(y,z)jk}_{1}+A^{(y,z)jk}_{2})\otimes\mathbb{1}\big{]}|\hskip
0.7pt\psi\rangle\Bigg{\\}}.\right.\end{split}$ (11)

The optimisation problem is then

$\begin{split}\max_{|\hskip
0.7pt\psi\rangle,\\{A^{x}_{a}\\},\\{B^{y}_{b}\\}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ &\left.W_{d,n}(|\hskip
0.7pt\psi\rangle,\\{A^{x}_{a}\\},\\{B^{y}_{b}\\})\right.\\\
\text{s.t.}\quad\quad\quad&\left.|\hskip
0.7pt\psi\rangle\in\mathbb{C}^{d}\otimes\mathbb{C}^{d},\leavevmode\nobreak\
\leavevmode\nobreak\ \langle\psi\hskip 0.7pt|\hskip
0.7pt\psi\rangle=1\right.\\\
&\left.A^{x}_{a},B^{y}_{b}\in\mathcal{L}_{\text{sa}}(\mathbb{C}^{d})\leavevmode\nobreak\
\leavevmode\nobreak\ \forall a,b,x,y\right.\\\ &\left.A^{x}_{a}\geq
0\leavevmode\nobreak\ \leavevmode\nobreak\ \forall a,x,\leavevmode\nobreak\
\leavevmode\nobreak\ B^{y}_{b}\geq 0\leavevmode\nobreak\ \leavevmode\nobreak\
\forall b,y\right.\\\ &\left.\sum_{a}A^{x}_{a}=\mathbb{1}\leavevmode\nobreak\
\leavevmode\nobreak\ \forall x,\leavevmode\nobreak\ \leavevmode\nobreak\
\sum_{b}B^{y}_{b}=\mathbb{1}\leavevmode\nobreak\ \leavevmode\nobreak\ \forall
y.\right.\end{split}$ (12)

In this work, we consider various approaches to solve the optimisation problem
(12). In particular, from Proposition 1 it follows that $n$ MUBs exist in
dimension $d$ if and only if the solution of the above optimisation problem is
$\binom{n}{2}\sqrt{d(d-1)}$. We will facilitate the problem using knowledge
about the optimal realisation of the Bell inequality from Ref. [9].

First of all, we notice that the value $\binom{n}{2}\sqrt{d(d-1)}$ can only be
achieved in dimension $d$ with the maximally entangled state [9]. Without loss
of generality, we therefore impose that $|\hskip 0.7pt\psi\rangle=|\hskip
0.7pt\phi^{+}_{d}\rangle=\frac{1}{\sqrt{d}}\sum_{j=1}^{d}|\hskip
0.7ptj\rangle\otimes|\hskip 0.7ptj\rangle$. We can then use the fact that for
any two operators $A$ and $B$ on $\mathbb{C}^{d}$ we have that
$\langle\phi^{+}_{d}\hskip 0.7pt|A\otimes B|\hskip
0.7pt\phi^{+}_{d}\rangle=\frac{1}{d}\operatorname{tr}(A^{T}B)$, where
$(.)^{T}$ is the transposition in the basis $\\{|\hskip 0.7ptj\rangle\\}$. As
a second simplification, we notice that in order to saturate the bound
$\binom{n}{2}\sqrt{d(d-1)}$ in dimension $d$, Alice’s measurement operators
$A^{x}_{1}$ and $A^{x}_{2}$, and all of Bob’s measurement operators must be
trace-1 [9]. For such operators we have that $\langle\phi^{+}_{d}\hskip
0.7pt|A^{x}_{1}\otimes\mathbb{1}|\hskip
0.7pt\phi^{+}_{d}\rangle=\langle\phi^{+}_{d}\hskip
0.7pt|A^{x}_{2}\otimes\mathbb{1}|\hskip 0.7pt\phi^{+}_{d}\rangle=\frac{1}{d}$.
The second term in Eq. (11) is then a constant, $-\binom{n}{2}\sqrt{d(d-1)}$,
and does not influence the optimisation problem. The simplified Bell
expression finally reads

$\begin{split}W^{+}_{d,n}&\left.(\\{A^{x}_{a}\\},\\{B^{y}_{b}\\})=\sum_{(y,z)\in\text{Pairs}[n]}\Bigg{\\{}\right.\\\
&\left.\frac{1}{d}\sum_{j,k}\operatorname{tr}\big{(}(A^{(y,z)jk}_{1}-A^{(y,z)jk}_{2})^{T}(B^{y}_{j}-B^{z}_{k})\big{)}\Bigg{\\}},\right.\end{split}$
(13)

and its maximum quantum value $W^{+}_{d,n}$ satisfies:

$W^{+}_{d,n}\leq W_{\text{MUB}}(d,n)\equiv
n(n-1)\sqrt{d(d-1)}\leavevmode\nobreak\ .$ (14)

The simplified optimisation problem becomes

$\begin{split}\max_{\\{A^{x}_{a}\\},\\{B^{y}_{b}\\}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
&\left.W^{+}_{d,n}(\\{A^{x}_{a}\\},\\{B^{y}_{b}\\})\right.\\\
\text{s.t.}\quad\quad&\left.A^{x}_{a},B^{y}_{b}\in\mathcal{L}_{\text{sa}}(\mathbb{C}^{d})\leavevmode\nobreak\
\leavevmode\nobreak\ \forall a,b,x,y\right.\\\ &\left.A^{x}_{a}\geq
0\leavevmode\nobreak\ \leavevmode\nobreak\ \forall a,x,\leavevmode\nobreak\
\leavevmode\nobreak\ B^{y}_{b}\geq 0\leavevmode\nobreak\ \leavevmode\nobreak\
\forall b,y\right.\\\ &\left.\operatorname{tr}A^{x}_{a}=1\leavevmode\nobreak\
\leavevmode\nobreak\ \forall a,x,\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname{tr}B^{y}_{b}=1\leavevmode\nobreak\ \leavevmode\nobreak\ \forall
y,b\right.\\\ &\left.\sum_{a}A^{x}_{a}=\mathbb{1}\leavevmode\nobreak\
\leavevmode\nobreak\ \forall x,\leavevmode\nobreak\ \leavevmode\nobreak\
\sum_{b}B^{y}_{b}=\mathbb{1}\leavevmode\nobreak\ \leavevmode\nobreak\ \forall
y.\right.\end{split}$ (15)

The optimal value of this optimisation problem is $n(n-1)\sqrt{d(d-1)}$ if and
only if $n$ MUBs exist in dimension $d$.

We may further simplify the optimisation problem. From Ref. [9] we know that
in the optimal realisation the $B^{y}_{b}$ operators are rank-1 projections,
$B^{y}_{j}=|\hskip 0.7ptb^{y}_{j}\rangle\langle b^{y}_{j}\hskip 0.7pt|$ (they
are projections onto the basis elements of MUBs). For such operators, we have
that

$(B^{y}_{j}-B^{z}_{k})^{3}=[1-\operatorname{tr}(B^{y}_{j}B^{z}_{k})](B^{y}_{j}-B^{z}_{k}),$
(16)

which implies that the spectrum of $B^{y}_{j}-B^{z}_{k}$ is contained in
$\\{0,\pm\lambda^{yz}_{jk}\\}$, where
$\lambda^{yz}_{jk}\equiv\sqrt{1-\operatorname{tr}(B^{y}_{j}B^{z}_{k})}=\sqrt{1-|\langle
b^{y}_{j}\hskip 0.7pt|\hskip 0.7ptb^{z}_{k}\rangle|^{2}}$. Moreover, we have
that $\operatorname{tr}[(B^{y}_{j}-B^{z}_{k})^{2}]=2(\lambda^{yz}_{jk})^{2}$
and $\operatorname{tr}(B^{y}_{j}-B^{z}_{k})=0$, and therefore
$B^{y}_{j}-B^{z}_{k}$ has one eigenvalue $\lambda^{yz}_{jk}$, one eigenvalue
$-\lambda^{yz}_{jk}$ and the rest of the eigenvalues are 0. Furthermore, in
the optimal realisation we have that $(A^{(y,z)jk}_{1})^{T}$ is the rank-1
projection onto the eigenspace of $B^{y}_{j}-B^{z}_{k}$ corresponding to
$\lambda^{yz}_{jk}$, and $(A^{(y,z)jk}_{2})^{T}$ is the rank-1 projection onto
the eigenspace of $B^{y}_{j}-B^{z}_{k}$ corresponding to $-\lambda^{yz}_{jk}$.
With these final simplifications the Bell expression reads

$\displaystyle W^{+B}_{d,n}(\\{B^{y}_{b}\\})$ $\displaystyle=$
$\displaystyle\frac{2}{d}\sum_{(y,z)\in\text{Pairs}[n]}\sum_{j,k}\sqrt{1-\operatorname{tr}(B^{y}_{j}B^{z}_{k})}$
(17) $\displaystyle=$
$\displaystyle\frac{2}{d}\sum_{(y,z)\in\text{Pairs}[n]}\sum_{j,k}\sqrt{1-|\langle
b_{j}^{y}|b_{k}^{z}\rangle|^{2}}$ (18)

and the corresponding optimisation problem is

$\begin{split}\max_{\\{B^{y}_{b}\\}}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ &\left.W^{+B}_{d,n}(\\{B^{y}_{b}\\})\right.\\\
\text{s.t.}\quad&\left.B^{y}_{b}\in\mathcal{L}_{\text{sa}}(\mathbb{C}^{d})\leavevmode\nobreak\
\leavevmode\nobreak\ \forall a,b,x,y\right.\\\
&\left.(B^{y}_{b})^{2}=B^{y}_{b}\leavevmode\nobreak\ \leavevmode\nobreak\
\forall b,y\right.\\\ &\left.\operatorname{tr}B^{y}_{b}=1\leavevmode\nobreak\
\leavevmode\nobreak\ \forall b,y\right.\\\
&\left.\sum_{b}B^{y}_{b}=\mathbb{1}\leavevmode\nobreak\ \leavevmode\nobreak\
\forall y.\right.\end{split}$ (19)

Note that projectivity and self-adjointness implies positive semidefiniteness,
and therefore positive semidefiniteness does not need to be imposed.

The optimal value of the optimisation problem (19), denoted as $W_{\rm
max}(d,n)$, is $n(n-1)\sqrt{d(d-1)}$ if and only if $n$ MUBs exist in
dimension $d$. In fact, Alice’s measurements have been completely removed from
the problem, and one may regard Eq. (19) as a purely geometrical problem: find
$n$ orthonormal bases, $\big{\\{}\\{|\hskip
0.7ptb_{j}^{y}\rangle\leavevmode\nobreak\ |\leavevmode\nobreak\
j\in[d]\\}\leavevmode\nobreak\ |\leavevmode\nobreak\ y\in[n]\big{\\}}$,
maximising the “MUB-ness measure” of Eq. (18).

In the following three sections we apply three numerical methods to solve the
optimisation problems (15) and (19) in order to numerically tackle Zauner’s
conjecture.

##  3 See-saw SDP

###  3.1 Methodology

One arrives at a relatively simple method of optimising problem (15) by
noticing that the objective function is bi-linear in the $A^{x}_{a}$ and
$B^{y}_{b}$ matrices. That is, if every $A^{x}_{a}$ is fixed, then the problem
simplifies to optimising a linear functional of the $B^{y}_{b}$ matrices with
a series of linear and positive semidefinite constraints. This is a standard
SDP, for which there exist efficient solvers.

The _see-saw_ optimisation technique starts with fixing the set of $A^{x}_{a}$
matrices satisfying the constraints of the problem (15), either with random
values or based on some prior knowledge. We then solve the problem for the
$B^{y}_{b}$ matrices, which is a standard SDP. Then, we fix the $B^{y}_{b}$
matrices to the optimum found, and solve the resulting SDP for the $A^{x}_{a}$
matrices, and so on. By repeating this process, the system eventually
converges to a stable result, i.e. the value of the objective function does
not change beyond a given threshold within a certain window of iterations (a
change of less than $10^{-9}$ for $10$ iterations in our implementation).
Although the see-saw method has never been proven to converge to the global
optimum, for our current problem it has never failed to converge within the
chosen precision if given sufficient time, always to the value expected (i.e.,
to $W_{\text{MUB}}(d,n)$ whenever it is known that $n$ MUBs exist in dimension
$d$). We implemented the see-saw algorithm with the help of the SDP solving
library MOSEK [24].

###  3.2 Results

The values obtained with the see-saw algorithm are shown in Table 1 (the
values displayed are $1-W_{d,n}/W_{\text{MUB}}(d,n)$ for easier comparison
across different $n$ and $d$, where $W_{d,n}$ is the result of the
optimisation). Notice that whenever $n$ MUBs exist in dimension $d$, the see-
saw method correctly converges to the MUB solution, and whenever it is known
that $n$ MUBs do not exist in dimension $d$, the method does indeed converge
to a value less than $W_{\text{MUB}}(d,n)$. For the unknown case of four MUBs
in dimension six, the see-saw method could not find four MUBs, supporting
Zauner’s conjecture. Furthermore, the optimal measurements found by the see-
saw method are numerically very close to the “four most distant bases” of Ref.
[21] (see also Section 6.3). Note, however, that the results simply mean that
the method could not find four MUBs in dimension six, but one cannot rule out
the possibility that they exist.

The method also has consistent convergence, for example, for $d=2$ and $n=2$,
all of $10000$ see-saw optimisations from random starts converged to the
correct value of $W_{\text{MUB}}(2,2)=2.82843$ up to five decimal places,
albeit finding a different set of optimum matrices. Similarly, for $10000$
optimisations for $d=2$ and $n=4$, all optimisations converged to the same
value of $16.72616$ (correctly signifying non-existence). We did not perform a
similar sized convergence analysis for the larger dimensions due to the time
required to optimise, however, several runs were always performed to certify a
minimal level of consistency.

All of the results for this method were obtained on a desktop PC with 8GB of
RAM using 4 cores, with times varying between milliseconds for the smallest
problem ($d=2$, $n=2$) and hours for the largest one ($d=6$, $n=4$). Since SDP
solvers are efficiently parallelised, this method offers good parallel
scaling, however, the memory requirement is the highest of all of our methods,
since it requires explicit storage and optimisation of the $A^{x}_{a}$
matrices.

$n$ $d$ | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---
2 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000
3 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000
4 | 0.01440 | 0.00000 | 0.00000 | 0.00000 | 0.00004
5 | - | 0.00391 | 0.00000 | 0.00000 | -
6 | - | - | 0.00186 | 0.00000 | -
7 | - | - | - | 0.00091 | -
Table 1: The values $1-W_{d,n}/W_{\text{MUB}}(d,n)$, where $W_{d,n}$ is the
result of the optimisation, obtained with the see-saw method after convergence
for various dimensions ($d$) and numbers of bases ($n$), to 5 decimal places.
The values depicted are consistently obtained by multiple runs. The values in
bold indicate that the value was at least $10^{-5}$ away from zero, meaning
that the method could not find $n$ MUBs in dimension $d$. The “-” symbol
indicates that we have not performed the optimisation. Note that in all the
known cases, the algorithm predicts correctly existence/nonexistence, and it
predicts that four MUBs do not exist in dimension six.

##  4 Non-linear SDP

###  4.1 Methodology

An alternative approach to the optimisation is to focus on problem (19), which
features only the $B^{y}_{b}$ matrices and thus contains fewer variables for a
reduction in the size of the search space as well as the memory required. The
downside, however, is that the objective function is now non-linear (not even
bi-linear) and thus many efficient solvers (i.e. for standard SDP systems) can
no longer be applied. To optimise this problem we adapt a method based on the
work by Yamashita et al. for optimising a non-linear SDP using a primal-dual
interior point method [25]. This method, assuming a few basic conditions
(discussed later), is guaranteed to converge to a Karush–Kuhn–Tucker (KKT)
point, a point satisfying a series of constraints known as the KKT conditions,
which are necessary for optimality [26].

These KKT conditions are only sufficient (imply a global minimum) in a subset
of cases, the main of which being that the problem is convex, which is
unfortunately not true in our case. An interesting property of our search-
space, however, is that by taking the derivative of our objective function it
can be shown that all local minima are global minima if the constraints are
met, thus implying that our search space is, in fact, a series of
discontinuous convex regions. This signifies that for the MUB-existence case
this method will always converge to MUBs, although no claim can be made for
the non-existence case.

In order to implement the method of Yamashita et al., we parametrise the
measurement operators $\\{B^{x}_{y}\\}$ by a real vector
$\mathbf{x}=(x_{i})_{i}$. To be able to deal with real numbers instead of
complex ones, we note that every self-adjoint matrix $B=B_{r}+\mathrm{i}B_{i}$
(where $B_{r}$ is real symmetric and $B_{i}$ is real anti-symmetric) can be
mapped to the real symmetric matrix $\hat{B}$ via

$B\mapsto\hat{B}=\begin{bmatrix}B_{r}&B_{i}\\\ -B_{i}&B_{r}\end{bmatrix}.$
(20)

It is easy to verify that $B\geq 0$ if and only if $\hat{B}\geq 0$, and
$\operatorname{tr}B=\frac{1}{2}\operatorname{tr}\hat{B}$.

We therefore define—in line with the method of Yamasitha et al.—a matrix
$X(\mathbf{x})=\sum_{i}C_{i}x_{i}+D$, which is a block diagonal matrix
containing the $\hat{B}^{y}_{b}$ matrices on its diagonal in such a way that
the linear constraints $\sum_{b}B^{y}_{b}=\mathbb{1}$ and
$\operatorname{tr}B^{y}_{b}=1$ of the optimisation problem (19) are already
enforced. The real parameters $x_{i}$ correspond to those elements of the
$\hat{B}^{y}_{b}$ matrices that are free after enforcing the linear
constraints. While the constraint $B^{y}_{b}\geq 0$ is superfluous for the
problem (19), we chose to include this in our optimisation problem, as this
constraint is heavily used in the method of Yamashita et al. With the
parametrisation above, this constraint is equivalent to $X(\mathbf{x})\geq 0$.
The last remaining constraint is projectivity, $(B^{y}_{b})^{2}=B^{y}_{b}$ for
all $y$ and $b$, which is equivalent to $X^{2}(\mathbf{x})=X(\mathbf{x})$. We
enforce this constraint through
$g(\mathbf{x})\equiv|\\!|X^{2}(\mathbf{x})-X(\mathbf{x})|\\!|_{F}^{2}=0$,
where $|\\!|.|\\!|_{F}$ is the Frobenius norm. Further, we denote the
objective function in terms of $\mathbf{x}$ by $W(\mathbf{x})$, suppressing
the $d,n$ indices whenever it does not lead to confusion.

The method requires introducing Lagrange multipliers (dual variables) for
every constraint. In our case, there is a single inequality constraint
$X(\mathbf{x})\geq 0$, to which we assign the dual variable $Z$, which is a
matrix with the same dimensions as $X(\mathbf{x})$. Furthermore, we have a
single equality constraint, $g(\mathbf{x})=0$, to which we assign the dual
variable $y$, which is a scalar. The resulting Lagrangian reads

$L(\mathbf{x},y,Z)=W(\mathbf{x})-yg(\mathbf{x})-\operatorname{tr}[Z^{T}X(\mathbf{x})].$
(21)

The algorithm for solving the optimisation problem is iterative, and each
iteration begins with the calculation of $G$, the Hessian of the Lagrangian.
In our case this can be quite expensive so we opt to use the alternative
update method also proposed in Ref. [25] based on the
Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm, which approximates the
Hessian without requiring the full calculation of the second derivatives. This
$G$ is then used to form a series of linear equations which we solve using a
stabilised bi-conjugate gradient method in order to obtain the update
directions for the primal ($\mathbf{x}$) and dual ($y,Z$) variables. Following
this, a simple line search is performed to find the optimum step size, then
the variables are updated. This process is then repeated until the following
barrier KKT conditions are met for some barrier parameter $\mu$:

$\displaystyle r(\mathbf{x},y,Z,\mu)\equiv\begin{pmatrix}\nabla
L(\mathbf{x},y,Z)\\\ g(\mathbf{x})\\\ X(\mathbf{x})Z-\mu
I\end{pmatrix}=\begin{pmatrix}0\\\ 0\\\ 0\end{pmatrix}$ (22)

Through repeated iterations of this method with values of $\mu$ converging
towards zero, this algorithm has been proven to always converge to a KKT point
of the system, assuming certain conditions. The first of these is that the
functions $W(\mathbf{x})$ and $g(\mathbf{x})$ are both twice continuously
differentiable, which in our case is true: $W(\mathbf{x})$ is simply a sum of
square roots of polynomials of $\mathbf{x}$, whilst $g(\mathbf{x})$ is a
vector of polynomials of $\mathbf{x}$. The second condition is that the vector
$\mathbf{x}$ must remain within a finite set during the optimisation, which
for us is true since infinite values are non-optimal for the MUB functional.
The third condition is that the matrices $C_{i}$ must be linearly independent,
which for us is true since they serve to place a single component of
$\mathbf{x}$ to one (diagonal) or two (off-diagonal) positions in the matrix
$X(\mathbf{x})$.

Similarly to the the see-saw method, there exists some non-determinism to this
method due to it starting from a randomised interior point. This is found by
first beginning with a random vector, then performing gradient descent using
the derivatives of the constraints until the vector satisfies all constraints
within some precision. This vector is then used as the starting point for the
iterative method. Although perhaps one could start with a “good guess” for the
vector and optimise from there, we decided to begin from a random point to
attempt to cover a wider search space.

###  4.2 Results

The values obtained through this algorithm are shown in Table 2. The optimal
values, as well as the optimal bases found, agree with the result of the see-
saw method up to high numerical precision. In particular, we find the same set
of four bases to numerically optimise the six dimensional case (see Section
6.3 for further details).

Regarding performance, a large amount of time can be saved if certain
parameters, such as the initial step size, are chosen correctly. Notably,
there appears to be no universally optimal set of parameters, and many of the
results are obtained through the manual tweaking of the parameters for the
specific system. All of the results for this method were obtained on a
standard desktop PC using 4 cores, with the times varying between milliseconds
for the smallest problem and an hour for the largest one. This method has a
significantly reduced memory cost compared to the see-saw method, since here
the $A^{x}_{a}$ matrices are not used, as well as in general offering faster
optimisation.

$n$ $d$ | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---
2 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000
3 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000
4 | 0.01440 | 0.00000 | 0.00000 | 0.00000 | 0.00004
5 | - | 0.00391 | 0.00000 | 0.00000 | -
6 | - | - | 0.00161 | 0.00000 | -
7 | - | - | - | 0.00091 | -
Table 2: The values $1-W_{d,n}/W_{\text{MUB}}(d,n)$, where $W_{d,n}$ is the
result of the optimisation, obtained at an approximate KKT point with
$|\\!|r(\mathbf{x},y,Z,\mu)|\\!|\equiv\sqrt{|\\!|\nabla
L(\mathbf{x},y,Z)|\\!|^{2}+|g(\mathbf{x})|^{2}+|\\!|X(\mathbf{x})Z-\mu
I|\\!|_{F}^{2}}\leq 10^{-5}$ for various dimensions ($d$) and numbers of bases
($n$), to 5 decimal places. The values depicted are consistently obtained by
multiple runs. The values in bold indicate that the value was at least
$10^{-5}$ away from zero, meaning that the method could not find $n$ MUBs in
dimension $d$. Note that in all the known cases, the algorithm predicts
correctly existence/nonexistence, and it predicts that four MUBs do not exist
in dimension six.

##  5 The MUB problem as a ground state problem

$n$ $d$ | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---
2 | 0. | 0. | 0. | 0. | 0. | 0. | 0. | 0. | 0.
3 | 0. | 0. | 0. | 0. | 0. | 0. | 0. | 0. | 0.
4 | 0.01440 | 0. | 0. | 0. | 0.00004215 | 0. | 0. | 0. | 0.000005938
5 | - | 0.003912 | 0. | 0. | - | 0. | 0. | 0. | -
6 | - | - | 0.001609 | 0. | - | 0. | 0. | - | -
7 | - | - | - | 0.0009073 | - | 0. | 0. | - | -
8 | - | - | - | - | - | 0. | 0. | - | -
Table 3: Monte Carlo results. Relative deviation from the MUB optimum in Eq.
(23) $1-W_{\rm max}/W_{\rm MUB}$. “$0.$” indicates that $n$ MUBs have been
found in dimension $d$ up to at least $10^{-10}$ precision. The positive
values in bold are the numerical maxima found over at least three independent
Monte Carlo simulations. For $d=10$, $n=4$, three out of ten independent
simulations found the indicated value; it remains possible that this does not
correspond yet to the true optimum.

###  5.1 Methodology

Finally, we apply an optimisation method inspired by statistical physics,
namely simulated annealing [27]. Given $n$ orthonormal bases
$\\{|b_{j}^{y}\rangle;j\in[d]\\}({y\in[n]})$ in a Hilbert space of dimension
$d$, our goal is to maximise the expression of Eq. (18), which we rewrite here
for completeness (we keep the $d,n$ dependence implicit throughout this
section, and keep the $+B$ superscript implicit in Eq. (18)):

$W[{\bf x}]=\frac{2}{d}\sum_{(y,z)\in{\rm
Pairs}[n]}\sum_{j,k=1}^{d}\sqrt{1-|\langle b_{j}^{y}|b_{k}^{z}\rangle|^{2}},$
(23)

where ${\bf x}=\\{|b_{j}^{y}\rangle\\}$ defines the collection of $n$ bases.
When the $n$ bases are mutually unbiased ($|\langle
b_{j}^{y}|b_{k}^{z}\rangle|^{2}=1/d$ for all $y,z\in\text{Pairs}[n]$ and all
$j,k\in[d]$), we have $W=W_{\rm MUB}=n(n-1)\sqrt{d(d-1)}$. In general, the
maximal value is $W_{\rm max}\leq W_{\rm MUB}$, with equality if and only if
$n$ MUBs exist in dimension $d$. Our strategy is then to optimise the bases by
maximising $W$ (or equivalently, minimising $-W$) via simulated annealing. As
detailed below, this amounts to parametrising the vectors $|b_{j}^{y}\rangle$
by some parameters ${\bf x}$, and regarding $-W({\bf x})$ as the energy of the
configuration ${\bf x}$. The ground state (namely, the lowest-energy
configuration) of $-W$ is found by sampling ${\bf x}$ with probability
$\propto\exp[\beta W({\bf x})]$, progressively ramping up the inverse
temperature $\beta$. In contrast to gradient-descent-based approaches, this
allows to explore a variety of local minima of $-W$ at nonzero temperature,
with the hope of finally converging to the global minimum when the temperature
approaches zero. As our results show, this hope is indeed confirmed by solid
evidence. If $n$ MUBs do exist, they are obtained at the end of the
optimisation, saturating the value $W_{\rm max}=W_{\rm MUB}$. Otherwise, the
algorithm converges to an optimum (presumably the global optimum) $W_{\rm
max}<W_{\rm MUB}$, supporting that $n$ MUBs do not exist in dimension $d$.

Parametrising the bases.– Each vector $|b_{j}^{y}\rangle$ is simply
parametrised by its decomposition in the canonical basis:
$|b_{j}^{y}\rangle=\sum_{i=1}^{d}U_{ji}^{y}|e_{i}\rangle$, where ${\bf
x}:=\\{U^{y}\\}_{y=1}^{n}$ are $n$ complex $d\times d$ unitary matrices. The
overlaps are then obtained as $\langle
b_{j}^{y}|b_{k}^{z}\rangle=\sum_{i=1}^{d}(U_{ji}^{y})^{*}U_{ki}^{z}$, where
$(.)^{\ast}$ is the complex conjugation.

Simulated annealing.– The basic idea of simulated annealing [27] is to
consider:

$\langle W\rangle_{\beta}=\int d{\bf x}\leavevmode\nobreak\ W({\bf
x})\leavevmode\nobreak\ \frac{e^{\beta W({\bf
x})}}{Z_{\beta}}\leavevmode\nobreak\ ,$ (24)

where $Z_{\beta}=\int d{\bf x}\leavevmode\nobreak\ e^{\beta W({\bf x})}$ is a
normalisation factor. Eq. (24) corresponds to an effective thermal average,
where the parameters ${\bf x}$ are sampled from a Gibbs distribution $e^{\beta
W({\bf x})}/Z_{\beta}$, in which $-W({\bf x})$ plays the role of the energy,
and $1/\beta$ is the temperature. We have that [27]:

$\langle W\rangle_{\beta}\underset{\beta\to\infty}{\longrightarrow}W_{\rm
max}\leavevmode\nobreak\ .$ (25)

In words, progressively ramping up the inverse temperature $\beta$, the
thermal average Eq. (24) converges towards the global maximum of the function
$W$. At each value of $\beta$, the Gibbs distribution in Eq. (24) is sampled
via a Markov-chain Monte Carlo algorithm [28].

Monte Carlo sampling.– The Markov chain is a list of samples $\\{{\bf
x}_{i}\\}_{i=1}^{N}$, generated in such a way that in the limit of infinitely
many samples, the average value of $W(\\{{\bf x}_{i}\\})$ over the $N$ samples
converges towards the exact average value, up to $O(1/\sqrt{N})$ corrections:

$\frac{1}{N}\sum_{i=1}^{N}W({\bf x}_{i})=\langle
W\rangle_{\beta}+O(1/\sqrt{N})\leavevmode\nobreak\ .$ (26)

In order to sample ${\bf x}$ according to the Gibbs distribution [Eq. (24)],
we implement a Metropolis algorithm [28]. That is, we start from an arbitrary
initial configuration ${\bf x}_{1}$, and then iterate:

  1. 1.

propose a new configuration ${\bf x}_{\rm new}$ (see below)

  2. 2.

compute the difference $\Delta=\beta[W({\bf x}_{\rm new})-W({\bf x}_{i})]$

  3. 3.

if $\Delta>0$, accept the move

  4. 4.

if $\Delta<0$, accept the move with probability $e^{\Delta}$

  5. 5.

if the move is accepted, update ${\bf x}_{i+1}={\bf x}_{\rm new}$; otherwise
${\bf x}_{i+1}={\bf x}_{i}$.

Implementation of the updates.– In our implementation, we ramp $\beta$
linearly from $\beta_{i}=1$ to $\beta_{f}\approx 10^{4}$–$10^{5}$ in $n_{\rm
steps}=10^{3}$ steps. For each value of $\beta$, we attempt $n_{\rm
attempts}=10^{5}$ Metropolis updates. Each move consists of selecting randomly
one basis among $n$, and rotate randomly all its elements. Specifically, the
moves are proposed as follows

  1. 1.

choose randomly and uniformly one of the bases $y\in[n]$

  2. 2.

draw $2d^{2}$ independent random numbers
$\\{(r_{jk},s_{jk})\\}_{(j,k)\in[d]^{2}}$, uniformly in the interval
$[-\epsilon,\epsilon]$ (see below for the choice of $\epsilon$)

  3. 3.

define $U^{y}_{\rm new}$ as $U_{jk}^{y}+r_{jk}+is_{jk}$

  4. 4.

make $U^{y}_{\rm new}$ unitary via the Gram–Schmidt procedure.

The parameter $\epsilon$, which defines the typical amplitude of the proposed
moves, is adapted throughout the algorithm in order to ensure a constant
acceptance rate. Intuitively, when the temperature is very high, large moves
involving a potentially large change in energy are required to efficiently
explore the parameter space. Progressively ramping down the temperature, the
bases start to stabilise in the vicinity of the maxima of $W$, and large moves
become often rejected by the Metropolis rule. On the other hand, if $\epsilon$
is very small, the moves will not efficiently explore the parameter space. As
a compromise, we adapt $\epsilon$ such that the acceptance rate $r_{\rm
accept}$ of the Metropolis update (that is, for each value of $\beta$, $r_{\rm
accept}$ is the number $N_{\rm accept}$ of accepted moves divided by the
number $N_{\rm attempts}$ of attempts), is kept between $r_{\min}=0.32$ and
$r_{\max}=0.48$. We initialise $\epsilon=1$ at the beginning of the
simulation, and whenever $r_{\rm accept}<r_{\min}$, we change $\epsilon$ to
$0.8\epsilon$. Similarly, whenever $r_{\rm accept}>r_{\max}$, we change
$\epsilon$ to $1.2\epsilon$.

Optimal bases.– As the goal of the simulation is not to accurately estimate
the thermal average $\langle W\rangle_{\beta}$ [Eq. (24)], but only to
efficiently find its global maximum, we do not carry out a detailed evaluation
of the error on $\langle W\rangle_{\beta}$ as estimated from our samples.
Throughout the simulation, we record the optimal bases ${\bf x}_{\rm opt}$
encountered, corresponding to the maximal value of $W_{\rm opt}$ found so far.
As a last step of the optimisation, we set $\beta=\infty$ and start a new
simulation starting from ${\bf x}_{\rm opt}$ as the initial configuration.
Effectively, this amounts to only accepting the moves that increase $W$, in
order to achieve as many digits of precision as needed for $W_{\rm opt}$, as
well as for the basis elements themselves.

Figure 1: Monte Carlo simulation for $d=6$, $n=4$. As a function of the
inverse temperature $\beta$, mean value of the Bell operator [cf. Eq. (24)].
Three independent simulations are plotted, and the dashed red line indicates
the maximal value found throughout the simulations (see text). Inset: Zoom of
the main panel for $1000<\beta<3000.$ Figure 2: Monte Carlo simulation for
$d=10$, $n=4$. As a function of the inverse temperature $\beta$, mean value of
the Bell operator [cf. Eq. (24)]. Ten independent simulations are plotted, and
the dashed red line indicates the maximal value found throughout the
simulations. Inset: Zoom of the main panel for $10000<\beta<20000.$

###  5.2 Results

The results of this optimisation are summarised in Table 3. For a meaningful
comparison of different dimensions $d$, and different number of bases $n$, we
indicate the relative deviation to the MUB optimum: $1-W_{\rm max}/W_{\rm
MUB}$.

Power-of-prime dimensions.– All dimensions $2\leq d\leq 9$ except $d=6$ are
powers of prime, for which the maximal number of MUBs is exactly $d+1$. Our
numerical simulations are consistent with this fact, and our variational
optimisation systematically reconstructs $n$ MUBs for all $n\leq d+1$ 111For
$d=8$ we actually stopped at $n=8$, and for $d=9$ at $n=5$, because of the
high numerical cost of these higher-dimensional instances.. For $d=2,3,4,5$,
we also performed the optimisation for $n=d+2$, where three independent
simulations gave the same optimum within numerical accuracy, also in agreement
with the see-saw and non-linear SDP methods, and with the fact that there
cannot exist $d+2$ MUBs in dimension $d$. The corresponding optimal solutions
are analysed analytically in Section 6 for $d=2$ and 3.

d=6.– In the case $d=6$, we do not find more than $n=3$ MUBs. For $n=4$, three
independent simulations gave consistently the optimum $W_{\rm opt}(6,4)=W_{\rm
MUB}(1-0.00004215)\approx 65.723938549\dots$ (within numerical accuracy). The
optimal bases found also coincide within numerical precision with those found
by the see-saw and non-linear SDP techniques of the previous sections (also
see Section 6.3 for a close analytical construction). The complete evolution
of $\langle W\rangle_{\beta}$ is illustrated in Fig. 1 for these three
simulations.

d=10.– The simulations for $d=10$ are reaching the limits of our current
implementation. We did find $n=3$ MUBs, but not $n=4$. As illustrated in Fig.
2 over ten independent simulations, three of them converged to the same
optimum $W_{\rm opt}(10,4)=W_{\rm MUB}(1-0.000005938)\approx
113.8413197\dots$. This value is reported in Table 3, and is our best estimate
for the true optimum222Using our open source code, Markus Grassl was able to
run the simulation for $d=10$ and $n=4$ on a computer cluster 2000 times. He
reported to us that the most frequent optimum in these runs agrees with our
findings, however, in two instances he found the slightly larger value of
$W_{\rm opt}(10,4)\approx 113.8414358$.. Even though some of our simulations
are trapped in local optima, our results support the conjecture that no more
than three MUBs exist in dimension 10.

##  6 Analytic constructions

In this section we describe analytic constructions that match the best bases
found numerically for the cases of four bases in dimension two, five bases in
dimension three, and four bases in dimension six.

###  6.1 Dimension two

One can parametrise any rank-1 qubit projection $B$ using the Bloch
representation

$B=\frac{1}{2}(\mathbb{1}+\vec{r}\cdot\vec{\sigma}),$ (27)

where $\vec{r}=(x,y,z)$ is a unit vector in $\mathbb{R}^{3}$,
$\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is a vector of the Pauli
matrices, and $\vec{r}\cdot\vec{\sigma}=x\sigma_{x}+y\sigma_{y}+z\sigma_{z}$.
Accordingly, for the case of four rank-1 projective measurements on a qubit,
we use the parametrisation
$B^{y}_{1}=\frac{1}{2}(\mathbb{1}+\vec{r}_{y}\cdot\vec{\sigma})$ for
$y\in\\{1,2,3,4\\}$, and $B^{y}_{2}=\mathbb{1}-B^{y}_{1}$. Consider then the
four measurements parametrised by the four vectors

$\begin{split}\vec{r}_{1}&\left.=\frac{1}{\sqrt{3}}(1,1,1)\right.\\\
\vec{r}_{2}&\left.=\frac{1}{\sqrt{3}}(1,-1,-1)\right.\\\
\vec{r}_{3}&\left.=\frac{1}{\sqrt{3}}(-1,1,-1)\right.\\\
\vec{r}_{4}&\left.=\frac{1}{\sqrt{3}}(-1,-1,1),\right.\end{split}$ (28)

defining a regular tetrahedron on the Bloch sphere. It is straightforward to
verify that by plugging these measurements into Eq. (17) we get

$W^{+B}_{2,4}(\\{B^{y}_{b}\\})=4(\sqrt{3}+\sqrt{6})\approx 16.7262,$ (29)

which agrees (after normalisation) with the values in Tables 1, 2, and 3 up to
numerical precision.

###  6.2 Dimension three

For the three-dimensional case, we parametrise each basis
$y\in\\{1,\ldots,5\\}$ by a unitary matrix $U_{y}$, whose columns correspond
to the basis vectors. Fixing the computational basis allows us to write the
first basis as $U_{1}=\mathbb{1}$. We then take the second basis to be the
Fourier basis

$U_{2}=\frac{1}{\sqrt{3}}\begin{bmatrix}1&1&1\\\
1&\omega_{3}&\omega_{3}^{2}\\\ 1&\omega_{3}^{2}&\omega_{3}\end{bmatrix},$ (30)

where $\omega_{3}=\mathrm{e}^{\frac{2\pi\mathrm{i}}{3}}$. We define the
remaining three bases as

$\begin{split}U_{3}&\left.=\operatorname{diag}(\omega_{6},\bar{\omega}_{6},1)U_{2}\right.\\\
U_{4}&\left.=\operatorname{diag}(1,\omega_{6},\bar{\omega}_{6})U_{2}\right.\\\
U_{5}&\left.=\operatorname{diag}(\bar{\omega}_{6},1,\omega_{6})U_{2},\right.\end{split}$
(31)

where $\omega_{6}=\mathrm{e}^{\frac{2\pi\mathrm{i}}{6}}$, $\bar{\omega}_{6}$
is its complex conjugate, and $\operatorname{diag}(x,y,z)$ is a diagonal
matrix of $x$, $y$, and $z$. Note that all the bases $U_{2}$, $U_{3}$,
$U_{4}$, and $U_{5}$ are unbiased to $U_{1}$. It is straightforward to verify
that by plugging the measurements $B^{y}_{b}$ corresponding to the bases
$U_{y}$ into Eq. (17), we get

$W^{+B}_{3,5}(\\{B^{y}_{b}\\})=8(\sqrt{2}+\sqrt{5}+\sqrt{6})\approx 48.7982,$
(32)

which agrees (after normalisation) with the values in Tables 1, 2, and 3 up to
numerical precision.

###  6.3 Dimension six

All three numerical methods converged to (numerically) the same set of four
bases in dimension six. Just like in the case of five bases in dimension
three, this set has the property that one basis is unbiased to the other three
bases. Upon closer inspection, one finds that these bases are numerically very
close to the “four most distant bases” of Ref. [21], which were found via
maximising the MUBness measure

$D^{2}\propto\sum_{(y,z)\in{\rm Pairs}[n]}\sum_{j,k}|\langle
b_{j}^{y}|b_{k}^{z}\rangle|^{2}(1-|\langle
b_{j}^{y}|b_{k}^{z}\rangle|^{2})\leavevmode\nobreak\ .$ (33)

One can again parametrise the bases by four unitary matrices such that the
first one is $U_{1}=\mathbb{1}$. The rest of the bases are based on Eq. (6) of
Ref. [21], in which the authors describe a family of three unitary matrices,
depending on two parameters, $\theta_{t}$ and $\theta_{x}$. The “optimal”
values of these parameters (optimality here originally means maximising the
distance measure of Ref. [20]) is then determined by finding the unique real
solution of Eq. (20) in Ref. [21], plugging it into Eq. (19) of Ref. [21] to
obtain $\theta_{t}$, and plugging it into Eq. (21) of Ref. [21] to obtain
$\theta_{x}$. These analytic values correspond to approximately
$(\theta_{x},\theta_{t})\approx(0.9852276,1.0093680)$.

If we substitute the resulting bases $\\{B^{y}_{b}\\}$ into our optimisation
problem in Eq. (19), we obtain (analytically)

$W_{6,4}^{+B}(\\{B^{y}_{b}\\})\approx 65.7239381$ (34)

Comparing with our optimum $\approx 65.7239385$, we observe that: 1) the
analytical optimum (for the bases of Ref. [21]) and our numerical optimum
differ after the eighth significant figure; and 2) consistently, our numerical
optimum is larger than the analytical one. This should not be a surprise, for
the two solutions optimise different MUBness measures, respectively Eq. (18)
and Eq. (33). As a further comparison, one can look at the overlaps
$\operatorname{tr}(B^{y}_{j}B^{z}_{k})=|\langle b^{y}_{j}\hskip 0.7pt|\hskip
0.7ptb^{z}_{k}\rangle|^{2}$ of the bases found. We indeed find that one of the
bases is unbiased to the other three, i.e.

$\operatorname{tr}(B^{1}_{j}B^{y}_{k})=|\langle b^{y}_{j}\hskip 0.7pt|\hskip
0.7ptb^{z}_{k}\rangle|^{2}=0.166667\quad\forall j,y,k.$ (35)

Apart from these, there are three different values of overlaps, whose
approximate values $0.124$, $0.181$, and $0.152$ agree with our numerical
findings up to two or three significant digits.

In conclusion, our three numerical methods, together with the approach of Ref.
[21] all appear to converge to essentially the same set of four bases in
dimension six. Small differences in the overlaps originate in the different
MUBness measures which are optimised. The results of these four independent
different numerical approaches may signify that Zauner’s conjecture is indeed
correct, although they cannot rule out the contrary.

##  7 Conclusions

We reformulated the existence problem of MUBs as an optimisation problem,
using a recently found family of Bell inequalities. We then applied three
numerical methods suitable for optimising Bell inequalities in order to tackle
the existence problem: see-saw SDP optimisation, non-linear SDP, and Monte
Carlo techniques. The results of all these numerical optimisations is in full
accordance with the known cases in dimensions $d=2,3,4,5$, where we find $d+1$
MUBs. Furthermore, whenever it is known that $n$ MUBs do not exist in a given
dimension $d$, all the different algorithms converge to the same set of bases
in all cases (with a slight difference between the see-saw method and the
other two methods for $d=4$, $n=6$), and these bases are not MUBs.

We applied our numerical techniques to the open case of four MUBs in dimension
six. All three algorithms suggest that there do not exist four MUBs in
dimension six, by converging to a Bell value strictly smaller than the
hypothetical MUB value. Moreover, the bases found by all three algorithms are
very close numerically with the “four most distant bases” in dimension six of
Ref. [21]. Hence, our findings provide further numerical evidence for Zauner’s
conjecture. In the next composite dimension, $d=10$, our Monte Carlo results
suggest that no more than $n=3$ MUBs exist.

It is important to point out that the numerical methods used in this work are
heuristic, i.e. there is no guarantee of convergence to the global optimum. As
such, heuristic numerics can never provide a rigorous proof of the non-
existence of MUBs (only that of existence, by explicitly finding MUBs). To
overcome this shortcoming, a plausible future direction towards a rigorous
numerical proof is using a variant of the Navascués–Pironio–Acín hierarchy of
SDPs [29] for maximising a Bell inequality in a fixed dimension. While such
numerical optimisation is computationally significantly more expensive than
those in our work, it provides certifiable upper bounds on Bell inequality
violations, and therefore could in principle be used to rigorously prove the
non-existence of MUBs.

Code availability. The numerical findings presented in this paper can be
reproduced using the codes made available on public repositories. For the see-
saw SDP, see https://github.com/Lumorti/seesaw. For the non-linear SDP, see
https://github.com/Lumorti/nonlinear. For the Monte Carlo simulations, see
https://github.com/mariaprat/mubs-montecarlo.

##  8 Acknowledgements

The authors thank Markus Grassl for providing a large set of data obtained by
running our code on a computer cluster. 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. 754510, the Government of Spain
(FIS2020-TRANQI, Severo Ochoa CEX2019-000910-S), Fundacio Cellex, Fundacio
Mir-Puig, Generalitat de Catalunya (CERCA, AGAUR SGR 1381 and QuantumCAT), the
ERC AdG CERQUTE and the AXA Chair in Quantum Information Science, the Agence
Nationale de la Recherche (Qu-DICE project ANR-PRC-CES47), the John Templeton
Foundation (Grant No. 61835), ICFO and SPIE under a Maria Yzuel Fellowship
Award.

## References

  * [1] I D Ivanovic.  Geometrical description of quantal state determination.  Journal of Physics A: Mathematical and General, 14(12):3241–3245, 1981.  doi:10.1088/0305-4470/14/12/019.
  * [2] G. Brassard C. H. Bennett.  Quantum cryptography: Public key distribution and coin tossing.  Proceedings of IEEE International Conference on Computers, Systems and Signal Processing (IEEE, 1984), 175:8, 1984.  doi:10.1016/j.tcs.2011.08.039.
  * [3] Artur K. Ekert.  Quantum cryptography based on Bell’s theorem.  Phys. Rev. Lett., 67:661–663, 1991.  doi:10.1103/PhysRevLett.67.661.
  * [4] Dagmar Bruß.  Optimal eavesdropping in quantum cryptography with six states.  Phys. Rev. Lett., 81:3018–3021, 1998.  doi:10.1103/PhysRevLett.81.3018.
  * [5] Armin Tavakoli, Alley Hameedi, Breno Marques, and Mohamed Bourennane.  Quantum random access codes using single $d$-level systems.  Phys. Rev. Lett., 114:170502, 2015.  doi:10.1103/PhysRevLett.114.170502.
  * [6] Máté Farkas and Jędrzej Kaniewski.  Self-testing mutually unbiased bases in the prepare-and-measure scenario.  Phys. Rev. A, 99:032316, 2019.  doi:10.1103/PhysRevA.99.032316.
  * [7] H. Bechmann-Pasquinucci and N. Gisin.  Bell inequality for qunits with binary measurements.  Quantum Info. Comput., 3(2):157–164, 2003.  doi:10.26421/QIC3.2-6.
  * [8] Jędrzej Kaniewski, Ivan Šupić, Jordi Tura, Flavio Baccari, Alexia Salavrakos, and Remigiusz Augusiak.  Maximal nonlocality from maximal entanglement and mutually unbiased bases, and self-testing of two-qutrit quantum systems.  Quantum, 3:198, 2019.  doi:10.22331/q-2019-10-24-198.
  * [9] Armin Tavakoli, Máté Farkas, Denis Rosset, Jean-Daniel Bancal, and Jędrzej Kaniewski.  Mutually unbiased bases and symmetric informationally complete measurements in Bell experiments.  Science Advances, 7(7):eabc3847, 2021.  doi:10.1126/sciadv.abc3847.
  * [10] Thomas Durt, Berthold-Georg Englert, Ingemar Bengtsson, and Karol Życzkowski.  On mutually unbiased bases.  International Journal of Quantum Information, 08(04):535–640, 2010\.  doi:10.1142/S0219749910006502.
  * [11] William K Wootters and Brian D Fields.  Optimal state-determination by mutually unbiased measurements.  Annals of Physics, 191(2):363–381, 1989.  doi:10.1016/0003-4916(89)90322-9.
  * [12] Paweł Wocjan and Thomas Beth.  New construction of mutually unbiased bases in square dimensions.  Quantum Info. Comput., 5(2):93–101, 2005.  doi:10.26421/QIC5.2-1.
  * [13] Mihály Weiner.  A gap for the maximum number of mutually unbiased bases.  Proc. Amer. Math. Soc., 141:1963–1969, 2013.  doi:10.1090/S0002-9939-2013-11487-5.
  * [14] Gerhard Zauner.  Quantendesigns: Grundzüge einer nichtkommutativen Designtheorie.  PhD thesis, 1999.
  * [15] P. Oscar Boykin, Meera Sitharam, Pham Huu Tiep, and Pawel Wocjan.  Mutually unbiased bases and orthogonal decompositions of Lie algebras.  Quantum Info. Comput., 7(4):371–382, 2007.  doi:10.26421/QIC7.4-6.
  * [16] Stephen Brierley and Stefan Weigert.  Constructing mutually unbiased bases in dimension six.  Phys. Rev. A, 79:052316, 2009.  doi:10.1103/PhysRevA.79.052316.
  * [17] Philippe Jaming, Máté Matolcsi, Péter Móra, Ferenc Szöllősi, and Mihály Weiner.  A generalized Pauli problem and an infinite family of MUB-triplets in dimension 6.  Journal of Physics A: Mathematical and Theoretical, 42(24):245305, may 2009.  doi:10.1088/1751-8113/42/24/245305.
  * [18] Gary McConnell, Harry Spencer, and Afaq Tahir.  Evidence for and against Zauner’s MUB conjecture in $\mathbb{C}^{6}$.  2021\.  doi:10.48550/arXiv.2103.08703.
  * [19] Sander Gribling and Sven Polak.  Mutually unbiased bases: polynomial optimization and symmetry.  2021\.  doi:10.48550/arXiv.2111.05698.
  * [20] Ingemar Bengtsson, Wojciech Bruzda, Åsa Ericsson, Jan-Åke Larsson, Wojciech Tadej, and Karol Życzkowski.  Mutually unbiased bases and Hadamard matrices of order six.  Journal of Mathematical Physics, 48(5):052106, 2007.  doi:10.1063/1.2716990.
  * [21] Philippe Raynal, Xin Lü, and Berthold-Georg Englert.  Mutually unbiased bases in six dimensions: The four most distant bases.  Phys. Rev. A, 83:062303, 2011.  doi:10.1103/PhysRevA.83.062303.
  * [22] Edgar A. Aguilar, Jakub J. Borkała, Piotr Mironowicz, and Marcin Pawłowski.  Connections between mutually unbiased bases and quantum random access codes.  Phys. Rev. Lett., 121:050501, 2018.  doi:10.1103/PhysRevLett.121.050501.
  * [23] Nicolas Brunner, Daniel Cavalcanti, Stefano Pironio, Valerio Scarani, and Stephanie Wehner.  Bell nonlocality.  Rev. Mod. Phys., 86:419–478, 2014.  doi:10.1103/RevModPhys.86.419.
  * [24] MOSEK ApS.  MOSEK Fusion API for C++ 9.2.49, 2021.  URL: https://docs.mosek.com/9.2/cxxfusion/index.html.
  * [25] Hiroshi Yamashita, Hiroshi Yabe, and Kouhei Harada.  A primal–dual interior point method for nonlinear semidefinite programming.  Mathematical programming, 135(1):89–121, 2012.  doi:10.1007/s10107-011-0449-z.
  * [26] Stephen Boyd and Lieven Vandenberghe.  Convex Optimization.  Cambridge University Press, 2004.  doi:10.1017/CBO9780511804441.
  * [27] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi.  Optimization by simulated annealing.  Science, 220(4598):671–680, 1983.  doi:10.1126/science.220.4598.671.
  * [28] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  The Journal of Chemical Physics, 21(6):1087–1092, 1953.  doi:10.1063/1.1699114.
  * [29] Miguel Navascués, Stefano Pironio, and Antonio Acín.  Bounding the set of quantum correlations.  Phys. Rev. Lett., 98:010401, 2007.  doi:10.1103/PhysRevLett.98.010401.


# Higher-order exceptional points using lossfree negative-index materials

Xin-Zhe Zhang1, Li-Ting Wu2, Ru-Zhi Luo1, and Jing Chen1,3
<EMAIL_ADDRESS>$1$ MOE Key Laboratory of Weak-Light Nonlinear Photonics,
School of Physics, Nankai University, Tianjin 300071, China
$2$ School of Information and Communication Engineering, Nanjing Institute of
Technology, Nanjing 211167, China
$3$ Collaborative Innovation Center of Extreme Optics, Shanxi University,
Taiyuan, Shanxi 030006, China

###### Abstract

Negative-index materials (NIMs) support optical anti-parity-time
(anti-$\mathcal{PT}$) symmetry even when they are lossless. Here we prove the
feasibility in achieving higher-order exceptional points (EPs) in lossfree
waveguide arrays by utilizing the anti-$\mathcal{PT}$ symmetry induced by NIM.
Numerical simulation about a third-order EP fits well with the coupled-mode
theory. A scheme of achieving fourth-order EPs is also discussed. This work
highlights the potential of lossfree NIMs in the study of non-Hermitian
optics.

123

###### pacs:

123

##  I Introduction

Exceptional points (EPs) refer to the singular degeneracies of non-Hermitian
wave/quantum systems R01 ; R02 ; R03 ; R04 ; R05 ; R06 ; R07 , where all the
eigenvalues and the eigenvectors of the effective Hamiltonian coalesce
simultaneously. Its novel topology enables interesting mode switching
behaviors when circling around it R06 ; R07 ; R08 ; R09 ; R10 . Furthermore,
the coalescent eigenfunction at EPs is very sensitive to tiny perturbation in
the effective Hamiltonian, so high-sensitivity applications R11 ; R12 ; R13 ;
R14  can be envisioned. The high sensitivity is also associated with stopped-
light effect and enhanced density of states, and can be utilized for other
interesting purposes such as in realizing coherent absorbers and lasers R15 ;
R16 ; R17 .

In most cases only the second-order EPs (EP2 for brief) are considered because
EPs are generally believed to exist mainly in open systems. A good example is
that in a parity-time ($\mathcal{PT}$) symmetric system that asks for delicate
balance among spatially distributed gain and loss R01 ; R02 . If the order of
EP is increased, it would provide more degrees of freedoms in artificially
designing the topology, and the sensitivity can be further enhanced.
Consequently, people have proposed many schemes in achieving higher-order EPs,
e.g. in various kinds of photonic crystals, microcavities, and waveguides
(WGs) R11 ; R12 ; R18 ; R19 ; R20 ; R21 ; R22 ; R23 ; R24 ; R25 ; R26 ; R27 ;
R28 ; R29 ; R30 . However, the presence of gain and loss still dramatically
hinders the transfer of EPs from a scientific curiosity to realistic
applications of our daily life. Ways to access high-order EPs without
resorting to gain and loss are thus desired. Such a target is, in principle,
achievable because non-Hermitian physics covers many miscellaneous categories
including but not limited to the $\mathcal{PT}$ symmetry. In addition, to
achieve a non-Hermitian Hamiltonian, besides introducing imaginary components
(gain and loss) to the diagonal elements, we can also just set the off-
diagonal elements unequal. The later method does not require gain and loss,
and the whole energy can be conserved. A good example is the recently
demonstrated anti-$\mathcal{PT}$ symmetry associated with lossfree negative-
index materials (NIMs) R31 ; R32 ; R33 ; R34 ; R35 ; R36 . Albeit the energy
in the lossfree NIM-dielectric system is conserved R35 ; R36 , due to the
backward propagation of field in NIM, the coupling between a NIM WG and an
ordinary dielectric WG can be modeled by an anti-$\mathcal{PT}$ Hamiltonian of
unequal off-diagonal elements. EP2s are shown to exist in the spectra of the
guided waves R35 ; R36 .

In this article, we would prove the feasibility in realizing an optical third-
order EP (EP3) by using the anti-$\mathcal{PT}$ symmetry induced by lossfree
NIMs. An effective Hamiltonian is developed by using the coupled-mode theory
in order to explain the existence of EP3. The transfer-matrix method (TMM) is
utilized to numerically calculate eigensolutions of the guided mode, prove the
existence of EP3, and reveal features of it. A scheme of realizing a fourth-
order EP (EP4) is also discussed. This work proves that NIMs and the
associated anti-$\mathcal{PT}$ symmetry have great potential in the study of
non-Hermitian optics in lossfree environments for various applicable purposes.

Figure 1:  (Color online) Schematic of the configuration under investigation,
which contains three straight lossfree WGs. WG1 is made of a NIM with
$\epsilon_{\rm NIM}=-3$ and $\mu_{\rm NIM}=-0.556$. WG2 and WG3 are made of
dielectrics.

This article is organized as follows. In Section 2.1 we firstly propose the
main concept of the coupled-WG structure and the physical mechanism of EP3 by
using coupled-mode theory. In Section 2.2 we provide numerical calculation and
analysis about the guided modes by using TMM. We show that the Hamiltonian
from the coupled-mode theory can explain main features of the results from
TMM, and proves the observed degeneracy is indeed an EP3. The analysis also
provides more detailed information about how the eigenmodes evolve near EP3.
Discussion about the importance of this study is provided in Section 3. We
also present a simple scheme of achieving an EP4 in this Section. Summary is
made at the end of this article.

##  II Theory and Analysis

###  II.1 Structure and Effective Hamiltonian

Let us consider the structure shown in Fig. 1. It contains three straight WGs
surrounded by air. All the media in this structure are lossfree. The lower two
WGs (WG2 and WG3) are made of dielectrics with $\epsilon>0$ and $\mu=1$. The
top WG1 is made of a NIM with $\epsilon_{\rm NIM}<0$ and $\mu_{\rm NIM}<0$.
Because NIM requires an intrinsic dispersion of $\partial(\epsilon_{\rm
NIM}\omega)/\partial\omega>0$ and $\partial(\mu_{\rm
NIM}\omega)/\partial\omega>0$ so as to give a positive energy density R31 ;
R32 ; R33 ; R36 , in this article we would keep the angular frequency
$\omega_{0}$ a constant, and test the variation of the wavevectors $\beta$ of
the eigenmodes versus a geometric parameter of the structure.

As been discussed in R35 ; R36 , because the total energy should be conserved
in this loss-free system, and the propagating directions of energy in the NIM
and dielectric WGs are opposite to each other, the coupling between the top
NIM WG1 and the adjacent dielectric WG2 is anti-$\mathcal{PT}$ symmetric R35 ;
R36 . As for WG2 and WG3, their interaction can be described by using a
Hermitian matrix. Only keeping the nearest-neighbor interaction, the coupled-
mode theory gives an effective Hamiltonian on the eigensolutions $\beta$ of
the guided mode

$\left[\begin{array}[]{ccc}k_{\rm WG1}&-A&0\\\ A&k_{\rm WG2}&B\\\ 0&B&k_{\rm
WG3}\\\ \end{array}\right]\left[\begin{array}[]{ccc}\psi_{1}\\\ \psi_{2}\\\
\psi_{3}\\\ \end{array}\right]=\beta\left[\begin{array}[]{ccc}\psi_{1}\\\
\psi_{2}\\\ \psi_{3}\\\ \end{array}\right],$ (1)

where $k_{\rm WG\it{i}}$ is the resonant wavevector of mode in separate WG$i$
($i=1,2,3$), and $\psi_{i}$ represents an associated field component of it.
Parameter $A$ represents the strength of the anti-$\mathcal{PT}$ coupling, and
$B$ is the Hermitian coupling strength between WG2 and WG3. Both $A$ and $B$
are real.

Equation (1) generally has three solutions. Here let us assume the resonances
in WG2 and WG3 are degenerated but different from that of WG1,

$k_{\rm WG1}=k_{0}+3\delta,k_{\rm WG2}=k_{\rm WG3}=k_{0},$ (2)

where the factor $3$ before $\delta$ is intentionally introduced in order to
make below analysis concise. Substituting them into Eq. (1) and assume

$y=\beta-k_{0}-\delta,$ (3)

the three solutions of $\beta$ can be found by solving

$y^{3}+C_{1}y+C_{0}=0,$ (4)

where the coefficients are given by

$C_{1}=A^{2}-B^{2}-3\delta^{2},C_{0}=(A^{2}+2B^{2}-2\delta^{2})\delta.$ (5)

Because $A$, $B$, $k_{0}$ and $\delta$ are all real, $C_{0}$ and $C_{1}$ are
also real valued.

Assuming the three solutions are $y_{1,2,3}$, Eq. (4) has many general
properties such as $y_{1}+y_{2}+y_{3}=0$ and $y_{1}y_{2}y_{3}=-C_{0}$. A
notable feature is that it can support an EP3 with three identical solutions
of $y=0$ when

$C_{1}=C_{0}=0.$ (6)

Equation (6) can be satisfied only when the conditions of

$\delta=0,A^{2}=B^{2}$ (7)

are met simultaneously. These are the existence conditions of EP3, the main
conclusion of this article.

Because parameters $A$ and $B$ are tunable by managing the distances between
adjacent WGs, to gain a deep insight about the formation of EP3 and the
variation of the associated eigenvectors, let us check how the eigensolutions
vary about EP3. Substituting $k_{WG1}=k_{WG2}=k_{WG3}=k_{0}$ into Eq. (1) we
can get

$(\beta-k_{0})^{3}+(\beta-k_{0})(A^{2}-B^{2})=0.$ (8)

Assuming all the solutions are real, we can sort them in ascending order. The
solution $\beta_{2}$ is between the other two and is given by

$\beta_{2}=k_{0},$ (9)

which is a constant and does not depend on the values of $A$ or $B$. The other
two solutions are given by

$\beta_{1,3}=k_{0}\mp\sqrt{B^{2}-A^{2}},$ (10)

which are complex (real) in the region of broken (exact) phase when
$B^{2}<A^{2}$ ($B^{2}>A^{2}$). As for the eigenvectors, from Eq. (1) we can
find

$\Psi_{2}=\frac{1}{\sqrt{A^{2}+B^{2}}}\left[B,0,-A\right]^{T},\Psi_{1,3}=\frac{1}{\sqrt{2}B}\left[A,\pm\sqrt{B^{2}-A^{2}},-B\right]^{T},$
(11)

respectively.

Once the condition of $A^{2}=B^{2}$ is satisfied, a degeneracy takes place,
where all the three eigensolutions coalesce together. This degeneracy is an
EP3, and the eigensolution and eigenvector are given by

$\beta_{\rm EP3}=k_{0},\Psi_{\rm EP3}=\frac{1}{\sqrt{2}}\left[1,0,-{\rm
sign}\left(\frac{B}{A}\right)\right]^{T},$ (12)

respectively, where the function sign($x$) returns 1 ($-1$) when $x>0$
$(x<0)$.

If the two conditions of Eq. (6) are not satisfied simultanously, e. g. when
$\delta\neq 0$, we could not access EP3. Now the system at most supports an
EP2. One solution of $\beta$ is always real, and the other two solutions are
determined by

$\Delta=\left(\frac{C_{0}}{2}\right)^{2}+\left(\frac{C_{1}}{3}\right)^{3}.$
(13)

When $\Delta=0$, the two solutions are identical and an EP2 is achieved. When
moving from $\Delta=0$, for example, by changing the distance between two
adjacent WGs so as to modify $A$ or $B$, the system would enter either the
exact or the broken non-Hermitian phase. If $\Delta>0$, two complex conjugate
solutions are achieved, and the system is within the broken phase. If
$\Delta<0$, two real solutions with different values are found, and the system
is within the exact phase. In the next section we would demonstrate this
phenomenon.

Figure 2:  (Color online) Variation of the eigensolutions $\beta$ versus the
distance $a$ by using TMM. (a,c) In general we can always get an EP2, where
the other dispersion curve is real and always persists. (b, d) By finely
tuning parameters of the structure we can get an EP3. Red circles are the best
fitted results by using Eq. (1).

###  II.2 Numerical Simulation and Analysis

Above theory about the eigensolutions and eigenvectors is based on the
effective Hamiltonian from the coupled-mode theory. When studying the guided
mode in coupled WGs we should still resort to methods based on Maxwell’s
equations. Here we analyze the transverse-electrical modes in the structure by
using TMM. In each layer the field is expressed as
$E_{y}=[E_{+}\exp(jk_{i}x)+E_{-}\exp(-jk_{i}x)]\exp(-j\beta z)$, where
$k_{i}^{2}+\beta^{2}=\epsilon_{i}\mu_{i}\omega_{0}^{2}/c^{2}$, $c$ is the
speed of light. The transmission/reflection properties of the structure are
summarized by $[E_{t},0]^{T}=\bm{M}[E_{0},E_{r}]^{T}$, where $E_{0}$, $E_{t}$
and $E_{r}$ are the incident, transmitted, and reflected fields, $\bm{M}$ is
the transfer matrix. At the top and bottom air media surrounding the WGs,
field is required to exponentially decay away ($\beta>\omega_{0}/c$) so that
out-going boundary conditions are met. Wave-guiding is satisfied by the self-
sustained condition of $M_{22}=0$ R37 , so that $E_{r}$ and $E_{t}$ are
nonzero even when no incidence is present ($E_{0}=0$). In the phase-broken
regions where $\beta$ is complex, similar process can be utilized as well by
scanning $\beta$ in the two dimensional space spanned by ${\rm Re}\\{\beta\\}$
and ${\rm Im}\\{\beta\\}$. As for the associated distribution of fields, they
can be calculated also by TMM when an eigensolution is found.

The geometric and optical parameters of the structure are set as follows. The
thicknesses of all WGs are 4 cm. The distance between WG2 and WG3 is set to be
6 cm so that $B$ is a constant here. The distance $a$ between WG1 and WG2,
which determines the value of $A$, is variable in our study. NIM is assumed to
be the documented one with $\epsilon_{\rm NIM}=1-\omega_{e}^{2}/\omega^{2}$
and $\mu_{\rm NIM}=1-F\omega^{2}/(\omega^{2}-\omega_{m}^{2})$, where
$\omega_{e}=2\pi\times 10$ GHz, $\omega_{m}=2\pi\times 4$ GHz, and $F=0.56$
R32 ; R33 . At $\omega_{0}=2\pi\times 5$ GHz (a free-space wavelength of 6 cm)
the dispersion gives $\epsilon_{\rm NIM}=-3$ and $\mu_{\rm NIM}=-0.556$.

Figs. 2(a) and 2(c) display the spectra of $\beta$ verses $a$ when WG2 and WG3
are made of the same dielectric of $\epsilon=4.504$. In this case $k_{\rm
WG2}=k_{\rm WG3}\neq k_{\rm WG1}$ so we could only get an EP2. From the curves
we can see the lower branch $\beta_{1}$ is always real, but the upper two
branches of $\beta_{2}$ and $\beta_{3}$ coalesce together in forming an EP2
when $a=7.09$ cm. When $a$ is smaller than 7.09 cm, the imaginary parts of
these two branches are no longer zero, and the non-Hermitian phase is broken.

In order to realize EP3, we should finely tune the parameters of the
structure. We find that when $\epsilon$ of the two dielectric WGs are modified
to 4.50465083, an EP3 is achieved, see Figs. 2(b) and 2(d). From Fig. 2(b) we
can see just this tiny $0.015\%$ tuning would sharply modify the dispersion
curves. Now the initial lower branch $\beta_{1}$ coalesce with the upper two
branches $\beta_{2,3}$ at a single point at $a_{c}=6.69$ cm and forms an EP3.
When $a$ becomes smaller, one branch of $\beta$ is always real, while the
other two become complex. The middle branch is always real and is almost a
constant against $a$. Note that due to the coalescence at EP3, it is
impossible to find the exact one-to-one correspondence among the three
branches of $\beta$ in the regions of broken and exact phases, so the complex
branches are still labeled in ascending order of $\rm Im\\{\beta\\}$.

Above numerical calculations are based on TMM, which is rigorous. Before
discussing characteristics of EP3 we should first check how they fit with the
coupled-mode theory of Eq. (1). To do it, we perform a best fitting about the
dispersion curves by numerical solving Eq. (1). We first fit the curves shown
in Figs. 2(b) and 2(d) because they are govern by $\delta=0$ and require less
fitting parameters. The value of $k_{0}$ is given by the numerical value of
$\beta_{2}=1.325144$ cm-1 obtained from TMM. As for A and B, following R36  we
assume

$A=B\exp\left(-\frac{a-a_{c}}{L_{d}}\right),$ (14)

where $a_{c}=6.69$ cm is the position of EP3, and $L_{d}$ is the decay length.
The best fitting then requires $B=1.22\times 10^{-3}$ cm-1 and $L_{d}=1.3$ cm.
From the red circles of Figs. 2(b) and 2(d) we can see the coupled-mode theory
predicts most features of the dispersion curves from TMM. The only minor
discrepancy is about $\rm Re\\{\beta\\}$ in the broken phase region, which
might be fixed by considering high-order dispersion effect.

We then try to fit the results shown in Figs. 2(a) and 2(c). Since the only
difference of them from Figs. 2(b) and 2(d) is that the resonances in WG2 and
WG3 are perturbed by the same amount, we would maintain all the parameters
about Figs. 2(b) and 2(d) unchanged, and only introduce a tiny perturbation to
$k_{\rm WG2,WG3}$, that $k_{\rm WG2}=k_{\rm WG3}=k_{0}+\Delta k$. A best
fitting is achieved when $\Delta k=-1.8\times 10^{-4}$ cm-1, as shown as the
red circles. We can see once again they fit well with these from TMM. It
confirms that the coupled-mode theory of Eq. (1) is an accurate model about
the guided modes in the coupled WGs.

Figure 3:  (Color online) Distributions of field $E_{y}$ at (a) EP2 and (b)
EP3 shown in Fig. 2.

Now we can try to analyze the formation of EP3 and the associated
characteristics of it. We first calculate the distributions of $E_{y}$ at EP2
and EP3 of Fig. 2, and show them in Fig. 3. We can see EP3 is sharply
difference from EP2 because the field inside WG2 is nearly zero, which is in
agreement with Eq. (12). The patterns of fields in WG1 and WG3 are similar in
these two scenarios.

Figure 4:  (Color online) Distributions of field $E_{y}$ at the chosen points
of Figs. 2(b) and Fig. 2(d) supporting EP3.

We then pay attention to the structure supporting EP3, and calculate the
distribution of $E_{y}$ at six points around EP3 in the dispersion curves of
Fig. 2(b) and 2(d). The result are shown in Fig. 4. For the three points
beyond EP3 ($a>a_{c}$), all the solutions of $\beta$ are real. As a
consequence, the field $E_{y}$ is real, see Figs. 4(a) to 4(c). A notable
feature is that for the point $p2$ in the middle branch $\beta_{2}$, the field
inside WG2 is very weak, similar to EP3. This phenomenon is also in agreement
with Eq. (11) .

When $a<a_{c}$, two branches of $\beta$ becomes complex. Now the field $E_{y}$
becomes complex, so in Figs. 4(d) to 4(f) we also plots the imaginary parts of
it. From the distribution of fields we can see in these broken phase scenario,
the patterns of the real part of $E_{y}$ are almost identical with each other
(also with that of EP3). The differences in the distribution of $E_{y}$ are
carried by the imaginary part. For the two complex-$\beta$ points, the
imaginary parts of $E_{y}$ are opposite with each other and is especially
strong around WG2. This phenomenon can be explained by using Eq. (11), that in
the region of broken phase, the complex amplitude of the basic vector
$\psi_{2}$ in $\Psi_{1,3}$ is purely complex and possesses a $\pi/2$ or
$3\pi/2$ phase difference from these of the basic vectors $\psi_{1,3}$.

To prove that at EP3 the three dispersion branches $\beta$ coalesce together,
we test whether their eigenvectors follow the prediction of Eq. (11). It can
be done by checking how the relative magnitude of field inside WG2 varies with
$a$. Here we assume that $\psi_{i}$ in the eigenvector represents $E_{y}$
inside WG$i$ ($i=1,2,3$), calculate the integral of $|E|^{2}$ inside each WG,
and find the ratio $\alpha$ of $|E|^{2}$ confined in WG2 by using

$\alpha=\frac{\int_{\rm WG2}|E_{y}|^{2}dx}{\int_{\rm
WG1}|E_{y}|^{2}dx+\int_{\rm WG2}|E_{y}|^{2}dx+\int_{\rm WG3}|E_{y}|^{2}dx}.$
(15)

Albeit this approach is very rough, e.g. because the fields outside WGs are
ignored, the variation of $\alpha$ versus $a$ shown in Fig. 5(a) agrees well
with the dependence of $|\psi_{2}|^{2}$ on $A$ given by Eq. (11). For modes in
the middle branch $\beta_{2}$, $\alpha$ is nearly zero, which is in agreement
with the null $\psi_{2}$ given by Eq. (12). As for the other two branches,
when $a$ is greater than $a_{c}$, the magnitude of $A$ is much smaller than
$B$, i.e. $|A|\ll|B|$. The eigenvectors are close to $[0,\pm 1,1]^{T}$, and
$\alpha$ approaches $50\%$. When $a$ decreases so that $A$ becomes comparable
with $B$, $\alpha$ decreases sharply to zero. At $a_{c}$ where $A=B$, as
expected, a coalescence takes place. Because the curves of $\alpha$ are
continuously varying, we can make the conclusion that this point at $a=a_{c}$
is not an accidental degeneracy of EP2 with other modes but a standard EP3.

Figure 5:  (Color online) (a) Ratio of field localized inside WG2, where the
$y$ axis is plotted logarithmic. (b) Group velocities $v_{g}$ of the three
branches. Results about complex $\beta$ are not shown because the associated
$\alpha$ and $v_{g}$ are not well defined.

Since EPs in $\mathcal{PT}$-symmetric WGs stop light R15 ; R16 , we also
calculate the group velocity $v_{g}$ of the three branches. Here $v_{g}$ is
given by the ratio of the Poynting vector $S_{z}=-\int E_{y}H_{x}^{*}dx/2$ to
the energy density $W=\int wdx$ via the formula of $v_{g}=S_{z}/W$. When
calculating the energy density we have adapted
$w=\epsilon_{0}\partial(\epsilon_{\rm
NIM}\omega)/\partial\omega|E|^{2}/4+\mu_{0}\partial(\mu_{\rm
NIM}\omega)/\partial\omega|H|^{2}/4$ in order to guarantee a positive energy.
At 5 GHz, the utilized dispersion of NIM gives $\partial(\epsilon_{\rm
NIM}\omega)/\partial\omega=5$ and $\partial(\mu_{\rm
NIM}\omega)/\partial\omega=4.975$.

Figure 5(b) shows the variation of $v_{g}$ versus $a$ of the three branches
$\beta$. When $a>a_{c}$ so that $A<B$, the $\beta_{2}$ branch possesses a
negative group velocity. The reason is that for modes in this branch, over
half of the field is localized in WG1 that is NIM and supports backward
propagation of field. As for the other two branches, they are generally
positive. When approaching $a_{c}$, all the branches coalesce together to
$v_{g}=0$, which implies that at $a_{c}$ the backward energy flux in the NIM
WG balances the forward one in the dielectric WGs and the surrounding air.
This phenomenon is also demonstrated in R36  about EPs in NIM-dielectric WGs.
Once again, Fig. 5 confirms that at $a_{c}$ an EP3 is achieved.

##  III Discussion

Up to now we have proved that it is feasible to achieve an EP3 in a lossfree
WG system containing NIM. This kind of EP3 is produced via the hybridization
of the anti-$\mathcal{PT}$ symmetric coupling in the NIM-dielectric pair and
the Hermitian coupling in the dielectric-dielectric pair. No loss or gain is
required. This work highlights the great potential of NIM in overcoming the
obstacles of ordinary non-Hermitian optics, and the possibilities of combining
anti-$\mathcal{PT}$, $\mathcal{PT}$, and Hermitian couplings for various
purposes. Albeit in this article we only consider NIM-dielectric-dielectric
configuration, EP3 can be achieved in the dielectric-NIM-NIM configuration as
well because the coupling in the NIM-NIM pair is also Hermitian.

Our work provides a useful route in designing lossfree systems in achieving
EPs with further increased orders. For example, here we can propose a
schematic structural design for an EP4 by using two dielectric and two NIM
WGs. Arranging these four WGs parallel in the dielectric-dielectric-NIM-NIM
order, and assuming their resonant wavevectors at $\omega_{0}$ are $k_{1}$,
$k_{2}$, $k_{1}$, and $k_{2}$, respectively, the effective Hamiltonian based
on the coupled-mode theory with nearest-neighbor interaction can be expressed
as

$\left[\begin{array}[]{cccc}k_{1}&A&0&0\\\ A&k_{2}&B&0\\\ 0&-B&k_{1}&C\\\
0&0&C&k_{2}\\\ \end{array}\right]\left[\begin{array}[]{cccc}\psi_{1}\\\
\psi_{2}\\\ \psi_{3}\\\ \psi_{4}\\\
\end{array}\right]=\beta\left[\begin{array}[]{cccc}\psi_{1}\\\ \psi_{2}\\\
\psi_{3}\\\ \psi_{4}\\\ \end{array}\right],$ (16)

where $A$, $B$, and $C$ are all real for simplification. Assuming

$k_{1,2}=k_{0}\pm\delta,$ (17)

it is then easy to prove that under the conditions of

$B=\pm(|A|+|C|),\delta^{2}=|AC|,$ (18)

all the four eigensolutions coalesce together in forming an EP4 at

$\beta_{\rm EP4}=k_{0}.$ (19)

The hierarchical construction of higher-order EPs by using above proposed
method deserves a further discussion.

##  IV Conclusion

In summary, here we show it is feasible to achieve an EP3 in loss-free WGs by
utilizing the anti-$\mathcal{PT}$ symmetry induced by NIM. TMM simulation
agrees well with the prediction of coupled-mode theory. A structure in
supporting EP4 is also designed. This work highlights the great potential of
loss-free NIM, together with the associated anti-$\mathcal{PT}$ symmetry, in
the study of non-Hermitian optics.

This work was supported by the Natural National Science Foundation of China
(NSFC) (12104227, 12274241), and the Scientific Research Foundation of Nanjing
Institute of Technology (YKJ202021).

## References

## References

  * (1) Feng L, El-Ganainy R and Ge L 2017 Nat. Photo. 11 75262
  * (2) El-Ganainy R, Makris K G, Khajavikhan M, Musslimani Z H, Rotter S and Christodoulides D N 2018 Nat. Phys. 14 11
  * (3) H. Cao and J. Wiersig, Rev. Mod. Phys. 87, 61 (2015).
  * (4) S. K. Ozdemir, S. Rotter, F. Nori, and L. Yang, Nat. Mater. 18, 783 (2019).
  * (5) M. A. Miri and A. Alù, Science 363, eaar7709 (2019).
  * (6) J. B. Khurgin, Y. Sebbag, E. Edrei, R. Zektzer, K. Shastri, U. Levy, and F. Monticone, Optica 8, 563 (2021).
  * (7) H. Z. Chen, T. Liu, H. Y. Luan, R. J. Liu, X. Y. Wang, X. F. Zhu, Y. B. Li, Z. M. Gu, S. J. Liang, H. Gao, L. Lu, L. Ge, S. Zhang, J. Zhu, and R. M. Ma, Nat. Phys. 16, 571 (2020).
  * (8) X. Shu, A. Li, G. Hu, J. Wang, A. Alu, and L. Chen, Nat. Comm. 13, 2123 (2022).
  * (9) H. Nasari, G. Lopez-Galmiche, H. E. Lopez-Aviles, A. Schumer, A. U. Hassan, Q. Zhong, S. Rotter, P. LiKamWa, D. N. Christodoulides, and M. Khajavikhan, Nature 605, 256 (2022).
  * (10) J. Doppler, A. A. Mailybaev, J. Bohm, U. Kuhl, A. Girschik, F. Libisch, T. J. Milburn, P. Rabl, N. Moiseyev, and S. Rotter, Nature 537, 76 (2016).
  * (11) Z. Lin, A. Pick, M. Loncar, and A. W. Rodriguez, Phys. Rev. Lett. 117, 107402 (2016).
  * (12) H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Nature 548, 187 (2017).
  * (13) W. Chen, S. K. Ozdemir, G. Zhao, J. Wiersig, and L. Yang, Nature 548, 192 (2017).
  * (14) J. Wiersig, Phys. Rev. A 93, 033809 (2016).
  * (15) T. Goldzak, A. A. Mailybaev, and N. Moiseyev, Phys. Rev. Lett. 120, 013901 (2018).
  * (16) X. Z. Zhang, R. Z. Luo, and J. Chen, Opt. Express 30, 38753 (2022).
  * (17) A. Pick, B. Zhen, O. D. Miller, C. W. Hsu, F. Hernandez, A. W. Rodriguez, M. Soljacic, and S. G. Johnson, Opt. Express 25, 012325 (2017).
  * (18) H. Jing, S. K. Ozdemir, H. Lu, and F. Nori, Sci. Rep. 7, 3386 (2017).
  * (19) J. Schnabel, H. Cartarius, J. Main, G. Wunner, and W. D. Heiss, Phys. Rev. A 95, 053868 (2019).
  * (20) L. Pan, S. Chen, and X. Cui, Phys. Rev. A 99, 011601(R) (2019).
  * (21) N. Habler and J. Scheuer, Phys. Rev. A 101, 043828 (2020).
  * (22) X. Wang, G. Guo, and J. Berakdar, Phys. Rev. Appl. 15, 034050 (2021).
  * (23) Y. L. Wu, P. J. Zhou, T. Li, W. S. Wan, and Y. Zou, Opt. Express 29, 6080 (2021).
  * (24) W. Xiong, Z. Li, Y. Song, J. Chen, G. Q. Zhang, and M. Wang, Phys. Rev. A 104, 063508 (2021).
  * (25) A. Novitsky, F. Morozko, D. Gao, L. Gao, A. Karabchevsky, and D. V. Novitsky, Phys. Rev. B 106, 195410 (2022).
  * (26) Q. Zhong, J. Kou, S. K. Ozdemir, and R. El-Ganainy, Phys. Rev. Lett. 125, 203602 (2020).
  * (27) J. Wiersig, Phys. Rev. A 106, 063526 (2022).
  * (28) H. Benisty, A. Lupu, and A. Degiron, Phys. Rev. A, 91, 053825 (2015)
  * (29) W. D. Heiss and G. Wunner, J. Phys. A: Math. Theor., 49, 495303 (2016).
  * (30) G. Demange and E. M. Graefe, J. Phys. A: Math. Theor. 45 025303 (2012).
  * (31) V. G. Veselago, Sov. Phys. Usp. 10, 509 (1968).
  * (32) A. C. Peacock and N. G. R. Broderick, Opt. Express 11, 2502 (2003).
  * (33) I. V. Shadrivov, A. A. Sukhorukov, and Y. S. Kivshar, Phys. Rev. E 67, 057602 (2003).
  * (34) A. Yariv, IEEE J. Quantum Electron. 9(9), 919–933 (1973).
  * (35) T. Mealy and F. Capolino, Phys. Rev. A 107, 012214 (2023).
  * (36) L. T. Wu, X. Z. Zhang, R. Z. Luo, and J. Chen, accepted by Opt. Express
  * (37) S. Longhi, Phys. Rev. A 82, 031801(R) (2010).