Optimal Antipodal Configuration of $2d$ Points
on a Sphere in $\mathbb{R}^{d}$ for Covering

Sergiy Borodachov

Department of Mathematics, Towson University, 8000 York Rd., Towson, MD, 21252

###### Abstract

We show that among antipodal $2d$-point configurations on the sphere $S^{d-1}$
in $\mathbb{R}^{d}$, the set of vertices of a regular cross-polytope inscribed
in $S^{d-1}$ uniquely solves the best-covering problem (this is new for $d\geq
5$) and the maximal polarization problem for potentials given by a function of
the distance squared with a positive and convex second derivative ($d\geq 3$).

Keywords: Antipodal spherical code, regular cross-polytope, best-covering
problem, polarization problem, convex polytope, radial projection, area
argument.

MSC 2020: 52A20, 52B05, 52B10, 52B12, 52C17.

##  1 Introduction

Let
$S^{d-1}=\\{(x_{1},\ldots,,x_{d})\in\mathbb{R}^{d}:x_{1}^{2}+\ldots+x_{d}^{2}=1\\}$,
$d\geq 2$, be the unit sphere in the Euclidean space $\mathbb{R}^{d}$. For a
given point configuration $\omega_{N}:=\\{{\bf x}_{1},\ldots,{\bf
x}_{N}\\}\subset S^{d-1}$, let

$\rho(\omega_{N},S^{d-1}):=\max\limits_{{\bf x}\in
S^{d-1}}\min\limits_{i=\overline{1,N}}\left|{\bf x}-{\bf x}_{i}\right|$

denote its mesh norm (relative to the sphere). The classical optimal covering
problem on the sphere is to minimize the mesh-norm.

######  Problem 1.1.

For a given $N\in\mathbb{N}$, find the quantity

$\rho_{N}(S^{d-1}):=\inf\limits_{\omega_{N}\subset
S^{d-1}}\rho(\omega_{N},S^{d-1})$ (1)

and optimal-covering $N$-point configurations on $S^{d-1}$; i.e.,
configurations $\omega^{\ast}_{N}\subset S^{d-1}$ attaining the infimum on the
right-hand side of (1).

Solution to Problem 1.1 is known on $S^{1}$ for every $N$ (equally spaced
points), on $S^{2}$ for $N=1-8$, $10$, $12$, and $14$, see the works by Fejes-
Tóth [16, 17, 18], Schütte [25], and Wimmer [29], on $S^{3}$ for $N=1-6$ and
$8$, see the works by Galiev [19], Böröczky and Wintsche [7], and Dalla,
Larman, Mani-Levitska, and Zong [14], and on $S^{d-1}$, $d>4$, for $1\leq
N\leq d+2$, see [19, 7]. For more information on optimal covering, see, among
others, books by Fejes-Tóth [16, 17], Rogers [24], and Böröczky [6].

The main goal of this paper is a further study of the optimal-covering
property of the set of vertices of a regular cross-polytope inscribed in
$S^{d-1}$, which is any $2d$-point configuration of the form
$\overline{\omega}_{2d}:=\\{\pm{\bf a}_{1},\ldots,\pm{\bf a}_{d}\\}$, where
$\\{{\bf a}_{1},\ldots,{\bf a}_{d}\\}$ is an orthonormal basis in
$\mathbb{R}^{d}$. The optimal covering property of $\overline{\omega}_{2d}$ is
known on $S^{d-1}$ for $d=2$ (a basic result), $d=3$ [16, 17], and $d=4$ [14].
It is a well-known open question for $d\geq 5$. We resolve it here for
antipodal configurations.

The proof in [14] for $d=4$ uses the area argument together with the Euler-
Poincaré formula and the fact that (cf. [5]) the largest ($d-1$)-dimensional
area of a spherical simplex inscribed in a given spherical cap of angular
radius less than $\pi/2$ is that of a regular simplex. If the number $\tau$ of
($d-1$)-dimensional simplices obtained after triangulating the facets of the
convex hull of a given $2d$-point configuration is less than or equal to
$2^{d}$, the area argument combined with the above mentioned fact from [5]
complete the proof. If $\tau>2^{d}$ then the Euler-Poincaré formula is used to
provide an efficient estimate for the number of one-dimensional edges of the
convex hull. Since the number of terms in the Euler-Poincaré formula grows
with $d$, such an estimate is not possible when $d\geq 5$ and the argument
fails. However, if one makes an additional assumption that the configurations
are antipodal, then it can be shown that $\tau=2^{d}$, and the area argument
can be used together with the above mentioned fact from [5] to complete the
proof, see Theorem 2.2.

The paper is structured as follows. We state our main result, Theorem 2.2, in
Section 2 and its consequences in Section 3. Section 4 contains the proof of
Theorem 2.2. In Section 5 we state and prove Theorem 5.1 on the minimum value
of the potential of a regular cross-polytope, which we use in Section 6 to
prove consequences of Theorem 2.2. The known proof of Lemma 4.1 is given in
Appendix (Section 7) for completeness.

##  2 Main result

Recall that a point configuration on $S^{d-1}$ is called antipodal if together
with a point ${\bf x}$ it contains $-{\bf x}$. We say that a given point set
in $\mathbb{R}^{d}$ is in general position if it is not contained in any
hyperplane (or, equivalently, in any $(d-1)$-dimensional affine subspace). It
will be convenient for us to recast the optimal-covering problem in terms of
dot products.

######  Remark 2.1.

A point configuration $\omega_{N}\subset S^{d-1}$ is a solution to Problem 1.1
if and only if $\omega_{N}$ maximizes the quantity

$\eta(\omega_{N},S^{d-1}):=\min\limits_{{\bf x}\in
S^{d-1}}\max\limits_{i=\overline{1,N}}{\bf x}\cdot{\bf x}_{i}$

over all $N$-point configurations $\omega_{N}=\\{{\bf x}_{1},\ldots,{\bf
x}_{N}\\}$ on $S^{d-1}$.

We allow configurations, where points may coincide. However, if an antipodal
configuration of $2d$ points on $S^{d-1}$ is in general position, its points
are pairwise distinct.

The main result of this paper is the following.

######  Theorem 2.2.

Let $d\geq 2$ and $\omega_{2d}$ be an antipodal configuration of $2d$ points
on $S^{d-1}$. Then

$\eta(\omega_{2d},S^{d-1})\leq\frac{1}{\sqrt{d}}.$ (2)

Equality in (2) holds if and only if $\omega_{2d}$ is the set of vertices of a
regular cross-polytope inscribed in $S^{d-1}$.

For $d=3$ and $4$, Theorem 2.2 is a special case of theorems proved in [16,
17] and in [14], respectively. For $d=2$, it is a basic result.

##  3 A consequence of Theorem 2.2 for polarization

The best-covering problem is a limiting case of the polarization problem
stated next. Let $f:[0,4]\to(-\infty,\infty]$ be a function finite and
continuous on $(0.4]$ such that $f(0)=\lim\limits_{t\to 0^{+}}f(t)$. We will
call $f$ a potential function. For a given point configuration
$\omega_{N}=\\{{\bf x}_{1},\ldots,{\bf x}_{N}\\}\subset S^{d-1}$, denote

$P_{f}(\omega_{N},S^{d-1}):=\min\limits_{{\bf x}\in
S^{d-1}}\sum\limits_{i=1}^{N}f\left(\left|{\bf x}-{\bf
x}_{i}\right|^{2}\right)$

and let

$\mathcal{P}_{f}(S^{d-1},N):=\sup\limits_{\omega_{N}\subset
S^{d-1}}P_{f}(\omega_{N},S^{d-1}).$ (3)

The (max-min) $N$-point polarization problem on the sphere is stated in the
following way.

######  Problem 3.1.

Find quantity (3) and $N$-point configurations $\omega_{N}^{\ast}$ on
$S^{d-1}$ that attain the supremum on the right-hand side of (3).

In the case of a sphere, solution to Problem 3.1 is known on the unit circle
$S^{1}$ for every $N\geq 1$ and $f(t)=t^{-s/2}$, $s>0$, and $f(t)=-t^{-s/2}$,
$-1\leq s<0$, which correspond to the Riesz potential, as well as for
$f(t)=\frac{1}{2}\ln\frac{1}{t}$, which corresponds to the logarithmic
potential, see the works by Stolarsky, Ambrus, Nikolov, Rafailov, Ball,
Erdélyi, Saff, Hardin, and Kendall [26, 1, 22, 2, 15, 20]. The solution is the
set of vertices of a regular $N$-gon inscribed in $S^{1}$. In fact, paper [20]
established this for arbitrary interactions given by a decreasing and convex
function of the geodesic distance on $S^{1}$. Optimal $N$-point configurations
for polarization on $S^{1}$ were also characterized by Bosuwan and Ruengrot
[12] for the Riesz potential with $s=-2,-4,\ldots,2-2N$, $N\geq 2$.

On the sphere $S^{d-1}$, $d\geq 3$, the solution to Problem 3.1 is known for
$1\leq N\leq d$ (basic result) and $N=d+1$, see the works by Su [28] and the
author [8] showing the optimality of a regular simplex. Also, for $N=2d$, the
optimality of a regular cross-polytope $\overline{\omega}_{2d}$ was shown
among centered configurations by Boyvalenkov, Dragnev, Hardin, Saff, and
Stoyanova [13]. A configuration $\omega_{2d}\subset S^{d-1}$ is called
centered if there is a point ${\bf y}\in S^{d-1}$ such that
$-\frac{1}{\sqrt{d}}\leq{\bf y}\cdot{\bf x}_{i}\leq\frac{1}{\sqrt{d}}$,
$i=1,\ldots,2d$.

Other settings of polarization problem were studied in [26, 21, 15, 12, 3, 4,
9, 13]. More extensive reviews on polarization (including the continuous
version and asymptotics) can be found, for example, in book [11].

Polarization problem can also be recast in terms of dot products. Let
$g(t):=f(2-2t)$, $t\in[-1,1]$. Then $f$ is a potential function if and only if
$g:[-1,1]\to(-\infty,\infty]$ is a function finite and continuous on $[-1,1)$
such that $g(1)=\lim\limits_{t\to 1^{-}}g(t)$. For every configuration
$\omega_{N}=\\{{\bf x}_{1},\ldots,{\bf x}_{N}\\}\subset S^{d-1}$, we have

$P^{g}(\omega_{N},S^{d-1}):=\min\limits_{{\bf x}\in
S^{d-1}}\sum\limits_{i=1}^{N}g({\bf x}\cdot{\bf x}_{i})=\min\limits_{{\bf
x}\in S^{d-1}}\sum\limits_{i=1}^{N}f\left(\left|{\bf x}-{\bf
x}_{i}\right|^{2}\right)=P_{f}(\omega_{N},S^{d-1}).$

Thus, a point configuration $\omega_{N}$ maximizes the quantity
$P_{f}(\omega_{N},S^{d-1})$ if and only if it maximizes the quantity
$P^{g}(\omega_{N},S^{d-1})$.

Inequality (2) in Theorem 2.2 is equivalent to the following statement.

######  Corollary 3.2.

Any antipodal $2d$-point configuration on $S^{d-1}$, $d\geq 2$, is centered.

Paper [13] proved universal bounds for polarization and used them to show that
regular cross-polytope $\overline{\omega}_{2d}$ is optimal for polarization
among all centered configurations on $S^{d-1}$. This result together with
Corollary 3.2 imply the following (the equality in (4) below follows from
Theorem 5.1).

######  Corollary 3.3.

Suppose $g:[-1,1]\to(-\infty,\infty]$ is a function continuous on $[-1,1)$
with $g(1)=\lim\limits_{t\to 1^{-}}g(t)$ and differentiable on $(-1,1)$ such
that $g^{\prime\prime}$ is non-negative and convex on $(-1,1)$. Suppose also
that $\omega_{2d}\subset S^{d-1}$, $d\geq 2$, is any antipodal $2d$-point
configuration. Then

$P^{g}(\omega_{2d},S^{d-1})\leq
P^{g}(\overline{\omega}_{2d},S^{d-1})=d\left(g\left(\frac{1}{\sqrt{d}}\right)+g\left(-\frac{1}{\sqrt{d}}\right)\right).$
(4)

Here we provide an elementary proof of Corollary 3.3 that avoids the use of
universal bounds for polarization and is valid for any centered configuration.
Concerning uniqueness of the optimal configuration, we have the following.

######  Corollary 3.4.

Under assumptions of Corollary 3.3, if $g^{\prime\prime}$ is strictly positive
on $(-1,1)$, then equality holds throughout (4) if and only if $\omega_{2d}$
is the set of vertices of a regular cross-polytope inscribed in $S^{d-1}$.

We remark that a similar uniqueness question for centered configurations
remains open. Observe that we do not require the potential function $g$ to be
monotone in Corollaries 3.3 and 3.4.

##  4 Proof of Theorem 2.2

We need the following auxiliary statement proved by Böröczky in [5] (see also
[6, Lemma 6.7.2]). Let $r:\mathbb{R}^{d}\setminus\\{{\bf 0}\\}\to S^{d-1}$,
$r({\bf x})={\bf x}/\left|{\bf x}\right|$, denote the radial projection onto
$S^{d-1}$. We will call a $(d-1)$-simplex a simplex in $\mathbb{R}^{d}$ with
$d$ vertices having dimension $(d-1)$.

######  Lemma 4.1.

Let $H$ be a hyperplane in $\mathbb{R}^{d}$, $d\geq 2$, with equation
$x_{d}=a$, where $0<a<1$. Let $Y$ be a non-degenerate $(d-1)$-simplex
inscribed in $H\cap S^{d-1}$. Then the $(d-1)$-dimensional volume of the
radial projection $r(Y)$ of $Y$ onto $S^{d-1}$ is the largest if and only if
$Y$ is a regular $(d-1)$-simplex inscribed in $H\cap S^{d-1}$.

For completeness, we provide the proof of Lemma 4.1 in the Appendix.

Let $D(\omega_{2d})$ denote the convex hull of a point configuration
$\omega_{2d}\subset S^{d-1}$. Points from $\omega_{2d}$ will also be called
vertices. Recall that a hyperplane ${\bf x}\cdot{\bf a}=\alpha$ in
$\mathbb{R}^{d}$ is called a supporting hyperplane for a convex body $B$ if
for every ${\bf z}\in B$, we have ${\bf z}\cdot{\bf a}\leq\alpha$, while for
some ${\bf y}\in B$, we have ${\bf y}\cdot{\bf a}=\alpha$.

######  Lemma 4.2.

Let $\omega_{2d}\subset S^{d-1}$, $d\geq 2$, be an antipodal configuration of
$2d$ points in general position. Then the interior of $D(\omega_{2d})$
contains the origin, and the boundary of $D(\omega_{2d})$ is the union of
$2^{d}$ $(d-1)$-simplices whose vertices are in $\omega_{2d}$, pairwise
intersections have $(d-1)$-dimensional measure $0$, and the hyperplane
containing each simplex does not pass through the origin.

###### Proof.

Pick an arbitrary subset of $d$ linearly independent vectors from
$\omega_{2d}$ and denote it by $\\{{\bf y}_{1},\ldots,{\bf y}_{d}\\}$. Such a
subset exists, since $\omega_{2d}$ is in general position. Let
$\Pi:=\\{-1,1\\}^{d}$. Then any of the $2^{d}$ sets
$T_{\boldsymbol{\sigma}}:=\\{\sigma_{1}{\bf y}_{1},\ldots,\sigma_{d}{\bf
y}_{d}\\}$, where $\boldsymbol{\sigma}=(\sigma_{1},\ldots,\sigma_{d})\in\Pi$,
is a linearly independent subset of $\omega_{2d}$ and, hence, is a set of
vertices of a $(d-1)$-simplex, which we will denote by
$F_{\boldsymbol{\sigma}}$. For every $\boldsymbol{\sigma}\in\Pi$, the set
$T_{\boldsymbol{\sigma}}$ is contained in a unique hyperplane
$H_{\boldsymbol{\sigma}}$. Then the set
$T_{-\boldsymbol{\sigma}}=-T_{\boldsymbol{\sigma}}$ is contained in the
hyperplane $-H_{\boldsymbol{\sigma}}$ with $-H_{\boldsymbol{\sigma}}\neq
H_{\boldsymbol{\sigma}}$, since $\omega_{2d}=T_{\boldsymbol{\sigma}}\cup
T_{-\boldsymbol{\sigma}}$ is in general position. Then $D(\omega_{2d})$ is
contained in the closed subset of $\mathbb{R}^{d}$ bounded by
$H_{\boldsymbol{\sigma}}$ and $-H_{\boldsymbol{\sigma}}$. Then the simplex
$F_{\boldsymbol{\sigma}}$ is contained in the boundary $\partial
D(\omega_{2d})$ of $D(\omega_{2d})$ for every $\boldsymbol{\sigma}\in\Pi$.

Since $\omega_{2d}$ is in general position, the origin is in the interior of
$D(\omega_{2d})$ as a point on the line segment joining two points in relative
interiors of opposite facets.

Let now ${\bf x}$ be any point in $\partial D(\omega_{2d})$. Then ${\bf
x}\neq{\bf 0}$ and there is a unique set of numbers
$\alpha_{1},\ldots,\alpha_{d}$ such that ${\bf x}=\alpha_{1}{\bf
y}_{1}+\ldots+\alpha_{d}{\bf y}_{d}$. For some
$\boldsymbol{\sigma}=(\sigma_{1},\ldots,\sigma_{d})\in\Pi$, we have ${\bf
x}=\beta_{1}\sigma_{1}{\bf y}_{1}+\ldots+\beta_{d}\sigma_{d}{\bf y}_{d}$,
where $\beta_{i}:=\left|\alpha_{i}\right|\geq 0$, $i=1,\ldots,d$. Let
$\beta:=\sum_{i=1}^{d}\beta_{i}$. Then $\beta>0$ and $(1/\beta){\bf x}\in
F_{\boldsymbol{\sigma}}\subset\partial D(\omega_{2d})$. If it were that
$\beta>1$, then, since $(1/\beta){\bf x}\in H_{\boldsymbol{\sigma}}$, the
origin and ${\bf x}$ would lie in different half-spaces relative to
$H_{\boldsymbol{\sigma}}$, which is a supporting hyperplane for
$D(\omega_{2d})$. Since the origin is in $D(\omega_{2d})$, we have ${\bf
x}\notin D(\omega_{2d})$; that is, ${\bf x}\notin\partial D(\omega_{2d})$. If
it were that $\beta<1$, then ${\bf x}$ would be in the relative interior of
the line segment joining $(1/\beta){\bf x}$ and ${\bf 0}$. Since
$(1/\beta){\bf x}$ is in $D(\omega_{2d})$ and ${\bf 0}$ is in its interior,
then ${\bf x}$ is also in the interior of $D(\omega_{2d})$; that is, ${\bf
x}\notin\partial D(\omega_{2d})$. This contradiction shows that $\beta=1$.
Then ${\bf x}\in F_{\boldsymbol{\sigma}}$.

Thus, $\partial
D(\omega_{2d})=\bigcup\limits_{\boldsymbol{\sigma}\in\Pi}F_{\boldsymbol{\sigma}}$.
Assume that the intersection of two simplices $F_{{\boldsymbol{\sigma}}}$ and
$F_{{\boldsymbol{\sigma}}^{\prime}}$, with
${\boldsymbol{\sigma}}=(\sigma_{1},\ldots,\sigma_{d})\neq{\boldsymbol{\sigma}}^{\prime}=(\sigma_{1}^{\prime},\ldots,\sigma_{d}^{\prime})$
is non-empty. For every point ${\bf z}$ in both simplices, there are numbers
$\beta_{1},\ldots,\beta_{d},\beta_{1}^{\prime},\ldots,\beta_{d}^{\prime}\geq
0$ such that $\sum_{i=1}^{d}\beta_{i}=\sum_{i=1}^{d}\beta_{i}^{\prime}=1$ and

${\bf z}=\beta_{1}\sigma_{1}{\bf y}_{1}+\ldots+\beta_{d}\sigma_{d}{\bf
y}_{d}=\beta_{1}^{\prime}\sigma_{1}^{\prime}{\bf
y}_{1}+\ldots+\beta_{d}^{\prime}\sigma_{d}^{\prime}{\bf y}_{d}.$

Since vectors ${\bf y}_{1},\ldots,{\bf y}_{d}$ are linearly independent, we
have $\beta_{i}\sigma_{i}=\beta_{i}^{\prime}\sigma_{i}^{\prime}$,
$i=1,\ldots,d$. Since $\sigma_{j}\neq\sigma_{j}^{\prime}$ for some $j$, we
have $\beta_{j}=-\beta_{j}^{\prime}$. In view of the non-negativity, we have
$\beta_{j}=\beta_{j}^{\prime}=0$. Therefore, ${\bf z}$ belongs to a lower-
dimensional face of $F_{{\boldsymbol{\sigma}}}$. Thus, ${\rm
Vol}_{d-1}(F_{{\boldsymbol{\sigma}}}\cap
F_{{\boldsymbol{\sigma}^{\prime}}})=0$. Furthermore, the hyperplane
$H_{\boldsymbol{\sigma}}$ containing $F_{\boldsymbol{\sigma}}$ does not
contain the origin, since if it did, vectors ${\bf y}_{1},\ldots,{\bf y}_{d}$
would be linearly dependent. ∎

Recall that an intersection $U$ of a convex polytope $P$ with its supporting
hyperplane, such that ${\rm dim}\ \\!U=d-1$, is called a facet of $P$. The
polytope $D(\omega_{2d})$ in Lemma 4.2 has $2^{d}$ facets.

###### Proof of Theorem 2.2.

We first verify that
$\eta(\overline{\omega}_{2d},S^{d-1})=\frac{1}{\sqrt{d}}$. Without loss of
generality, we can assume that
$\overline{\omega}_{2d}=\widetilde{\omega}_{2d}:=\\{\pm{\bf
e}_{1},\ldots,\pm{\bf e}_{d}\\}$, where ${\bf e}_{1},\ldots,{\bf e}_{d}$ are
standard basis vectors in $\mathbb{R}^{d}$. Then for every vector ${\bf
x}=(x_{1},\ldots,x_{d})\in S^{d-1}$, we have
$\max\limits_{i=\overline{1,d}}{\bf x}\cdot(\pm{\bf
e}_{i})=\max\limits_{i=\overline{1,d}}\left|x_{i}\right|\geq\frac{1}{\sqrt{d}}$,
since $\sum_{i=1}^{d}\left|x_{i}\right|^{2}=1$, with equality occuring
whenever each coordinate of ${\bf x}$ is $\pm\frac{1}{\sqrt{d}}$.

Assume to the contrary that there is an antipodal configuration
$\omega_{2d}\subset S^{d-1}$ such that $\eta(\omega_{2d},S^{d-1})>1/\sqrt{d}$.
Then $D(\omega_{2d})$ contains ${\bf 0}$ in its interior (if it did not, we
would have $\eta(\omega_{2d},S^{d-1})\leq 0<1/\sqrt{d}$). Furthermore,
$\omega_{2d}$ is in general position and, hence, its points are pairwise
distinct. Polytope $D(\omega_{2d})$ contains the sphere $S$ of radius
$1/\sqrt{d}$ centered at the origin. If it didn’t, then any point ${\bf a}\in
S\setminus D(\omega_{2d})$ would be strictly separated from $D(\omega_{2d})$
by some hyperplane $L=\\{{\bf x}:{\bf x}\cdot{\bf v}=c\\}$, where we can take
$\left|{\bf v}\right|=1$ and ${\bf v}\cdot{\bf a}>c$. Then for every ${\bf
x}_{i}\in\omega_{2d}$, we would have ${\bf x}_{i}\cdot{\bf v}<c<{\bf
a}\cdot{\bf v}\leq\left|{\bf a}\right|\left|{\bf v}\right|=1/\sqrt{d}$ which
would contradict the contrary assumption.

By Lemma 4.2, boundary of $D(\omega_{2d})$ is the union of $(d-1)$-simplices
(which we denote by $F_{1},\ldots,F_{2^{d}}$) with vertices in $\omega_{2d}$.
Their pairwise intersections have $(d-1)$-dimensional volume $0$, and there
are $a_{i}>0$ and ${\bf z}_{i}\in S^{d-1}$, $i=1,\ldots,2^{d}$, such that the
hyperplane $H_{i}:=\\{{\bf x}:{\bf x}\cdot{\bf z}_{i}=a_{i}\\}$ contains
$F_{i}$. Denote by $k$ an index such that the radial projection $r(F_{k})$
onto $S^{d-1}$ of the simplex $F_{k}$ has $(d-1)$-dimensional volume at least
$1/2^{d}$ of the $(d-1)$-dimensional volume of $S^{d-1}$. If $V$ denotes a
regular $(d-1)$-simplex inscribed in $H_{k}\cap S^{d-1}$, then in view of
Lemma 4.1, we have

$2^{-d}{\rm Vol}_{d-1}(S^{d-1})\leq{\rm Vol}_{d-1}(r(F_{k}))\leq{\rm
Vol}_{d-1}(r(V)).$

Since $H_{k}$ does not contain the origin (by Lemma 4.2), there are no
antipodal pairs among the vertices of $F_{k}$. Then the remaining $d$ points
from $\omega_{2d}$ are contained in the hyperplane ${\bf x}\cdot{\bf
z}_{k}=-a_{k}$. We have $a_{k}=\max\limits_{i=\overline{1,2d}}{\bf
z}_{k}\cdot{\bf x}_{i}\geq\eta(\omega_{2d},S^{d-1})>1/\sqrt{d}$. The
$(d-1)$-dimensional volume (denoted by $\nu$) of the radial projection $r(W)$
of a regular simplex $W$ incribed in $M\cap S^{d-1}$, where $M=\\{{\bf x}:{\bf
x}\cdot{\bf z}_{k}=1/\sqrt{d}\\}$, will be strictly larger than ${\rm
Vol}_{d-1}(r(V))$. This is because $a_{k}>1/\sqrt{d}$ and the radius of the
intersection of a hyperplane perpendicular to ${\bf z}_{k}$ with $S^{d-1}$
decreases as the hyperplane moves further away from the origin. Since $r(W)$
is the same as the radial projection of any facet of
$D(\overline{\omega}_{2d})$ and there are exactly $2^{d}$ facets, we have

$2^{-d}{\rm Vol}_{d-1}(S^{d-1})\leq{\rm Vol}_{d-1}(r(V))<{\rm
Vol}_{d-1}(r(W))=2^{-d}{\rm Vol}_{d-1}(S^{d-1}).$

This contradiction proves (2).

To complete the proof of Theorem 2.2, assume that equality holds in
$\eqref{3}$. Then the origin is in the interior of $D(\omega_{2d})$. Assume to
the contrary that there is a facet $Q$ of $D(\omega_{2d})$ such that the
hyperplane containing it is at a distance strictly greater than $1/\sqrt{d}$
from the origin. Then by Lemma 4.1, the radial projection of $Q$ onto
$S^{d-1}$ has $(d-1)$-dimensional volume strictly less than $\nu$. Equality in
(2) implies that every other facet $Y$ of $D(\omega_{2d})$ is contained in a
hyperplane whose distance to the origin is at least $1/\sqrt{d}$. By Lemma
4.1, ${\rm Vol}_{d-1}(r(Y))\leq\nu$ and, by Lemma 4.2, $D(\omega_{2d})$ has
$2^{d}$ facets. Then the radial projection onto $S^{d-1}$ of the whole
boundary of $D(\omega_{2d})$ has $(d-1)$-dimensional volume strictly less than
$2^{d}\nu={\rm Vol}_{d-1}(S^{d-1})$. This contradiction shows that the
hyperplane containing each facet of $D(\omega_{2d})$ is at a distance exactly
$1/\sqrt{d}$ from the origin. Each facet of $D(\omega_{2d})$ must be a regular
simplex. If some facet $U$ were not, by Lemma 4.1, we would have ${\rm
Vol}_{d-1}(r(U))<\nu$ while for any other facet $J$, we would have ${\rm
Vol}_{d-1}(r(J))\leq\nu$ leading to a similar contradiction. Let ${\bf
a}_{1},\ldots,{\bf a}_{d}$ be the vertices of one of the facets of
$D(\omega_{2d})$. Then $\omega_{2d}=\\{\pm{\bf a}_{1},\ldots,\pm{\bf
a}_{d}\\}$. It is not difficult to verify that ${\bf a}_{i}\cdot{\bf
a}_{j}=0$, $i\neq j$. Then $\\{{\bf a}_{1},\ldots,{\bf a}_{d}\\}$ is an
orthonormal basis in $\mathbb{R}^{d}$ and $\omega_{2d}$ is the set of vertices
of a regular cross-polytope inscribed in $S^{d-1}$. ∎

##  5 Minimum of the potential of a regular cross-polytope

In order to show the optimality of $\overline{\omega}_{2d}$ for Problem 3.1,
one will need to know the quantity $P_{f}(\overline{\omega}_{2d},S^{d-1})$ by
locating the absolute minima of the potential of $\overline{\omega}_{2d}$ on
the sphere $S^{d-1}$. This was done earlier for the vertices of a regular
$N$-gon inscribed in $S^{1}$ and for the vertices of a regular simplex, cross-
polytope, and cube inscribed in $S^{d-1}$ for Riesz potential functions
($s\neq 0$) and their horizontal translations, see the works by Stolarsky,
Nikolov, and Rafailov [26, 27, 22, 23]. For general potentials, this has been
recently done in [8] for vertices of a regular simplex. Below, we extend one
of the results of [27, 23] by finding the absolute minima of the potential of
$\overline{\omega}_{2d}$ for potential functions $g$ that have a convex second
derivative. Our proof is different from the one in [27, 23] and uses
polynomial interpolation and convexity of $g^{\prime\prime}$. Without loss of
generality, we can assume in this section that
$\overline{\omega}_{2d}=\widetilde{\omega}_{2d}=\\{\pm{\bf
e}_{1},\ldots,\pm{\bf e}_{d}\\}$, where ${\bf e}_{1},\ldots,{\bf e}_{d}$ is
the standard basis in $\mathbb{R}^{d}$.

######  Theorem 5.1.

Let $d\geq 2$ and $g:[-1,1]\to(-\infty,\infty]$ be a function continuous on
$[-1,1)$ and differentiable on $(-1,1)$ such that $g(1)=\lim\limits_{t\to
1^{-}}g(t)$ and $g^{\prime\prime}$ is convex on $(-1,1)$. Then the potential

$p^{g}(\widetilde{\omega}_{2d},{\bf x}):=\sum\limits_{{\bf
y}\in\widetilde{\omega}_{2d}}g({\bf x}\cdot{\bf y})$

achieves its absolute minimum over $S^{d-1}$ at any point of $S^{d-1}$ whose
every coordinate is $1/\sqrt{d}$ or $-1/\sqrt{d}$ (these points are vertices
of the cube dual to $\widetilde{\omega}_{2d}$). Furthermore,

$P^{g}(\widetilde{\omega}_{2d},S^{d-1})=d\left(g\left(\frac{1}{\sqrt{d}}\right)+g\left(-\frac{1}{\sqrt{d}}\right)\right).$

We also remark that in the upcoming paper [10], we obtain the locations of
absolute minima of the potential of any configuration on $S^{d-1}$ which is a
tight spherical design of an even strength or a ($2m-1$)-design contained in
the union of $m$ parallel hyperplanes. Such is, for example,
$\overline{\omega}_{2d}$ for $m=2$.

###### Proof of Theorem 5.1.

Let $h(t):=g(t)+g(-t)$. Then $h$ is even and $h^{\prime\prime}$ is convex on
$(-1,1)$. Let $p$ be the Hermite interpolating polynomial for $h$ at points
$t_{1}=-1/\sqrt{d}$ and $t_{2}=1/\sqrt{d}$. Then the even polynomial
$(p(t)+p(-t))/2$ is also Hermite for $h$ at $t_{1}$ and $t_{2}$. By
uniqueness, $p$ must be even. Since $p$ has degree at most $3$, it has the
form $p(t)=at^{2}+b$. Furthermore,
$\frac{a}{d}+b=p(1/\sqrt{d})=h\left(1/\sqrt{d}\right)$.

We also have, $h(t)\geq p(t)$, $t\in[-1,1]$. Indeed, assume to the contrary
that $v(t):=h(t)-p(t)$ is negative for some $t=t_{0}\in(-1,1)$. Note that
$v(t_{1})=v(t_{2})=v^{\prime}(t_{1})=v^{\prime}(t_{2})=0$. Then the Mean value
theorem and the Rolle’s theorem imply that there are points
$-1<\tau_{1}<\tau_{2}<\tau_{3}<1$ such that $v^{\prime\prime}(\tau_{1})<0$ and
$v^{\prime\prime}(\tau_{2})=v^{\prime\prime}(\tau_{3})=0$ (if $t_{0}<t_{1}$)
or $v^{\prime\prime}(\tau_{1})<0$, $v^{\prime\prime}(\tau_{2})>0$, and
$v^{\prime\prime}(\tau_{3})<0$ (if $t_{1}<t_{0}<t_{2}$) or
$v^{\prime\prime}(\tau_{1})=v^{\prime\prime}(\tau_{2})=0$ and
$v^{\prime\prime}(\tau_{3})<0$ (if $t_{0}>t_{2}$). None of these cases is
possible, since $h^{\prime\prime}$ (and, hence, $v^{\prime\prime}$) is convex
on $(-1,1)$. Thus, $v$ is non-negative on $(-1,1)$; that is, $h(t)\geq p(t)$,
$t\in(-1,1)$. We extend this inequality to the endpoints by passing to the
limit.

Let ${\bf x}=(x_{1},\ldots,x_{d})$ be an arbitrary point on $S^{d-1}$. Then

$\begin{split}p^{g}(\widetilde{\omega}_{2d},{\bf
x})&=\sum\limits_{i=1}^{d}\left(g(x_{i})+g(-x_{i})\right)=\sum\limits_{i=1}^{d}h(x_{i})\geq\sum\limits_{i=1}^{d}p(x_{i})=\sum\limits_{i=1}^{d}(ax_{i}^{2}+b)\\\
&=a+d\cdot b=d\cdot
h\left(\frac{1}{\sqrt{d}}\right)=d\left(g\left(\frac{1}{\sqrt{d}}\right)+g\left(-\frac{1}{\sqrt{d}}\right)\right)=p^{g}(\widetilde{\omega}_{2d},{\bf
x}^{\ast}),\end{split}$

where ${\bf x}^{\ast}$ is any point on $S^{d-1}$ whose every coordinate is
$1/\sqrt{d}$ or $-1/\sqrt{d}$. ∎

##  6 Proof of Corollaries 3.3 and 3.4

###### Proof of Corollary 3.3.

We choose an arbitrary centered configuration $\omega_{2d}=\\{{\bf
y}_{1},\ldots,{\bf y}_{2d}\\}\subset S^{d-1}$ and let ${\bf x}^{\ast}\in
S^{d-1}$ be a point such that

$-\frac{1}{\sqrt{d}}\leq{\bf x}^{\ast}\cdot{\bf y}_{i}\leq\frac{1}{\sqrt{d}},\
\ \ i=1,\ldots,2d.$ (5)

In fact, $-\frac{1}{\sqrt{d}}\leq\pm{\bf x}^{\ast}\cdot{\bf
y}_{i}\leq\frac{1}{\sqrt{d}}$ for all $i$. Let $-t_{1}\leq\ldots\leq-
t_{2d}\leq t_{2d}\leq\ldots\leq t_{1}$ denote all the dot products that ${\bf
x}^{\ast}$ and $-{\bf x}^{\ast}$ form with points ${\bf y}_{i}$. Since
$g^{\prime\prime}$ is non-negative on $(-1,1)$, $g^{\prime}$ is non-decreasing
on $(-1,1)$, and, hence, the function $h(t)=g(t)+g(-t)$ is non-decreasing on
$[0,1)$. Since $t_{1},\ldots,t_{2d}\in[0,1/\sqrt{d}]$, we have

$\begin{split}P^{g}&(\omega_{2d},S^{d-1})\leq\frac{1}{2}\sum\limits_{i=1}^{2d}g({\bf
x}^{\ast}\cdot{\bf y}_{i})+\frac{1}{2}\sum\limits_{i=1}^{2d}g(-{\bf
x}^{\ast}\cdot{\bf
y}_{i})=\frac{1}{2}\sum\limits_{i=1}^{2d}\left(g(t_{i})+g(-t_{i})\right)\\\
&=\frac{1}{2}\sum\limits_{i=1}^{2d}h(t_{i})\leq d\cdot
h\left(\frac{1}{\sqrt{d}}\right)=d\left(g\left(\frac{1}{\sqrt{d}}\right)+g\left(-\frac{1}{\sqrt{d}}\right)\right)=P^{g}(\overline{\omega}_{2d},S^{d-1}),\end{split}$
(6)

where the last equality in (6) holds in view of Theorem 5.1. This proves (4)
for any centered configuration $\omega_{2d}$. In view of Corollary 3.2, we
have Corollary 3.3. ∎

###### Proof of Corollary 3.4.

Assume that equality holds throughout (4) for a given antipodal configuration
$\omega_{2d}\subset S^{d-1}$. Since $\omega_{2d}$ is antipodal, ${\bf
x}^{\ast}$ can be chosen in the beginning of the proof of Corollary 3.3 with
the additional property that $\max\limits_{i=\overline{1,2d}}{\bf
x}^{\ast}\cdot{\bf y}_{i}=\eta(\omega_{2d},S^{d-1})$ (inequalities (5) will
then hold in view of Theorem 2.2). Since equality now holds throughout (6), we
have $\sum_{i=1}^{2d}h(t_{i})=2d\cdot h(1/\sqrt{d})$. Since $0\leq t_{i}\leq
1/\sqrt{d}$, $i=1,\ldots,2d$ and $h$ is strictly increasing on $[0,1)$, we
have $t_{i}=1/\sqrt{d}$, $i=1,\ldots,2d$. Then ${\bf x}^{\ast}$ forms only dot
products $1/\sqrt{d}$ and $-1/\sqrt{d}$ with points of $\omega_{2d}$. Thus,
$\eta(\omega_{2d},S^{d-1})=\max\limits_{i=\overline{1,2d}}{\bf
x}^{\ast}\cdot{\bf y}_{i}=1/\sqrt{d}$. By the uniqueness part of Theorem 2.2,
configuration $\omega_{2d}$ is the set of vertices of a regular cross-polytope
inscribed in $S^{d-1}$. ∎

##  7 Appendix. Proof of Lemma 4.1.

The assertion of Lemma 4.1 is trivial for $d=2$. Therefore, we assume that
$d\geq 3$. Let $Y$ be a $(d-1)$-simplex inscribed in $H\cap S^{d-1}$ whose
radial projection $r(Y)$ onto $S^{d-1}$ has the largest $(d-1)$-dimensional
volume. Let ${\bf y}_{1},\ldots,{\bf y}_{d}$ be the vertices of $Y$. Let
$\Omega(Y)$ be the intersection of the closed unit ball $B^{d}$ centered at
the origin with the cone, denoted by ${\rm cone}\\{{\bf y}_{1},\ldots,{\bf
y}_{d}\\}$, constructed as the convex hull of the rays starting at ${\bf 0}$
and passing through each vertex of $Y$. Then ${\rm Vol}_{d-1}(r(Y))=d\cdot{\rm
Vol}_{d}(\Omega(Y))$.

Assume to the contrary that $Y$ is not a regular simplex. Then there is a
vertex of $Y$ such that some two edges stemming out of it have non-equal
lengths. Without loss of generality, we can assume that $\left|{\bf
y}_{3}-{\bf y}_{1}\right|\neq\left|{\bf y}_{3}-{\bf y}_{2}\right|$. Let $L$ be
the hyperplane that is the perpendicular bisector for the line segment with
endpoints ${\bf y}_{1}$ and ${\bf y}_{2}$. Observe that ${\bf 0}\in L$. Denote
by ${\bf y}_{i}^{\prime}$ the orthogonal projection of the point ${\bf y}_{i}$
onto $L$, $i=3,\ldots,d$, and let $Y^{\prime}$ be the convex hull of $\\{{\bf
y}_{1},{\bf y}_{2},{\bf y}_{3}^{\prime},\ldots,{\bf y}_{d}^{\prime}\\}$. Since
${\bf y}_{1},{\bf y}_{2},{\bf y}_{3}^{\prime},\ldots,{\bf y}_{d}^{\prime}\in
H\cap B^{d}$, we have $Y^{\prime}\subset H\cap B^{d}$ and
$r(Y^{\prime})\subset C:=\\{(x_{1},\ldots,x_{d})\in S^{d-1}:x_{d}\geq a\\}$.
Also, $r(Y^{\prime})={\rm cone}\\{{\bf y}_{1},{\bf y}_{2},{\bf
y}_{3}^{\prime},\ldots,{\bf y}_{d}^{\prime}\\}\cap S^{d-1}$.

Recall that a Steiner symmetrization of a set $A\subset\mathbb{R}^{d}$
relative to a hyperplane $L$ is the set

${\rm
St}(A):=\bigcup\limits_{\ell\in\mathcal{Q}}\left(\frac{1}{2}(A\cap\ell)+\frac{1}{2}(\widetilde{A}\cap\ell)\right),$

where $\mathcal{Q}$ is the set of all lines passing through points of $A$ that
are perpendicular to hyperplane $L$ and $\widetilde{A}$ is the reflection of
$A$ with respect to $L$. The Steiner symmetrization of ${\rm cone}\\{{\bf
y}_{1},\ldots,{\bf y}_{d}\\}$ is contained in ${\rm cone}\\{{\bf y}_{1},{\bf
y}_{2},{\bf y}_{3}^{\prime},\ldots,{\bf y}_{d}^{\prime}\\}$. Since $L$ passes
through ${\bf 0}$, whenever $A,D\subset\mathbb{R}^{d}$ are such that ${\rm
St}(A)\subset D$, we have ${\rm St}(A\cap B^{d})\subset D\cap B^{d}$.
Therefore, ${\rm St}(\Omega(Y))\subset\Omega(Y^{\prime})$, where
$\Omega(Y^{\prime}):={\rm cone}\\{{\bf y}_{1},{\bf y}_{2},{\bf
y}_{3}^{\prime},\ldots,{\bf y}_{d}^{\prime}\\}\cap B^{d}$. Since Steiner
symmetrization preserves the volume, we have

$\begin{split}{\rm Vol}_{d-1}(r(Y))&={d}{{\rm Vol}_{d}(\Omega(Y))}={d}{{\rm
Vol}_{d}[{\rm St}(\Omega(Y))]}\\\ &\leq{d}{{\rm
Vol}_{d}(\Omega(Y^{\prime}))}={\rm Vol}_{d-1}(r(Y^{\prime})).\end{split}$

Since the $(d-1)$-dimensional volume of $r(Y)$ is positive, so is the one of
$r(Y^{\prime})$ and, hence of $Y^{\prime}$. Then $Y^{\prime}$ is a
$(d-1)$-simplex. Since ${\bf y}_{3}\notin L$, its projection ${\bf
y}^{\prime}_{3}$ is an interior point of $H\cap B^{d}$. Then there is a point
${\bf y}_{3}^{\prime\prime}\in H\cap B^{d}$ into which we can move ${\bf
y}^{\prime}_{3}$ so that the new simplex $Y^{\prime\prime}$ with vertices
${\bf y}_{1},{\bf y}_{2},{\bf y}_{3}^{\prime\prime},{\bf
y}_{4}^{\prime}\ldots,{\bf y}_{d}^{\prime}$ contains $Y^{\prime}$ together
with some open set that lies outside of $Y^{\prime}$. Then we will have ${\rm
Vol}_{d-1}(r(Y^{\prime}))<{\rm Vol}_{d-1}(r(Y^{\prime\prime}))$ with
$Y^{\prime\prime}\subset H\cap B^{d}$. The simplex $Y^{\prime\prime}$ is
contained in some $(d-1)$-simplex $Y^{\prime\prime\prime}$ inscribed in $H\cap
S^{d-1}$. Then ${\rm Vol}_{d-1}(r(Y))<{\rm
Vol}_{d-1}(r(Y^{\prime\prime\prime}))$. This contradicts the assumption that
$r(Y)$ has the largest volume over all $(d-1)$-simplices $Y$ inscribed in
$H\cap S^{d-1}$.

Thus, if ${\rm Vol}_{d-1}(r(Y))$ is maximal, then $Y$ must be a regular
simplex. At the same time, for all regular $(d-1)$-simplices $V$ inscribed in
$H\cap S^{d-1}$, the volumes ${\rm Vol}_{d-1}(r(V))$ are the same. Then ${\rm
Vol}_{d-1}(r(V))$ is maximal.

## References

  * [1] Ambrus, Analytic and probabilistic problems in discrete geometry, Ph.D. Thesis, University College London, London, 2009.
  * [2] G. Ambrus, K.M. Ball, T. Erdélyi, Chebyshev constants for the unit circle, Bull. Lond. Math. Soc. 45 (2013), no. 2, 236–248.
  * [3] G. Ambrus, S. Nietert, Polarization, sign sequences and isotropic vector systems, Pacific J. Math. 303 (2019), no. 2, 385–399.
  * [4] L. Bétermin, M. Faulhuber, S. Steinerberger, A variational principle for Gaussian lattice sums, https://arxiv.org/pdf/2110.06008.pdf.
  * [5] K. Böröczky Jr., On an extremum property of the regular simplex in $S^{d}$. In: Intuitive Geometry, K. Böröczky, G. Fejes Tóth (eds.), Colloq. Math. Soc. J’anos Bolyai 48, Elsevier, 117–121, 1987.
  * [6] K. Böröczky, Jr. Finite packing and covering, Cambridge University Press, 2004.
  * [7] K. Böröczky, Jr., G. Wintsche, Covering the sphere by equal spherical balls, Discrete Comput. Geom. Algorithms and Combinatorics 25 (Springer, Berlin, 2003), 235–251.
  * [8] S.V. Borodachov, Polarization problem on a higher-dimensional sphere for a simplex, Discrete and Computational Geometry 67 (2022), no. 2, 525–542.
  * [9] S.V. Borodachov, Min-max polarization for certain classes of sharp configurations on the sphere (submitted), https://arxiv.org/abs/2203.13756.
  * [10] S.V. Borodachov, Absolute minima of potentials of a certain class of spherical designs (in preparation, presented at the Workshop ”Optimal Point Configurations on Manifolds”, ESI, Vienna, January 17–21, 2022. https://www.youtube.com/watch?v=L-szPTFMsX8).
  * [11] S.V. Borodachov, D.P. Hardin, E.B. Saff, Discrete Energy on Rectifiable Sets, Springer Monographs in Mathematics, 2019.
  * [12] N. Bosuwan, P. Ruengrot, Constant Riesz potentials on a circle in a plane with an application to polarization optimality problems, ScienceAsia 43 (2017), 267–274.
  * [13] P.G. Boyvalenkov, P.D. Dragnev, D.P. Hardin, E.B. Saff, M.M. Stoyanova, On polarization of spherical codes and designs (submitted), https://arxiv.org/abs/2207.08807.
  * [14] L. Dalla, D.G. Larman, P. Mani-Levitska, C. Zong, The blocking numbers of convex bodies. Discret. Comput. Geom. 24 (2000), no. 2–3, 267–277. The Branko Grünbaum birthday issue.
  * [15] T. Erdélyi, E.B. Saff, Riesz polarization inequalities in higher dimensions, J. Approx. Theory 171 (2013), 128–147.
  * [16] L. Fejes Tóth, Lagerungen in der Ebene, auf der Kugel und im Raum, (German) Die Grundlehren der mathematischen Wissenschaften in Einzeldarstellungen mit besonderer Berücksichtigung der Anwendungsgebiete, Band LXV. Springer-Verlag, Berlin-Göttingen-Heidelberg, 1953.
  * [17] L. Fejes Tóth, Regular figures, A Pergamon Press Book The Macmillan Company, New York, 1964.
  * [18] L. Fejes Tóth, Kreisüberdeckungen der Sphäre, Stud. Sci. Math. Hung. 4 (1969), 225–247.
  * [19] S.I. Galiev, Multiple packings and coverings of a sphere, Diskret. Mat. 8 (1996), no. 3, 148–160.
  * [20] D.P. Hardin, A.P. Kendall, E.B. Saff, Polarization optimality of equally spaced points on the circle for discrete potentials, Discrete Comput. Geom. 50 (2013), no. 1, 236–243.
  * [21] H.L. Montgomery, Minimal theta functions, Glasgow Math. J. 30 (1988), no. 1, 75–85.
  * [22] N. Nikolov, R. Rafailov, On the sum of powered distances to certain sets of points on the circle, Pacific J. Math. 253 (2011), no. 1, 157–168.
  * [23] N. Nikolov, R. Rafailov, On extremums of sums of powered distances to a finite set of points, Geom. Dedicata 167 (2013), 69–89.
  * [24] C.A. Rogers, Packing and covering, Cambridge University Press, 2008.
  * [25] K. Schütte, Überdeckungen der Kugel mit höchstens acht Kreisen, Math. Ann. 129 (1955), 181–186.
  * [26] K.B. Stolarsky, The sum of the distances to certain pointsets on the unit circle, Pacific J. Math. 59 (1975), no. 1, 241–251.
  * [27] K.B. Stolarsky, The sum of the distances to $N$ points on a sphere, Pacific J. Math. 57 (1975), no. 2, 563–573.
  * [28] Y. Su, Discrete minimal energy on flat tori and four-point maximal polarization on $S^{2}$, Ph.D. Thesis, Vanderbilt University, Nashville, TN, 2015.
  * [29] L. Wimmer, Covering the sphere with equal circles, Discrete and Computational Geometry 57 (2017), no. 4, 763–781.


algorithm[ALG.]

††thanks: Correspondence should be addressed to
<EMAIL_ADDRESS>
# Sensorless wavefront correction in two-photon microscopy across different
turbidity scales

Maximilian Sohmen *  — Institute of Biomedical Physics, Medical University of
Innsbruck, Müllerstraße 44, 6020 Innsbruck, Austria  Molly A. May *  —
Institute of Biomedical Physics, Medical University of Innsbruck, Müllerstraße
44, 6020 Innsbruck, Austria  Nicolas Barré *  — Institute of Biomedical
Physics, Medical University of Innsbruck, Müllerstraße 44, 6020 Innsbruck,
Austria  Monika Ritsch-Marte *  — Institute of Biomedical Physics, Medical
University of Innsbruck, Müllerstraße 44, 6020 Innsbruck, Austria  Alexander
Jesacher *  — Institute of Biomedical Physics, Medical University of
Innsbruck, Müllerstraße 44, 6020 Innsbruck, Austria

###### Abstract

Adaptive optics (AO) is a powerful tool to increase the imaging depth of
multiphoton scanning microscopes. For highly scattering tissues, sensorless
wavefront correction techniques exhibit robust performance and present a
straight-forward implementation of AO. However, for many applications such as
live-tissue imaging, the speed of aberration correction remains a critical
bottleneck. Dynamic Adaptive Scattering compensation Holography (DASH) – a
fast-converging sensorless AO technique introduced recently for scatter
compensation in nonlinear scanning microscopy – addresses this issue. DASH has
been targeted at highly turbid media, but to-date it has remained an open
question how it performs for mild turbidity, where limitations imposed by
phase-only wavefront shaping are expected to impede its convergence. In this
work, we study the performance of DASH across different turbidity regimes, in
simulation as well as experiments. We further provide a direct comparison
between DASH and a novel, modified version of the Continuous Sequential
Algorithm (CSA) which we call Amplified CSA (a-CSA).

##  I Introduction

Dynamic wavefront correction is a powerful approach for extending the imaging
depth of nonlinear microscopy in scattering tissues such as the mouse brain.
However, to make the benefits of this technology accessible to an even broader
range of applications, some remaining limitations need to be overcome. A
particular problem is posed by the speed at which an aberration-compensating
pattern can be retrieved, as imaging into live specimens requires to outpace
the decorrelation time imposed by the constantly changing tissue.
To this end, we have recently introduced “Dynamic Adaptive Scattering
compensation Holography” (DASH), a fast-converging indirect wavefront sensing
algorithm for scatter correction in nonlinear scanning microscopy [1]. DASH
employs a programmable phase-diffractive element to shape two beams
simultaneously: a test beam, whose wavefront is varied in each step to explore
possible signal improvements, and a corrected beam, whose wavefront is
continuously improved using the information gained from interference with the
test beam. DASH employs phase-only wavefront shaping, which bears the
advantage of power efficiency. On the other hand, discarding amplitude
information introduces errors in the resulting wavefronts. In this work we
investigate the impact of these errors on the correction performance of DASH
for different regimes of turbidity. Furthermore, we compare DASH to an
alternative approach, which does not suffer from the “phase-only” restriction:
a modified version of the Continuous Sequential Algorithm (CSA) [2], which
follows a pixel-by-pixel testing approach. Our modification to CSA consists in
amplifying the intensity of the tested pixel relative to all other pixels,
which increases signal contrast. This amplification allows the application of
CSA in situations with low signal-to-noise ratio (SNR), as often encountered
in practice. We refer to our modified CSA algorithm as _Amplified_ Continuous
Sequential Algorithm, or a-CSA.

This work is structured as follows: In Section II we discuss how turbidity can
be quantified, how a scattering medium of specific turbidity can be modelled
numerically, and how in a two-photon imaging experiment a scattering volume of
‘tunable’ mean free path can be emulated by a phase-mask displayed on a
spatial light modulator (SLM). In Section III, the principles behind DASH and
a-CSA are explained. In Section IV, we detail the implementation of a-CSA and
DASH in our numerical simulations (Section IV.2), our experiment (Sections
IV.3), and present a comparison of algorithm performance for two-photon
excited fluorescence (TPEF) imaging of a homogeneous ‘dye-slide’ sample in
different regimes of turbidity (Section IV.4) as well as mouse brain (Section
IV.5). In Section IV.6, implications of a low SNR on algorithm performance are
discussed.

##  II Quantifying turbidity

It is essential for the present work to define what we mean when speaking of
“low” or “high” turbidity. The scattering properties of materials and tissues
are often quantified using the scattering mean free path $l_{s}$, i.e., the
expectation value of a photon’s free travelling path before it is scattered.
This is defined via the Beer-Lambert law,
$|U_{0}(L)|^{2}=|U_{0}(0)|^{2}~{}\mathrm{e}^{-L/l_{s}}$, where
$|U_{0}(L)|^{2}$ represents the intensity of the _non-scattered_ light after
travelling a distance $L$, and $|U_{0}(0)|^{2}$ the incident light intensity.
The transport mean free path $l_{t}$ takes scattering anisotropy into account:
$l_{t}=l_{s}/(1-g)$, where $g=\langle\cos\theta\rangle$ is the expectation
value of the cosine of the scattering angle $\theta$. For instance, in a
material which predominantly scatters into the forward direction (causing
small scattering angles), $l_{t}$ is much larger than $l_{s}$. Conversely, in
an isotropic scatterer $l_{t}=l_{s}$. Typical values of $l_{t}$ for brain
tissue range between a few tens to hundreds of micrometers [3, 4, 5].

Our goal is to model the effect of a (in general _three-dimensional_ , 3D)
scattering medium on a light field propagating in positive $z$-direction by a
_two-dimensional_ (2D) phase mask, located at axial position
$z=z_{\text{scat}}$, with transmission function
$\exp\left(\mathrm{i}\Phi(\boldsymbol{\rho})\right)$. Here,
$\Phi(\boldsymbol{\rho})$ denotes the scattering-related phase shifts
experienced by a field point at the 2D lateral coordinate $\boldsymbol{\rho}$.
Starting from this ansatz, we note the relation between
$\Phi(\boldsymbol{\rho})$ and the overlap integral (OI) of scattered field
$U_{\text{scat}}$ and incident field $U_{0}$:

$\displaystyle\text{OI}\left[U_{\text{scat}},U_{0}\right]$
$\displaystyle=\left|\int
U_{\text{scat}}(\boldsymbol{\rho})\,U^{*}_{0}(\boldsymbol{\rho})\,\mathrm{d}\boldsymbol{\rho}\right|^{2}$
$\displaystyle=\left|\int\left|U_{0}(\boldsymbol{\rho})\right|^{2}\mathrm{e}^{\mathrm{i}\Phi(\boldsymbol{\rho})}\,\mathrm{d}\boldsymbol{\rho}\right|^{2},$
(1)

where, assuming that all photons are scattered into the forward direction and
absence of absorption, the fields must be normalised such that
$\int|U_{\text{scat}}|^{2}\,\mathrm{d}\boldsymbol{\rho}=1=\int|U_{0}|^{2}\,\mathrm{d}\boldsymbol{\rho}$.
The equality (1) is most intuitive if the integral is evaluated in the plane
of the 2D scattering mask, but for freely propagated fields the OI in fact
stays constant in all transverse planes at $z\geq z_{\text{scat}}$.

Equivalently, the OI can be interpreted as the unscattered intensity fraction,

$\text{OI}[U_{0}(L),U_{0}(0)]=\mathrm{e}^{-L/l_{t}}.$ (2)

For instance, if $U_{0}(L)$ and $U_{0}(0)$ are orthogonal ($\text{OI}=0$),
every photon must have scattered at least once. Conversely, one can assume
that no photon has scattered if $U_{0}(L)=U_{0}(0)\Rightarrow\text{OI}=1$
(neglecting the rare case of successive scattering events that exactly
compensate).
Combining Eqs. 1–2 we can quantify a computed phase mask in terms of the
corresponding “thickness” expressed in units of the transport mean free path
$l_{t}$:111We note that the relation (3) is consistent with the considerations
made in Ref. [6] (see Eq. 4 therein), which lead to the derivation of the
scattering-phase theorem.

$L/l_{t}=-\ln(\text{OI})=-\ln\left(\left|\int\left|U_{0}(\boldsymbol{\rho})\right|^{2}\mathrm{e}^{\mathrm{i}\Phi(\boldsymbol{\rho})}~{}\mathrm{d}\boldsymbol{\rho}\right|^{2}\right)$
(3)

For the case of dominant forward scattering and negligible absorption, this
relation allows us to compute a 2D phase mask $\Phi(\boldsymbol{\rho})$ that
leads to a speckle pattern in the object plane which is in many ways similar
to that of a voluminous 3D scatter medium of the same transport mean free path
$l_{t}$. In the experiments described later in this work, we will exploit this
fact to simulate different regimes of turbidity by displaying computed 2D
scatter masks of specific $l_{t}$ on an SLM. Of course the equivalence between
a 3D and a 2D scatterer – even if they exhibit the same $l_{t}$ – does not
encompass all physical properties; for instance, the isoplanatic patch (i.e.,
the ‘corrected field of view’) obtained through an AO wavefront correction
will be smaller for a 3D than for a 2D scatterer. However, concerning the
aspects studied in this work (e.g., the algorithm convergence at a single
field point), a 3D and a 2D scatterer of same $l_{t}$ can be regarded as
equivalent.

##  III Approaches for sensorless AO in nonlinear microscopy

Most indirect (or sensorless) AO schemes construct an aberration compensation
phase mask from measurements of the TPEF signal for many different test
patterns or “modes” $M_{n}$ displayed on an SLM, $n\in\\{1,\ldots,N\\}$
denoting the mode index. These test patterns can be pre-defined, for instance
as a set of Zernike [7] or Hadamard modes [8, 9], or directly derived from the
signal of previous test patterns, such as in genetic algorithms [10]. The
phase patterns are commonly imaged into the back focal plane (BFP) of the
objective lens, where – by the Fourier transform property of the lens – the
test patterns have a homogeneous effect on the point-spread function (PSF)
over the entire focal plane. However, there exist also sample-conjugated
schemes, where the SLM is imaged directly into the sample space [11, 12, 13,
14].

Figure 1: Sensorless AO methods for nonlinear microscopy. a-CSA uses a
superpixel basis, DASH a plane-wave basis of test modes on a fluorescent
target. Respective examples of a test mode (bottom) and the compensation
pattern at the corresponding step are shown on the left and right. For a-CSA
the order of tested modes is by increasing distance to the pupil center, for
DASH the order is by increasing spatial frequency. During phase-stepping, the
relative phase $\phi_{p}$ of the test mode with respect to the compensation
pattern is varied (see main text).

A particular case is the “traditional adaptive optics” regime, where the given
phase aberration $\Phi(\boldsymbol{\rho})$ can be approximated by a weighted
sum of the available test modes,
$\Phi(\boldsymbol{\rho})=\sum_{n=1}^{N}a_{n}\,M_{n}(\boldsymbol{\rho})$, and
each $a_{n}$ has a rather small phase magnitude (up to about 1 rad RMS over
the pupil). Here, $\boldsymbol{\rho}$ is the 2D coordinate vector in the BFP.
In this case, the TPEF signal $S$ generated in the focus has a smooth
dependence on the mode magnitudes $a_{n}$, allowing to approximate it by a
simple function such as a multi-dimensional Lorentzian or parabola, e.g.,
$S\propto 1-\sum_{n}\sum_{m}\alpha_{nm}\,a_{n}\,a_{m}$. In many cases, the
cross-talk matrix $\alpha$ can be diagonalized by choice of an adequate mode
basis [15]. Then, $N+1$ measurements can be sufficient to characterize the
phase aberration [16], although usually $2N+1$ or even $3N$ measurements are
taken. It is important to note that this particular case does not necessarily
coincide with low turbidity (i.e., a small value of $L/l_{t}$), since a large
number of modes, even if their individual magnitudes are small, can still sum
up to a large total aberration.
Outside the traditional regime, in what is often called the “scattering”
regime, it is usually required to take many more measurements to find a
suitable corrective phase mask. Typically, in the scattering regime the
complexity of the aberration is beyond the capabilities of the correction
device, meaning that the number of scattering modes is larger than the number
of correctable ones. Additional reasons which can hinder full aberration
compensation include an SNR too small to measure the contributions of less
significant modes, or simply lack of time to measure all of the (many)
scattering modes. Nevertheless, it has been shown that even correcting only a
limited number of scattering modes can suffice to restore a diffraction-
limited focus, which in this context can be understood as an intensity-
enhanced speckle [17].
There exist several approaches for sensorless AO which can operate in the
scattering regime [2, 18, 10, 19, 20, 1, 21]. Two of them, which we have
identified as fast converging in numerical simulations, are studied in this
work and outlined in Fig. 1.

The first approach, known as continuous sequential algorithm (CSA) [2],
operates on a single-pixel basis: sequentially, each pixel’s phase is adjusted
to maximize the signal. Although a basis of single pixels appears advantageous
due to its intrinsic orthogonality, the approach is usually impractical
because the interference contrast from single-pixel phase modulation is
typically buried in noise. This inspired alternative approaches, such as
hybrid methods using larger “superpixels”, each of them featuring an internal
spatial phase pattern [20, 22]. However, we point out that the SNR problem of
CSA can also be remedied by simply amplifying the intensity in the test
superpixel, as this increases the interference contrast. This variant, which
we call a-CSA (“amplified CSA”), requires the laser power on each superpixel
to be controllable. Then, as we show in Section IV.3, a-CSA can work in
experiments where only a comparably small number of modes needs to be
corrected.

The mode basis of the second approach, DASH, is a set of plane waves which
are, in a sense, Fourier-related to the single-pixel basis. Importantly, for
this approach the SNR is uncritical, as the power fraction contained in each
plane wave can be adapted directly via the SLM hologram. This makes this
method more practical than CSA. However, to shape test beam and corrected beam
without introducing artefacts, the SLM would be required to manipulate both,
the spatial phase and amplitude distribution of the incident laser beam, just
as for a-CSA. In DASH, the amplitude part is disregarded and only the phase
part is included in the hologram, ensuring high power efficiency but
inevitably leading to errors in the generated wavefronts and undesired
diffraction orders. An example of the effects of phase-only shaping is given
in Fig. 2. The left image shows a target arrangement of 256 spots with random
phases we wish to create in the focal plane using an SLM in the Fourier plane.
Discarding the amplitude part of the synthetic hologram (a superposition of
the corresponding 256 plane waves) leads to the image on the right. The
standard deviation of the phases at the target sites is about 0.25 rad,
corresponding to about 4% of the wavelength, and the relative amplitude error
at the target sites is about 20%. Additionally, we see that weak “ghost spots”
appear outside the target square.

In line with our earlier works [1, 14], we will in the following for both
a-CSA and DASH denote by $f$ the fraction of the total pupil intensity which
is contained in the test mode.

Figure 2: Artefacts related to phase-only light shaping. Field modulus of a
target pattern (left) and the simulated pattern reconstructed by phase-only
wavefront shaping in the Fourier plane (right), illustrating introduced
artefacts.

##  IV Comparison of algorithm performance

In this section we compare the performance of DASH and a-CSA over different
scales of turbidity in numerical simulations and experiment.

###  IV.1 General recipe

The procedure for our systematic comparison between DASH and a-CSA is as
follows.

  1. 1.

Both algorithms are initialised with an empty correction mask.

  2. 2.

We calculate a scattering phase mask that emulates a certain degree of
turbidity and is kept constant throughout each algorithm run. In the
simulations, this scattering mask is simply added to the BFP-conjugate
wavefront during image propagation; in the experiment, it is displayed on our
SLM.

  3. 3.

We imprint a test mode on the excitation beam and carry out a phase-stepping
procedure, i.e., measure the TPEF signal generated in the object plane for $P$
different global phase offsets $\phi_{p}=2\pi\,p/P$ (with $p=1\ldots P$)
applied to the test mode.222In the experiment we typically set $P=5$; in the
simulations, especially for the noise-free case, we find that $P$ can be
reduced to $3$ for faster computation without loss of performance. From the
phase-stepping procedure, we extract the parameters for which the test mode
maximizes the TPEF signal.

  4. 4.

We directly update our correction mask by including the tested mode with the
retrieved optimal parameters. In the simulations, this correction mask is
added to the wavefront incident on the scatterer; in the dye-slide experiment,
it is displayed on our SLM, on top of the scattering mask.

  5. 5.

We now go back to step 3, test the next mode, update the correction mask, and
so forth, until the correction mask contains the full number of test modes
with their optimal contributions.

  6. 6.

Once the full set of modes has been tested, we can start again from the first
test mode for another correction run. We typically run 3 full iterations,
which has shown to ensure algorithm convergence in most cases.

Detailed algorithm descriptions are provided in the Supplementary
Material.333This also includes a minor modification to the DASH algorithm
compared to Ref. [1].

###  IV.2 Numerical model

In our numerical model, we consider aberrations defined on a square grid of
64${\times}$64 pixels accounting for $N^{2}_{\text{scat}}={}$4096 scattering
modes. The aberration phase masks are calculated by adding cosines of all
spatial frequencies supported by the grid size with uniformly distributed
random phases. The amplitudes of the cosines follow a Gaussian weighting with
standard deviation $\sigma$, such that cosines with higher spatial frequencies
are increasingly attenuated, as expected for most scattering scenarios found
in nature. By varying $\sigma$ we can hence tune the spatial frequency content
of the model scatterer. For details on our specific implementation of
calculating the scatter mask, the reader is referred to our Supplementary
Material (Algorithm 4).

The simulated SLM for (square) correction patterns features 32${\times}$32
pixels. Both the scattering mask as well as the simulated SLM are located in
the BFP, in Fourier relation to the object plane, and have a side length
measuring $2k_{0}\,\text{NA}$, where $k_{0}$ is the vacuum wavevector of light
and $\text{NA}={}$0.7 the numerical aperture. Each set of modes is tested
repeatedly for 3 iterations. In our experiments, the SLM illumination does not
feature an ideal, flat-top intensity profile, but exhibits a (weak) Gaussian
shape of $\mathrm{e}^{-1/2}$ distance approximately equal to the pupil radius.
To reflect the experimental situation, also in our simulations we assume a
light intensity distribution across the BFP which is symmetrically Gaussian
and of corresponding width.
The fluorescent sample is assumed to be a homogeneous, absorption-free
fluorescent layer of $d_{s}=10$ µm thickness.444NB that a thin 3D fluorescent
sample layer as assumed here leads to smaller signal enhancements and slower
algorithm convergence than an (infinitely thin) 2D layer. In the limit of a
homogeneous 3D sample volume with thickness much larger than the Rayleigh
range of the focused beam, the second-order nonlinearity of TPEF imaging is
insufficient for providing any signal enhancement at all [23]. We model this
sample volume by grid points on $N_{s}$ distinct 2D planes with axial
interspaces of $d_{s}/(N_{s}-1)$, and the nominal focus plane located at its
center. Normally, we set $N_{s}=6$ layers (interspaced by 2 µm), since a
larger number of planes does not notably improve the accuracy of results and
we have checked that the exact choice of $N_{s}$ is uncritical for our
conclusions. To simulate the PMT counts, we simply propagate the light field
into the $N_{s}$ planes, calculate the squared intensity on each grid point,
and take the sum over all points.

In Section IV.4, we show simulation results at high-SNR conditions, for which
photon shot noise can be neglected. For scenarios with low SNR, as discussed
in Section IV.6, we simulate shot noise by varying the summed PMT counts
according to Poisson statistics.

###  IV.3 Experimental procedure

Figure 3: Experimental setup. SLM = spatial light modulator; PMT =
photomultiplier tube. Inset: Exemplary test phase mask for a-CSA. Every
superpixel except the one under test is superposed by a shallow blazed grating
to diffract the desired fraction of light intensity off the optical axis where
it is blocked by an iris aperture. This enables amplitude tunability while
using a standard, phase-only SLM.

We compare DASH and a-CSA experimentally in a home-built TPEF microscope
featuring a phase-only SLM (Meadowlark HSP1920-500-1200-HSP8,
1920${\times}$1152 pixels of side length 9.2 µm) located in a BFP-conjugate
plane of the excitation path. This SLM serves two purposes in parallel: First,
by displaying a ‘scattering’ phase mask of defined transport mean free path
$l_{t}$ (see Section II) it allows to _emulate_ the effect of a scattering
medium in the light path. Second, by running our sensorless wavefront
correction schemes and displaying the retrieved phase compensation patterns
_on top_ of the given scattering mask, we can test and compare the algorithm
performances.
A sketch of the setup is shown in Fig. 3. For excitation, we use a mode-locked
Ti:sapphire laser (MaiTai DeepSee, Spectra Physics) with emission maximum set
to 800 nm wavelength. The epi-TPEF is collected by a water-dipping objective
(Olympus XLUMPLFLN20XW, $\text{NA}=1$) and directed towards a photomultiplier
tube (PMT, Hamamatsu H10769A-40) using a dichroic beam splitter (AHF
Analysentechnik, HC665LP) and an additional emission filter in front of the
PMT (AHF, ET680SP-2P8). In our measurements, to operate with a more
homogeneous pupil illumination, we artificially reduce the NA to about 0.7
using an aperture in a pupil-conjugate plane (not shown in Fig. 3).

For our systematic comparison, the sample consists of a thin layer of
rhodamine solution sandwiched between a glass slide and a coverslip of 1 mm
and 170 µm thickness, respectively. Before each algorithm run, we perform an
initial precorrection run to compensate for optical-system wavefront
distortions such as the spherical aberrations introduced by the coverslip.
Starting from this precorrection, the wavefront correction algorithms
afterwards only have to correct for the artificial scatterer displayed on our
SLM, ensuring compatibility with the numerically simulated correction runs.

For DASH, we typically set the intensity in the test mode to 30% of the total
intensity (i.e., $f=0.3$, see Supplementary Material). For a-CSA, we amplify
the relative power in the test-superpixel by superposing _all other_ a-CSA-
superpixels (consisting of many physical SLM pixels) with a blazed grating of
defined modulation depth [24]. Thus, for all except the test-superpixel, power
is diffracted away from the optical axis and cut by an aperture as indicated
in Fig. 3. It is important that the mean phase of the blazed grating is nought
to prevent an effect on the zero-order beam, which carries the phase
mask.555Note that by encoding our phase masks in the zero order and dumping
power in higher orders we avoid dispersion of the femtosecond laser pulses.
Our blazed gratings have a period of 4 SLM-pixels and a modulation depth of
approximately $\pi$ rad, resulting in $(1-\beta)\sim 50$% of the laser power
being dumped (measured value) and a power fraction $f\approx 1/(\beta N^{2})$
in the case of $N{\times}N$ superpixels. Naturally, dumping power on the iris
also decreases the total TPEF signal, wherefore we have to compensate by
increasing the total laser power. Since in practice the available power is
necessarily finite, this approach for amplitude modulation is only applicable
in regimes of turbidity which do not require too many (i.e., too small) a-CSA-
superpixels.

For compatibility, we conduct our experiments with the same number of test
modes as in the simulations for the different regimes, i.e. 16, 64, and 256
modes for low, medium, and high turbidity, respectively. Our SLM holograms
measure 560$\times$560 SLM-pixels and entirely fill the reduced objective
pupil. For a-CSA this means that we form 4$\times$4 square superpixels (each
140$\times$140 SLM-pixels) for the low turbidity, 8$\times$8 superpixels (each
70$\times$70 SLM-pixels) for the medium turbidity, and 16$\times$16
superpixels (each 35$\times$35 SLM-pixels) for the high turbidity scenario.
Analogous to the simulation, DASH and a-CSA are executed for three full
iterations. The results are summarized in Fig. 4.

###  IV.4 Results systematic comparison

For our systematic comparison, we study three different scenarios (A, B, C)
concerned with increasing levels of turbidity.

In Scenario A we study low turbidity, with an effective scatterer thickness of
a single transport mean free path, $L/l_{t}=1$, and a spatial frequency
distribution of the scatterer chosen accordingly narrow, $\sigma=1$. An
example of such a scatter mask is shown on the top left of Fig. 4 (A). In the
focal plane, such a mild scatterer typically leads to an intensity
distribution which is still spatially confined and only moderately deviates
from the aberration-free focus, as shown in Fig. 4 (A), bottom left. Here, the
intensity scale is normalized to the peak intensity of the aberration-free
focus, such that the maximum value (typically around 0.4 for Scenario A)
equals the respective Strehl ratio. In this low-turbidity case it is expected
that only a small number of modes is needed for adequate compensation,
wherefore we correct $4{\times}4=16$ modes. The plots in the center and right
column of Fig. 4 (A) show the TPEF signal enhancement simulated numerically
(center) and measured experimentally (right), respectively, for a-CSA (blue)
and DASH (orange). Specifically, we plot the TPEF signal measured when the
established compensation pattern is applied at the respective measurement
number.666Note that these curves are often similar, but in general not
identical to the signal collected during an algorithm run, as the latter
results from an interference between a partly compensated and a test beam.
In the numerical simulations (Fig. 4 (A), center), we observe that for a-CSA
the signal initially increases fast, but then flattens off from the second
iteration onwards. DASH, in contrast, exhibits a less monotonic signal
evolution than a-CSA, first increasing more slowly, but ultimately surpassing
the final signal level of a-CSA. We attribute the signal flattening and lesser
total performance of a-CSA to the fact that its compensation patterns are
inherently displayed at the native resolution of the test-mode basis (i.e.,
4$\times$4 superpixels in Scenario A) and therfore do not exploit the full
resolution supported by the SLM (32$\times$32 pixels in the simulation). This
coarser spatial discretization leads to increased diffraction losses compared
to DASH (which always exploits the full SLM resolution, regardless of the
number of tested modes) during the course of the optimisation. Additionally,
we observe that for DASH the steepest changes in signal usually occur at the
beginning of each iteration. This is due to the fact that we intentionally
sort our test modes in ascending order with regard to the angle between their
propagation direction and the optical axis. Since small scattering angles are
typically more dominant in fluorescence imaging settings, this tends to speed
up the algorithm convergence. Both methods achieve comparable final Strehl
ratios around 0.75 and 0.85 for a-CSA and DASH, respectively (see Supplemental
Material).
Our experimental measurements (Fig. 4 (A), right) show the same general
trends. When comparing absolute values of signal enhancement, it is important
to keep in mind the critical dependence on the sample thickness: in our
simulations, e.g., we observe for increasing the thickness as 0 $\rightarrow$
4 µm ${\rightarrow}$ 10 µm (at constant fluorophore density) a tendency of
decreasing final enhancements of 2.14 $\rightarrow$ 1.32 $\rightarrow$ 1.29 in
case of a-CSA and 2.9 $\rightarrow$ 1.8 $\rightarrow$ 1.5 in case of DASH.
Given that in our experiments the sample layer thickness is controllable and
measurable only with limited accuracy, it seems fair to claim good agreement
between the observed values in simulation and experiment.

Figure 4: Scatter correction in different turbidity regimes. The three
scenarios A–C correspond to an increasing degree of scattering with (A)
$L/l_{t}=1,\sigma=1$, (B) $L/l_{t}=3,\sigma=3$, and (C) $L/l_{t}=5,\sigma=5$,
respectively. The left column shows examples of the phase aberration applied
in the objective pupil and the resulting TPEF intensity distribution in the
focal plane. The plots in the center and right column show how the TPEF signal
improves during correction in simulation (center) and experiment (right).
Here, solid curves represents the mean value over 5 repeated runs in the
simulation or experiment, each initialised with a different random scatterer,
and the ribbons represent the respective standard error of the mean. Kinks in
the experimental a-CSA curves (blue asterisks) are caused by the circular,
underfilled pupil (see main text). Additional information, such as retrieved
SLM compensation masks and resulting focal spots, is provided in the
Supplementary Material.

In Scenario B, we assume medium turbidity with $L/l_{t}=3$ and an intermediate
contribution of modes of higher spatial frequency, $\sigma=3$. We increase the
number of correctable modes to 64 to account for this greater spatial
frequency content. Before correction, the typical Strehl ratios in this
scenario are on the order of 5% (Fig. 4 (B), bottom left).
Our numerical simulations (Fig. 4 (B), center) indicate that again, a-CSA
initially improves faster than DASH, but levels off at a lower final signal.
The total signal enhancement achievable by both correction algorithms is
higher than for the low-turbidity case; the final Strehl ratios are around 0.4
for a-CSA and 0.6 for DASH. We suspect that the initial delay of DASH in
comparison to a-CSA is related to the appearance of higher diffraction orders
(“ghost foci” as shown in Fig. 2) stemming from the phase-only field shaping
of DASH. This is one of several reasons which lead us to believe that the use
of a complex field-shaping technique might bear great potential for
improvement of algorithm performance in the future.
In our experimental measurements (Fig. 4 (B), right), we observe similar
behaviour. Initially, a-CSA improves faster than DASH, but ultimately DASH
delivers higher signal enhancement. For a-CSA, we observe a kink in the signal
enhancement at the transition between two iterations (marked by the blue
asterisk in Fig. 4), which is not present in the simulations. We attribute
this kink to the circular pupil in our experiment which cuts power from
superpixels located near the corners of the square SLM pattern. Since we step
through the a-CSA superpixels in order of their distance to the pupil center,
this effect is most dominant towards the end of each iteration.

In Scenario C, we assume high turbidity, with $L/l_{t}=5$ and $\sigma=5$,
where without correction typical Strehl ratios are on the order of 1%. In this
scenario, we correct 256 modes.
In the numerical simulations (Fig. 4 (C), center), we observe that the initial
speed advantage of a-CSA compared to DASH decreases, but the final
enhancements achievable with both algorithms are even higher than in Scenario
B. Comparing between the two algorithms, our simulations again deliver a
better performance for DASH than for a-CSA. Final Strehl ratios are around
0.35 for a-CSA and 0.45 for DASH.
These trends are well supported by our experimental data (Fig. 4 (C), right).

Graphical animations of our simulated correction runs are provided as GIF
files (see Supplementary Material). These animations show the phase patterns
displayed on the SLM during the correction algorithm as well as the evolution
of the focal plane intensity distribution.

###  IV.5 Experimental comparison for a biological sample

Both DASH and a-CSA can offer striking quality improvements for imaging
through turbid media, such as in microscopy of layers deep inside tissue. To
demonstrate this, we use our TPEF scanning microscope to image microglia
expressing green fluorescent protein (CX3Cr-1GFP, cf. Ref. [1]), located 200
µm deep inside the corpus striatum of a mouse brain slice fixed via perfusion.
For this, we adjust our wavelength to the excitation maximum at 900 nm,
retrieve a precorrection for optical system aberrations by focussing on
microglia directly below the coverslip, and subsequently move the objective
focus mechanically 200 µm down into the brain tissue.

Figure 5 (A) shows an example image of a microglia cell in the deep layer,
recorded with only the precorrection applied. Light scattering in the brain
tissue above leads to rather low contrast between structures in Fig. 5 (A), as
illustrated by the plot of signal intensity along the black dashed line, shown
in the inset. Figures 5 (B, C, D) show the same microglia cell after running
one out of three AO correction algorithms, respectively, each for 3 iterations
of 256 modes. Figure 5 (B) shows the cell after application of the correction
mask shown the upper left inset, which has been obtained by performing a-CSA
with power-overhead in the test mode ($f\sim 2/256$) on top of the
precorrection. The improvement in signal is on the order of a factor of 2 to 3
across the cell body, and processes extending from the cell into its
surrounding are starting to become visible, especially when a logarithmic
color map is applied (upper right inset). In Fig. 5 (C), the cell is shown
after performing regular CSA (i.e., without power overhead in the test mode,
$f=1/256$) on top of the precorrection. The measured modulation SNR was
insufficent for algorithm convergence; as is apparent, the signal quality did
not improve compared to the precorrection alone (Fig. 5 A), or even became
slightly worse. Figure 5 (D) shows the cell after performing DASH (3
iterations of 256 modes, $f=0.3$) on top of the precorrection. The signal
intensity across the cell body is increased by a factor of about 5; the
contrast between structures is clearly enhanced, allowing to distinguish
processes extending from the cell body into the surrounding tissue.

Note that for taking the images of Fig. 5 we have started from the
precorrection of Fig. 5 (A) simply for illustration purposes, to disentangle
the correction of (mild) optical system aberrations from the correction of
scattering inside the brain tissue, and for compatibility with the dye-slide
simulations and experiments. This, however, is not an experimental necessity;
typically, images of a similar quality as shown in Fig. 5 (B, D) can also be
obtained by performing DASH or a-CSA without any precorrection.

Figure 5: Image enhancement for a biological sample. Murine microglia
expressing green fluorescent protein, $200$ µm deep inside fixed brain tissue.
Image (A) is recorded using only a precorrection for optical system
aberrations. (B)–(D) show the same cell after precorrection plus a subsequent
correction (3 iterations of 256 modes) for scattering inside the tissue: (B) =
DASH, (C) = regular CSA, (D) = a-CSA. Insets show the TPEF along the white,
dashed line; the size of the scalebar corresponds to 10 µm.

###  IV.6 Low-SNR scenarios

The practicality of every sensorless wavefront correction scheme depends on
its robustness with respect to noise. It has already been shown that DASH
performs well in this regard in comparison to alternative methods [1]. Let us
now compare performances at low-SNR conditions. The key for all such methods
is being able to discern the TPEF signal modulation – caused by phase-stepping
of the test beam – on top of the noise floor, since this modulation carries
the relevant information. For this, the test mode must contain a significant
percentage of the total light power.

In the following, we use numerical simulations to compare which SNR each of
the three methods DASH,777For DASH we slightly increase the power in the test
beam compared to before ($f=0.3\rightarrow 0.35$), which has proven helpful in
low-SNR situations [1]. a-CSA, and regular CSA888CSA without controlling the
power in individual pixels, as in Ref. [2]. requires to operate successfully.
To this aim, we simulate shot noise in the PMT readout and then repeatedly
execute the methods starting from a decreasing initial signal level. We stop
repeating once the signal enhancement $\eta$ achieved by a method (after 3
full iterations) drops below a certain threshold $\eta_{\text{th}}$. We set
this threshold to $\eta_{\text{th}}=0.75\left(\eta_{\text{max}}-1\right)+1$,
where $\eta_{\text{max}}$ is the enhancement for the noise-free case (i.e.,
simulations as in Section IV.4). The exact choice of the threshold prefactor
(0.75 in our case) is largely arbitrary and uncritical for our conclusions. We
assume the same sample and investigate the same three turbidity scenarios A,
B, C as discussed above.

Our simulations show that for weak scattering (A) and correction of $N^{2}=16$
modes, DASH signal enhancement crosses threshold when the photon level (before
correction) exceeds around 100 photons. a-CSA delivers comparable performance
when the test-superpixel contains about 10 times more laser power than any of
the other superpixels, corresponding to around 40% of the total power. Regular
CSA only crosses threshold for signal levels from around 1000 detected photons
onwards, i.e., about 10 times higher than required for DASH.
For medium turbidity (B) and $N^{2}=64$ modes, DASH crosses threshold at about
100 photons collected per measurement (before correction). Comparable
performance can be achieved using a-CSA if the test-superpixel amplification
factor is around 50, again representing roughly 40% of the total laser power.
Regular CSA crosses threshold at around 3 k photons.
Finally, for high turbidity (C) and $N^{2}=256$ modes, DASH crosses threshold
at around 400 counted photons (before correction). Again, comparable
performance can be achieved using a-CSA with a test-superpixel amplification
factor around 50. Regular CSA, in contrast, requires about 15 k photons.

These results are summarized in Table 1. Of course, in practice the required
signal levels will also depend on the nature of the sample, wherefore these
numbers can only serve as a rough orientation. For example, increasing the
sample thickness will make it increasingly difficult to obtain a successful
correction.

Nevertheless, from our results we can draw three main conclusions. First, as
expected, regular CSA requires a much higher SNR than DASH to operate
successfully. Achieving a higher SNR requires sending more optical power into
the sample volume, making regular CSA unfavourable, e.g., for imaging of
fragile specimens. Second, it is possible to successfully operate a-CSA at the
same SNR conditions as DASH, if a sufficient amplification of the test-
superpixel can be provided. However, it needs to be stressed that if the
amplification is realized using a superpixel method (as in Section IV.3),
discarding light from other pixels, this high performance of a-CSA comes at
the price of wasting much optical power. For instance, if the hologram
features $N{\times}N$ superpixels and the test-superpixel is supposed to
contain a fraction $f$ of the total pupil intensity, the incident laser power
must be increased by a factor of $(N^{2}-1)f/(1-f)$ to keep the total optical
power in the objective pupil constant. Even for a moderate number of modes of
64, i.e., $N=8$, achieving $f=0.3$ requires a total power increase by a factor
of 27. This steep scaling with $N$ in practice quickly becomes unfavorable for
application of a-CSA in deep tissue imaging, where typically many correctable
modes are needed. Third, for DASH, in contrast, such problems do not arise,
since the power in the test mode is independent of $N$ and can be tuned
conveniently through the depth of the corresponding grating on the SLM.

Table 1: Photon counts (before correction) required at low SNR for comparable performance as in the noise-free case. Regular CSA requires much higher signal levels than DASH; a-CSA requires substantial power in the test-superpixel (see main text). Turbidity | low (A) | medium (B) | high (C)
---|---|---|---
DASH | 100 | 100 | 400
CSA | 1 k | 3 k | 15 k

##  V Discussion and Summary

In this work, we have devised a method to compute a 2D phase mask whose effect
on a light field is in many ways equivalent to that of a voluminous scatter
medium of corresponding transport mean free path. This enables software-
controlled systematic studies where the degree of scattering is controlled
precisely, tunable over a large range, and manual exchange of physical scatter
materials is not needed.

We have used this method to investigate the performance of DASH for different
scales of turbidity, ranging from the “traditional AO” regime to the regime of
multiply-scattered photons, and compared it to an alternative approach which
we call Amplified Continuous Sequential Algorithm (a-CSA). a-CSA operates on a
square-superpixel mode basis and is shown to work also in low-SNR situations
(unlike regular CSA) due to an increased relative power in the test mode. In
practice, this amplification of the test superpixel can be achieved by
attenuating all other superpixels, e.g., by adding phase gratings of defined
modulation amplitude. While easy to implement, this approach works only for a
limited number of correctable modes due to its inefficient use of light power,
and only for a-CSA superpixels consisting of a sufficient number of physical
SLM pixels.

We found in both numerical simulations and experiments that for low turbidity
DASH initially improves more slowly than a-CSA, pointing to effects of the
phase-only light shaping principle behind DASH. This initial-speed
disadvantage of DASH disappears for increasing turbidity. Furthermore, we
observed that independent of the degree of turbidity, from the end of the
second iteration onwards, DASH outperforms a-CSA. We attribute this
performance advantage to the fact that DASH can always exploit the full
resolution of the SLM, hence minimizes diffraction losses compared to
superpixel methods such as a-CSA.
As an overall tendency, we observed that the signal enhancement achievable
using DASH or a-CSA grows for increasing degrees of turbidity, highlighting
the potential benefits these algorithms promise for nonlinear imaging in
highly scattering environments.

We have illustrated the practical improvements these methods can yield for
two-photon microscopy by the example of GFP-microglia 200 µm deep inside mouse
brain tissue. Furthermore, we emphasize that a-CSA and DASH were performed
using identical hardware on the same experimental setup. Therefore, depending
on the task, it is possible to execute either (or a combination) of both
routines on a pure software level.

A critical challenge for wavefront correction in highly scattering samples
remains the size of the corrected field of view (“isoplanatic patch”, IP). For
fixed brain tissue, as shown in Fig. 5, IP diameters are on the order of 20 to
30 µm, whereas live brain tissue tends to scatter photons at larger angles,
decreasing the IP size to just a few µm [14]. Several strategies have been
proposed to increase IPs, including multi-conjugate AO [25, 26, 27] or the
application of individual corrections for many sample points [13, 14, 21].
Whether a-CSA and DASH may benefit from such strategies in terms of the IP
size is an interesting question for future studies.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

All authors contributed to conceiving the experiment. M.S. conducted the
experiments. M.S. and A.J. wrote the manuscript. A.J. performed the numerical
simulations. N.B. contributed to Section II and the Supplementary Materials.

## Funding

The authors acknowledge funding from the Austrian Science Fund (FWF) under the
grants no. P32146-N36, M 3060-NBL, I3984.

## Acknowledgements

We would like to thank Kai K. Kummer and Jeiny Luna Choconta, Institute for
Physiology (Innsbruck Med. Univ.), for preparing the biological samples.

## Data Availability Statement

The datasets are available on request.

## References

  * May  _et al._ [2021a] M. A. May, N. Barré, K. K. Kummer, M. Kress, M. Ritsch-Marte, and A. Jesacher, Fast holographic scattering compensation for deep tissue biological imaging, Nature Communications 12, 1 (2021a).
  * Vellekoop and Mosk [2008] I. M. Vellekoop and A. Mosk, Phase control algorithms for focusing light through turbid media, Optics communications 281, 3071 (2008).
  * Oheim  _et al._ [2001] M. Oheim, E. Beaurepaire, E. Chaigneau, J. Mertz, and S. Charpak, Two-photon microscopy in brain tissue: parameters influencing the imaging depth, Journal of neuroscience methods 111, 29 (2001).
  * Chaigneau  _et al._ [2011] E. Chaigneau, A. J. Wright, S. P. Poland, J. M. Girkin, and R. A. Silver, Impact of wavefront distortion and scattering on 2-photon microscopy in mammalian brain tissue, Opt. Express 19, 22755 (2011).
  * Jacques [2013] S. L. Jacques, Optical properties of biological tissues: a review, Physics in Medicine & Biology 58, R37 (2013).
  * Wang  _et al._ [2011] Z. Wang, H. Ding, and G. Popescu, Scattering-phase theorem, Optics letters 36, 1215 (2011).
  * Booth  _et al._ [2002] M. J. Booth, M. A. A. Neil, and T. Wilson, New modal wave-front sensor: application to adaptive confocal fluorescence microscopy and two-photon excitation fluorescence microscopy, J. Opt. Soc. Am. A 19, 2112 (2002).
  * Conkey  _et al._ [2012a] D. B. Conkey, A. M. Caravaca-Aguirre, and R. Piestun, High-speed scattering medium characterization with application to focusing light through turbid media, Optics express 20, 1733 (2012a).
  * Blochet  _et al._ [2017] B. Blochet, L. Bourdieu, and S. Gigan, Focusing light through dynamical samples using fast continuous wavefront optimization, Opt. Lett. 42, 4994 (2017).
  * Conkey  _et al._ [2012b] D. B. Conkey, A. N. Brown, A. M. Caravaca-Aguirre, and R. Piestun, Genetic algorithm optimization for focusing through turbid media in noisy environments, Optics express 20, 4840 (2012b).
  * Paudel  _et al._ [2015] H. P. Paudel, J. Taranto, J. Mertz, and T. Bifano, Axial range of conjugate adaptive optics in two-photon microscopy, Optics Express 23, 20849 (2015).
  * Tao  _et al._ [2017] X. Tao, T. Lam, B. Zhu, Q. Li, M. R. Reinig, and J. Kubby, Three-dimensional focusing through scattering media using conjugate adaptive optics with remote focusing (CAORF), Opt. Express 25, 10368 (2017).
  * Papadopoulos  _et al._ [2020] I. N. Papadopoulos, J.-S. Jouhanneau, N. Takahashi, D. Kaplan, M. Larkum, J. Poulet, and B. Judkewitz, Dynamic conjugate F-SHARP microscopy, Light: Science & Applications 9, 1 (2020).
  * May  _et al._ [2021b] M. A. May, K. K. Kummer, M.-L. Edenhofer, J. L. Choconta, M. Kress, M. Ritsch-Marte, and A. Jesacher, Simultaneous scattering compensation at multiple points in multi-photon microscopy, Biomedical Optics Express 12, 7377 (2021b).
  * Kubby [2013] J. A. Kubby,  _Adaptive optics for biological imaging_ (CRC press, 2013).
  * Booth [2006] M. J. Booth, Wave front sensor-less adaptive optics: a model-based approach using sphere packings, Opt. Express 14, 1339 (2006).
  * Vellekoop and Mosk [2007] I. M. Vellekoop and A. P. Mosk, Focusing coherent light through opaque strongly scattering media, Opt. Lett. 32, 2309 (2007).
  * Tang  _et al._ [2012] J. Tang, R. N. Germain, and M. Cui, Superpenetration optical microscopy by iterative multiphoton adaptive compensation technique, Proceedings of the National Academy of Sciences 109, 8434 (2012).
  * Papadopoulos  _et al._ [2017] I. N. Papadopoulos, J.-S. Jouhanneau, J. F. Poulet, and B. Judkewitz, Scattering compensation by focus scanning holographic aberration probing (F-SHARP), Nature Photonics 11, 116 (2017).
  * Rodríguez  _et al._ [2021] C. Rodríguez, A. Chen, J. A. Rivera, M. A. Mohr, Y. Liang, R. G. Natan, W. Sun, D. E. Milkie, T. G. Bifano, X. Chen, _et al._ , An adaptive optics module for deep tissue multiphoton imaging in vivo, Nature Methods 18, 1259 (2021).
  * Blochet  _et al._ [2021] B. Blochet, W. Akemann, S. Gigan, and L. Bourdieu, Fast wavefront shaping for two-photon brain imaging with large field of view correction, bioRxiv (2021).
  * Hu  _et al._ [2020] Q. Hu, J. Wang, J. Antonello, M. Hailstone, M. Wincott, R. Turcotte, D. Gala, and M. Booth, A universal framework for microscope sensorless adaptive optics: Generalized aberration representations, APL Photonics 5, 100801 (2020).
  * Katz  _et al._ [2014] O. Katz, P. Heidmann, M. Fink, and S. Gigan, Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations, Nature photonics 8, 784 (2014).
  * Bagnoud and Zuegel [2004] V. Bagnoud and J. D. Zuegel, Independent phase and amplitude control of a laser beam by use of a single-phase-only spatial light modulator, Optics letters 29, 295 (2004).
  * Kam  _et al._ [2007] Z. Kam, P. Kner, D. Agard, and J. W. Sedat, Modelling the application of adaptive optics to wide-field microscope live imaging, Journal of microscopy 226, 33 (2007).
  * Simmonds and Booth [2013] R. D. Simmonds and M. J. Booth, Modelling of multi-conjugate adaptive optics for spatially variant aberrations in microscopy, Journal of Optics 15, 094010 (2013).
  * Wu and Cui [2015] T.-W. Wu and M. Cui, Numerical study of multi-conjugate large area wavefront correction for deep tissue microscopy, Optics Express 23, 7463 (2015).

## Supplementary Material

###  .1 Algorithms for indirect wavefront sensing

In this Section, we provide the algorithms for DASH and a-CSA in pseudo-code.
Algorithm .1 describes general steps that both approaches have in common. The
function wavefront_sensing takes as inputs a set of test modes $M$ and a
correction pattern $C$, initialized by arrays of ones. The algorithm loops $I$
times through the entire set of $N$ modes. For each mode $M_{n}$, the phase
stepping procedure phase_step is executed, which returns an array of complex
numbers $U_{\text{SLM}}$ which carries the amplitude and phase pattern for a
particular phase step of the excitation beam. measure_signal measures a single
TPEF signal $S_{p}$ for the pattern $U_{\text{SLM}}$ displayed on the SLM
corresponding to a particular phase step $p=1,\ldots,P$. From the set of $P$
TPEF measurements per mode the complex-valued coefficient $a_{n}$ is
calculated, containing magnitude and phase estimates of the mode $M_{n}$.
Finally, the correction beam $C$ is updated by the function update using the
new information ($a_{n},M_{n}$), improving the wavefront correction.

DASH and a-CSA differ in how $U_{\text{SLM}}$ and the corrected beam part $C$
are updated. The respective sub-functions are outlined in algorithm .1 for
a-CSA and in algorithm .1 for DASH.

The algorithms contain a few variables that are not specifically declared as
inputs and can be viewed as “global”, for instance the number of iterations
($I$), modes ($N$) or phase steps ($P$). A particular variable that needs
explanation is the scalar value $f$, which ranges between 0 and 1 and
determines the power fraction contained in the test mode.

Compared to our earlier implementation of DASH [1], the Algorithm .1 outlined
below features a minor modification: the normalization of the corrected beam
($|C_{i,n}|$ in Eq. 1 of Ref. [1]) which erases the amplitude information from
the corrected field $C_{i,n}$ is replaced by a different scalar normalization
factor, which can lead to slightly better performance.

General algorithm

$\mathbf{Inputs{:}}$ $M$, a list of $N$ 2D real-valued input modes of size
$N_{\text{SLM}}\times N_{\text{SLM}}$; $C$, a 2D array of size
$N_{\text{SLM}}\times N_{\text{SLM}}$, initialized with ones.

procedure wavefront_sensing($M$, $C$)

for $i=1\,\mathbf{to}\,I$ do

for $n=1\,\mathbf{to}\,N$ do

for $p=1\,\mathbf{to}\,P$ do

$U_{SLM}\leftarrow\textsc{phase\\_step}(C,M_{n},p)$

$S_{p}\leftarrow\textsc{measure\\_signal}(U_{\text{SLM}})$

end for

$a_{n}\leftarrow\frac{1}{P}\sum^{P}_{p=1}\sqrt{S_{p}}~{}\exp(-\text{i}\frac{2\pi}{P}p)$

$C\leftarrow\textsc{update}(C,M_{n},a_{n})$

end for

end for

return $C$

end procedure

Specific functions for a-CSA

procedure phase_step($C$, $M_{n}$, $p$)

$U_{SLM}\leftarrow\sqrt{1-f}~{}(1-M_{n})~{}C+\sqrt{f}~{}M_{n}~{}\exp\left(\text{i}\frac{2\pi}{P}p\right)$

return $U_{\text{SLM}}$

end procedure

procedure update($C$, $M_{n}$, $a_{n}$)

$C\leftarrow
C~{}(1-M_{n})+M_{n}~{}\exp{\left(-\text{i}~{}M_{n}~{}\mathrm{angle}(a_{n})\right)}$

return $C$

end procedure

Specific functions for DASH

procedure phase_step($C$, $M_{n}$, $p$)

$U_{\text{SLM}}\leftarrow\sqrt{1-f}~{}C+\sqrt{f}\exp\left(\text{i}(M_{n}+\frac{2\pi}{P}p)\right)$

return $\exp{\left(\text{i}~{}\mathrm{angle}(U_{\text{SLM}})\right)}$

end procedure

procedure update($C$, $M_{n}$, $a_{n}$)

$C\leftarrow C+a_{n}^{*}~{}\exp{\left(\text{i}~{}M_{n}\right)}$

$\gamma\leftarrow\sqrt{\frac{1}{N_{\text{SLM}}^{2}}\sum|C|^{2}}$

$C\leftarrow\frac{C}{\gamma}$ $\triangleright$ Normalization

return $C$

end procedure

###  .2 Construction of scatter mask in the numerical simulation and
experiment

Algorithm .2 describes the algorithmic procedure of constructing the
scattering phase mask displayed at the SLM. The function make_scatterer takes
the following inputs: the desired spatial frequency content characterized by
$\sigma$, the side length of the phase mask $N_{\text{scat}}$, and the desired
RMS phase magnitude of the scatterer $a_{\text{scat}}$. The pixel indices are
given by $x$ and $y$. The function rand creates equally distributed random
numbers in the interval $[0,1[$. fft represents the fast Fourier transform and
norm normalizes a real-valued array to an RMS value of $a_{\text{scat}}$.

Scatterer calculation

procedure make_scatterer($\sigma$, $N_{\text{scat}}$, $a_{\text{scat}}$)

$W\leftarrow\exp{\left(-(x^{2}+y^{2})/(2\sigma^{2})\right)}$

$A\leftarrow\exp\left(\text{i}~{}2\pi~{}\textsc{rand}(N_{\text{scat}},N_{\text{scat}})\right)\cdot
W$

$A[1]\leftarrow 0$ $\triangleright$ discard DC term

$\Phi_{\text{scat}}\leftarrow\mathrm{real}\left(\textsc{fft}(A)\right)$

$\Phi_{\text{scat}}\leftarrow\textsc{norm}(\Phi_{\text{scat}},a_{\text{scat}})$

return $\Phi_{\text{scat}}$

end procedure

###  .3 Animations of DASH and CSA correction procedures

Figure 6: Exemplary frame from the GIF animation for scenario B. Focus quality
and phase correction pattern after optimization of the 10th mode in the first
optimization iteration.

We provide three GIF animations visualizing the progress of correction for
DASH and CSA999Note that for the noise-free case as discussed here
(SNR${}\rightarrow\infty$) the performance of regular and a-CSA is identical,
wherefore we use the more general label “CSA”. for the three scenarios A, B
and C as discussed in the main document. The animations show the evolution of
focus quality (left column) and corresponding correction phase masks over
three measurement iterations. The scalebar for the intensity images is
normalized to the maximum intensity of an aberration-free focus. An exemplary
frame of the animation for scenario B is shown in Fig. 6.