# Critical points of a perturbed Otha-Kawasaki functional

Matteo Rizzi

(Date)

###### Abstract

In the paper, we consider a small perturbation of the Otha-Kawasaki functional
and we construct at least four critical points close to suitable translations
of the Schwarz P surface with fixed volume.

###### Contents

  1. 1 Introduction
  2. 2 The proof of Theorem : Lyapunov-Schmidt reduction
    1. 2.1 The volume constraint
    2. 2.2 The auxiliary equation
    3. 2.3 The bifurcation equation
  3. 3 Solving the auxiliary equation
    1. 3.1 The linear problem
    2. 3.2 The proof of Proposition : a fixed point argument
  4. 4 Solving the bifurcation equation.
  5. 5 Appendix

##  1 Introduction

A diblock copolymer is a complex molecule where chains of two different kinds
of monomers, say A and B, are grafted togheter. Diblock copolymer melts are
large collections of diblock copolymers. The experiments show that, above a
certain temperature, these melts behave like fluids, that is the monomers are
mixed in a disordered way, while below this critical temperature phase
separation is observed. Some common periodic structures observed in
experiments are spheres, cylinders, gyroids and lamellae (see figure 1).

Figure 1: The most commonly observed periodic structures are spheres,
cylinders, gyroids and lamellae

These patterns can be found by minimizing some energy. It looks reasonable to
describe the phenomenon through an energy given by the sum of the perimeter,
that forces the separation surfaces to be minimal, plus some nonlocal term
that keeps trace of the long-range interactions between monomers. More
explicitly, one can take the functional

$\displaystyle\mathcal{E}(u):=\frac{1}{2}\int_{\Omega}|\nabla
u|dx+\gamma\int_{\Omega}\int_{\Omega}G(x,y)(u(x)-m)(u(y)-m)dxdy$ (1)

as an energy. Here $\Omega$ is a bounded domain of $\mathbb{R}^{3}$, that can
be seen as the container where the diblock copolimer melt is confined, $u$ is
a bounded variation function in $\Omega$ with values in $\\{\pm 1\\}$ (for
instance, we can assume that $u(x)=1$ if there are only monomers of type A at
$x$, $u(x)=-1$ if there are only monomers of type B at $x$),
$\int_{\Omega}|\nabla u|dx$ is its total variation, or equivalently the
perimeter of the set $\\{x\in\Omega:u(x)=1\\}$, $G$ is the Green’s function of
$-\Delta$ on $\Omega$, that is the disrtibutional solution to

$\displaystyle\begin{cases}-\Delta_{x}G(x,y)=\delta_{y}(x)-\frac{1}{|\Omega|}&\text{in
}\Omega\\\ \partial_{\nu(x)}G(x,y)=0&\text{on }\partial\Omega.\end{cases}$

$G$ turns out to be the sum of the Green’s function of $-\Delta$ over
$\mathbb{R}^{3}$ and a regular part $R(x,y)$, namely

$\displaystyle G(x,y)=\frac{c}{|x-y|}+R(x,y),$

(see [26]). $\gamma\geq 0$ is a parameter depending on the material, that we
will assume to be small.

This energy appears as the $\Gamma$-limit as $\varepsilon\to 0$ of the
approximating functionals

$\displaystyle\mathcal{E}_{\varepsilon}(u)=\frac{\varepsilon}{2}\int_{\Omega}|\nabla
u|^{2}dx+\frac{1}{\varepsilon}\int_{\Omega}\frac{(1-u^{2})^{2}}{4}dx$
$\displaystyle+\frac{16\gamma}{3}\int_{\Omega}\int_{\Omega}G(x,y)(u(x)-m)(u(y)-m)dxdy,$

introduced by Otha and Kawasaki (see [2, 6, 7, 8]).

In a more geometric way our functional is given by

$\displaystyle
J_{\gamma}(E):=P_{\Omega}(E)+\gamma\int_{\Omega}\int_{\Omega}G(x,y)(u_{E}(x)-m)(u_{E}(y)-m)dxdy$
(2)

where

$\displaystyle E:=\\{x\in\Omega:u(x)=1\\},$

so that $u_{E}=\chi_{E}-\chi_{\Omega\backslash E}$. The first variation of
$J_{\gamma}$ is given by

$\displaystyle
J_{\gamma}^{{}^{\prime}}(E)[\varphi]=\int_{\Sigma}(H_{\Sigma}(x)+4\gamma
v_{E}(x))\varphi(x)d\sigma(x),$ (3)

while its second variation is given by

$\displaystyle
J_{\gamma}^{{}^{\prime\prime}}(E)[\varphi]=\int_{\Sigma}L\varphi(x)\varphi(x)d\sigma(x),$
(4)

where

$\displaystyle
L\varphi=-\Delta_{\Sigma}\varphi-|A|^{2}\varphi+8\gamma\int_{\Sigma}G(\cdotp,y)\varphi(y)d\sigma(y)+4\gamma\partial_{\nu}v\varphi.$
(5)

Here $\varphi$ is in the space

$\displaystyle W:=\bigg{\\{}w\in
H^{1}(\Sigma):\int_{\Sigma}w(x)\nu_{i}(x)d\sigma(x)=0,$ $\displaystyle\text{
}1\leq i\leq 3\bigg{\\}},$ (6)

$\Sigma:=\partial E$ and

$\displaystyle v_{E}(x):=\int_{T^{3}}G(x,y)(u_{E}(y)-m)dy$ (7)

is the unique solution to the problem

$\displaystyle\begin{cases}-\Delta v_{E}=u_{E}-m&\text{in $T^{3}$}\\\
\int_{T^{3}}v_{E}dx=0.\end{cases}$ (8)

For an explicit computation of the first and the second variation, see for
instance [9]. In the sequel, $\Omega$ will always be the $3$-dimensional torus
$T^{3}$, that is the quotient of the cube $[0,1]^{3}$ by the equivalence
relation that identifies the opposite faces. It is known that $J_{\gamma}$ is
translation invariant, that is $J_{\gamma}(E+\xi)=J_{\gamma}(E)$, for any
$\xi\in T^{3}$ (see [2],[9]), thus, once we find a critical point of it, any
translation in $T^{3}$ is still critical.

There are several results in the literature about critical points of this
functional. For instance, an interesting problem is to understand whether all
global minimizers are periodic, like the patterns described above (spheres,
cylinders, gyroids and lamellae, see Figure 1). This is known to be true in
dimension one (see [22]), but the problem is still open in higher dimension.
We refer to [1, 31] for further results. Some other authors, such as Ren and
Wei [26, 27, 28, 29, 30], constructed explicit examples of stable periodic
local minimizers, that is with positive second variation. Moreover, Acerbi
Fusco and Morini [2] showed that any stable critical point is actually a local
minimizer with respect to small $L^{1}$ perturbations.

Here we add a small linear perturbation that corresponds to an external force
$f$ applied to the system, that can be taken to be
$C^{0,1}_{loc}(\mathbb{R}^{3})$ and periodic, with triple period $1$. The
energy becomes

$\displaystyle
I_{\gamma}(E):=J_{\gamma}(E)+\gamma\int_{\Omega}f(x)u_{E}(x)dx.$ (9)

The additional linear term breakes the translation invariance. We will
construct at least four critical points $F_{j}$ of $I_{\gamma}$, $1\leq j\leq
4$, for $\gamma$ small enough, that are close to suitable translations of the
Schwarz’ P surface $\Sigma$ (see figure 2), under the volume constraint

$\displaystyle\mathcal{L}_{3}(F_{j})=\mathcal{L}_{3}(E),$ (10)

where $E$ is the interior of $\Sigma$.

######  Remark 1.

The Schwartz P surface can be seen as a periodic surface in $\mathbb{R}^{3}$,
with triple period $1$. Moreover, it divides the Torus into two components, an
interior and an exterior. In the sequel, $E$ will denote the interior part.

Figure 2: Schwarz’ P surface

We will use a technique based on a finite dimensional Lyapunov-Schmidt
reduction (see [4], Chapter $2.2$), and on the Lusternik-Schnirelman theory
(see [3], Chapter $9$) for the multiplicity.

For $0<\alpha<1$ and for any integer $k\geq 0$, we introduce the Hölder spaces

$\displaystyle C^{k,\alpha}_{s}(\overline{\Sigma}):=\\{w\in
C^{k,\alpha}(\overline{\Sigma}):w(x)=w(T_{j}x),1\leq j\leq 3\\},$ (11)

where $T_{j}$ are the reflections defined by

$\displaystyle T_{1}(x_{1},x_{2},x_{3})=(-x_{1},x_{2},x_{3})$ $\displaystyle
T_{2}(x_{1},x_{2},x_{3})=(x_{1},-x_{2},x_{3})$ $\displaystyle
T_{3}(x_{1},x_{2},x_{3})=(x_{1},x_{2},-x_{3}).$

Here it is understood that we have put the origin in the centre of the cube
(see Figure 2), in such a way that these spaces consist of functions that
respect the simmetries of $\Sigma$, that is the simmetries with respect to the
coordinate planes $\\{x_{j}=0\\}$, $1\leq j\leq 3$. We endow these spaces with
the norm

$\displaystyle||w||_{C^{k,\alpha}(\Sigma)}=\sum_{j=0}^{k}||\nabla^{j}w||_{L^{\infty}(\Sigma)}+\sup_{x\neq
y}\sup_{|\beta|=k}\frac{|\partial_{\beta}w(x)-\partial_{\beta}w(y)|}{d(x,y)^{\alpha}},$
(12)

where $d$ is the geodesics distance on $\Sigma$.

######  Theorem 1.

Let $I_{\gamma}$ be defined as in (9) and $\nu(x)$ be the outward-pointing
unit normal to the Schwarz P surface $\Sigma$. Then there exists
$\gamma_{0}>0$ such that, for any $0<\gamma<\gamma_{0}$, there exist
$\xi_{j}\in T^{3}$, $1\leq j\leq 4$, and $w_{\gamma,j}\in
C^{2,\alpha}_{s}(\overline{\Sigma})$, with

$\displaystyle||w_{\gamma,j}||_{C^{2,\alpha}(\Sigma)}\leq c\gamma,$ (13)

such that the sets $F_{j}$ defined as the interior of

$\displaystyle\Gamma_{j}:=\\{x+\xi_{j}+\nu(x)w_{\gamma,j}(x):x\in\Sigma\\}$
(14)

are critical points of $I_{\gamma}$ under the volume constraint

$\displaystyle\mathcal{L}_{3}(F_{j})=\mathcal{L}_{3}(E).$ (15)

######  Remark 2.

(i) If we take $f\equiv 0$, we find a unique critical point $F$, that is the
interior of

$\displaystyle\Gamma:=\\{x+\nu(x)w_{\gamma}(x):x\in\partial E\\},$ (16)

where $w_{\gamma}$ is a small correction, namely
$||w_{\gamma}||_{C^{2,\alpha}(\Sigma)}\leq c\gamma$, found by means of the
implicit function Theorem (see Remark 4). Then any translation $F+\xi$ is
still a critical point of $J_{\gamma}$. A similar result was proved by
Cristoferi (see [11], Theorem $4.18$), who constructed a critical point of
$J_{\gamma}$ close to any smooth periodic strictly stable constant mean
curvature surface.

(ii) We stated the theorem in the case of $I_{\gamma}$ for simplicity. The
same proof should yield existence and multiplicity results also for regular
nonlinear perturbations and different coefficients in the nonlocal and forcing
terms.

A similar result was obtained by Bonacini and Cristoferi [5], who studied a
nonlocal version of the isoperimetric problem, that is they considered a small
nonlocal perturbation of the perimeter and showed that the unique minimizers
$F$ under the volume constraint $\mathcal{L}_{N}(F)=m$ are the balls, provided
$m$ is small enough. The critical points we construct here are not necessarily
stable, since we apply the Lusternik-Schnirelmann theory (see [3], chapter
$9$).

A crucial tool in the proof is nondegeneracy up to translations of the Jacobi
operator of the Schwarz P surface. In [25], Ross showed that the Schwarz P
surface is a critical point of the area and it is volume preserving stable,
that is it the second variation of the area is non-negative on any normal
variation with zero average. More precisely, setting $I_{0}:=P_{\Omega}$, we
have

$\displaystyle
I_{0}^{{}^{\prime\prime}}(E)(\varphi,\varphi)=\int_{\Sigma}|\nabla_{\Sigma}\varphi|^{2}-|A|^{2}\varphi^{2}d\sigma\geq
0$ (17)

for any $\varphi\in H^{1}(\Sigma)$ satisfying

$\displaystyle\int_{\Sigma}\varphi d\sigma=0,$ (18)

(see Theorem $1$ of [25]). Let $\nu(x)$ denote the exterior unit normal to
$\Sigma$ at $x$. Since $I_{0}$ is translation invariant, then
$\nu_{i}(x):=(\nu(x),e_{i})$ are Jacobi fields of $\Sigma$, that is they
satisfy

$\displaystyle-\Delta_{\Sigma}\nu_{i}-|A|^{2}\nu_{i}=0$ $\displaystyle\text{in
}\Sigma,$ (19)

(see [2],[9]). Moreover, Grosse-Brauckmann and Wohlgemuth showed in ([18])
that $\Sigma$ is nondegenerate up to translations, that is there are no other
nontrivial Jabobi fields. In other words

$\displaystyle Ker(I_{0}^{{}^{\prime\prime}}(E))=span\\{\nu_{i}\\}_{1\leq
i\leq 3}.$ (20)

######  Remark 3.

Let us observe that the $\nu_{i}$’s are linearly independent. In fact, if not,
there would exist a constant vector $b=(b_{1},b_{2},b_{3})\neq 0$ such that
$0=(b,\nu(x))$ for any $x\in\Sigma$, but this contradicts the geometry of
$\Sigma$.

We note that the $\nu_{i}$’s have zero average, since

$\displaystyle\int_{\Sigma}\nu_{i}(x)d\sigma(x)=\int_{T^{3}}\text{div}e_{i}=0.$
(21)

In addition, we decompose $H^{1}(\Sigma)$ into the orthogonal sum

$\displaystyle H^{1}(\Sigma)=span\\{\nu_{i}\\}_{1\leq i\leq 3}+W,$ (22)

(see (6) for the definition of $W$), and we define

$\displaystyle W^{0}:=\bigg{\\{}w\in
W:\int_{\Sigma}w(x)d\sigma(x)=0\bigg{\\}}.$ (23)

The above discussion can be rephrased by saying that

$\displaystyle\int_{\Sigma}|\nabla_{\Sigma}w|^{2}-|A|^{2}w^{2}d\sigma\geq
c||w||_{H^{1}(\Sigma)}^{2}$ $\displaystyle\text{for any }w\in W^{0}.$ (24)

Aknowledgments The author is supported by the PRIN project Variational and
perturbative aspects of nonlinear differential problems. The author is also
particularly grateful to F. Mahmoudi for his precious collaboration.

##  2 The proof of Theorem 1: Lyapunov-Schmidt reduction

We need to find at least four sets $F$ of the form (14) and a Lagrange
multiplier $\lambda\in\mathbb{R}$ such that

$\displaystyle H_{\partial F}(y)+4\gamma v_{F}(y)+\gamma f(y)=\lambda$
$\displaystyle\text{ }\forall y\in\partial F,$ (25)

or equivalently

$\displaystyle I^{{}^{\prime}}_{\gamma}(F)=\lambda.$ (26)

Exploiting the variational nature of the problem and the fact that
$H_{\Sigma}=0$, equation (25) is equivalent to

$\displaystyle\lambda=4\gamma v_{E}(x)+Lw(x)+Q(w)(x)+\gamma f(y),$
$\displaystyle\forall x\in\Sigma,$ (27)

where $y$ is seen as a function of $x$ depending on the parameter $\xi$,
namely $y=x+\xi+w(x)\nu(x)$, and

$\displaystyle
Q(w):=J^{{}^{\prime}}_{\gamma}(F)-J^{{}^{\prime}}_{\gamma}(E)-J^{{}^{\prime\prime}}_{\gamma}(E)w.$
(28)

Writing

$\displaystyle Lw=-\Delta_{\Sigma}w-|A|^{2}w+\gamma\tilde{L}w,$ (29)

where

$\displaystyle\tilde{L}w=8\int_{\Sigma}G(\cdotp,\zeta)w(\zeta)d\sigma(\zeta)+4\partial_{\nu}v_{E}w,$
(30)

we can see that (27) is equivalent to

$\displaystyle-\Delta_{\Sigma}w-|A|^{2}w=\lambda+\mathcal{F}(\gamma,\xi,w),$
(31)

where the nonlinear functional $\mathcal{F}$ is given by

$\displaystyle\mathcal{F}(\gamma,\xi,w)(x)=-4\gamma
v_{E}(x)-\gamma\tilde{L}w(x)-Q(w)(x)-\gamma f(y),$ $\displaystyle\forall
x\in\Sigma.$ (32)

The unknowns are the function $w$, $\xi\in T^{3}$ and $\lambda\in\mathbb{R}$.

###  2.1 The volume constraint

Now we will consider the relation between the volume of $F$ and $w$. In order
to do so, we point out that there exists a global parametrization

$\displaystyle\phi:Y\to\Sigma,$ (33)

defined on an open set $Y\in\mathbb{R}^{2}$ (see [14], section $3$), that
induces a change of coordinates on a neighbourhood of $\Sigma$ given by

$\displaystyle
X(\text{y}_{1},\text{y}_{2},z):=\phi(\text{y}_{1},\text{y}_{2})+z\nu(\text{y}_{1},\text{y}_{2}),$
(34)

where, with an abuse of notation, $\nu(\text{y}_{1},\text{y}_{2})$ is the
outward-pointing unit normal to $\Sigma$ at $\phi(\text{y}_{1},\text{y}_{2})$.
The volume of $F$ is given by

$\displaystyle\mathcal{L}_{3}(F)=\mathcal{L}_{3}(E)+\int_{Y}d\text{y}\int_{0}^{w(\text{y})}\det
JX(\text{y},z)dz,$

where $JX$ is the Jacobian of $X$. We expand

$\displaystyle\det JX(\text{y},z)=\det
JX(\text{y},0)+zA(\text{y})+z^{2}B(\text{y}),$

thus we get

$\displaystyle\mathcal{L}_{3}(F)=\mathcal{L}_{3}(E)+\int_{Y}d\text{y}\int_{0}^{w(\text{y})}\bigg{(}\det
JX(\text{y},0)+zA(\text{y})+z^{2}B(\text{y})\bigg{)}dz$
$\displaystyle=\mathcal{L}_{3}(E)+\int_{Y}\det
JX(\text{y},0)w(\text{y})d\text{y}+\int_{Y}\bigg{(}\frac{1}{2}w^{2}(\text{y})A(\text{y})+\frac{1}{3}w^{3}(\text{y})B(\text{y})\bigg{)}d\text{y}.$

Since $\det
JX(\text{y},0)=(\nu(\text{y}),\partial_{\text{y}_{1}}\phi\times\partial_{\text{y}_{2}}\phi)\neq
0$ for any y$\in Y$,

$\displaystyle\mathcal{L}_{3}(F)=\mathcal{L}_{3}(E)+\int_{\Sigma}w(x)d\sigma(x)+\int_{\Sigma}\tilde{Q}(x,w(x))d\sigma(x),$
(35)

where

$\displaystyle\tilde{Q}(x,w)=\frac{1}{\det
JX(x)}\bigg{(}\frac{1}{2}w^{2}(x)A(x)+\frac{1}{3}w^{3}(x)B(x)\bigg{)}.$ (36)

Therefore the volume constraint is equivalent to an equation of the form

$\displaystyle\int_{\Sigma}w(x)dx=-\int_{\Sigma}\tilde{Q}(x,w(x))d\sigma(x).$
(37)

###  2.2 The auxiliary equation

The aim is to solve (31) under the volume constraint (37). However, since, by
(20) and (24), the Jacobi operator $-\Delta_{\Sigma}-|A|^{2}$ is non
degenerate up to translations, we can actually solve the system

$\displaystyle\begin{cases}-\Delta_{\Sigma}w-|A|^{2}w=\lambda+P\mathcal{F}(\gamma,\xi,w)&\text{in
}\Sigma\\\ \partial_{n}w=0&\text{on }\partial\Sigma,\\\
\int_{\Sigma}wd\sigma=-\int_{\Sigma}\tilde{Q}(x,w(x))d\sigma(x),\end{cases}$
(38)

where $P:L^{2}(\Sigma)\to\tilde{W}$ is the projection onto the space

$\displaystyle\tilde{W}:=\bigg{\\{}\varphi\in
L^{2}(\Sigma):\int_{\Sigma}\varphi(x)\nu_{i}(x)d\sigma(x)=0,$ $\displaystyle
1\leq i\leq 3\bigg{\\}},$ (39)

$\partial_{n}w:=(\nabla_{\Sigma}w,n)$ and $n$ is the outward pointing unit
normal to $\partial\Sigma$ in $\Sigma$. This will be done by a fixed point
argument in the following Proposition, proved in section $3$.

######  Proposition 1.

For any $\xi\in T^{3}$ and for any $\gamma$ sufficiently small, there exists a
unique solution $(w_{\gamma,\xi},\lambda_{\gamma,\xi})\in
C^{2,\alpha}_{s}(\overline{\Sigma})\times\mathbb{R}$ to problem (38)
satisfying

$\displaystyle||w_{\gamma,\xi}||_{C^{2,\alpha}(\Sigma)}+|\lambda_{\gamma,\xi}|\leq
C\gamma,$ (40) $\displaystyle\int_{\Sigma}w(x)\nu_{i}(x)d\sigma(x)=0,$
$\displaystyle\text{ }1\leq i\leq 3,$ (41)

for some constant $C>0$. Moreover, the solution is Lipschitz continuous with
respect to the parameter $\xi$, that is

$\displaystyle||w_{\gamma,\xi_{1}}-w_{\gamma,\xi_{2}}||_{C^{2,\alpha}(\Sigma)}+|\lambda_{\gamma,\xi_{1}}-\lambda_{\gamma,\xi_{1}}|\leq
C\gamma|\xi_{1}-\xi_{2}|,$ $\displaystyle\forall\xi_{1},\xi_{2}\in T^{3}.$
(42)

######  Remark 4.

If we take $f\equiv 0$, in order to get the right correction $w$, we just
solve (38) for $\xi=0$, due to the translation invariance of $J_{\gamma}$ (see
Remark 2). We do not need the Lyapunov-Schmidt reduction.

###  2.3 The bifurcation equation

In order to conclude the proof of Theorem 1, we have to find at least four
points $\xi\in T^{3}$ such that
$(Id-P)\mathcal{F}(\gamma,\xi,w_{\gamma,\xi})(x)=0$, or equivalently

$\displaystyle\int_{\Sigma}\mathcal{F}(\gamma,\xi,w_{\gamma,\xi})(x)\nu_{i}(x)d\sigma(x)=0,$
(43)

for $i=1,2,3$.

Since $\partial_{n}w_{\gamma,\xi}=0$ on $\partial\Sigma$ and the same is true
for the $\nu_{i}$’s, an integration by parts yields

$\displaystyle\int_{\Sigma}(-\Delta_{\Sigma}w_{\gamma,\xi}-|A|^{2}w_{\gamma,\xi})(x)\nu_{i}(x)d\sigma(x)=0,$

for $i=1,2,3$, thus by (38) we can see that $w$ solves

$\displaystyle P(I_{\gamma}^{{}^{\prime}}(F)-\lambda)=0,$ (44)

or equivalently

$\displaystyle
I_{\gamma}^{{}^{\prime}}(F)-\lambda=\sum_{i=1}^{3}A_{i,\gamma,\xi}\nu_{i}.$
(45)

Since, by construction,

$\displaystyle\mathcal{F}(\gamma,\xi,w_{\gamma,\xi})=-I_{\gamma}^{{}^{\prime}}(F)-\Delta_{\Sigma}w_{\gamma,\xi}-|A|^{2}w_{\gamma,\xi}=$
$\displaystyle-
I_{\gamma}^{{}^{\prime}}(F)+\lambda+P\mathcal{F}(\gamma,\xi,w)=-\sum_{i=1}^{3}A_{i,\gamma,\xi}\nu_{i}+P\mathcal{F}(\gamma,\xi,w),$

and (45) holds, we can see that (43) is equivalent to

$\displaystyle A_{i,\gamma,\xi}=0$ $\displaystyle\text{for }i=1,2,3.$ (46)

Equation (46) is solvable thanks to the Lusternik-Schnirelmann theory and the
compactness of the Torus. We recall that the Torus $T^{3}$ has category $4$
(see [3], example $9.4$, (iii)).

######  Proposition 2.

Equation (46) is satisfied if $\xi$ is a critical point of the function
$\Phi_{\gamma}:T^{3}\to\mathbb{R}$ defined by

$\displaystyle\Phi_{\gamma}(\xi):=I_{\gamma}(F),$ (47)

where $F$ is the interior of

$\displaystyle\Gamma:=\\{x+\xi+w_{\gamma,\xi}(x)\nu(x):x\in\Sigma\\}.$

The proof of Proposition 2 will be carried out in Section $4$. It is possible
to see that $\Phi_{\gamma}$ actually admits at least $4$ critical points, due
to Theorem $9.10$ of [3] applied to $I_{\gamma}$, with $M=T^{3}$. The
compactness of the torus $T^{3}$ is crucial, since it guarantees that
$I_{\gamma}$ is bounded from below on $M$ and the Palais-Smale condition is
satisfied.

##  3 Solving the auxiliary equation

The aim of this section is to prove Proposition 1. First, in Section $3,1$, we
will treat the corresponding linear problem, then, in Section $3,2$, we will
solve problem (38) by a fixed point argument.

###  3.1 The linear problem

######  Proposition 3.

Let $a\in\mathbb{R}$ and $\varphi\in C^{0,\alpha}_{s}(\overline{\Sigma})$ be
such that

$\displaystyle\int_{\Sigma}\varphi\nu_{i}d\sigma=0$ $\displaystyle\text{for
}i=1,2,3.$ (48)

Then there exists a unique solution $(w,\lambda)=\Psi(\varphi,a)\in
C^{2,\alpha}_{s}(\overline{\Sigma})\times\mathbb{R}$ to the problem

$\displaystyle\begin{cases}-\Delta_{\Sigma}w-|A|^{2}w=\lambda+\varphi&\text{in
}\Sigma\\\ \partial_{n}w=0&\text{on }\partial\Sigma\\\
\int_{\Sigma}w\nu_{i}d\sigma=0&\text{for }1\leq i\leq 3,\\\
\int_{\Sigma}wd\sigma=a.\end{cases}$ (49)

Moreover, we have the stimate

$\displaystyle||w||_{C^{2,\alpha}(\Sigma)}+|\lambda|\leq
c(||\varphi||_{C^{0,\alpha}(\Sigma)}+|a|).$ (50)

######  Remark 5.

Since the $\nu_{i}$’s are linearly independent (see Remark 3), then the matrix

$\displaystyle L_{ki}:=\int_{\Sigma}\nu_{k}\nu_{i}d\sigma$ (51)

is invertible (for a detailed proof, see the appendix).

###### Proof.

Step (i): existence and uniqueness.
First we look for a weak solution $w\in W$. We write any $w\in W$ as

$\displaystyle w=w_{0}+\frac{1}{|\Sigma|}\int_{\Sigma}wd\sigma,$

with $w_{0}\in\ W^{0}$. The linear problem can be rephrased as follows

$\displaystyle\begin{cases}-\Delta_{\Sigma}w_{0}-|A|^{2}w_{0}=\lambda+\varphi+|A|^{2}\frac{a}{|\Sigma|}&\text{in
}\Sigma\\\ \int_{\Sigma}w_{0}=0.\end{cases}$ (52)

We note that the right-hand side of (52) is orthogonal to $\nu_{i}$, for
$i=1,2,3$, due to the fact that

$\displaystyle\int_{\Sigma}|A|^{2}\frac{a}{|\Sigma|}\nu_{i}(x)d\sigma=\int_{\Sigma}\bigg{(}\Delta_{\Sigma}\nu_{i}+|A|^{2}\nu_{i}\bigg{)}(x)\frac{a}{|\Sigma|}d\sigma=0,$
(53)

since $\partial_{n}\nu_{i}=0$ on $\partial\Sigma$, and

$\displaystyle\int_{\Sigma}\nu_{i}(x)d\sigma(x)=\int_{E}\text{div}(e_{1})dx=0.$
(54)

In addition, the norm defined by

$\displaystyle||w||=\int_{\Sigma}|\nabla_{\Sigma}w|^{2}-|A|^{2}w^{2}$ (55)

is equivalent to the $H^{1}(\Sigma)$-norm on $W^{0}$, thus the functional

$\displaystyle
I(w)=\int_{\Sigma}|\nabla_{\Sigma}w|^{2}-|A|^{2}w^{2}d\sigma-\int_{\Sigma}\bigg{(}\varphi+|A|^{2}\frac{a}{|\Sigma|}\bigg{)}wd\sigma$

is bounded from below by

$\displaystyle I(w)\geq
c||w||_{H^{1}(\Sigma)}^{2}-||\varphi||_{L^{2}(\Sigma)}||w||_{H^{1}(\Sigma)},$
(56)

on $W^{0}$, hence it is coercive on it. Moreover, this functional is also
w.l.s.c. and strictly convex on $W^{0}$, therefore any minimizing sequence
$w_{k}\in W^{0}$ weakly converges, up to subsequence, to the unique minimizer
$w_{0}\in W^{0}$, which satisfies the Euler-Lagrange equation

$\displaystyle\int_{\partial\Sigma}\partial_{n}wvd\sigma_{\partial\Sigma}+\int_{\Sigma}(-\Delta_{\Sigma}w_{0}-|A|^{2}w_{0})vd\sigma=$
$\displaystyle\lambda\int_{\Sigma}vd\sigma+\sum_{i=1}^{3}\beta_{i}\int_{\Sigma}\nu_{i}vd\sigma+\int_{\Sigma}\varphi
vd\sigma+\int_{\Sigma}|A|^{2}\frac{a}{|\Sigma|}vd\sigma,$ (57)

for any $v\in H^{1}(\Sigma)$, for some Lagrange multipliers
$\lambda,\beta_{i}\in\mathbb{R}$. Since $\varphi\in
C^{0,\alpha}(\overline{\Sigma})$, then $w\in C^{2,\alpha}(\overline{\Sigma})$
(see for instance [24]). Taking the test functions $v\in C^{1}_{c}(\Sigma)$,
we can see that $w$ satsfies

$\displaystyle-\Delta_{\Sigma}w_{0}-|A|^{2}w_{0}=\lambda+\sum_{i=1}^{3}\beta_{i}\nu_{i}+\varphi+|A|^{2}\frac{a}{|\Sigma|}$
$\displaystyle\text{in }\Sigma,$

in the classical sense. Taking now $v\in C^{1}(\overline{\Sigma})$, we can see
that the Neumann boundary condition is satisfied in the classical sense too.
Moreover, $w$ respects the required simmetries because of the symmetries of
the laplacian and uniqueness. Taking $\nu_{j}$ as a test function in (57),
using (54), (48), (53) the Neumann boundary condition and the fact that
$\partial_{n}\nu_{i}=0$ on $\partial\Sigma$, we get

$\displaystyle\sum_{i=1}^{3}\beta_{i}\int_{\Sigma}\nu_{i}\nu_{j}d\sigma=0,$

therefore by Remark 5, $\beta_{i}=0$.

Step (ii): Regularity estimates.
Multiplying (52) by $w_{0}$, integrating by parts and using (24), the Neumann
boundary conditions and Hölder’s inequality, we can see that

$\displaystyle
c||w_{0}||_{H^{1}(\Sigma)}^{2}\leq\int_{\Sigma}|\nabla_{\Sigma}w_{0}|^{2}-|A|^{2}w_{0}^{2}d\sigma=\int_{\Sigma}\varphi
w_{0}d\sigma+\frac{a}{|\Sigma|}\int_{\Sigma}|A|^{2}w_{0}d\sigma\leq$
$\displaystyle||w_{0}||_{L^{2}(\Sigma)}(||\varphi||_{L^{2}(\Sigma)}+\tilde{c}|a|)\leq||w_{0}||_{H^{1}(\Sigma)}(||\varphi||_{L^{2}(\Sigma)}+\tilde{c}|a|).$

Since $||w||_{H^{1}(\Sigma)}^{2}=||w_{0}||_{H^{1}(\Sigma)}^{2}+a^{2}$, then

$\displaystyle||w||_{H^{1}(\Sigma)}\leq c(||\varphi||_{L^{2}(\Sigma)}+|a|).$

In order to estimate $\lambda$, we integrate (49) and we get

$\displaystyle\lambda|\Sigma|+\int_{\Sigma}\varphi
d\sigma=-\int_{\Sigma}|A|^{2}wd\sigma,$

since, by the Neumann boundary conditions,

$\displaystyle\int_{\Sigma}\Delta_{\Sigma}wd\sigma=\int_{\partial\Sigma}\partial_{n}wd\sigma_{\partial\Sigma}=0,$
(58)

thus

$\displaystyle|\lambda|\leq
c(||\varphi||_{L^{2}(\Sigma)}+||w||_{L^{2}(\Sigma)}).$

To sum up, we have the estimate

$\displaystyle|\lambda|+||w||_{H^{1}(\Sigma)}\leq
c(||\varphi||_{L^{2}(\Sigma)}+|a|),$ (59)

In order to get the estimate with respect to the norms we are interested in,
we point out that, by the Sobolev embeddings

$\displaystyle||w||_{L^{\infty}(B_{\delta}(x))}\leq
c||w||_{W^{2,2}(B_{\delta}(x))}\leq
c(||w||_{L^{2}(B_{2\delta}(x))}+||\varphi+\lambda||_{L^{2}(B_{2\delta}(x))}+|a|)\leq$
$\displaystyle c(||\varphi||_{L^{2}(\Sigma)}+|a|),$

for any $\delta>0$ small but fixed and $x\in\Sigma$ such that
$d(x,\partial\Sigma)>\delta$ (here, $B_{\delta}(x)$ is the geodesic ball of
radius $\delta$ centered at $x$ in $\Sigma$). In particular,

$\displaystyle||w||_{L^{\infty}(\Sigma)}\leq
c(||\varphi||_{L^{\infty}(\Sigma)}+|a|).$

By the Hölder’s regularity estimates, we conclude that,

$\displaystyle||w||_{C^{2,\alpha}(\Sigma)}\leq
c(||w||_{L^{\infty}(\Sigma)}+||\varphi+\lambda||_{C^{0,\alpha}(\Sigma)})\leq
c(||\varphi||_{C^{0,\alpha}(\Sigma)}+|a|),$

(see [15], Chapter $6$, Theorem $6.30$). Since the same is true for
$|\lambda|$, the proof is over. ∎

###  3.2 The proof of Proposition 1: a fixed point argument

Now we are ready to show existence, uniqueness and Lipschitz continuity with
respect to $\xi$ of the solution $(w,\lambda)$ to (38).

Step (i): Existence and uniqueness.

We solve our problem by a fixed point argument. In fact the map

$\displaystyle
T(w,\lambda)=\Psi\bigg{(}P\mathcal{F}(\gamma,\xi,w),-\int_{\Sigma}\tilde{Q}(w)\bigg{)}$

is a contraction on the product $B\times\Lambda$, where
$\Lambda=(-C\gamma,C\gamma)$ and

$\displaystyle B:=\\{w\in W\cap
C^{2,\alpha}_{s}(\overline{\Sigma}):\partial_{n}w=0\text{ on $\partial\Sigma$,
}||w||_{C^{2,\alpha}(\Sigma)}<C\gamma\\},$ (60)

provided $C$ is large enough. In fact

$\displaystyle||\mathcal{F}(\gamma,\xi,w)||_{C^{0,\alpha}(\Sigma)}\leq\gamma(4||v_{E}||_{C^{2,\alpha}(\Sigma)}+||f||_{C^{0,\alpha}(\Sigma)})+c\gamma||w||_{C^{2,\alpha}(\Sigma)}\leq$
$\displaystyle\gamma(4||v_{E}||_{C^{2,\alpha}(\Sigma)}+||f||_{C^{0,\alpha}(\Sigma)})+cC\gamma^{2}<C\gamma$

provided $C>2(4||v_{E}||_{C^{2,\alpha}(\Sigma)}+||f||_{C^{0,\alpha}(\Sigma)})$
and $\gamma$ is small enough. Similarly, we can see that
$\mathcal{F}(\gamma,\xi,w)$ is Lipschitz continuous in $w$ with Lipschitz
constant of order $\gamma$.

In addition, the second component fulfills

$\displaystyle\bigg{|}\int_{\Sigma}\tilde{Q}(w)\bigg{|}\leq
c||w||_{C^{2,\alpha}(\Sigma)}^{2}\leq cC^{2}\gamma^{2}<C\gamma$

if $\gamma$ is small enough, and the same is true for the Lipschitz constant.

Lipschitz continuity with respect to $\xi$.

In order to prove (42), we point out that, if we set $w_{i}:=w_{\xi_{i}}$ and
$y_{i}:=x+\xi_{i}+w_{i}(x)\nu(x)$, for $i=1,2$,

$\displaystyle||f(y_{1})-f(y_{2})||_{C^{0,\alpha}(\Sigma)}\leq
c(|\xi_{1}-\xi_{2}|+||w_{1}-w_{2}||_{C^{2,\alpha}(\Sigma)})$
$\displaystyle||\tilde{L}w_{1}-\tilde{L}w_{2}||_{C^{0,\alpha}(\Sigma)}\leq
c||w_{1}-w_{2}||_{C^{2,\alpha}(\Sigma)}$

and

$\displaystyle||Q(w_{1})-Q(w_{2})||_{C^{0,\alpha}(\Sigma)}\leq
c(||w_{1}||_{C^{2,\alpha}(\Sigma)}+||w_{2}||_{C^{2,\alpha}(\Sigma)})||w_{1}-w_{2}||_{C^{2,\alpha}(\Sigma)}\leq$
$\displaystyle cC\gamma||w_{1}-w_{2}||_{C^{2,\alpha}(\Sigma)}.$

Similarly, we can show that

$\displaystyle\bigg{|}\int_{\Sigma}(\tilde{Q}(w_{1})-\tilde{Q}(w_{2}))d\sigma\bigg{|}\leq
c\gamma||w_{1}-w_{2}||_{C^{2,\alpha}(\Sigma)},$

thus, applying $\Psi$,

$\displaystyle|\lambda_{1}-\lambda_{2}|+||w_{1}-w_{2}||_{C^{2,\alpha}(\Sigma)}\leq
c\gamma(||w_{1}-w_{2}||_{C^{2,\alpha}(\Sigma)}+|\xi_{1}-\xi_{2}|).$

In conclusion, for $\gamma$ small enough,

$\displaystyle|\lambda_{1}-\lambda_{2}|+\frac{1}{2}||w_{1}-w_{2}||_{C^{2,\alpha}(\Sigma)}\leq
c\gamma|\xi_{1}-\xi_{2}|.$

##  4 Solving the bifurcation equation.

The parametrization $\phi:Y\to\Sigma$ of $\Sigma$ introduced in (33) induces a
parmetrization $\beta:Y\to\Gamma:=\partial F$ given by

$\displaystyle\beta(\text{y}_{1},\text{y}_{2}):=\phi(\text{y}_{1},\text{y}_{2})+\xi+w_{\gamma,\xi}(\text{y}_{1},\text{y}_{2})\nu(\text{y}_{1},\text{y}_{2}).$
(61)

The volume element can be expressed in terms of $\phi$ in this way

$\displaystyle|\beta_{\text{y}_{1}}\times\beta_{\text{y}_{2}}|=|\phi_{\text{y}_{1}}\times\phi_{\text{y}_{2}}|+L^{1}_{\xi}w_{\xi}+Q^{1}_{\xi}w_{\xi},$

where $L^{1}_{\xi}$ depends linearly on $w_{\xi}$ and on its gradient and
$Q^{1}_{\xi}$ is quadratic in the same quantites. More precisely, they satisfy
the estimates

$\displaystyle\begin{cases}|L^{1}_{\xi}w|\leq c||w||_{C^{2,\alpha}(\Sigma)}\\\
|Q^{1}_{\xi}(w)|\leq c||w||_{C^{2,\alpha}(\Sigma)}^{2}.\end{cases}$ (62)

Using the Taylor expansion of the function $\frac{1}{1+s}$, we can show that
the outward-pointing unit normal to $\Gamma$ is

$\displaystyle\nu_{\Gamma}=\frac{\beta_{\text{y}_{1}}\times\beta_{\text{y}_{2}}}{|\beta_{\text{y}_{1}}\times\beta_{\text{y}_{2}}|}=\frac{\phi_{\text{y}_{1}}\times\phi_{\text{y}_{2}}}{|\phi_{\text{y}_{1}}\times\phi_{\text{y}_{2}}|}+\tilde{L}^{1}_{\xi}w_{\xi}+\tilde{Q}^{1}_{\xi}w_{\xi}=$
(63) $\displaystyle\nu+\tilde{L}^{1}_{\xi}w_{\xi}+\tilde{Q}^{1}_{\xi}w_{\xi},$

with $\tilde{L}^{1}_{\xi}$ and $\tilde{Q}^{1}_{\xi}$ satisfying (62).

Now we point out that, if $\xi$ is a critical point of $\Phi_{\gamma}$, then

$\displaystyle\partial_{\xi_{i}}\Phi_{\gamma}(\xi)=0.$ (64)

We will rephrase this fact in a more convenient way, that will be more
suitable for the forthcoming computations. We define the one-parameter family
of diffeomorphisms

$\displaystyle y_{t}:Y\to\mathbb{R}^{3}$

by

$\displaystyle
y_{t}(\text{y}_{1},\text{y}_{2}):=\phi(\text{y}_{1},\text{y}_{2})+\xi+te_{i}+w_{\gamma,\xi+te_{i}}(\text{y}_{1},\text{y}_{2})\nu(\text{y}_{1},\text{y}_{2}),$
(65)

for $i=1,2,3$; $\Gamma_{t}:=y_{t}(Y)$ is the image of $y_{t}$. By
construction, $\Gamma_{t}$ is actually a submanifold of $T^{3}$ and
$\Gamma_{0}=\Gamma$. In terms of $\Gamma_{t}$, condition (64) is equivalent to

$\displaystyle\frac{d}{dt}I_{\gamma}(\Gamma_{t})|_{t=0}=0.$ (66)

By a result of Fall and Mahmoudi (see [12]),

$\displaystyle
0=\frac{d}{dt}I_{\gamma}(\Gamma_{t})|_{t=0}=\int_{\Gamma}(H_{\Gamma}+4\gamma
v_{F}+f)(\zeta,\nu_{\Gamma})d\sigma_{\Gamma}+\frac{1}{|\partial\Gamma|}\int_{\partial\Gamma}(\zeta,\nu^{\Gamma}_{\partial\Gamma})ds,$
(67)

where

$\displaystyle\zeta=\frac{d}{dt}y_{t}(x)|_{t=0}=e_{i}+\partial_{\xi_{i}}w_{\xi}\nu.$
(68)

and $\nu^{\Gamma}_{\partial\Gamma}$ is the unit normal to $\partial\Gamma$ in
$\Gamma$. The boundary term vanishes by periodicity and by the symmetries of
the problem. Using the parametrization $\beta$ of $\Gamma$ and expansions (4)
and (4), the latter relation becomes

$\displaystyle\int_{Y}\bigg{\\{}(H_{\Gamma}+4\gamma
v_{F}+f)(\beta(\text{y}_{1},\text{y}_{2}))$
$\displaystyle(e_{i}+\partial_{\xi_{i}}w_{\xi}\nu,\nu+\tilde{L}^{1}_{\xi}w_{\xi}+\tilde{Q}^{1}_{\xi}w_{\xi})$
$\displaystyle(|\phi_{x}\times\phi_{y}|+L^{1}_{\xi}w_{\xi}+Q^{1}_{\xi}w_{\xi})\bigg{\\}}d\text{y}_{1}d\text{y}_{2}=0.$

By the auxiliary equation, we know that

$\displaystyle(H_{\Gamma}+4\gamma
v_{F}+f)(\beta(\text{y}_{1},\text{y}_{2}))=\sum_{k=1}^{3}A_{k,\gamma,\xi}\nu_{k}(\text{y}_{1},\text{y}_{2})+\lambda,$
(69)

thus

$\displaystyle\sum_{k=1}^{3}A_{k,\gamma,\xi}\bigg{(}\int_{\Sigma}\nu_{k}\nu_{i}d\sigma+b_{ki}\bigg{)}+\lambda\int_{\Gamma}(\zeta,\nu_{\Gamma})d\sigma_{\Gamma}=0,$
$\displaystyle\text{for }i=1,2,3,$ (70)

with $b_{ki}=O(\gamma)$. Moreover, once again by [12], we know that

$\displaystyle\frac{d}{dt}\mathcal{L}_{3}(F_{t})=\int_{\Gamma}(\zeta,\nu_{F})d\sigma_{\Gamma},$

hence, by the volume constraint,

$\displaystyle\int_{\Gamma}(\zeta,\nu_{F})d\sigma_{\Gamma}=0,$

thus we get

$\displaystyle\sum_{k=1}^{3}A_{k,\gamma,\xi}\bigg{(}\int_{\Sigma}\nu_{k}\nu_{i}d\sigma+b_{ki}\bigg{)}=0,$
$\displaystyle\text{for }i=1,2,3.$ (71)

Since the matrix $L_{ki}$ is invertible (see Remark 5) and the coefficients
$b_{ki}$ are small, the matrix $L_{ki}+b_{ki}$ is invertible too, therefore
$A_{k,\gamma,\xi}=0$ for $k=1,2,3$.

##  5 Appendix

Proof of Remark 5

We argue by contradiction. If the statement were not true, there would exist a
vector $c=(c_{1},c_{2},c_{3})\neq 0$ such that $Lc=0$, or equivalently

$\displaystyle\sum_{j=1}^{3}\bigg{(}\int_{\Sigma}\nu_{i}(x)\nu_{j}(x)d\sigma(x)\bigg{)}c_{j}=0.$
(72)

Furthermore, writing $\nu_{i}$ as a linear combination of an orthonormal basis
$\\{e_{i}\\}_{1\leq i\leq 3}$ of $span\\{\nu_{i}\\}_{1\leq i\leq 3}$, namely

$\displaystyle\nu_{i}(x)=\sum_{k=1}^{3}\nu_{ik}e_{k}(x),$

where

$\displaystyle\nu_{ik}:=\int_{\Sigma}\nu_{i}(x)e_{k}(x)d\sigma(x),$

we can see that, setting $a_{k}:=\sum_{j=1}^{3}\nu_{jk}c_{j}$, (72) is
equivalent to

$\displaystyle
0=\sum_{k=1}^{3}\bigg{(}\int_{\Sigma}\nu_{i}(x)e_{k}(x)d\sigma(x)\bigg{)}a_{k}=\int_{\Sigma}\nu_{i}(x)a(x)d\sigma(x)$

with $a(x):=\sum_{k=1}^{3}a_{k}e_{k}(x)\in span\\{e_{i}\\}_{1\leq i\leq
3}=span\\{\nu_{i}\\}_{1\leq i\leq 3}$, so in particular $a\equiv 0$. On the
other hand, $a_{k}=0$ for any $k$ is equivalent to

$\displaystyle\int_{\Sigma}c(x)e_{k}(x)d\sigma(x)=0$ $\displaystyle\text{
}\forall 1\leq k\leq 3,$

with $c(x)=\sum_{j=1}^{3}c_{j}\nu_{j}(x)$. Thus $c\equiv 0$, that is $c_{j}=0$
for any $1\leq j\leq 3$, a contradiction.

## References

  * [1] G. Alberti, R. Choksi, F.Otto _Uniform energy distribution for an isoperimetric problem with long-range interactions_ J. Amer. Math. Soc. 22 (2009), no. 2, 569-605.
  * [2] E. Acerbi, N. Fusco, M. Morini _Minimality via second variation for a nonlocal isoperimetric problem_ Comm. Math. Phys. 322 (2013), no.2, 515-557.
  * [3] A. Ambrosetti, A. Malchiodi _Nonlinear analysis and semilinear elliptic problems_.
  * [4] A. Ambrosetti, A. Malchiodi _Perturbation methods and semilinear elliptic problems on $\mathbb{R}^{N}$_.
  * [5] M. Bonacini, R. Cristoferi _Local and global minimality results for a nonlocal isoperimetric problem on $\mathbb{R}^{N}$_ SIAM J. Math. Anal. 46, (2014), no. 4, 2310-2349.
  * [6] R. Choksi, M. A. Peletier, _On the phase diagram for microphase separation of diblock copolymers: an approach via a nonlocal Cahn-Hilliard functional_ SIAM J. Appl. Math. 69 (2009), no. 6, 1712-1738.
  * [7] R. Choksi, M. A. Peletier, _Small volume fraction limit of the diblock copolymer problem: I. Sharp-interface functional_ SIAM J. Math. Anal. 42 (2010), no. 3, 1334-1370.
  * [8] R. Choksi, M. A. Peletier, _Small volume-fraction limit of the diblock copolymer problem: II. Diffuse-interface functional_ SIAM J. Math. Anal. 43 (2011), no. 2, 739-763.
  * [9] R. Choksi, P. Sternberg _On the first and second variations of a nonlocal isoperimetric problem_ J. Reine Angew. Math. 611, (2007), 75-108.
  * [10] M. Cicalese, E. Spadaro _Droplet minimizers of an isopertimetric problem with long range interactions_ Comm. Pure Appl. Math. 66 (2013), 1298-1333.
  * [11] R. Cristoferi _On periodic critical points and local minimizers of the Ohta-Kawasaki functional_ , 2015, preprint.
  * [12] M.M. Fall, F. Mahmoudi _Hypersurfaces with free boundary and large constant mean curvature: concentration along submanifolds_ Ann. Sc. Norm. Super. Pisa Cl. Sci. (5) 7 (2008), no. 3, 407-446.
  * [13] G. Gamov _Mass defect curve and nuclear constitution_ Proc. R. Soc. Lond. Ser. A 126 (1930), 632-644.
  * [14] P.J.F. Gandy, J. Klinowski _Exact computation of the triply periodic Schwarz P minimal surface_ , Chemical Physics Letters 322 (2000), 579-586.
  * [15] D. Gilbarg, N. S. Trudinger _Elliptic partial Differential Equations of Second Order_ Classics in Mathematics, Springer-Verlag, Berlin (2001), reprint of the 1998 edition.
  * [16] D. Goldman, C. Muratov, S. Serfaty _The $\Gamma$-limit of the two-dimensional Otha-Kawasaki energy. II. Droplet arrangement via the renormalized energy_ Arch. Ration. Mech. Anal., 212 (2014), 445-501.
  * [17] H. Groemer _Geometric Applications of Fourier Series and Spherical Harmonics_ Enciclopedia Math. Appl. 61, Cambridge University Press, Cambridge, UK (1996).
  * [18] K. Grosse Brauckmann, M. Wohlgemuth _The gyroid is embedded and has constant mean curvature companions_ Calc. Var. Partial Differential Equations 4 (1996), no. 6, 499-523.
  * [19] V. Julin _Isoperimetric problem with a Coulombic repulsive term_ Indiana Univ. Math. J., to appear.
  * [20] V. Julin, G. Pisante _Minimality via Second Variation for Microphase Separation of Diblock Copolymer Melts_ preprint.
  * [21] M. Morini, P. Sternberg _Cascade of minimizers for a nonlocal isoperimetric problem in thin domains_ SIAM Journal on Mathematical Analysis, 46 (2014), 2033-2051.
  * [22] S. Müller, _Singular perturbations as a selection criterion for periodic minimizing sequences_ Calc. Var. Part. Diff. Eq. 1 (1993), no. 2, 169-204.
  * [23] C. B. Muratov _Theory of domain patterns in systems with long-range interacrions of Coulomb type_ Phys. Rev. E (3), 66 (2002), 1-25.
  * [24] G. Nardi _Schauder estimate for solutions of Poisson’s equation with Neumann boundary condition_ preprint, 2015.
  * [25] M. Ross, _Schwarz’ P and D surfaces are stable_ Differential geom. Appl. 2 (1992), no. 2, 179-195.
  * [26] X. Ren, J. Wei _Many droplet pettern in the cylindrical phase of diblock copolymer morphology_ Rev. Math. Phys. 19 (2007), no. 8, 879-921.
  * [27] X. Ren, J. Wei _Stability of spot and ring solutions of the diblock copolymer equation_ J. Math. Phys. 45 (2004), no. 11, 4106-4133.
  * [28] X. Ren, J. Wei _Wriggled lamellar solutions and their stability in the diblock copolymer problem_ SIAM J. Math. Anal. 37 (2005), no. 2, 455-489.
  * [29] X. Ren, J. Wei _Many droplet pattern in the cylindrical phase of diblock copolymer morphology_ Rev. Math. Phys. 19, (2007), no. 8, 879-921.
  * [30] X. Ren, J. Wei _Spherical solutions to a nonlocal free boundary problem from diblock copolymer morphology_ SIAM J. Math. Anal. 39 (2008), no. 5, 1497-1535.
  * [31] E. N. Spadaro _Uniform energy and density distribution: diblock copolymers’ functional_ Interfaces Free Bound. 11 (2009), no. 3, 447-474.
  * [32] P. Sternberg, I. Topaloglu _On the global minimizers of a nonlocal isoperimetric problem in two dimensions_ Interfaces Free Bound., 13 (2011), 155-169.
  * [33] I. Topaloglu _On a nonlocal isoperimetric problem on the two-sphere_ Comm. Pure Appl. Anal., 12 (2013), 597-620.


# Phase Memory of Orbital Angular Momentum in Multiple Scattering Environment

Igor Meglinski  Correspondence<EMAIL_ADDRESS>College of Engineering
and Physical Sciences, Aston University, Birmingham, B4 7ET, UK  Ivan
Lopushenko<EMAIL_ADDRESS>Optoelectronics and Measurement
Techniques, University of Oulu, P.O. Box 4500, Oulu, FI-90014, Finland  Anton
Sdobnov<EMAIL_ADDRESS>Optoelectronics and Measurement Techniques,
University of Oulu, P.O. Box 4500, Oulu, FI-90014, Finland  Alexander Bykov
<EMAIL_ADDRESS>Optoelectronics and Measurement Techniques, University
of Oulu, P.O. Box 4500, Oulu, FI-90014, Finland

###### Abstract

Recent advancements in wavefront shaping techniques have facilitated the study
of complex structured light’s propagation with orbital angular momentum (OAM)
within various media. The introduction of a spiral phase modulation to the
Laguerre–Gaussian (LG) beam during its paraxial propagation is facilitated by
the negative gradient of the medium’s refractive index’s temporal change,
resulting in an accelerated retardation in OAM twist. This approach attains
remarkable sensitivity to even the slightest variations in the medium’s
refractive index ($\sim 10^{-6}$). The phase memory of OAM is revealed as the
ability of twisted light preserving initial helical phase even propagating
through the turbid tissue-like multiple scattering medium. The results confirm
fascinating opportunities of the exploiting OAM light in biomedical
applications, e.g. such as non-invasive trans-cutaneous glucose diagnosis and
optical communication through biological tissues and other optically dense
media.

Keywords: orbital angular momentum (OAM), Laguerre–Gaussian (LG) beam,
refractive index, multiple scattering, phase memory

While using polarised light for the studies of different properties of matter
in astrophysics, material science, and biomedicine has already a long history
[1], the shaped light possessing Orbital Angular Momentum (OAM) [2] has been
added to the potential practical toolkit only recently [3]. The shaped OAM
light plays an emerging role in both classical and quantum science, and offers
fascinating opportunities for exploring new fundamental ideas, as well as for
being used in practical applications [4]. The research in the field of shaped
light with OAM achieved high recognition in generation and characterization of
the exotic vector laser beams [5], improving telecommunication technologies
[6], optical trapping of cells [7] and micro-particles [8]. Despite the
technical challenges, the shaped OAM light is emerged as one of the most
exciting front-lines of contemporary research, promising to defy the
deficiencies of current optical sensing techniques [9]. OAM-based twisted
light holds great potential for various biomedical applications in the field
of biological tissue diagnosis. Apart from the precise optical manipulation
and sorting of biological particles and cells utilizing optical optical
tweezers [10, 11] the emerging biomedical applications based on OAM-based
twisted light include imaging and microscopy [12], generation twisted light
from typical on-chip devices [13], optical communication in biological media
[14] and other [9, 15].

We explore the prediction capacities of OAM of Laguerre-Gaussian (LG) beams by
analyzing their helical wavefront change along propagation through a tissue-
like medium. In frame of the paraxial approximation LG beam is defined as [16,
17, 18]:

$\displaystyle
LG_{p}^{\ell}(\rho,\phi,z)=\sqrt{\frac{2p!}{\pi(\lvert\ell\rvert+p)!{w}^{2}(z)}}\left[\frac{\rho\sqrt{2}}{{w}(z)}\right]^{\lvert\ell\rvert}\times$
(1) $\displaystyle\times
L_{p}^{\lvert\ell\rvert}\left[\frac{2\rho^{2}}{{w}^{2}(z)}\right]\exp\left[-\frac{\rho^{2}}{{w}^{2}(z)}\right]\exp[i(2p+\lvert\ell\rvert+1)\times$
$\displaystyle\times\arctan(z/z_{R})]\exp\left[\frac{-ik\rho^{2}z}{2(z^{2}+z_{R}^{2})}\right]\exp[-i\ell\varphi]\exp[-ikz].$

Here, $k=2\pi/\lambda$, $\lambda$ is the wavelength of laser radiation,
${w}(z)={w}(0)\sqrt{1+(z/z_{R})^{2}}$, $z_{R}=\pi{w}^{2}(0)/\lambda$, ${w}(0)$
corresponds to the zero-order Gaussian beam waist and is adjusted to fit
experimental image, $\\{\rho,\phi,z\\}$ represents the cylindrical coordinates
system utilized for characterization of beam propagation along $z$ axis in
terms of radial ($\rho$) and angular ($\phi$) coordinates.

Phase evolution of the LG beam along its propagation in a medium is defined as
[18]:

$\displaystyle\Psi(\rho,\phi,z)=\arg\left(LG_{p}^{\ell}(\rho,\phi,z)\right)=$
(2) $\displaystyle=\frac{-k\rho^{2}z}{2(z^{2}+z_{R}^{2})}-\ell\phi-kz+G(z),$

where $G(z)=(2p+\lvert\ell\rvert+1)\arctan(z/z_{R})$ corresponds to the Gouy
phase [19]. The helical phase of an LG beam refers to the phase front that
winds around the beam’s axis, is a result of the azimuthal phase term in the
beam’s electric field expression, which corresponds to the topological charge
$(\ell)$ and radial index $(p)$.

Figure 1:  (a) Schematic representation of the central part of the Mach-
Zehnder interferometer experiment: the LG beam, carrying OAM imparted by a
spatial light modulator (SLM), traverses a cuvette containing the probed
medium; interference between the transmitted LG beam and a reference plane
wave (an initially expanded Gaussian beam) is then analyzed for intensity
and/or retrieved phase distributions. (b) Spiral modulation of the intensity
of $LG^{5}_{0}$ beam observed experimentally (left) during the gradual
increase of the medium’s refractive index, alongside the prediction from
theoretical modeling [20, 21] (right); Blue line corresponds to accelerated
retardation of OAM twist along LG beam propagation in the medium. The inset
(left) depicts relative phase twist ($\Theta=\dfrac{\Psi}{\ell}$) of OAM of
$LG^{5}_{0}$ beam during gradual increase of medium refractive index ($\Delta
n=3.69\times 10^{-4}$). (c= Relative phase twist of OAM of $LG^{3}_{0}$
(circles) and $LG^{5}_{0}$ (squares) beams during the gradual increase of
medium refractive index within a range of $2\times 10^{-5}$; Highlighted
colored areas ($\Delta n=5\times 10^{-6}$) feature corresponding twist of OAM
for $LG^{3}_{0}$ and $LG^{5}_{0}$ beams, respectively.

The LG beam retains its helical structure (LABEL:eq:LGbeamPhase) while
propagating through the transparent medium, preserving the phase composition
$\Psi(\rho,\phi,z)$ associated with OAM and the characteristic helix-like
wavefront (Fig. 1). Here and below we consider, so-called, scalar LG beams
[22] with homogeneous annular polarisation distribution ($LG^{3}_{0}$ and
$LG^{5}_{0}$). The observed an accelerated retardation in the OAM rotation
along with relative twist of phase ($\Theta=\dfrac{\Psi}{\ell}$) of the LG
beam (see Fig. 1-b), as it propagates through the medium with the negative
gradient of the temporal change of the refractive index, are due to spatial
dispersion of Poynting vector trajectories, defined as [18, 23]:

$\displaystyle
L(r_{0},\varphi_{0},\zeta_{s})=\int_{0}^{\zeta_{s}}\sqrt{\left(\frac{dr}{d\zeta}\right)^{2}+r^{2}\left(\frac{d\varphi}{d\zeta}\right)^{2}}d\zeta,$
(3) $\displaystyle
r(\zeta)=r_{0}\sqrt{1+4\zeta^{2}},\quad\varphi(\zeta)=\frac{\ell}{2r_{0}^{2}}\arctan(2\zeta)+\varphi_{0}.$

Figure 2:  (a) Phase distribution (speckle patterns) observed experimentally
for the $LG^{3}_{0}$ beam propagated through the low $(d/l^{*}=2)$ scattering
(left) and multiple $(d/l^{*}=9.6)$ scattering (right) media. The axial
annular zone (embossed by contours) corresponds to the $LG^{3}_{0}$ beam as if
it were passing through a medium devoid of scattering. (b) Phase variations
manifest at the single speckle grain within a deliberately chosen sector of
the $LG^{3}_{0}$ axial annular area for low-scattering (indicated by black
circles) and multiple-scattering (represented by red squares) media. The
observed phase changes are contingent upon the initial phase configuration
($-3\pi/10\leq\Psi\leq 3\pi/10$) established at SLM (schematically shown in
inset). (c) The ensuing alterations in the phase mapping of the speckle
pattern within the designated areas ($150\times 150~{}\mu m$) highlighted in
the $LG^{3}_{0}$ axial annular domain (see Fig. 2-a), aligning with the
prescribed initial phase configuration established at the SLM
($-3\pi/10\leq\Psi\leq 3\pi/10$). The upper and lower rows present,
respectively, scenarios for low and multiple scattering environment; scale bar
corresponds to $150~{}\mu m$. Video 1 presents the speckle pattern phase
dynamics corresponding to the phase twist at the SLM.

Here, $L(r_{0},\varphi_{0},\zeta_{s})$ is the length of Poynting vector
trajectory along LG beam propagation; $r_{0},\varphi_{0}$ define the radial
and angular coordinates of the trajectory at the starting point ($z=0$);
$r=\dfrac{\rho}{{w}(0)}$ and ${\zeta=\dfrac{z}{k{w}^{2}(0)}}$ are the
dimensionless cylindrical coordinates; $\zeta_{s}$ defines the point to which
trajectory length is evaluated; $r(\zeta)$, $\varphi(\zeta)$ are the
corresponding coordinates at the trajectory ($0\leq\zeta\leq\zeta_{s}$), the
direction of which ($\mathbf{p}=\\{p_{\rho},p_{\phi},p_{z}\\}$) defined as
[16, 18]:

$\displaystyle p_{\rho}=\frac{\omega k\rho
z}{\left(z_{R}^{2}+z^{2}\right)}\lvert LG^{\ell}_{p}(\rho,\phi,z)\rvert^{2},$
(4) $\displaystyle p_{\phi}=\frac{\omega l}{\rho}\lvert
LG^{\ell}_{p}(\rho,\phi,z)\rvert^{2},$ $\displaystyle p_{z}=\omega k\lvert
LG^{\ell}_{p}(\rho,\phi,z)\rvert^{2}\text{, }$

where $\omega$ is the angular frequency of the light.

In LG beam the Poynting vector trajectories follow the spiral paths around the
optical axis of the beam [18, 23], as presented in Fig. 1-a. Gradual increase
of the medium’s refractive index gives a proportional elongation of these
spiral-like trajectories and increased length deviation between them. The
increment of the phase component
$\left(\frac{-k\rho^{2}z}{2(z^{2}+z_{R}^{2})}-kz\right)$ along the longer
trajectories leads to the emergence of a gradient in the LG beam phase
(LABEL:eq:LGbeamPhase) that in turn manifests itself as a twisting effect in
its transverse distribution (see Fig. 1-b). Accordingly, when the
characteristic size of the medium (e.g. the thickness of the medium) is much
larger than the wavelength ($d\gg\lambda$) even a minor variation in the
medium refractive index leads to a significant alteration of the transverse
spatial phase distribution of the LG beam (see Fig. 1-c). A prominent twist of
OAM is observed experimentally up to minuscule changes of the medium
refractive index ($\Delta n=10^{-6}$ (see Fig. 1-c).

A gradual decrease in the medium’s refractive index induces a spiral
modulation of the LG beam, observed via interference with a plane wave (an
expanded Gaussian beam), as an accelerated retardation of the OAM twist (see
Fig. 1-b). The experimental results obtained align closely with the
theoretical predictions as presented in Fig. 1-b. In case of the LG beams with
lower topological charge ($\ell=3$), the Poynting vector trajectories tend to
less tightly wound around the optical axis [18, 23]. Therefore, the OAM twist
for $LG^{3}_{0}$ beams with lower topological charge exhibits higher accuracy
in predicting even negligible changes in the refractive index compared to the
$LG^{5}_{0}$ beam (see Fig. 1-c).

While the twist of OAM of LG beam along propagation through a transparent
medium (see Fig. 1) primarily arises due to gradual increase of the medium
refractive index, causing a corresponding enhancement of the lengths of spiral
the Poynting vector trajectories (3), in turbid medium, the scattering of the
LG beam leads to the speckle pattern formation. Arising from the superposition
of partial components of the helical wavefront, the speckle interference
pattern manifests spatially varying intensity and phase distributions. This
occurrence contributes to the disruption of the LG beam helical structure
composition in low-scattering media and its complete degradation in multiple
scattering media (as presented in Fig.A1, see Appendix A).

Disperse turbid media exhibiting low or multiple scattering of light are
characterised by their optical depth ($\sim d/l^{*}$), which serves as a
quantitative measure of the extent of attenuation and scattering strength
experienced by light along its propagation through. Here, $d$ is the thickness
of the scattering medium and $l^{*}$ is the transport mean free path. A
commonly used guideline is that if $d/l^{*}\sim 10$ or larger the medium is
considered as multiple or diffuse scattering [24], whereas for the single and
an intermediate (‘snake-like photons’) scattering $d/l^{*}$ is lower ($\sim
3-6$).

In a low scattering medium ($d/l^{*}\sim 2$), the LG beam propagates with
minimal disruption, enabling it to maintain its initial OAM state and
doughnut-like spatial intensity profile [21, 25], as well as the helical phase
structure (Fig. 2-a). The multiple scattering ($d/l^{*}\sim 10$) results in a
diffusive spread of the LG beam’s intensity profile and the destruction of the
helical phase front, leading to creation of complex speckle pattern (see Fig.
2-a). Nevertheless, despite the influence of strong diffuse scattering, the
phase speckle pattern maintains the modulation of the initial phase of the LG
beam, representing a manifestation of the OAM’s memory in multiple scattering.

The preservation of the helical phase structure by LG beam along propagation
through the scattering medium (see Fig. 2-b and Fig. 2-c) is attributed due to
the similarities in the deviations of the Poynting vector trajectories around
the beam axis causing due to rotational symmetry.

Influenced by the scattering medium, the Poynting vector trajectories
significantly transform in relation to their original spiral paths (Fig. 3-a),
resulting in different partial components of the helical wavefront of the LG
beam experiencing varying degrees of phase distortions.

Due to the rotational symmetry, the deviations in the Poynting vector
trajectories tend to be similar in magnitudes and directions around the beam
axis. Therefore, despite the substantial distortions of phase induced by
multiple scattering, the overall preservation of OAM of the LG beam is clearly
observed in the phase dynamics of the speckle pattern, corresponding to the
modulation of the initial phase at the SLM (see Fig. 2-c and Video 1).
Remarkably, during propagation through the multiple scattering medium
($d/l^{*}\sim 10$), the phase memory of OAM is distinctly pronounced within
the axial annular region of the LG beam, as embossed in Fig. 2-a (see also
Video 1), with no discernible presence in the central and external areas
beyond the annular region.

When light passes through a turbid tissue-like medium, even with a low
scattering ($d/l^{*}<2$), the wavefront undergoes rapid deformation, resulting
in a distinct speckle pattern. Notably, while traversing a more scattering
medium, polarization scrambling occurs on a different length scale [26],
leading to the adoption of a chaotic or disordered nature in the polarization
state of light.

In other words, the initial orientation of the electric field vector
representing the polarization of light undergoes unpredictable changes due to
interactions with scattering elements in the medium. Polarization scrambling
results in a loss of the initial alignment, leading to a more random
distribution of polarization states. The degree of polarization of light is
defined as the measure of the alignment of the electric field vectors within
the light wave [27]:

$\displaystyle
P=\dfrac{2L}{l_{s}}\mathrm{sinh}\left(\dfrac{l_{s}}{\xi_{L}}\right)\exp\left(-\dfrac{L}{\xi_{L}}\right)\text{,
}$ (5)

where $\xi_{L}$ is the characteristic depolarization length
$\left(\xi_{L}=l_{s}/\left(\sqrt{3\mathrm{ln}\frac{10}{7}}\right)\right)$,
$l_{s}$ is the elastic mean free path
($l_{s}=l^{*}(1-\langle\cos\theta\rangle)$), $\theta$ is the scattering angle,
and $L$ corresponds to the Poynting vector trajectories (3) of LG beam (see
Fig.3-a).

In disordered scattering media, such as biological tissues and turbid tissue-
like materials, light encounters multiple scattering events with particles or
structures [28], leading to the formation of branched transport channels [29].
In a similar fashion, the preservation mechanism within the axial annular
region is associated with the exponential localization of the transmission
eigenchannels, aligning with the shortest optical pathway. [30].

Figure 3: (a) The Poynting vector trajectories, schematically depicted in
Cartesian coordinates, for the scalar linearly $x$-polarized $LG^{3}_{0}$ beam
as if it were passing from the SLM to the axial annular zone (embossed by
contours) at the detector through a medium devoid of scattering. (b) Radial
distribution of the degree of polarization ($P$) and relative phase shift
($\Delta\Psi$) occurring along the open channel as propagated through a
multiple scattering medium ($d/l^{*}\sim 10$) with the characteristic length
of depolarization ($\xi_{L}$). The colourful background represents
longitudinal profile of the intensity of $LG^{3}_{0}$ beam along the radial
component ($\rho$).

Fig. 3-b shows gradual transformation of polarization degree along the radial
component ($\rho$) of the $LG^{3}_{0}$ beam, corresponding to the spread of
the pathways ($L/\xi_{L}$) within the medium from SLM to the detecting area as
presented in Fig. 3-a. In other words, upon propagation of the scalar linearly
$x$-polarized $LG^{3}_{0}$ beam through multiple scattering medium
($d/l^{*}\sim 10$) the polarization remains effectively along the
eigenchannel, whereas the relative phase distortion
$\Delta\Psi=\dfrac{\Psi|_{L/\xi_{L}}-\Psi|_{z/\xi_{L}}}{2\pi}$ becomes minimal
within the axial annular zone (see Fig. 3-b).

Eventually, upon propagation through a turbid tissue-like scattering medium
($d/l^{*}\sim 10$), the phase memory of Orbital Angular Momentum (OAM)
exhibits spatial variability when examined along $\rho$ — the radial position
(schematically shown in Fig. 3-a). The overall helical phase structure
associated with the OAM remains preserved, clearly observed in the single
speckle grains within an arbitrarily chosen area within the annular region of
the LG beam (see Fig. 2-a). However, in the areas outside of the doughnut-like
ring, the phase memory effect disappears. Conversely, in a low-scattering
environment ($d/l^{*}\sim 2$), the phase memory is more consistently observed
in all regions (see Fig.B1 in Appensix B, and Video 2).

Computational analyses corroborate these experimental findings, affirming that
the phase memory of OAM in LG beams propagated through turbid tissue-like
scattering media is not absolute. It is susceptible to various factors,
including scattering, absorption, and anisotropy of scattering, potentially
leading to the degradation or alteration of the OAM’s phase content.

The obtained results offer fascinating opportunities of the exploiting OAM
light in biomedical applications, e.g. such as non-invasive trans-cutaneous
glucose diagnosis and optical communication through biological tissues and
other disperse multiple scattering materials.

## Optical setup and processing

The experimental configuration employed for measurements is illustrated in Fig
C1 in Appendix C. In summary, a Gaussian light beam at $640~{}nm$ is divided
into sample and reference beams using a polarizing beam splitter. A spatial
light modulator (PLUTO-2-NIR-011, Holoeye, Germany) is utilized in the
reference arm to generate various LG beams carrying OAM. These OAM-carrying
beams traverse the sample – scattering environment (refer to Appendix D for
details). Simultaneously, the reference Gaussian beam undergoes expansion
through a series of lenses and is directed to a beam splitter, where it
interferes with the sample beam. Subsequently, the interference pattern is
recorded using a charge-coupled device (CCD) camera. To establish a comparison
between experimental data and theoretical simulations, the methodology
outlined in Appendix E is employed.

The utilized setup facilitates the acquisition of interference patterns in
both on-axis [31] and off-axis [32] regimes. The determination of relative OAM
twist in the on-axis regime involves assessing the change in polar angle of a
petal caused by its rotation around the LG beam center. Conversely, in the
off-axis regime, the LG beam phase is retrieved through a fast Fourier
transform approach [32]. A detailed elucidation of the experimental data
processing methodology is presented in Appendix F and G, respectively.

###### Acknowledgements.

Authors acknowledge the support from the Leverhulme Trust and The Royal
Society (Ref. no.: APX111232 APEX Awards 2021), UKKi UK-Israel innovation
researcher mobility and Academy of Finland (grant projects 325097, 351068).

## References

  * Mischenko  _et al._ [2006] M. Mischenko, L. Travis, and A. Lacis,  _Multiple Scattering of Light by Particles, Radiative Transfer and Coherent Backscattering_ (Cambridge University, Cambridge, 2006).
  * Torres and Torner [2011] J. Torres and L. Torner,  _Twisted Photons: Applications of Light with Orbital Angular Momentum_ (Wiley-VCH Verlag GmbH, Germany, 2011).
  * Rubinsztein-Dunlop and et al. [2016] H. Rubinsztein-Dunlop and et al., Roadmap on structured light, J. Opt. 19, 013001 (2016).
  * Shen and et al. [2019] Y. Shen and et al., Optical vortices 30 years on: OAM manipulation from topological charge to multiple singularities, Light Sci. Appl. 8, 90 (2019).
  * Rego and et al. [2019] L. Rego and et al., Generation of extreme-ultraviolet beams with time-varying orbital angular momentum, Science 364, eaaw9486 (2019).
  * Wang  _et al._ [2012] J. Wang, J. Yang, I. Fazal, and et al., Terabit free-space data transmission employing orbital angular momentum multiplexing, Nat. Photon. 6, 488–496 (2012).
  * Zhu  _et al._ [2020] R. Zhu, T. Avsievich, A. Popov, and I. Meglinski, Optical tweezers in the studies of red blood cells, Cells 9, 545 (2020).
  * Padgett and Bowman [2012] M. Padgett and R. Bowman, Tweezers with a twist, Nat. Photon. 5, 343 (2012).
  * Weng and Pan [2023] Y. Weng and Z. Pan, Orbital angular momentum based sensing and their applications: A review, J. Light. Technol. 41, 2007 (2023).
  * Bustamante  _et al._ [2021] C. J. Bustamante, Y. R. Chemla, S. Liu, and M. D. Wang, Optical tweezers in single-molecule biophysics, Nat. Rev. Methods Primers 1, 25 (2021).
  * Avsievich  _et al._ [2020] T. Avsievich, R. Zhu, A. Popov, A. Bykov, and I. Meglinski, The advancement of blood cell research by optical tweezers, Rev. Phys. 5, 100043 (2020).
  * Ritsch-Marte [2017] M. Ritsch-Marte, Orbital angular momentum light in microscopy, Philos. Trans. A Math. Phys. Eng. Sci. 375, 20150437 (2017).
  * Yang  _et al._ [2021] H. Yang, Z. Xie, H. He, Q. Zhang, and X. Yuan, A perspective on twisted light from on-chip devices, APL Photonics 6, 110901 (2021).
  * Perez  _et al._ [2022] N. Perez, D. Preece, R. Wilson, R. Wilson, and A. Bezryadina, Conservation of orbital angular momentum and polarization through biological waveguides, Sci. Rep. 12, 14144 (2022).
  * Chen  _et al._ [2021] J. Chen, C. Wan, and Q. Zhan, Engineering photonic angular momentum with structured light: a review, Adv. Photonics 3, 064001 (2021).
  * Allen  _et al._ [1999] L. Allen, M. J. Padgett, and M. Babiker,  _Progress in Optics_ , edited by E. Wolf, Vol. XXXIX (Elsevier, North-Holland, Amsterdam, 1999) pp. 291–372.
  * Rosales-Guzmán  _et al._ [2018] C. Rosales-Guzmán, B. Ndagano, and A. Forbes, A review of complex vector light fields and their applications, J. Opt. 20, 123001 (2018).
  * Berry and McDonald [2008] M. V. Berry and K. T. McDonald, Exact and geometrical optics energy trajectories in twisted beams, J. Opt. A: Pure Appl. Opt. 10, 035005 (2008).
  * Andrews and Babiker [2012] D. Andrews and M. Babiker,  _The Angular Momentum of Light_ (Cambridge University Press, Cambridge, 2012).
  * Lopushenko  _et al._ [2023] I. Lopushenko, O. Sieryi, A. Bykov, and I. Meglinski, Exploring the evolution of circular polarized light backscattered from turbid tissue-like disperse medium utilizing generalized Monte Carlo modeling approach with a combined use of Jones and Stokes–Mueller formalisms, J. Biomed. Opt. 29, 052913 (2023).
  * Doronin  _et al._ [2019] A. Doronin, N. Vera, J. P. Staforelli, P. Coelho, and I. Meglinski, Propagation of cylindrical vector laser beams in turbid tissue-like scattering media, Photonics 6, 56 (2019).
  * Rosales-Guzmán  _et al._ [2018] C. Rosales-Guzmán, B. Ndagano, and A. Forbes, A review of complex vector light fields and their applications, J. Opt. 20, 123001 (2018).
  * Bliokh  _et al._ [2013] K. Y. Bliokh, A. Y. Bekshaev, A. G. Kofman, and F. Nori, Photon trajectories, anomalous velocities and weak measurements: a classical interpretation, New J. Phys. 15, 073022 (2013).
  * Yoo  _et al._ [1990] K. M. Yoo, F. Liu, and R. R. Alfano, When does the diffusion approximation fail to describe photon transport in random media?, Phys. Rev. Lett. 64, 2647 (1990).
  * Gianani  _et al._ [2020] I. Gianani, A. Suprano, T. Giordani, N. Spagnolo, F. Sciarrino, D. Gorpas, V. Ntziachristos, K. Pinker, N. Biton, J. Kupferman, and S. Arnon, Transmission of vector vortex beams in dispersive media, Adv. Photon. 2, 036003 (2020).
  * de Aguiar  _et al._ [2017] H. B. de Aguiar, S. Gigan, and S. Brasselet, Polarization recovery through scattering media, Sci. Adv. 3, e1600743 (2017).
  * Bicout  _et al._ [1994] D. Bicout, C. Brosseau, A. S. Martinez, and J. M. Schmitt, Depolarization of multiply scattered waves by spherical diffusers: Influence of the size parameter, Phys. Rev. E 49, 1767 (1994).
  * Jacques [2013] S. L. Jacques, Optical properties of biological tissues: a review, Phys. Med. Biol. 58, R37 (2013).
  * Brandstötter  _et al._ [2019] A. Brandstötter, A. Girschik, P. Ambichl, and S. Rotter, Shaping the branched flow of light through disordered media, Proc. Natl. Acad. Sci. U.S.A. 116, 13260 (2019).
  * Yılmaz  _et al._ [2019] H. Yılmaz, C. Hsu, A. Yamilov, and H. Cao, Transverse localization of transmission eigenchannels, Nat. Photonics 13, 352 (2019).
  * Cui  _et al._ [2019] S. Cui, B. Xu, S. Luo, H. Xu, Z. Cai, Z. Luo, J. Pu, and S. Chávez-Cerda, Determining topological charge based on an improved fizeau interferometer, Opt. Exp. 27, 12774 (2019).
  * Vayalamkuzhi  _et al._ [2021] P. Vayalamkuzhi, S. Bhattacharya, _et al._ , Transform-based phase retrieval techniques from a single off-axis interferogram, Appl. Opt. 60, 5523 (2021).
  * Kumar and Nishchal [2019] P. Kumar and N. K. Nishchal, Modified mach–zehnder interferometer for determining the high-order topological charge of laguerre–gaussian vortex beams, J. Opt. Soc. Am. A 36, 1447 (2019).
  * Jiménez Riobóo  _et al._ [2009] R. Jiménez Riobóo, M. Philipp, M. Ramos, and J.-K. Krüger, Concentration and temperature dependence of the refractive index of ethanol-water mixtures: Influence of intermolecular interactions, The European Physical Journal E 30, 19 (2009).
  * Wrobel  _et al._ [2015] M. S. Wrobel, A. P. Popov, A. V. Bykov, M. Kinnunen, M. Jedrzejewska-Szczerska, and V. V. Tuchin, Measurements of fundamental properties of homogeneous tissue phantoms, J. Biomed. Opt. 20, 045004 (2015).
  * Otsu [1979] N. Otsu, A threshold selection method from gray-level histograms, IEEE Trans. Syst. Man Cybern. 9, 62 (1979).
  * Gonzalez  _et al._ [2004] R. Gonzalez, R. Woods, and S. Eddins,  _Digital Image Processing Using MATLAB_ (Pearson Prentice Hall, New Jersey, 2004).

APPENDIX A: Addition to Figure 2

Figure A1:  Experimentally observed spatial intensity distributions for
$LG^{3}_{0}$ and $LG^{5}_{0}$ beams propagated through a transparent medium
(a) and (d), as well as through the low $(z/l^{*}=2)$ and multiple
$(z/l^{*}=9.6)$ scattering environments, respectively, (b) and (e) and (c) and
(f). Scale bar is equal to $750~{}\mu m$.

APPENDIX B: Addition to Video 2

Figure B1:  Phase distribution (speckle patterns) observed experimentally for
the $LG^{3}_{0}$ beam propagated through the low $(d/l^{*}=2)$ scattering (a)
and multiple $(d/l^{*}=9.6)$ scattering (b) media. The axial annular zone
(embossed by contours) corresponds to the $LG^{3}_{0}$ beam as if it were
passing through a medium devoid of scattering. The resulting changes in the
phase distribution of speckle patterns observed within specified areas (size:
$150\times 150~{}\mu m$) selected at the center (left), axial annular region
(center), and outside (right) of the $LG^{3}_{0}$ beam profile (as indicated
in (a) and (b), presented for both low scattering (c) and multiple scattering
(d) environments. Respectively, the phase variations observed at the single
speckle grain within center (left), axial annular region (center), and outside
(right) of the $LG^{3}_{0}$ beam both low (e) and multiple (f) scattering. The
observed phase changes are contingent upon the initial phase configuration
($-3\pi/10\leq\Psi\leq 3\pi/10$) established at SLM see Video 2.

## APPENDIX C: Experimental setup

The modified Mach–Zehnder-based interferometer [33] is used to examine an
evolution of OAM of the LG beams propagated through the multiple scattering
environment (Fig. C1).

A coherent Gaussian light beam ($40~{}mW$, BioRay laser diode, Coherent, USA
with coherence length $>20~{}cm$) emitting at $640~{}nm$ serves as the light
source. To clear the optical mode, the laser beam is focused into a single-
mode optical fiber (Thorlabs, USA) F. The resulting output beam is collimated
using a beam collimator (Thorlabs, USA) BC to obtain a Gaussian beam with a
$1.6~{}mm$ waist diameter. Additionally, a polarizer (Thorlabs, USA) P is used
after the collimator to achieve horizontal linear polarization. The Gaussian
beam is then split into sample and reference beams using a polarizing beam
splitter (Thorlabs, USA) PBS.

Figure C1: LD – laser diode. P – polarizer. FM – fiber mount. BC – beam
colimator. SLM – spatial light modulator. M1, M2 – mirrors. L1, L2, L3, L4 –
lenses. PH – pinhole. PBS – polarizing beam splitter. HWP – half wave plate. S
– cuvette filled with the sample liquid. BS – beam splitter. NF – neutral
filter. O – objective. CCD – camera. The detailed description of the optical
setup is presented in the main article text.

The sample beam illuminates the phase-only spatial light modulator
(PLUTO-2-NIR-011, Holoeye, Germany) SLM, operating in reflective regime. To
produce LG beams with different moments, the corresponding forked diffraction
patterns are generated on the SLM. The diffracted light from the SLM is
directed using a set of mirrors M to lens L3 ($f=45~{}mm$, Thorlabs, USA).
This lens is used to focus the first-order diffraction through a pinhole
(Thorlabs, USA). The LG beam is then re-collimated using lens L4 ($f=45~{}mm$,
Thorlabs, USA). Finally, the sample beam passes through the sample S. The
reference beam goes through a half-wave plate (Thorlabs, USA) to control the
polarization orientation of the Gaussian beam. This reference beam is expanded
by lenses L1 ($f=30~{}mm$ Thorlabs, USA) and L2 ($f=70~{}mm$, Thorlabs, USA),
and directed to the beam splitter (Thorlabs, USA) BS, where the expanded
Gaussian beam interferes with the LG beam. The interference pattern is then
registered using a CMOS camera (DCC3240M, $1280\times 1024$, Thorlabs, USA)
CAM in combination with an objective (10$\times$, Nikon, Japan) O. The waist
beam diameter for $LG^{3}_{0}$ and $LG^{5}_{0}$ beams at the detector is $2.7$
and $3~{}mm$, respectively.

The setup enables the acquisition of interference patterns in both on-axis
[31] and off-axis [32] regimes. In the on-axis scenario, the definition of
relative OAM twist involves tracking the change in the polar angle of a petal
caused by its rotation around the center of the LG beam. Meanwhile, in the
off-axis regime, the retrieval of the LG beam phase is accomplished through a
fast Fourier transform approach [32].

## APPENDIX D: Samples: scattering environment

The selection of the ethanol-water solution as an optically transparent medium
for experimental inquiry is underpinned by careful consideration. The
quantification of the relative OAM twist, stemming from temperature-dependent
changes in the refractive index of the solution, is systematically conducted.
A precisely measured volume of $5~{}ml$ of an ethanol-water solution,
characterized by a water concentration of $51.89~{}mol\%$, is carefully
introduced into a glass cuvette. The cuvette, with a total thickness of
$5.6~{}mm$, features glass walls, each having a thickness of $1~{}mm$.

The investigation thoroughly addresses the intricate thermal dependencies
governing the refractive index of the ethanol-water amalgamation, as
comprehensively detailed in the authoritative work [34]. Subsequently, the
cuvette containing the liquid sample undergoes refrigeration until reaching a
temperature of $7.8\pm 0.5^{\circ}C$. After the cooling phase, the cuvette is
reintegrated into the experimental apparatus and subjected to a $25~{}minute$
duration at room temperature ($21^{\circ}C$). This deliberate protocol induces
modulations in the refractive index of the sample liquid during the heating
process.

The nuanced alterations in interference patterns, reflective of variations in
the refractive index, are meticulously captured by a high-speed camera
throughout the entirety of the 25-minute heating interval. The recording is
conducted at a frame rate of 10 frames per second, each frame having a
$1~{}msec$ exposure time. These experimental procedures adhere to rigorous
standards to ensure the fidelity and reliability of the acquired data.

Throughout the experimental duration, the sample temperature is vigilantly
monitored at 10-second intervals utilizing a FLIR C5 thermal camera from
Teledyne FLIR, USA.

The measurements are conducted in a low scattering ($d/l^{*}=2$) phantom with
a thickness of $1~{}mm$ ($\mu_{s}=10~{}mm^{-1}$) and the multiple-scattering
one ($d/l^{*}\sim 10$) with a thickness of $8~{}mm$ ($\mu_{s}=6~{}mm^{-1}$).
The meticulous methodology employed in the preparation of the phantoms is
comprehensively expounded upon in the work by Wrobel et al., as documented in
[35].

Figure F1: Schematic presentation of (a) Main steps and principles for
determining the relative twist of petals in the on-axis regime; (b) Main steps
and principles for the phase retrieval procedure in the off-axis regime.

## APPENDIX E: Computational modelling

Expressions (1)–(LABEL:eq:LGbeamPhase) describe intensity and phase
distributions of the LG beam in the medium with $n=1$. To assess the
conversion of OAM in the tissue-like scattering medium ($n\neq 1$) a large set
(e.g. $N_{ph}\sim 10^{9}$) of LG beam photon trajectories with starting points
$\rho_{0_{i}},\phi_{0_{i}},i\in[1...N_{ph}]$ is generated [21]. Azimuths phase
for each trajectory can be easily calculated as $-\ell\phi$ as follows from
(LABEL:eq:LGbeamPhase). The length $L_{i}$ of each trajectory is estimated
according to (3). Within cuvette (or within the sample) we estimate trajectory
lengths as $\Delta
L_{i}(r_{0_{i}},\varphi_{0_{i}})=L(r_{0_{i}},\varphi_{0_{i}},\zeta_{2})-L(r_{0_{i}},\varphi_{0_{i}},\zeta_{1})$,
where $\zeta_{2}-\zeta_{1}$ corresponds to cuvette (sample) thickness in
dimensionless units introduced in (3). Due to the non-zero contrast between
refractive indices of free space, glass and interior of the cuvette (sample)
the trajectories are additionally refracted according to Snell’s law, slightly
increasing $\Delta L_{i}$ values.

We estimate cuvette (sample) influence on the LG beam propagation via phase
retardation caused by the increasing path length of light within medium:

$\Delta\Psi_{i}=\dfrac{2\pi n\Delta L_{i}}{\lambda}.$ (6)

Here, $\Delta\Psi_{i}$ corresponds to the phase retardation along $i$’th
trajectory. By comparing the phase before cuvette $\Psi|_{L_{i}}$ and phase
after cuvette $\Psi|_{L_{i}+\Delta L_{i}}=\Psi|_{L_{i}}+\Delta\Psi_{i}$ for
different values of $n$, we obtain phase patterns that are twisted
differently, as seen on Fig. 1.

Upon computation of the intensity and phase beyond the cuvette (sample), it
becomes feasible to derive the interference pattern of LG beams with a plane
wave, yielding the characteristic petal pattern on the screen with
$\lvert\ell\rvert$ petals. These petals undergo a twist in correspondence to
variations in the refractive index $n$ within the cuvette interior. The
quantification of this twist is elucidated in the ensuing procedure, and the
obtained results exhibit notable concordance with experimental measurements,
as illustrated in Fig. 1.

The proposed computational methodology facilitates the assessment of
pathlengths $\Delta L_{i}$ in scenarios where the cuvette (sample) interior
manifests turbidity. In such instances, the trajectories of LG beams are
construed in the context of Monte Carlo photons [21], undergoing multiple
scattering events dictated by the medium’s scattering coefficient $\mu_{s}$,
absorption coefficient $\mu_{a}$, and anisotropy of scattering parameter $g$
($g=\langle\cos\theta\rangle$).

## APPENDIX F: Characterization of twist of the LG beam propagated through the
medium

The quantification of the relative OAM twist in the on-axis regime is
established by defining the change in polar angle of a petal due to its
rotation around the center of the LG beam. To ascertain this twist, a polar
coordinate system is introduced, with the origin situated at the center of the
LG beam for each recorded experimental image. Subsequently, a binary
representation is generated for each image by binarizing the intensity values,
with those exceeding a defined threshold set to one, and all other values set
to zero. The luminance threshold is determined utilizing the Otsu approach
[36], and can be manually adjusted either for an individual image or across a
series of images to ensure clarity in visualizing each petal on the binary
representation.

The identification of pixels corresponding to the outer boundaries of each
illuminated region in the binary image is achieved through the utilization of
the Moore-Neighbor tracing algorithm, adapted to adhere to Jacob’s stopping
criteria [37]. The preeminent bright spots are associated with the
interference pattern petals of the LG beam, contingent upon the selection of
an optimal contrast threshold. All other bright spots are deemed artifacts,
and their respective boundary pixels are systematically excluded from
subsequent data analysis.

By discerning the boundary pixels, individual polygons are meticulously
constructed for each petal, preserving their inherent geometric
characteristics. For every polygon, the polar coordinates of its centroid are
calculated. The derivation of final relative twist values is accomplished by
subtracting the angular coordinates of the centroids between different images,
considering the spatial orientation of the petals and ensuring a seamless
$0$-to-$2\pi$ transition. A schematic representation delineating the
underlying principles and key steps of this method is presented in Fig. F1-a.

## APPENDIX G: Phase retrieval

In the off-axis regime, the retrieval of the LG beam phase is accomplished
through a straightforward signal processing approach, as elucidated in the
work by Vayalamkuzhi et al. [32]. The schematic illustrating the principal
steps of this procedure is presented in Fig. F1-b. Initially, a fast Fourier
transform is employed on the detected interference pattern with the primary
objective of extracting the frequency spectrum. Subsequently, a judicious
selection is made to isolate the pertinent frequency spectrum corresponding to
the LG beam, demarcated by the red square in the FFT image.