# Optimal Transport with Defective Cost Functions with Applications to the
Lens Refractor Problem

Axel G. R. Turnquist  Department of Mathematics, University of Texas at
Austin, Austin, TX, 78712<EMAIL_ADDRESS>
###### Abstract.

We define and discuss the properties of a class of cost functions on the
sphere which we term defective cost functions. We then discuss how to extend
these definitions and some properties to cost functions defined on Euclidean
space and on surfaces embedded in Euclidean space. Some important properties
of defective cost functions are that they result in Optimal Transport mappings
which map to points along geodesics, have a nonzero mixed Hessian term, among
other important properties. We also compute the cost-sectional curvature for a
broad class of cost functions, to verify and some known examples of cost
functions and easily prove positive cost-sectional curvature for some new cost
functions. Finally, we discuss how we can construct a regularity theory for
defective cost functions by satisfying the Ma-Trudinger-Wang (MTW) conditions
on an appropriately defined domain. As we develop the regularity theory of
defective cost functions, we discuss how the results apply to a particular
instance of the far-field lens refractor problem.

##  1\. Introduction

Recently, much work has been done on deriving PDE formulations of freeform
optics problems with lenses and reflectors, see for example the work in [19,
17, 18, 1, 7], resulting in Optimal Transport PDE and, more generally,
generated Jacobian equations posed on the unit sphere
$\mathbb{S}^{2}\subset\mathbb{R}^{3}$. The general problem is to find, given
source and target light intensities modeled as probability measures
$\mu\subset\mathbb{S}^{2}$ and $\nu\subset\mathbb{S}^{2}$, the function that
determines the shape of a lens (or reflector) system that achieves the desired
light intensity output pattern. The PDE formulations for some of these optics
problems have an important interpretation in terms of Optimal Transport. On
the theoretical side, some interest in the menagerie of Optimal Transport
problems that arise from these optics problems stems from the rather “exotic”
nature of these cost functions in their Optimal Transport interpretation. For
optical problems concerning redirecting directional light to a desired far-
field intensity, the it is natural for the PDE to be posed on the sphere. Some
cost functions on the sphere which appear “exotic” actually satisfy some very
natural conditions and thus lead to well-behaved Optimal Transport mappings. A
regularity theory has been built for “well-behaved” cost functions in Theorem
4.1 of Loeper [11]. The cost functions we address in this manuscript do not
fit into the theory established by Loeper, but instead exhibit “defects” that
do not allow for mass to be transported too far on the unit sphere, hence we
deem these cost functions “defective”. However, not all is lost, and defective
cost functions can satisfy the Ma-Trudinger-Wang (MTW) conditions (see [12])
on appropriate domains.

We will pay special attention to a problem known as the far-field lens
refractor problem. The setup is as follows: directional light gets redirected
through a refractive medium to a desired directional target intensity. The
refractive medium could have 1) refractive index less than or 2) greater than
the surrounding medium, resulting in different cost functions in their Optimal
Transport PDE formulation. These cases are very different on the PDE level,
mainly due to the fact that in one case the cost-sectional curvature is
positive and in the other, negative. However, both cost functions exhibit the
same problem in that they do not allow for mass to be transported too far.
Most of the early work on the refractor problem, especially proving the
existence of weak solutions of the problem and computing the cost-sectional
curvature was done in [7] and further related optical problems (that do not
have Optimal Transport PDE formulations) were developed by Gutierrez and
coauthors in a series of works, for example the near-field refractor problem
work in [6].

The goal of this manuscript, is twofold. First, and foremost, is a greater
exploration of the specific PDE arising in the far-field lens refractor
problem and its solvability and regularity, with an eye for designing future
numerical discretizations which ideally should take into account the
solvability condition and yield provably convergent solutions even when the
solution of the PDE are known to be “weak” viscosity solutions. Secondly, our
goal is to develop a theory defective cost functions to facilitate the
analysis of such Optimal Transport cost functions arising in optics problem.
We can characterize many cost functions by in which direction the mass
transports and how far the mass is “allowed” to travel given some structure on
the cost function. Defective cost functions will transport mass along
geodesics and will not allow mass to be transported beyond a certain distance.
Thus, this manuscript is also an exploration into the structure of such cost
functions and thus many of the results in this paper do not exclusively apply
to the cost function arising from the lens refractor problem, but apply in
greater generality. This allows for the construction of a regularity theory
that does not rely on arguments involving ellipsoids or hyperboloids of
revolution, which was the approach taken by [7] and [9], but rather follows in
the footsteps of the line of work by Loeper in [10], and especially the paper
[11], which allows for the theory to apply to a wide variety of cost
functions. As a byproduct of our work, for a large class of cost functions, we
have discovered more fundamental conditions which will allow for Theorem 4.1
of Loeper’s work [11] to apply, as opposed to posing the MTW conditions as an
assumption in that theorem. This allows, then, for a much easier way for a
researcher working on applications to verify the hypotheses of Theorem 4.1 of
Loeper for many cost functions, since the computations are done in full in
this manuscript. An example of using the theorems in this manuscript to verify
the MTW conditions is shown in Section 4.1.

In this paper, our goal is to analyze solvability and regularity of such PDE
for a class of cost functions which we will call defective cost functions. We
desire to prove that smooth enough data lead to smooth solutions. In Section
2, we review the setup of the lens refractor problem, which serves as an
example to motivate the types of cost functions we will be considering.
Section 3 contains the core computations and definitions of this paper. We
motivate the definition of defective cost functions via a computation that
shows a solvability condition in the lens refractor problem. We then show an
example where the lens refractor problem cannot be solved. We then introduce
the conditions required to have cost functions leading to exponential-type
maps for Optimal Transport problems on $\mathbb{S}^{2}\subset\mathbb{R}^{3}$,
$\mathbb{R}^{d}$ and, surfaces $M\subset\mathbb{R}^{3}$. We introduce a
function $\beta$, which characterizes how far the mapping moves along a
geodesic. We then introduce the definition of a defective cost function, which
is a cost function that changes concavity as a function of the Riemannian
distance. We show that this is the condition for the lens refractor problem
that leads to the solvability condition observed earlier. Using these
definitions and some useful formulas, we show how defective cost functions
have nonzero mixed Hessian determinant on the sphere $\mathbb{S}^{2}$ and
$\mathbb{R}^{d}$. We end this section with a computation of the cost-sectional
curvature for defective cost functions on $\mathbb{S}^{2}$ and
$\mathbb{R}^{d}$ and apply the results to some known examples, confirming
previous works in the literature. In Section 4, we first show how defective
cost functions satisfy the MTW conditions, with the caveat that they must
also, independently, satisfy the cost-sectional curvature condition. We then
show a computation of power cost functions on the unit sphere and how they
satisfy every condition in Theorem 4.1 of Loeper [11], except for the strict
cost-sectional curvature condition, unless we are dealing with the squared
geodesic cost. We then show how one can control the distance the Optimal
Transport mapping moves mass by imposing some controls on the source and
target density. This allows us to build a $C^{\infty}$ regularity theory for
defective cost functions following the outline of the regularity theory in
Loeper [11]. In Section 5 we summarize the contributions of the paper.

##  2\. Preliminaries

In this section, we first introduce the far-field lens refractor problem, in
order to introduce a known physical application where defective cost functions
arise. The definition of defective cost functions appears in Definition 7. We
then show a simple example where the lens refractor problem cannot be solved
and argue in simple terms why this should make sense given the physical setup
of the problem.

###  2.1. Far-Field Lens Refractor Problem

####  2.1.1. Interior and Exterior Media Separated By Outer Lens Edge

The far-field refractor problem is described in [7] as follows. Light
radiating from a point in an inner medium with refractive index $n_{I}$ goes
in the direction $x$ and then passes through a lens with refractive index
$n_{O}$, whose interior edge is a sphere and whose outer edge is described by
the function $\mathcal{L}=u_{1}(x)x$ and enters the outer medium also with
refractive index $n_{O}$ in the direction $y$. Light radiating from the origin
has a radial intensity pattern $f(x)$. After passing through the lens, the
far-field intensity pattern is given by $g(y)$, see Figure 1. The inverse
problem, as in Section LABEL:sec:lens1 is to solve for the function $u_{1}(x)$
given $f$ and $g$.

Figure 1. Light with radial intensity $f(x)$ travels in the direction $x$ in a
region of refractive index $n_{I}$, enters the lens with refractive index
$n_{O}$, whose inside edge is spherical and outer edge is given by the
function $\mathcal{L}=u_{1}(x)x$, then exits the lens into a region of
refractive index $n_{O}$ in the direction $y$ making a far-field intensity
pattern $g(y)$.

For reference, the index of refraction of a vacuum is $1$, water is $1.333$,
glass is $1.52$, and diamond is $2.417$. These representative values will be
used throughout the paper. The setup of the lens refractor problem allows for
us to define two quantities, $n=n_{I}/n_{O}$, and $\kappa=n_{O}/n_{I}$. We
choose to use different notation, since the PDE formulations of each problem
will be different. We will refer to the case where $n=n_{I}/n_{O}>1$ as lens
refractor problem I and the case with $\kappa>1$ as lens refractor problem II.

####  2.1.2. Potential and Cost Functions

From [14], the derivation of the PDE for the lens refractor problem I shows
that the cost function is of the following form:

(1) $c(x,y)=-\log\left(n-x\cdot y\right).$

Such a cost function may seem a bit unusual in Optimal Transport contexts, but
it very closely resembles the cost function for the perhaps better known
reflector antenna problem $c(x,y)=-\log(1-x\cdot y)$, see [17, 18]. Some clear
differences are, however, that the cost in (1) is Lipschitz, whereas the cost
$c(x,y)=-\log(1-x\cdot y)$ for the reflector antenna problem is not. Also,
when the magnitude of the gradient of the potential function is zero, the
mapping for (1) satisfies $T(x)=x$, whereas for the reflector antenna we get
$T(x)=-x$. It should become clear as we analyze the lens problem in greater
detail, then, that it does not make sense to take the limit as $n\rightarrow
1$ since we cannot talk about the lens problem “converging” in any sense to
the reflector problem. In fact, in some sense the reflector antenna cost is
less problematic, even though it is not Lipschitz, and the associated
regularity theory for the solution of the resulting Optimal Transport PDE was
addressed in [11].

From [7], the derivation of the PDE for the lens refractor problem II shows
that the cost function is of the following form:

(2) $c(x,y)=\log\left(\kappa x\cdot y-1\right).$

We can immediately notice a potential cause of concern. Since $\kappa>1$, it
is possible on the sphere to have $x\cdot y=1/\kappa$. That means that the
cost function in Equation (2) is not Lipschitz. Moreover, there is an issue if
in our Optimal Transport formulation, we have source and target densities
which require that mass transports further than a distance
$\arccos(1/\kappa)$. Surprisingly, as will be seen in Section 3.6, reflector
antenna problem II is actually the less problematic of the two, even though at
first glance its issues seem more apparent.

####  2.1.3. PDE Formulation

From the conservation of light intensity, for any subset
$A\subset\mathbb{S}^{2}$ where $T(A)\subset\mathbb{S}^{2}$ we get

(3) $\int_{A}f(x)dx=\int_{T(A)}g(y)dy.$

Solving for the mapping using the law of refraction (given the shape of the
outer edge of the lens $u_{1}(x)$) and using the change of variables
$u(x)=\log u_{1}(x)$ for the lens refractor problem I and the change of
variables $u(x)=-\log u_{1}(x)$ for the lens refractor problem II, the
following PDE is derived

(4)
$\det\left(D^{2}u(x)+D^{2}_{xx}c(x,y)|_{y=T(x)}\right)=\frac{\left|D^{2}_{xy}c(x,y)|_{y=T(x)}\right|f(x)}{g(T(x))},$

where

(5) $\nabla u(x)=-\nabla_{x}c(x,y)|_{y=T(x)}.$

##  3\. Computations, Defective Cost Functions, Mixed Hessian Determinants,
and Cost-Sectional Curvature

###  3.1. Explicit Form of the Mapping for the Lens Reflector Problem

For $x,y\in\mathbb{S}^{2}$, we initiate our study by performing an explicit
computation of the mapping in terms of the gradient of the potential function
$\nabla u$, which we label as $p\in T_{x}\mathbb{S}^{2}$. This will inform us
about when the PDE (4), and hence the far-field lens refractor problem, will
be unsolvable. We solve for the mapping $T=y$ via the equation
$\nabla_{\mathbb{S}^{2},x}c(x,y)=-p$. Let $\hat{q}=x\times\hat{p}$.

####  3.1.1. Computation of Mapping and Solvability Condition for the Lens
Refractor Problem I

For $c(x,y)=-\log\left(n-x\cdot y\right)$, we get, in local tangent
coordinates $(\hat{p},\hat{q})$:

(6) $\left(\left\|p\right\|,0\right)=\left(\frac{-y\cdot\hat{p}}{n-y\cdot
x},\frac{-y\cdot\hat{q}}{n-y\cdot x}\right).$

Since $n-x\cdot y\neq 0$, we find that $y\cdot\hat{q}=0$. Thus, along with the
constraint $(y\cdot x)^{2}+(y\cdot\hat{p})^{2}=1$, we arrive at the following
expression for the mapping $T$:

(7) $T(x,p)=x\
\frac{\sqrt{1+(1-n^{2})\left\|p\right\|^{2}}+n\left\|p\right\|^{2}}{1+\left\|p\right\|^{2}}+p\
\frac{\sqrt{1+(1-n^{2})\left\|p\right\|^{2}}-n}{1+\left\|p\right\|^{2}}.$

Note that $p=0\implies T(x,0)=x$. Since the shape of the reflector satisfies
$\nabla u(x)=\nabla u_{1}(x)e^{u_{1}(x)}$ then $p=\nabla u_{1}(x)=0$ implies
the lens is “flat” with respect to the canonical spherical metric, which, of
course, makes sense physically because where the lens is flat the light is not
redirected but simply passes straight through.

An important observation from Equation (7) is that $\left\|p\right\|$ cannot
exceed a certain threshold without making the mapping become complex valued.
This consequently puts a hard constraint on how far the mapping $T(x)$ can
move mass at a point $x$. We designate the quantity $p^{*}=1/\sqrt{n^{2}-1}$
and any associated vector with magnitude $p^{*}$ as $p^{*}\hat{p}$ and at this
value the square roots vanish in Equation (7). Therefore,

(8)
$T(x,p^{*}\hat{p})=\frac{n}{1+(p^{*})^{2}}\left(x(p^{*})^{2}-p^{*}\hat{p}\right),$

(9) $=\frac{x}{n}-\frac{\sqrt{n^{2}-1}}{n}\hat{p}.$

Thus,

(10) $T(x,p^{*}\hat{p})\cdot x=\frac{1}{n},$

and thus denoting $y=T$, we see that we get the following solvability
condition for the case when $n=n_{O}/n_{I}>1$:

(11) $x\cdot y\geq\frac{1}{n}.$

The requirement in Equation (11) has been noted in previous works, such as
[7]. Clearly, if this condition is not satisfied, we cannot find an Optimal
Transport mapping that solves the PDE (4) . From the results in [7], we
actually see that this solvability condition is both necessary and sufficient
for the existence of solutions of the PDE (4) (in a weak sense). In Section
3.4, we will show how we work the solvability condition into the definition of
a defective cost function.

####  3.1.2. Computation of Mapping and Solvability Condition for the Lens
Refractor Problem II

For the cost function $c(x,y)=\log\left(\kappa x\cdot y-1\right)$, we
immediately note that for the cost function to be evaluated we require that

(12) $x\cdot y\geq 1/\kappa.$

This condition also appears when we solve for the mapping $T$. Following the
same procedure as above, we derive the following expression for the mapping
$T$:

(13)
$T(x,p)=x\frac{\sqrt{\kappa^{2}+(\kappa^{2}-1)\left\|p\right\|^{2}}+\left\|p\right\|^{2}}{\kappa(1+\left\|p\right\|^{2})}+p\frac{\sqrt{\kappa^{2}+(\kappa^{2}-1)\left\|p\right\|^{2}}-1}{\kappa(1+\left\|p\right\|^{2})}.$

By taking the limit $\left\|p\right\|\rightarrow 0$, we get $T(x)=x$. By
taking the limit, $\left\|p\right\|\rightarrow\infty$, we get

(14)
$\lim_{\left\|p\right\|\rightarrow\infty}T(x,p)=\frac{1}{\kappa}x+\hat{p}\frac{\sqrt{\kappa^{2}-1}}{\kappa},$

and thus $x\cdot y=\frac{1}{\kappa}$. Thus, we do not have any restrictions on
the magnitude of $\left\|p\right\|$, but mass cannot be transported beyond a
certain distance no matter what the potential function looks like.

######  Remark 1.

What we will show later in this section is that the solvability condition from
Equation (11) and the solvability condition from Equation (12), while they
look very similar, arise from two different properties of the cost function
$c(x,y)$, which we will encapsulate in the definition of defective cost
function in Section 3.4. Both cost functions can be written as
$c(x,y)=G(d_{\mathbb{S}^{2}}(x,y))$. The condition from Equation (11) arises
because $G^{\prime\prime}(z)=0$ for $z=\arccos(1/n)$. The condition from
Equation (12) arises because $G^{\prime}(z)=\infty$ for $z=\arccos(1/\kappa)$.

###  3.2. A Simple Example Showing Ill-Posedness

It is now apparent that for the far-field refractor problem we must have that
the Optimal Transport mapping cannot move mass beyond a certain distance. This
puts a hard constraint on the allowable source and target masses. Here, we
give an example demonstrating a source and target intensities that do not
satisfy the solvability condition. Depicted schematically in Figure 2, no
single lens of can take the mass in red to the mass in blue, no matter how
large the refractive index is. Physically, this makes sense that it would be
unsolvable since the maximum angle at which the light rays can be deflected by
a lens is a right angle, by taking the limit in the refractor problem I, i.e.
$\lim_{n\rightarrow\infty}\arccos(1/n)=\pi/2$ from Equation (11) or the limit
in the refractor problem II, i.e.
$\lim_{\kappa\rightarrow\infty}\arccos(1/\kappa)=\pi/2$ from Equation (12).
However, the light intensity in Figure 2 must transport more than halfway
around the sphere.

Figure 2. Schematic of source (red) and target (blue) distributions, in which
the Optimal Transport PDE for any lens refractor problem is ill-posed with the
cost function (1)

Explicit examples for a continuous mapping $T$ which are unsolvable via a lens
system can easily be furnished. Suppose we fix $x_{0}\in\mathbb{S}^{2}$, then
for any $x\in\mathbb{S}^{2}$, we define the source density:

(15) $f(x)=\frac{1}{\beta}e^{\frac{-\frac{1}{2}d^{2}(x,x_{0})}{\sigma^{2}}},$

where $\beta$ is the normalization parameter:

(16) $\beta=2\pi\int_{0}^{\pi}e^{-\frac{1}{2}(\theta/\sigma)^{2}}\sin\theta
d\theta,$

and a target density given by:

(17) $g(x)=\frac{1}{\beta}e^{\frac{-\frac{1}{2}d^{2}(x,-x_{0})}{\sigma^{2}}}.$

Choose $\sigma$ to be very small, for example $\sigma=0.01$. Then,
$\beta\approx 0.0006283$. In order to solve the lens refractor problem, we
must satisfy the conservation of energy equation (3). We choose a Borel set
$B_{x_{0}}(5\sigma)$ to be the open geodesic ball of radius $5\sigma$ centered
around $x_{0}$. Then:

(18)
$\frac{2\pi}{\beta}\int_{0}^{0.05}e^{-\frac{1}{2}(\theta/\sigma)^{2}}\sin\theta
d\theta\approx 0.999992.$

That is $99.9992\%$ of the source mass is contained within a radius $0.05$ of
the point $x_{0}$. Likewise, for the target mass, we have $99.9992\%$ of the
target mass is contained within a radius $0.05$ of the point $-x_{0}$. Thus,
the condition that the mapping be mass preserving requires that the mapping
$T$ satisfy:

(19)
$0.999992\approx\int_{B_{x_{0}}(0.05)}f(x)dx=\int_{T(B_{x_{0}}(0.05))}g(y)dy.$

By the fact that the mapping $T$ is continuous, the non-empty open set
$B_{x_{0}}(0.05)$ is mapped to the non-empty open set $T(B_{x_{0}}(0.05))$.
Since only $0.0008\%$ of the mass of $g$ is located outside the ball
$B_{-x_{0}}(0.05)$, we must have that $T(B_{x_{0}}(0.05))\cap
B_{-x_{0}}(0.05)\neq\emptyset$ and the intersection is an open set. Thus, the
continuous mapping $T$ satisfies $d_{\mathbb{S}^{2}}(x,T(x))\geq\pi-0.1$ in a
non-negligible set and violates the fact that mass is not allowed to move too
far. The result from [7] shows, furthermore, that we cannot even expect weak
solutions if the mass distributions require mass to move too far.

###  3.3. Cost Functions Leading to Exponential-Type Maps

For the squared geodesic cost function
$c(x,y)=\frac{1}{2}d_{\mathbb{S}^{2}}(x,y)^{2}$, a classical result of McCann
[13] showed that the mapping is of the form $T(x,p)=\text{exp}_{x}(p)$. For
other cost functions, like the one presented in this paper,
$c(x,y)=-\log(n-x\cdot y)$, $n>1$, is it true that we can express the map as
$T(x,p)=\text{exp}_{x}(\hat{p}\beta(\left\|p\right\|))$ for some function
$\beta:\mathbb{R}^{+}\rightarrow\mathbb{R}$? This leads to the following
definition:

######  Definition 2.

We will say that a cost function leads to a map of exponential type if the
Optimal Transport mapping can be expressed in the following way:

(20) $T(x)=\text{exp}_{x}\left(\frac{\nabla u(x)}{\left\|\nabla
u(x)\right\|}\beta(\left\|\nabla u(x)\right\|)\right)$

where $u(x)$ is the potential function and $\beta$ is a continuous function on
its appropriate domain.

For lens refractor problems I and II, the cost functions lead to maps of
exponential type. Since, Equation (7) can be expressed in the following way:

(21) $T(x,p)=xR_{1}(\left\|p\right\|)+pR_{2}(\left\|p\right\|),$

with
$R_{1}(\left\|p\right\|)^{2}+\left\|p\right\|^{2}R_{2}(\left\|p\right\|)^{2}=1$,
we see that $T(x,p)$ is on the great circle passing through $x$ in the
direction $p$. By the formula in Equation (7), the length of the geodesic is
$\arccos
R_{1}(\left\|p\right\|)=\arccos\left(\frac{\sqrt{1+(1-n^{2})\left\|p\right\|^{2}}+n\gamma^{2}}{1+\gamma^{2}}\right)=\beta(\left\|p\right\|)$.
Thus, Equation (7) can be rewritten as:

(22)
$T(x,p)=\text{exp}_{x}\left(\hat{p}\arccos\left(\frac{\sqrt{1+(1-n^{2})\left\|p\right\|^{2}}+n\left\|p\right\|^{2}}{1+\left\|p\right\|^{2}}\right)\right).$

For the values $n=1.333$, $n=1.52$, and $n=2.417$, typical values for the
refractive indices of water, glass, and diamond, the $\beta$ functions are
plotted in Figure 3. Here we note that the “worst” behavior we can expect from
such functions is the appearance of a cusp, which is seen in the figure.
Staying away from $p^{*}$, the $\beta$ functions are actually Lipschitz, see
Lemma 8 below.

Figure 3. The function $\beta$ plotted as a function of $\left\|p\right\|$.
Black indicates $n=1.333$, approximately the refractive index of water, red
indicates $n=1.52$, approximately the refractive index of glass, and blue
indicates $n=2.417$, approximately the refractive index of diamond.

For the cost function $c(x,y)=\frac{1}{\kappa}\log(\kappa x\cdot y-1)$, we get
the $\beta$ function:

(23)
$\beta(\left\|p\right\|)=\arccos\left(\frac{\sqrt{\kappa^{2}+(\kappa^{2}-1)\left\|p\right\|^{2}}+\left\|p\right\|^{2}}{\kappa(1+\kappa^{2}\left\|p\right\|^{2})}\right).$

For the values $n=1.333$, $n=1.52$, and $n=2.417$, the $\beta$ functions are
plotted in Figure 4. The $\beta$ functions are monotonically increasing, but
flatten out as $\left\|p\right\|\rightarrow\infty$.

Figure 4. The function $\beta$ plotted as a function of $\left\|p\right\|$.
Black indicates $n=1.333$, red indicates $n=1.52$, and blue indicates
$n=2.417$.

Thus, we say that $c(x,y)=-\log(n-x\cdot y)$ and $c(x,y)=\log(\kappa x\cdot
y-1)$ are cost functions leading to maps of exponential type. The following
theorem shows when this is possible and when $\beta$ is monotone and satisfies
$\beta(0)=0$.

######  Theorem 3.

Suppose the cost function on the unit sphere can be written as a function of
the dot product between $x$ and $y$ and (consequently) the distance between
$x$ and $y$, i.e. $c(x,y)=F(x\cdot y)=G(d_{\mathbb{S}^{2}}(x,y))$ and that $G$
is twice differentiable and $G^{\prime\prime}(z)$ is strictly positive or
strictly negative on an open interval $(0,b)$ for $z\in(0,b)$ and $b\leq\pi$
and $F$ is differentiable and $F^{\prime}(\zeta)\neq 0$ for
$\zeta\in(\cos(b),1)$ and $\lim_{\zeta\rightarrow 1}F^{\prime}(\zeta)\neq 0$.
Then, the Optimal Transport mapping arising from this cost function can be
written as

(24) $T(x,p)=\text{exp}_{x}\left(\hat{p}\beta(\left\|p\right\|)\right),$

where $\beta$ is monotonically increasing and $\beta(0)=0$ and
$\left\|p\right\|\in[0,\left|G^{\prime}(b)\right|)$.

###### Proof.

Suppose that the cost function can be written as $c(x,y)=F(x\cdot y)$.
Equation (5), i.e. $p=-\nabla_{x}F(x,y)$, for some tangent vector $p$ implies
that we have $(\left\|p\right\|,0)=\left(-F^{\prime}(x\cdot
y)y\cdot\hat{p},-F^{\prime}(x\cdot y)y\cdot\hat{q}\right)$ expressed in the
tangent coordinates $(\hat{p},\hat{q})$ for $\hat{q}\in T_{x}\mathbb{S}^{2}$
where $\hat{q}\cdot\hat{p}=0$. This then immediately shows that
$0=y\cdot\hat{q}F^{\prime}(x\cdot y)$. Thus, if $F^{\prime}(x\cdot y)\neq 0$
for $x\cdot y\in(\cos b,1)$, then $y\cdot\hat{q}=0$ (when $x\cdot y=1$ we also
have $y\cdot\hat{q}=0$ by taking the limit) and therefore, the mapping does
not go in the $\hat{q}$ direction and therefore, $T$ is of the form
$T(x,p)=R_{1}x+R_{2}p$ and hence $T$ lies on the geodesic connecting $x$ to
$T$. Thus, the equation $\left\|p\right\|=-F^{\prime}(x\cdot y)y\cdot\hat{p}$
tells us that as $\left\|p\right\|\rightarrow 0$, we have either
$y\cdot\hat{p}\rightarrow 0$ or $F^{\prime}(x\cdot y)\rightarrow 0$. If
$y\cdot\hat{p}\rightarrow 0$, since $y$ is on the great circle including $x$
in the direction $p$, then either $x\cdot y\rightarrow 1$ or $x\cdot
y=\rightarrow-1$. Likewise, if $F^{\prime}(x\cdot y)\rightarrow 0$, then since
$F^{\prime}(x\cdot y)\neq 0$ for $x\cdot y\in(\cos b,1)$, then $x\cdot
y\rightarrow 1$. Therefore, as $\left\|p\right\|\rightarrow 0$, we find that
$\beta(0)=0$ and also $T$ can be written generally as:

(25) $T(x,p)=\text{exp}_{x}\left(\hat{p}\beta(\left\|p\right\|)\right),$

where $\beta(\left\|p\right\|)=\arccos R_{1}(\left\|p\right\|)$ and exists for
$\left\|p\right\|\in[0,\left|G^{\prime}(b)\right|)$. Furthermore, if
$F^{\prime}(x\cdot y)<0$ on $(\cos b,1)$, then $y\cdot\hat{p}>0$, and
likewise, if $F^{\prime}(x\cdot y)>0$ on $(\cos b,1)$, then $y\cdot\hat{p}<0$.

Since $\cos\left(d_{\mathbb{S}^{2}}(x,y)\right)=x\cdot y$, we have that the
cost function can be written as a function of the geodesic distance on the
sphere $c(x,y)=F(x\cdot y)=G(d_{\mathbb{S}^{2}}(x,y))$. Letting
$z=d_{\mathbb{S}^{2}}(x,y)$, we have that:

(26) $G^{\prime}(z)=-F^{\prime}(x\cdot y)\sin(z),$

since $\cos(z)=x\cdot y$. Furthermore, since $y=(x\cdot
y)x+(y\cdot\hat{p})\hat{p}$, we get that $(y\cdot\hat{p})^{2}=\sin^{2}(z)$. If
$F^{\prime}(x\cdot y)<0$, then $y\cdot\hat{p}=\sin(z)$ and if
$F^{\prime}(x\cdot y)>0$, then $y\cdot\hat{p}=-\sin(z)$.

Therefore, we find that

(27) $\displaystyle\left\|p\right\|$ $\displaystyle=-F^{\prime}(x\cdot
y)y\cdot\hat{p}=G^{\prime}(z),$ if $\displaystyle F^{\prime}(x\cdot y)<0$ (28)
$\displaystyle\left\|p\right\|$ $\displaystyle=-G^{\prime}(z),$ if
$\displaystyle F^{\prime}(x\cdot y)>0.$

Thus, we get that $\beta$ and $\pm G^{\prime}$ are inverses of each other,
depending on the sign of $F^{\prime}(x\cdot y)$. Taking a derivative, we get

(29) $\displaystyle\beta^{\prime}(\left\|p\right\|)$
$\displaystyle=1/G^{\prime\prime}(z),$ if $\displaystyle F^{\prime}(x\cdot
y)<0$ (30) $\displaystyle\beta^{\prime}(\left\|p\right\|)$
$\displaystyle=-1/G^{\prime\prime}(z),$ if $\displaystyle F^{\prime}(x\cdot
y)>0.$

Therefore, If $G^{\prime\prime}(z)$ is strictly positive or negative, we see
that $\beta(\left\|p\right\|)$ is strictly monotone and differentiable.

∎

######  Remark 4.

The preceding theorem can be modified to accommodate the case $\beta(0)=\pi$,
as is the case with the logarithmic cost function $c(x,y)=-\log(1-x\cdot y)$,
mutatis mutandis. However, for the sake of exposition, in this paper we focus
on the case where $\beta(0)=0$, that is, when the potential function is flat
with respect to the induced metric on the sphere, the Optimal Transport
mapping satisfies $T(x)=x$.

The proof of Theorem 3 yielded three important relations. The first comes from
the definition of $\beta$ and the exponential map and by denoting
$z=d_{\mathbb{S}^{2}}(x,y)$:

(31) $z=\beta(\left\|p\right\|).$

The second relation is contained in Equations (27) and (28) and the third is
contained in Equations (29) and (30). Thus, we see that $\beta$ and $G$ are
related very closely.

We saw that for the cost functions for the far-field lens refractor problem
led to maps of exponential type. It should also be remarked here that the
squared geodesic cost function
$c(x,y)=d_{\mathbb{S}^{2}}(x,y)^{2}=\arccos(x\cdot y)^{2}$ and the cost
function arising from the reflector antenna:
$c(x,y)=-\log\left\|x-y\right\|=-\log(1-x\cdot y)$ lead to maps of exponential
type, where $\beta(\left\|p\right\|)$ is monotonically increasing and
monotonically decreasing, respectively. For a real-world example of a cost
function on the sphere which is does not lead to a map of exponential type, we
refer to the point-to-point reflector problem, whose setup and derivation are
contained in [19]. The cost function is:

(32) $c(x,y)=\log\left(\frac{1}{(T^{2}-L^{2})^{2}}-\frac{1-x\cdot
y}{2(T^{2}-L^{2})(T-L(x\cdot\hat{e}))(T-L(y\cdot\hat{e}))}\right),\ \
x,y\in\mathbb{S}^{2},$

where $T,L$ are parameters arising in the optical setup satisfying $T\geq L$.
$T$ is the length of the optical path and $L$ is the distance between the
points. There is a special direction in the problem which is the direction
from the source point to the target point, denoted by $\hat{e}$. Due to the
presence of the term $x\cdot\hat{e}$ arising from this special direction
$\hat{e}$, the cost function cannot be written as simply a function of $x\cdot
y$. For this reason, in computing the mapping for this kind of “exotic” cost
function, we observe that $T$ has the more general form:

(33) $T(x,p)=R_{1}x+R_{2}p+R_{3}q,$

where w.l.o.g. $q=x\times p$. The presence of the nonzero $R_{3}$ term due to
the cost (32) indicates that $T(x,p)$ does not travel along the great circle
defined by $x$ and $p$. Thus, we have an example of a cost function which
leads to a mapping $T$ which is not of exponential type.

######  Remark 5.

As we just saw, for the sphere, cost functions that are functions of the
Riemannian distance between $x$ and $y$ can be written as functions of the dot
product $x\cdot y$ and this simplifies the exposition and computations
greatly. This is the only manifold where this is true. The natural extension,
then, of costs leading to maps of exponential type in Euclidean space
$c(x,y):\Omega\times\Omega^{\prime}\rightarrow\mathbb{R}$, where
$\Omega,\Omega^{\prime}\subset\mathbb{R}^{n}$ are of the form
$c(x,y)=H(d(x,y))$, where $d(x,y)=\left\|x-y\right\|$. We may write such cost
functions as $c(x,y)=H(d(x,y))=J(\frac{1}{2}d(x,y)^{2})$ (where now the
function $J$ will play the role of the earlier function $F$) and thus,
$p=-\nabla_{x}c(x,y)$ results in the following equation:

(34) $p=J^{\prime}\left(\frac{1}{2}d(x,y)^{2}\right)(y-x),$

which indicates that $y-x=\frac{p}{J^{\prime}}$, which we see implies that the
mapping is along a geodesic from $x$ in the direction of $p$. The
generalization, then, of the condition in Theorem 3 that
$F^{\prime}(\zeta)\neq 0$ is given in the Euclidean case by
$J^{\prime}(\zeta)\neq 0$. Furthermore, we can find an expression for
$\beta^{\prime}\left(\left\|p\right\|\right)$. By taking the magnitude of
Equation (34) and using the definition of $J$ and $H$, we get:

(35)
$\left\|p\right\|=J^{\prime}\left(\frac{1}{2}d(x,y)^{2}\right)d(x,y)=H^{\prime}(d(x,y)),$

and thus, by letting $z=d(x,y)$ and assuming that $H^{\prime\prime}(z)\neq 0$,
we can invert this and also show that:

(36) $\beta^{\prime}(\left\|p\right\|)=\frac{1}{H^{\prime\prime}(z)},$

which is a result that held in the case of the sphere as well, see Equation
(29). This shows that the condition $H^{\prime\prime}(z)\neq 0$ is also
important in the Euclidean case. Furthermore, we have
$z=\beta(\left\|p\right\|)$ and $\left\|p\right\|=H^{\prime}(z)$ so $\beta$
and $H^{\prime}$ are inverses.

######  Remark 6.

The further generalization to Riemannian manifolds can be given as follows. If
$c(x,y)=H(d_{M}(x,y))=J(\frac{1}{2}d_{M}(x,y)^{2})$, then since

(37)
$p=-\nabla_{x}c(x,y)=-J^{\prime}\left(\frac{1}{2}d_{M}(x,y)^{2}\right)\nabla_{x}\left(\frac{1}{2}d_{M}(x,y)^{2}\right),$

we get:

(38)
$\frac{p}{J^{\prime}\left(\frac{1}{2}d_{M}(x,y)^{2}\right)}=-\nabla_{x}\left(\frac{1}{2}d_{M}(x,y)^{2}\right),$

and thus, by the classical results of McCann [13], we get

(39)
$y=\text{exp}_{x}\left(\frac{p}{J^{\prime}\left(\frac{1}{2}d_{M}(x,y)^{2}\right)}\right).$

Thus, we see generally that if the cost function can be written as a function
of the Riemannian distance on a manifold $M$, we get maps of exponential type.
The generalization, then, of the condition in Theorem 3 that
$F^{\prime}(\zeta)\neq 0$ is given in the Riemannian case by
$J^{\prime}(\zeta)\neq 0$. Since the distance from $x$ to $y$ is given by the
magnitude of the tangent vector in the argument of the exponential map, we
can, as in the Euclidean case, derive the relation:

(40)
$\left\|p\right\|=J^{\prime}\left(\frac{1}{2}d_{M}(x,y)^{2}\right)d_{M}(x,y)=H^{\prime}(d_{M}(x,y)),$

and thus, assuming $H^{\prime\prime}(z)\neq 0$, we derive Equation (36) for
the Riemannian manifold case. As in the Euclidean case, we have
$z=\beta(\left\|p\right\|)$ and $\left\|p\right\|=H^{\prime}(z)$ so $\beta$
and $H^{\prime}$ are inverses of each other.

###  3.4. Defective Cost Functions

For some cost functions of exponential type, the mass is unable to transport
beyond a certain distance. We observed this with the cost functions
$c(x,y)=-\log(n-x\cdot y)$ and $c(x,y)=\log(\kappa x\cdot y-1)$ from the far-
field refractor problem. Now, we identify properties of the cost function
which lead to this phenomenon and call such cost functions defective.

######  Definition 7.

A cost function $c(x,y):\mathbb{S}^{2}\times\mathbb{S}^{2}\rightarrow 0$,
which can be written $c(x,y)=F(x\cdot y)=G(d_{\mathbb{S}^{2}}(x,y))$, is
defective if $G^{\prime\prime}(z)=0$ or $G^{\prime}(z)=\infty$ for some
$z\in(0,\pi)$ and $G^{\prime\prime}(z)$ is either strictly positive or
negative for $z\in[0,z^{*})$, where $z^{*}\in(0,\pi)$ is the smallest value
for which $G^{\prime\prime}(z)=0$ or $G^{\prime}(z)=\infty$ and
$F^{\prime}(\zeta)\neq 0$ on the interval $(\cos z^{*},1)$ with
$\lim_{\zeta\rightarrow 1}F^{\prime}(\zeta)\neq 0$. If necessary in the
discussion, we will refer to cost functions for which
$G^{\prime\prime}(z^{*})=0$ as defective cost functions of the first type, and
cost functions for which $G^{\prime}(z^{*})=\infty$ as defective cost
functions of the second type. We also denote $p^{*}$ as the value
$p^{*}=\left|G^{\prime}(z^{*})\right|$.

Defective cost functions lead to maps of exponential type provided that the
distance transported is less than $z^{*}$. By the definition, defective cost
functions of the first type are those for which the concavity of the cost
function as a function of the Riemannian distance on $\mathbb{S}^{2}$ changes
for some $z\in(0,\pi)$. This change of concavity of the function was a noted
issue going back to work by Gangbo and McCann [5].

Defective cost functions of the second type are where $G^{\prime}(z)$ diverges
to infinity as $\left\|p\right\|\rightarrow\infty$. This implies by
differentiating, that $G^{\prime\prime}(z^{*})=\infty$. Then, the equations
$G^{\prime}(z)=\left\|p\right\|$ and
$G^{\prime\prime}(z)=1/\beta^{\prime}(\left\|p\right\|)$ imply that since
$G^{\prime\prime}(z)\neq 0$, we can invert $\left\|p\right\|=G^{\prime}(z)$ to
get $z=\beta(\left\|p\right\|)$ and thus
$\lim_{\left\|p\right\|\rightarrow\infty}\beta(\left\|p\right\|)=z^{*}$ and
$\lim_{\left\|p\right\|\rightarrow\infty}\beta^{\prime}(\left\|p\right\|)=0$.
So, we see that defective cost functions $c(x,y)=G(z)$ restricted to the
domain $z\in[0,z^{*})$ are “well-behaved”, but do not allow for mass to be
transported beyond a distance $z^{*}$.

The cost functions $c(x,y)=d_{\mathbb{S}^{2}}(x,y)^{2}$ and
$c(x,y)=-\log(1-x\cdot y)$ are not defective costs, whereas
$c(x,y)=-\log(n-x\cdot y)$ is a defective cost function of the first type
because $G^{\prime\prime}(z^{*})=0$, for $z^{*}=\arccos(1/n)$. The cost
function $c(x,y)=\frac{1}{\kappa}\log(\kappa x\cdot y-1)$ is a defective cost
function of the second type, because
$G^{\prime\prime}(\arccos(1/\kappa))=\infty$. For the cost
$c(x,y)=-\log(n-x\cdot y)=-\log(n-\cos(z))$, we get that
$G^{\prime\prime}(z)=0$ when $z=\arccos(1/n)$. From Equation (11), we see that
this value of $z$ is exactly the farthest the mapping may transport while
staying real-valued. For the cost function $c(x,y)=-\log(n-x\cdot
y)+\log(n+1)$, we have plotted the cost function as a function of geodesic
distance and denoted the inflection points $G^{\prime\prime}(z)=0$ with
circles, see Figure 5.

Figure 5. Cost function $c(x,y)=-\log(n-\cos(z))+\log(n+1)$ plotted as a
function of $z$. Black indicates $n=1.333$, approximately the refractive index
of water, red indicates $n=1.52$, approximately the refractive index of glass,
and blue indicates $n=2.417$, approximately the refractive index of diamond.
Below the inflection points, the cost functions are concave.

A defective cost function automatically fails to satisfy the injectivity for
each $x\in\mathbb{S}^{2}$ of the map $y\rightarrow-\nabla_{x}c(x,y)$ over the
set $\mathbb{S}^{2}\setminus\left\\{-x\right\\}$, which is one of the MTW
conditions presented in [11]. It should be noted also that the cost functions
in Theorem 4.1 of [11] lead to maps of exponential type that are not
defective.

For $\left\|p\right\|\leq c<p^{*}$, we also have the important fact that the
$\beta$ functions, defined in Definition 2 are Lipschitz for defective cost
functions below the value $p^{*}$.

######  Lemma 8.

Suppose that there exists a $z^{*}$ such that $G^{\prime\prime}(z^{*})=0$.
Denote by $p^{*}$ the value such that $\beta(p^{*})=z^{*}$. The $\beta$
function arising from a defective cost function $c(x,y)=F(x\cdot
y)=G\left(d_{\mathbb{S}^{2}}(x,y)\right)$, where $F^{\prime}(\zeta)\neq 0$, is
Lipschitz on the interval $[0,c]$, where $c<p^{*}$.

###### Proof.

This is an immediate consequence of Equation (29) or Equation (30). For
$z<z^{*}$, we see that the derivative of $\beta$ is bounded. ∎

The important relation in Equations (29) and (30) explains why the $\beta$
function for the lens refractor problem I has a cusp and why the $\beta$
function for the lens refractor problem II flattens out. Therefore, if we have
a $\beta$ function with a cusp, then it arises from a defective cost function
of the first type, and if the $\beta$ function limits to a horizontal
asymptote, where the asymptote is at a value strictly less than $\pi$, then it
arises from a defective cost function of the second type. Both are types of
cost functions which do not allow the mass to be transported a distance beyond
a fixed value strictly less than $\pi$.

######  Remark 9.

It should be clear how to extend the definition of a defective cost function
in Euclidean space. If we have a cost function
$c(x,y)=H(\left\|x-y\right\|)=J\left(\frac{1}{2}\left\|x-y\right\|^{2}\right)$,
then denote $z^{*}$ to be the smallest value for which $H^{\prime\prime}(z)=0$
or $H^{\prime}(z)=\infty$. Then, let $F$ satisfy $F^{\prime}(\zeta)\neq 0$ for
$\zeta<\sqrt{2z^{*}}$ and $H^{\prime\prime}(z)$ is either positive or negative
for $z\in[0,z^{*})$. We call such a $c$ a defective cost function.

######  Remark 10.

In the case of more general Riemannian manifolds, if we have a cost function
$c(x,y)=H(d_{M}(x,y))=J\left(\frac{1}{2}d_{M}(x,y)^{2}\right)$, then denote
$z^{*}$ to be the smallest value for which $H^{\prime\prime}(z)=0$ or
$H^{\prime}(z)=\infty$. Then, let $F$ satisfy $F^{\prime}(\zeta)\neq 0$ for
$\zeta<\sqrt{2z^{*}}$ and $H^{\prime\prime}(z)$ is either positive or negative
for $z\in[0,z^{*})$. We call such a $c$ a defective cost function. Note that
in this case, $z^{*}$ will necessarily be less than the injectivity radius.

####  3.4.1. Using Euclidean Distance for Computations on Manifolds

It is perhaps of interest to the computational Optimal Transport community to
investigate the possibility of using the ambient Euclidean distance in a cost
function, i.e. $c(x,y)=\frac{1}{2}d_{\mathbb{R}^{3}}(x,y)^{2}$ to solve the
Optimal Transport problem on a connected, orientable surface
$M\subset\mathbb{R}^{3}$. It is tempting to use the Euclidean distance is very
cheap to compute, whereas the cost function
$c(x,y)=\frac{1}{2}d_{M}(x,y)^{2}$, by contrast, is very expensive, since
pairwise geodesic distances on $M$ must be computed for every point. Here, we
demonstrate that the cost function
$c(x,y)=\frac{1}{2}d_{\mathbb{R}^{3}}(x,y)^{2}=H(d_{M}(x,y))$, while monotonic
(as a function of the Riemannian distance) for many manifolds, is defective
for $M=\mathbb{S}^{2}$. We have the cost function:

(41)
$c(x,y)=\frac{1}{2}d_{\mathbb{R}^{3}}(x,y)^{2}=2\sin^{2}\left(\frac{d_{M}(x,y)}{2}\right),$

so $G(z)=2\sin^{2}\left(\frac{z}{2}\right)$ and thus
$G^{\prime\prime}(z)=\cos(z)=0$ for $z=\pi/2$. Therefore, the cost function
$c(x,y)=\frac{1}{2}d_{\mathbb{R}^{3}}(x,y)^{2}$ is defective. We see what
happens when we solve for the mapping. Since for this cost function we have:

(42) $p=-\nabla_{x}c(x,y)=y-x,$

Taking the dot product of both sides with respect to $p$ and $q$,
respectively, yields:

(43) $\displaystyle\left\|p\right\|$ $\displaystyle=y\cdot p,$ (44)
$\displaystyle 0$ $\displaystyle=y\cdot q.$

Using the fact that $y=(y\cdot p)\hat{p}+(y\cdot x)x$, we get:

(45) $y=p+\sqrt{1-\left\|p\right\|^{2}}x.$

We thus we that $\left\|p\right\|\leq 1$. The case $\left\|p\right\|=1$
corresponds to $y=p$, which is the case where
$d_{\mathbb{S}^{2}}(x,y)=\frac{\pi}{2}$. Thus, we see that for this cost
function, mass is not allowed to move more than a distance of $\pi/2$.

We can also demonstrate the problem with some other surfaces, where the
computation of the Riemannian distance on the manifold is not trivial
computationally, for example the case of the $2$-torus
$\mathbb{T}^{2}\subset\mathbb{R}^{3}$ equipped with the standard Riemannian
metric. While the cost function
$c(x,y)=\frac{1}{2}d_{\mathbb{R}^{3}}(x,y)^{2}$ cannot be strictly written as
a simply a function of the Riemannian distance on the manifold $M$, it shares
the same property of defective cost functions in that it cannot transport mass
beyond a certain distance. Although we have the $2$-torus in mind, it should
be very clear how this argument extends to other surfaces and manifolds
embedded in higher dimensions. The mapping satisfies:

(46) $p=-\nabla_{x}c(x,y)=y-x.$

Therefore, taking the dot product with respect to $p$ and $q$, respectively,
we get:

(47) $\displaystyle\left\|p\right\|$ $\displaystyle=y\cdot p-x\cdot p,$ (48)
$\displaystyle 0$ $\displaystyle=y\cdot q-x\cdot q.$

Therefore, we have:

(49) $y=\left(\left\|p\right\|+x\cdot p\right)\hat{p}+\left(x\cdot
q\right)\hat{q}+(y\cdot\hat{n})\hat{n}$

Thus, the resulting mapping arising from Optimal Transport on the torus with
the squared Euclidean cost satisfies:

(50) $y=x+p+(y\cdot\hat{n}-x\cdot\hat{n})\hat{n}.$

In general, directly solving for $y\cdot\hat{n}$ yields a very complicated
formula even when a level-set formulation of the surface is known. For the
$2$-torus, the computation of $y\cdot\hat{n}$ in general requires one to solve
a quartic equation in $y\cdot\hat{n}$. However, a simple geometric argument
renders the need to compute $y\cdot\hat{n}$ unnecessary. The problem of
solving for the mapping $y$ is equivalent to the problem of finding a
particular point of intersection of the parametrized line
$y=x+p+\alpha\hat{n}$ and the torus, since $y\in M$. In general, we can easily
find $p$ which cannot lead to an Optimal Transport mapping as follows. Take
the tangent plane $\mathcal{T}_{x}$ and orthogonally project the set
$\mathbb{T}^{2}\subset\mathbb{R}^{3}$ onto $\mathcal{T}_{x}$. That is, define
$S:\mathbb{R}^{3}\rightarrow\mathcal{T}_{x}$ by
$S(z)=(1-\hat{n}\hat{n}^{T})z$. Then, denote $T^{\prime}=S(\mathbb{T}^{2})$.
For a fixed $x\in M$, a value of $p$ that is not allowed is one for which
$x+p\notin T^{\prime}$.

Even more generally, if $c(x,y)=f(d^{2}_{\mathbb{R}^{3}}(x,y))$, then the
problem of finding $y$, then we have:

(51) $p=-\nabla_{x}c(x,y)=f^{\prime}(d^{2}_{\mathbb{R}^{3}}(x,y))(y-x).$

Note that, by taking the magnitude of both sides, we get:

(52)
$\frac{\left\|p\right\|}{f^{\prime}(d^{2}_{\mathbb{R}^{3}}(x,y))}=d_{\mathbb{R}^{3}}(x,y).$

Therefore, taking the dot product of Equation (51) with respect to $p$ and
$q$, respectively, we get:

(53)
$\displaystyle\frac{\left\|p\right\|}{f^{\prime}(d^{2}_{\mathbb{R}^{3}}(x,y))}$
$\displaystyle=y\cdot p-x\cdot p,$ (54) $\displaystyle 0$
$\displaystyle=y\cdot q-x\cdot q.$

Therefore, we have:

(55) $\displaystyle y$
$\displaystyle=\left(\frac{\left\|p\right\|}{f^{\prime}(d^{2}_{\mathbb{R}^{3}}(x,y))}+x\cdot
p\right)\hat{p}+\left(x\cdot q\right)\hat{q}+(y\cdot\hat{n})\hat{n}$ (56)
$\displaystyle y$ $\displaystyle=\left(\beta(\left\|p\right\|)+x\cdot
p\right)\hat{p}+\left(x\cdot q\right)\hat{q}+(y\cdot\hat{n})\hat{n}$

Thus, the resulting mapping arising from Optimal Transport on the torus with
the squared Euclidean cost satisfies:

(57)
$y=\frac{\beta(\left\|p\right\|)}{\left\|p\right\|}p+x+(y\cdot\hat{n}-x\cdot\hat{n})\hat{n},$

where $z=\beta(\left\|p\right\|)$ which is found by inverting the equation
$\left\|p\right\|=zf^{\prime}(z^{2})$, when it is invertible.

Thus, finding the mapping $y$ is equivalent to the problem of finding the
closest point (to $x$) of intersection of the line
$y=\frac{\beta(\left\|p\right\|)}{\left\|p\right\|}p+x+\alpha\hat{n}$ with the
torus. Thus, we can see that no matter how $\beta$ increases as a function of
$\left\|p\right\|$, it is not possible for the mapping to reach around to the
“other side” of the torus. Thus, we see that all such cost functions built
from Euclidean distances are defective. This, of course, does not apply to the
case of the sphere (on which the logarithmic cost function $c(x,y)=-2\log
d_{\mathbb{R}^{3}}(x,y)$, for example, is not defective), since in that case,
$x\cdot p=x\cdot q=0$.

In the case where, a priori, the mapping is known to be close to identity, a
Euclidean cost can be used as an expedient for rapid computations, for example
using the Sinkhorn iterations proposed by Cuturi [2]. It should be emphasized
that this is not simply a claim that asymptotically the Euclidean distance
approximates the Riemannian distance. Rather, the Euclidean distance can be
used as a substitute distance up to a certain point.

###  3.5. Derivation of the Mixed Hessian Term
$\left|D^{2}_{xy}c(x,y)|_{y=T(x)}\right|$

The mixed Hessian term $\left|D^{2}_{xy}c(x,y)|_{y=T(x)}\right|$ appearing in
Equation (4) is important to compute to check the MTW conditions as well as
for numerical discretizations. See, for example, the paper [8] for an example
of a numerical discretization that uses the explicit form of the mixed
Hessian. In this subsection, we derive various expressions for computing the
mixed Hessian for cost functions leading to maps of exponential type. The
expressions we derive, which are Equations (67), (71) and (75), are for
mappings of exponential type, which arise from cost functions depending solely
on terms involving the dot product $x\cdot y$, see Section 3.3, however, the
derivation can be easily generalized to other, more general, cost functions
also arising in optics applications and will be explored in future work.

Here we derive a formula for the mixed Hessian via an integral formulation.
Using the notation established in Equation (21), define
$R(p)=pR_{2}(\left\|p\right\|)$. Then, for a region $E\subset
T_{x}\mathbb{S}^{2}$ define $U=R(E)=\left\\{pR_{2}(\left\|p\right\|)|p\in
E\right\\}$ and unit vectors $\hat{u}_{1},\hat{u}_{2}$, where
$\hat{u}_{1}=\hat{p}$ and $\hat{u}_{2}=x\times\hat{p}$. For any vector $v\in
T_{x}\mathbb{S}^{2}$, this defines a coordinate system on
$T_{x}\mathbb{S}^{2}$, i.e. $v=(u_{1},u_{2})$, where $u_{1}=v\cdot\hat{u}_{1}$
and $u_{2}=v\cdot\hat{u}_{2}$. Then, define the region
$T(x,E)=\left\\{\sqrt{1-u_{1}^{2}-u_{2}^{2}}|(u_{1},u_{2})\in
U\right\\}\subset\mathbb{S}^{2}$, see Figure 6.

Figure 6. Change in area formula from tangent coordinates $p$ to coordinates
on the sphere $T(x,p)=xR_{1}(\left\|p\right\|)+pR_{2}(\left\|p\right\|)$ via
the coordinates $(u_{1},u_{2})$ of the orthogonal projection of $T(x,p)$ onto
the tangent plane $T_{x}\mathbb{S}^{2}$.

Since the mixed Hessian satisfies the formula

(58) $\left|\det D^{2}_{xy}c\right|=\frac{1}{\left|\det D_{p}T\right|},$

we see that if we compute the quantity $\left|\det D_{p}T\right|$, then we can
compute the determinant of the mixed Hessian. We see that the quantity
$\left|\det D_{p}T\right|$ can be computed via computing the area of the
region $T(x,E)$. From Figure 6, we see:

(59) $\int_{T(x,E)}dS=\int_{E}\left|\det D_{p}T\right|dp.$

Also,

(60) $\int_{T(x,E)}dS=\int_{U}\frac{1}{\sqrt{1-\left\|u\right\|^{2}}}du,$

(61) $=\int_{E}\frac{1}{\sqrt{1-\left|R(p)\right|^{2}}}\left|\det\nabla
R(p)\right|dp.$

And therefore, by Equations (58), (59) and (61), we get the following
expression for the mixed Hessian:

(62) $\left|\det
D^{2}_{xy}c\right|=\frac{\sqrt{1-\left|R(p)\right|^{2}}}{\left|\det\nabla
R(p)\right|}.$

Since $R(p)=pR_{2}(\left\|p\right\|)$, we compute:

(63) $\displaystyle(\nabla R)_{11}$
$\displaystyle=\frac{p_{1}^{2}}{\left\|p\right\|}R_{2}^{\prime}(\left\|p\right\|)+R_{2}(\left\|p\right\|),$
(64) $\displaystyle(\nabla R)_{12}$ $\displaystyle=(\nabla
R)_{21}=\frac{p_{1}p_{2}}{\left\|p\right\|}R_{2}^{\prime}(\left\|p\right\|),$
(65) $\displaystyle(\nabla R)_{22}$
$\displaystyle=\frac{p_{2}^{2}}{\left\|p\right\|}R_{2}^{\prime}(\left\|p\right\|)+R_{2}(\left\|p\right\|),$

and thus,

(66) $\det\left(\nabla
R\right)=\left\|p\right\|R_{2}(\left\|p\right\|)R_{2}^{\prime}(\left\|p\right\|)+R_{2}(\left\|p\right\|)^{2},$

and hence

(67) $\left|\det
D^{2}_{xy}c\right|=\frac{\sqrt{1-\left\|p\right\|^{2}R_{2}(\left\|p\right\|)^{2}}}{\left|\left\|p\right\|R_{2}(\left\|p\right\|)R_{2}^{\prime}(\left\|p\right\|)+R_{2}(\left\|p\right\|)^{2}\right|}.$

We find that there is a more convenient expression for the determinant of the
mixed Hessian in terms of $R_{1}(\left\|p\right\|)$. First, we use the fact
that:

(68)
$R_{1}(\left\|p\right\|)^{2}+\left\|p\right\|^{2}R_{2}(\left\|p\right\|)^{2}=1,$

and thus, we get
$\left|R_{1}(\left\|p\right\|)\right|=\sqrt{1-\left\|p\right\|^{2}R_{2}(\left\|p\right\|)^{2}}$
and, also,

(69)
$R_{1}(\left\|p\right\|)R_{1}^{\prime}(\left\|p\right\|)+\left\|p\right\|R_{2}(\left\|p\right\|)^{2}+\left\|p\right\|^{2}R_{2}(\left\|p\right\|)R_{2}^{\prime}(\left\|p\right\|)=0,$

and thus,

(70)
$\left|\left\|p\right\|R_{2}(\left\|p\right\|)R_{2}^{\prime}(\left\|p\right\|)+R_{2}(\left\|p\right\|)^{2}\right|=\frac{\left|R_{1}(\left\|p\right\|)R_{1}^{\prime}(\left\|p\right\|)\right|}{\left\|p\right\|}.$

Thus, in terms of $R_{1}$ we get the very simple expression:

(71) $\left|\det
D^{2}_{xy}c\right|=\frac{\left\|p\right\|}{\left|R_{1}^{\prime}(\left\|p\right\|)\right|}.$

Using Equation (67), we derive a formula for the mixed Hessian term in terms
of the function $\beta$ and $F^{\prime}(x\cdot y)$ and
$G^{\prime\prime}(d_{\mathbb{S}^{2}}(x,y))$.

From the definition $T(x,p)=xR_{1}(\left\|p\right\|)+pR_{2}(\left\|p\right\|)$
and the fact that $\beta(\left\|p\right\|)=\arccos(R_{1})$, we get that
$\beta(\left\|p\right\|)=\arcsin(\left\|p\right\|R_{2}(\left\|p\right\|))$.
Therefore, since

(72)
$\beta^{\prime}(\left\|p\right\|)=\frac{R_{2}(\left\|p\right\|)+xR_{2}^{\prime}(\left\|p\right\|)}{\sqrt{1-\left\|p\right\|^{2}R_{2}(\left\|p\right\|)^{2}}},$

we get,

(73)
$\left|\beta^{\prime}(\left\|p\right\|)\right|=\frac{1}{\left|R_{2}(\left\|p\right\|)\right|}\left|D_{p}T\right|,$

from (67). Therefore,

(74)
$\left|D_{p}T\right|=\left|R_{2}(\left\|p\right\|)\right|\left|\beta^{\prime}(\left\|p\right\|)\right|=\frac{\sin\beta(\left\|p\right\|)}{\left\|p\right\|}\left|\beta^{\prime}(\left\|p\right\|)\right|,$

and thus, we have the following expression for the mixed Hessian in terms of
the $\beta$ function:

(75)
$\left|D^{2}_{xy}c(x,y)\right|=\frac{\left\|p\right\|}{\sin\left(\beta(\left\|p\right\|)\right)\left|\beta^{\prime}(\left\|p\right\|)\right|}$

This expression, while cumbersome to compute by hand, proffers a much simpler
route for checking the MTW conditions for all defective cost functions.
Provided that we satisfy $z\in[0,z^{*})$, the mixed Hessian is strictly
positive.

######  Lemma 11.

For a defective cost function $c(x,y)$, when $z$ satisfies $z<z^{*}$, we have

(76) $\left|D^{2}_{xy}c(x,y)\right|>0.$

###### Proof.

Let $z=d_{\mathbb{S}^{2}}(x,y)$. Then, we derive the important equality:

(77)
$\left|D^{2}_{xy}c(x,y)\right|=\frac{\left\|p\right\|}{\sin\left(\beta(\left\|p\right\|)\right)\left|\beta^{\prime}(\left\|p\right\|)\right|}=\frac{\left|G^{\prime\prime}(z)\right|}{\left|R_{2}(\left\|p\right\|)\right|}=\left|G^{\prime\prime}(z)\right|\left|F^{\prime}(x\cdot
y)\right|,$

since $y\cdot\hat{p}=\frac{\left\|p\right\|}{F^{\prime}(x\cdot
y)}=\left\|p\right\|R_{2}\left(\left\|p\right\|\right)$, therefore,
$R_{2}(\left\|p\right\|)=\frac{1}{F^{\prime}(x\cdot y)}$. Since
$\frac{1}{\left|\beta^{\prime}(\left\|p\right\|)\right|}=\left|G^{\prime\prime}(z)\right|\neq
0$ for $z\in[0,z^{*})$, and $F^{\prime}(\zeta)\neq 0$ for $\zeta\in(\cos
z^{*},1)$ and $\lim_{\zeta\rightarrow 1}F^{\prime}(\zeta)\neq 0$, we get that
$\left|D^{2}_{xy}c(x,y)\right|>0$. ∎

######  Remark 12.

A similar computation can be done for the Euclidean case for defective cost
functions, see Remark 5 for the background and assumptions, i.e.
$J^{\prime}(\zeta)\neq 0$ and $H^{\prime\prime}(z)\neq 0$. It can be shown
that for the $\mathbb{R}^{d}$ for such cost functions, we get:

(78)
$\left|D_{p}T\right|=\beta^{\prime}(\left\|p\right\|)\left(\frac{\beta(\left\|p\right\|)}{\left\|p\right\|}\right)^{d}.$

Let $\zeta=\frac{1}{2}\left\|x-y\right\|^{2}$. By Equation (35) we get
$J^{\prime}(\zeta)z=\left\|p\right\|$ implies
$J^{\prime}(\zeta)=\left\|p\right\|/\beta(\left\|p\right\|)$. Using Equation
(36), we get:

(79)
$\left|D^{2}_{xy}c(x,y)\right|=\frac{1}{\beta^{\prime}(\left\|p\right\|)}\left(\frac{\left\|p\right\|}{\beta(\left\|p\right\|)}\right)^{d}=\left|H^{\prime\prime}(z)\right|\left|J^{\prime}(\zeta)\right|^{d},$

which is the natural equivalent of Equation (77) and shows that the mixed
Hessian term is nonzero for cost functions in Euclidean space such that
$J^{\prime}(\zeta)\neq 0$ and $H^{\prime\prime}(z)\neq 0$.

###  3.6. Cost-Sectional Curvature for Defective Cost Functions

In this subsection, we present simple general formulas for checking the
positivity of the cost-sectional curvature for cost functions of the type
$c(x,y)=F(x\cdot y)$ on the sphere and
$c(x,y)=J\left(\frac{1}{2}\left\|x-y\right\|^{2}\right)$ on subsets of
Euclidean space. These formulas, interestingly, will only depend just the
first- and second-order derivatives of $F$ and $J$, respectively, (and an
expression of $x\cdot y$ in terms of $\left\|p\right\|$) even though the cost-
sectional curvature tensor requires taking four derivatives.

We will then use the formulas we derive to check the cost-sectional curvature
for various cost functions, including the cost from the lens reflector problem
I: $c(x,y)=-\log(n-x\cdot y)$ and the cost from the lens refractor problem II:
$c(x,y)=\log(\kappa x\cdot y-1)$. We will confirm that the cost function for
the lens reflector problem I does not satisfy the positive cost-sectional
curvature condition, but does for the lens reflector problem II, as was found
in [7]. However, we emphasize that the formulas we derive will be valid for
any defective cost function and as such extend beyond the lens refractor
problem.

The cost-sectional curvature condition comes from defining the following
fourth-order tensor, which was first defined in [12], but we will use the
formula from [10]. Note: there is a minus sign in front of the cost function
and that the derivatives with respect to $x$ need to be taken with respect to
the metric on either $\mathbb{S}^{2}$ or $\mathbb{R}^{d}$, respectively. Here,
we define the cost-sectional curvature tensor:

######  Definition 13.

On the domain that the map $y\mapsto-\nabla_{x}c(x,y)$ is injective,
$\left|D^{2}_{xy}c(x,y)\right|\neq 0$ and $c(x,y)\in C^{4}$, we define the map

(80)
$\mathfrak{G}_{c}(\xi,\eta)=D^{4}_{p_{k},p_{l},x_{i},x_{j}}\left[(x,p)\rightarrow-c(x,y)(p)|_{y=T(x)}\right]\xi_{i}\xi_{j}\eta_{k}\eta_{l}.$

An important condition to check for the purposes of regularity theory is
positive sectional curvature, usually known as the Aw condition, i.e. that on
an appropriate domain for all $x\in M$, $y\in M$, $\xi\in T_{x}M$, $\eta\in
T_{x}M$, $\xi\perp\eta$ that:

(81) $\mathfrak{G}_{c}(\xi,\eta)\geq 0.$

The Aw condition is an important condition to check such that smooth source
and target mass density functions which are bounded away from zero and
infinity lead to Optimal Transport mapping and potential functions that are
smooth as well. Without the Aw condition, it is possible to have smooth smooth
source and target mass density functions which are bounded away from zero and
infinity and possibly have an Optimal Transport mapping that is not even
continuous, see [10] and also [4].

A more stringent condition is known as the As condition, which is often the
condition used in deriving more precise regularity results, as in [11]. It
requires that there exist a constant $C_{0}>0$ such that

(82) $\mathfrak{G}_{c}(\xi,\eta)\geq
C_{0}\left|\xi\right|^{2}\left|\eta\right|.$

####  3.6.1. Computation of the cost-sectional curvature for $c(x,y)=F(x\cdot
y)$ on the sphere $\mathbb{S}^{2}$

To compute the Hessian on $\mathbb{S}^{2}$, we first compute the $3\times 3$
matrix Hessian for the cost function $c(x,y)=F(x\cdot y)$ using Euclidean
derivatives, since our metric on the sphere is induced by the surrounding
Euclidean space. Given a function $f:\mathbb{R}^{3}\rightarrow\mathbb{R}$, we
compute the Hessian as follows:

(83) $\nabla^{2}_{xx}f(x)=D^{2}_{xx}f(x)-(D_{x}f(x)\cdot x)\text{Id},$

where $D$ are the standard Euclidean (ambient) derivatives. We will thus
compute the term $\nabla^{2}_{xx}c(x,y)|_{y=T}$. We get:

(84) $D_{x}F(x\cdot y)\cdot x=F^{\prime}(x\cdot y)x\cdot
y=-\frac{\left\|p\right\|}{\sin\beta(\left\|p\right\|)}\cos\beta(\left\|p\right\|)=-\left\|p\right\|\cot\beta(\left\|p\right\|).$

We also compute:

(85) $\frac{\partial^{2}}{\partial x_{i}x_{j}}F(x\cdot
y)=F^{\prime\prime}(x\cdot y)y_{i}y_{j}.$

Thus, we get:

(86) $\nabla^{2}_{xx}F(x\cdot y)=F^{\prime\prime}(x\cdot
y)\begin{pmatrix}y_{1}^{2}&y_{1}y_{2}&y_{1}y_{3}\\\
y_{1}y_{2}&y_{2}^{2}&y_{2}y_{3}\\\
y_{1}y_{3}&y_{2}y_{3}&y_{3}^{2}\end{pmatrix}-F^{\prime}(x\cdot y)(x\cdot
y)\begin{pmatrix}1&0&0\\\ 0&1&0\\\ 0&0&1\end{pmatrix}.$

Recall that $y_{i}=x_{i}(x\cdot y)+\frac{p_{i}}{F^{\prime}(x\cdot y)}$. This
Hessian can be entirely expressed as a function of $x,p$, however, for
simplicity, we express in terms of $x,\zeta,p$, where $\zeta=x\cdot
y=\cos\beta(\left\|p\right\|)$:

(87)
$\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}=F^{\prime\prime}(\zeta)\left(x_{k}\zeta+\frac{p_{k}}{F^{\prime}(\zeta)}\right)\left(x_{l}\zeta+\frac{p_{l}}{F^{\prime}(\zeta)}\right)-\zeta
F^{\prime}(\zeta)\delta_{kl},$

Since $\xi$ and $p$ can be related via $\zeta=\cos\beta(\left\|p\right\|)$, we
relabel Equation (87) by $f_{1}(p)=F^{\prime\prime}(\zeta)\zeta^{2}$,
$f_{2}(p)=F^{\prime\prime}(\zeta)/(F^{\prime}(\zeta))^{2}$,
$f_{3}(p)=\zeta/F^{\prime}(\zeta)$, and $f_{4}(p)=-\zeta F^{\prime}(\zeta)$
and get

(88)
$\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}=f_{1}(p)x_{k}x_{l}+f_{2}(p)p_{k}p_{l}+f_{3}(p)(x_{k}p_{l}+x_{l}p_{k})+f_{4}(p)\delta_{kl}.$

With this simplification of notation, we proceed and compute:

(89)
$D_{p_{i}}\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}=D_{p_{i}}f_{1}(p)x_{k}x_{l}+D_{p_{i}}f_{2}(p)p_{k}p_{l}+f_{2}(p)(\delta_{ik}p_{l}+\\\
\delta_{il}p_{k})+D_{p_{i}}f_{3}(p)(x_{k}p_{l}+x_{l}p_{k})+f_{3}(p)(x_{k}\delta_{il}+x_{l}\delta_{ik})+D_{p_{i}}f_{4}(p)\delta_{kl}.$

Taking another derivative, we get:

(90)
$D_{p_{i}p_{j}}\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}=D_{p_{i}p_{j}}f_{1}(p)x_{k}x_{l}+D_{p_{i}p_{j}}f_{2}(p)p_{k}p_{l}+D_{p_{i}}f_{2}(p)(\delta_{jk}p_{l}+\delta_{jl}p_{k})+\\\
D_{p_{j}}f_{2}(p)(\delta_{ik}p_{l}+\delta_{il}p_{k})+f_{2}(p)(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})+D_{p_{i}p_{j}}f_{3}(p)(x_{k}p_{l}+x_{l}p_{k})+\\\
D_{p_{i}}f_{3}(p)(x_{k}\delta_{jl}+x_{l}\delta_{jk})+D_{p_{j}}f_{3}(p)(x_{k}\delta_{il}+x_{l}\delta_{ik})+D_{p_{i}p_{j}}f_{4}(p)\delta_{kl}.$

Now, we compute the following:

(91)
$\sum_{i,j,k,l}D_{p_{i}p_{j}}\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}\xi_{i}\xi_{j}\eta_{k}\eta_{l}=\left\langle
D^{2}_{pp}f_{1}(p)\xi,\xi\right\rangle(x\cdot\eta)^{2}+\left\langle
D^{2}_{pp}f_{2}(p)\xi,\xi\right\rangle(p\cdot\eta)^{2}+\\\
4(D_{p}f_{2}(p)\cdot\xi)(\xi\cdot\eta)(p\cdot\eta)+f_{2}(p)(\xi\cdot\eta)^{2}+2\left\langle
D^{2}_{pp}f_{3}(p)\xi,\xi\right\rangle(x\cdot\eta)(p\cdot\eta)+2(D_{p}f_{3}(p)\cdot\xi)(x\cdot\eta)(\xi\cdot\eta)+\\\
2(D_{p}f_{3}(p)\cdot\xi)(\xi\cdot\eta)(x\cdot\eta)+\left\langle
D^{2}_{pp}f_{4}(p)\xi,\xi\right\rangle\left|\eta\right|^{2}.$

We now use the fact that $\xi\perp\eta=0$ and $x\cdot\eta=0$:

(92)
$\sum_{i,j,k,l}D_{p_{i}p_{j}}\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}\xi_{i}\xi_{j}\eta_{k}\eta_{l}=\left\langle
D^{2}_{pp}f_{2}(p)\xi,\xi\right\rangle(p\cdot\eta)^{2}+\left\langle
D^{2}_{pp}f_{4}(p)\xi,\xi\right\rangle\left|\eta\right|^{2}.$

This shows that in order to satisfy condition As, we must check two
conditions, (1) the function $f_{4}(p)$ must be strictly negative definite
(strictly concave as a function of $\left\|p\right\|$) and (2)
$f_{2}(p)+f_{4}(p)$ must be strictly negative definite (strictly concave as a
function of $\left\|p\right\|$). In order to satisfy condition Aw, we simply
weaken the strictly negative definite condition to be just simply negative
definite.

We can use these formulas to check the cost-sectional curvature conditions Aw
and As for various cost functions on the sphere:

  1. (1)

For the far-field refractor problem I, we replace
$f_{2}(p)=F^{\prime\prime}(x\cdot y)/(F^{\prime}(x\cdot y))^{2}$ and
$f_{4}(p)=-(x\cdot y)F^{\prime}(x\cdot y)$ back in, we can now check the lens
problem:

(93) $f_{4}(p)=-(x\cdot y)F^{\prime}(x\cdot y)=-\frac{x\cdot y}{n-x\cdot y},$

and

(94) $f_{2}(p)=F^{\prime\prime}(x\cdot y)/(F^{\prime}(x\cdot y))^{2}=-1.$

Recall from Equation (7) that

(95) $x\cdot
y=\frac{\sqrt{1+(1-n^{2})\left\|p\right\|^{2}}+n\left\|p\right\|^{2}}{1+\left\|p\right\|^{2}}.$

The result, for different $n$ is shown in Figure 7. All are convex. This means
that the lens problem does not satisfy condition Aw, corroborating the result
in [7].

Figure 7. The function $f_{4}(p)$ plotted as a function of $\left\|p\right\|$.
Black indicates $n=1.333$, approximately the refractive index of water, red
indicates $n=1.52$, approximately the refractive index of glass, and blue
indicates $n=2.417$, approximately the refractive index of diamond.

  2. (2)

For the far-field refractor problem II, we compute:

(96) $f_{4}(p)=-(x\cdot y)F^{\prime}(x\cdot y)=-\frac{\kappa x\cdot
y}{\kappa(x\cdot y)-1},$

and

(97) $f_{2}(p)=F^{\prime\prime}(x\cdot y)/(F^{\prime}(x\cdot y))^{2}=-1.$

Then, we use Equation (13) to get:

(98) $x\cdot
y=\frac{\sqrt{\kappa^{2}+(\kappa^{2}-1)\left\|p\right\|^{2}}+\left\|p\right\|^{2}}{\kappa(1+\left\|p\right\|^{2})}.$

The result, in Figure 8 shows that $f_{4}(p)$ is concave in $\left\|p\right\|$
and thus satisfies condition Aw. In order to check As, we need to compute the
second derivative of the function in Equation (96) with respect to
$\left\|p\right\|$ using Equation (98). This is a straightforward exercise,
and a similar computation was done in [7], where it was shown that condition
As does not hold as $\left\|p\right\|\rightarrow\infty$. However, if we
restrict $\left\|p\right\|<p^{*}$, then $f_{4}(p)$ is strictly convex and thus
As can be satisfied.

Figure 8. The function $f_{4}(p)$ plotted as a function of $\left\|p\right\|$.
Black indicates $\kappa=1.333$, red indicates $\kappa=1.52$, and blue
indicates $\kappa=2.417$.

  3. (3)

For the squared geodesic cost $c(x,y)=\frac{1}{2}d_{\mathbb{S}^{2}}(x,y)^{2}$,
we have that $F(x\cdot y)=\frac{1}{2}\arccos(x\cdot y)^{2}$. Therefore, we
get:

(99) $f_{4}(p)=-(x\cdot y)F^{\prime}(x\cdot y)=\frac{(x\cdot y)\arccos(x\cdot
y)}{\sqrt{1-(x\cdot
y)^{2}}}=\frac{\left\|p\right\|\cos\left\|p\right\|}{\sqrt{1-\cos^{2}\left\|p\right\|}},$

which is concave in $\left\|p\right\|$. Also,

(100) $f_{2}(p)=F^{\prime\prime}(x\cdot y)/(F^{\prime}(x\cdot
y))^{2}=\frac{1}{(\arccos(x\cdot y))^{2}}-\frac{x\cdot y}{\arccos(x\cdot
y)\sqrt{1-(x\cdot y)^{2}}},$

which is convex in $\left\|p\right\|$, but the sum $f_{2}(p)+f_{4}(p)$ is
concave in $\left\|p\right\|$ as expected. Checking the condition As is a
straightforward exercise (take two derivatives), and the results from [11]
show that the condition As is satisfied for the squared geodesic cost
function.

####  3.6.2. Computation of the cost-sectional curvature for $c(x,y)=F(x\cdot
y)$ for Euclidean space $\mathbb{R}^{d}$

We can easily check the cost-sectional curvature conditions Aw and As in
Euclidean space. This can be contrasted with the previous computation on the
sphere in Section 3.6.1. Denote $\zeta=\frac{1}{2}\left\|x-y\right\|^{2}$. We
get:

(101) $\frac{\partial}{\partial
x_{i}}c(x,y)=J^{\prime}\left(\zeta\right)(x_{i}-y_{i}),$

and thus

(102) $\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}c(x,y)=J^{\prime\prime}(\zeta)(x_{i}-y_{i})(x_{j}-y_{j}),$

and

(103) $\frac{\partial^{2}}{\partial
x_{i}^{2}}c(x,y)=J^{\prime\prime}(\zeta)(x_{i}-y_{i})^{2}+J^{\prime}(\zeta),$

and thus,

(104)
$\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}=J^{\prime}(\zeta)\delta_{kl}+\\\
J^{\prime\prime}(\zeta)\begin{pmatrix}(x_{1}-y_{1})^{2}&(x_{1}-y_{1})(x_{2}-y_{2})&(x_{1}-y_{1})(x_{3}-y_{3})\\\
(x_{1}-y_{1})(x_{2}-y_{2})&(x_{2}-y_{2})^{2}&(x_{2}-y_{2})(x_{3}-y_{3})\\\
(x_{1}-y_{1})(x_{3}-y_{3})&(x_{2}-y_{2})(x_{3}-y_{3})&(x_{3}-y_{3})^{2}\end{pmatrix}.$

Since $x-y=-\frac{p}{J^{\prime}(\zeta)}$, we get:

(105)
$\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}=J^{\prime}(\zeta)\delta_{kl}-\frac{d}{d\zeta}\left(\frac{1}{J^{\prime}(\zeta)}\right)p_{k}p_{l}.$

Since we can relate $p$ and $\zeta$ via the equation
$\left\|p\right\|=J^{\prime}(\zeta)\sqrt{2\zeta}$, we may rewrite Equation
(105) as:

(106)
$\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}=f_{1}(p)\delta_{kl}+f_{2}(p)p_{k}p_{l},$

where $f_{1}(p)=J^{\prime}(\zeta)$ and
$f_{2}(p)=\frac{J^{\prime\prime}(\zeta)}{J^{\prime}(\zeta)^{2}}$. Thus, we can
compute:

(107)
$D_{p_{i},p_{j}}\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}=D_{p_{i}p_{j}}f_{1}(p)\delta_{kl}+D_{p_{i}p_{j}}f_{2}(p)p_{k}p_{l}+\\\
D_{p_{i}}f_{2}(p)(\delta_{jk}p_{l}+\delta_{jl}p_{k})+D_{p_{j}}f_{2}(p)(\delta_{ik}p_{l}+\delta_{il}p_{k})+f_{2}(p)(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}).$

Thus, we compute

(108)
$\sum_{i,j,k,l}D_{p_{i},p_{j}}\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}\xi_{i}\xi_{j}\eta_{k}\eta_{l}=\left\langle
D_{pp}f_{1}(p)\xi,\xi\right\rangle\left|\eta\right|^{2}+\left\langle
D_{pp}f_{2}(p)\xi,\xi\right\rangle(p\cdot\eta)^{2}+\\\
4(D_{p}f_{2}(p)\cdot\xi)(p\cdot\eta)(\xi\cdot\eta)+2f_{2}(p)(\xi\cdot\eta)^{2}.$

Since $\xi\perp\eta$, we get:

(109)
$\sum_{i,j,k,l}D_{p_{i},p_{j}}\left(\nabla^{2}_{xx}c(x,y)|_{y=T}\right)_{kl}\xi_{i}\xi_{j}\eta_{k}\eta_{l}=\left\langle
D_{pp}f_{1}(p)\xi,\xi\right\rangle\left|\eta\right|^{2}+\left\langle
D_{pp}f_{2}(p)\xi,\xi\right\rangle(p\cdot\eta)^{2}.$

As in the case of the sphere in Section 3.6.1, in order for condition As to
hold, it is necessary for two properties to hold: (1) that $f_{1}$ is strictly
concave in $\left\|p\right\|$ and that (2) $f_{1}+f_{2}$ is strictly concave
in $\left\|p\right\|$. For property Aw to hold, we weaken strict concavity to
simply concavity. We have

(110)
$f_{1}(p)=J^{\prime}(\zeta)=\frac{\left\|p\right\|}{\beta(\left\|p\right\|)},$

and

(111)
$f_{2}(p)=\frac{J^{\prime\prime}(\zeta)}{J^{\prime}(\zeta)^{2}}=\frac{1}{\left\|p\right\|}\left(\frac{1}{\left\|p\right\|\beta^{\prime}(\left\|p\right\|)}-\frac{1}{\beta(\left\|p\right\|)}\right).$

As an example, for the squared Euclidean distance cost,
$J(\zeta)=\zeta=\frac{1}{2}\left\|p\right\|^{2}$, for which
$J^{\prime}(\zeta)=1$, which is concave in $\left\|p\right\|$. Also,
$J^{\prime\prime}(\zeta)=0$, so we see that we satisfy the cost-sectional
curvature condition Aw, but As does not hold, which is corroborated by [15].

######  Remark 14.

Some explicit computations for more general Riemannian manifolds have been
done in, for example, [4] where there is an interesting example where for the
squared geodesic cost function $c(x,y)=\frac{1}{2}d_{M}(x,y)^{2}$, the cost-
sectional curvature is not positive for an ellipsoid of revolution. For
general cost functions and general Riemannian manifolds, the formula for the
cost-sectional curvature tensor is formidable, see, for example the survey
paper by Villani [16].

##  4\. Regularity and Solvability

###  4.1. Ma-Trudinger-Wang Conditions for Defective Cost Functions

It was proved, in [7] that there exist weak solutions to the lens refractor
problem I & II, provided that, most importantly, the conditions in Equations
(11) and (12) were met, respectively. In [9], it was proved that smooth
solutions exist for the lens refractor problem II, with an argument using
supporting ellipsoids and hyperboloids, respectively. Here, instead, we are
inspired to take the route taken in [11] and prove a regularity result for the
lens refractor problem II, but, more generally, for any defective cost
function that satisfies the MTW conditions.

We begin by stating the MTW conditions, formulated originally in [12], but we
focus on the Riemannian generalization as stated in [11]. Given a compact
domain $D\subset\mathbb{S}^{2}\times\mathbb{S}^{2}$, denote by
$\pi_{1}:\mathbb{S}^{2}\times\mathbb{S}^{2}\mapsto\mathbb{S}^{2}$, the
projection $\pi_{1}(x,y)=x$ and its inverse
$\pi_{1}^{-1}(x)=\left\\{x\right\\}\times\mathbb{S}^{2}$. For any
$x\in\pi_{1}(D)$, we denote by $D_{x}$ the set $D\cap\pi_{1}^{-1}(x)$. Then,
we introduce the following conditions:

######  Hypothesis 15.

  * $\mathbf{A0}$

The cost function $c$ belongs to $C^{4}(D)$.

  * $\mathbf{A1}$

For all $x\in\pi_{1}(D)$, the map $y\rightarrow-\nabla_{x}c(x,y)$ is injective
on $D_{x}$.

  * $\mathbf{A2}$

The cost function $c$ satisfies $\det D^{2}_{xy}c\neq 0$ for all $(x,y)$ in
$D$.

  * $\mathbf{Aw}$

The cost-sectional curvature is non-negative on $D$. That is, for all
$(x,y)\in D$, for all $\xi,\eta\in T_{x}\mathbb{S}^{2}$, $\xi\perp\eta$,

(112) $\mathfrak{G}_{c}(x,y)(\xi,\eta)\geq 0$

  * $\mathbf{As}$

The cost-sectional curvature is uniformly positive on $D$. That is, for all
$(x,y)\in D$, for all $\xi,\eta\in T_{x}\mathbb{S}^{2}$, $\xi\perp\eta$,

(113) $\mathfrak{G}_{c}(x,y)(\xi,\eta)\geq
C_{0}\left|\xi\right|^{2}\left|\eta\right|^{2}$

Denoting $z^{*}=\left\\{\min z:G^{\prime\prime}(z)=0\ \text{or}\
G^{\prime\prime}(z)=\infty,z\in(0,\pi)\right\\}$, we now choose the subdomain
$D$ as follows. For each $x\in\mathbb{S}^{2}$, we denote the corresponding
geodesic ball $B_{x}(z^{*})$. Then, let $0<\gamma<z^{*}$

(114) $D_{\gamma}=\cup_{x}\left\\{x\right\\}\times B_{x}(\gamma)$

Now that we have defined $D_{\gamma}$, the computations of Section 3 allow for
us to verify the MTW conditions. Assuming that $F\in C^{4}([\cos\gamma,1])$,
$G\in C^{4}([0,\gamma])$, then, $c$ satisfies $\mathbf{A0}$ on $D_{\gamma}$.
As long as $c$ satisfies the hypotheses in Theorem 3, then $\mathbf{A1}$ is
satisfied on $D_{\gamma}$. By the computation in Equation (77), if we have
$F^{\prime}(\zeta)\neq 0$ and $G^{\prime\prime}(z)\neq 0,\infty$, then $c$
satisfies $\mathbf{A2}$ on $D_{\gamma}$. As shown in Section 3.6, if
$f_{4}(p)$ is concave and $f_{4}(p)+f_{2}(p)$ is concave for
$\left\|p\right\|<p^{*}$, then $\mathbf{Aw}$ is satisfied for $c$ on
$D_{\gamma}$. The condition $\mathbf{As}$ can then be checked by taking two
derivatives of $f_{2}(p)$ and $f_{4}(p)$ with respect to $\left\|p\right\|$.

######  Remark 16.

Based on Remarks 5, 9, 12, and the work in Section 3.6.2, the MTW conditions
can be checked for a cost function on Euclidean space in an analogous
subdomain $D_{\gamma}\subset\mathbb{R}^{d}\times\mathbb{R}^{d}$.

####  4.1.1. Example of Computations for Cost Functions in Theorem 4.1 of
Loeper [11]: Power Costs on the Unit Sphere
$c(x,y)=\frac{1}{p}d_{\mathbb{S}^{2}}(x,y)^{p}$

As remarked in Section 1, an important result of the work in this paper is to
allow for the MTW conditions to be checked easily for a wide class of cost
functions, not just defective cost functions. Note that if $z^{*}\geq\pi$,
then our cost functions naturally fit into Theorem 4.1 of Loeper [11], which
are the cost functions that can be checked with the preexisting regularity
theory of Loeper. The benefit is that in Section 3 we have identified simple
conditions on the cost function which allow for the hypotheses of Theorem 4.1
of Loeper (including the MTW conditions) to be verified. As an illustration of
this, we verify the MTW conditions for the power cost functions on the sphere:
$c(x,y)=\frac{1}{s}d_{\mathbb{S}^{2}}(x,y)^{s}=\frac{1}{s}\arccos(x\cdot
y)^{s}$.

  1. (1)

First, we check that the power cost functions satisfy some preliminary
conditions from Theorem 4.1 in Loeper [11]. That is, $G(z)$ is smooth and
strictly increasing with $G^{\prime}(0)=0$. We immediately verify that
$G(z)=\frac{1}{s}z^{s}$ is smooth and strictly increasing for $s>0$ and
$G^{\prime}(z)=z^{s-1}$, so that means that if $s>1$ we have that
$G^{\prime}(0)=0$. From now on assume that $s>1$.

  2. (2)

Now, we verify the MTW conditions on
$D=\mathbb{S}^{2}\times\mathbb{S}^{2}\setminus\text{antidiag}$ using the
formulas from this manuscript. Clearly, $\mathbf{A0}$ is satisfied on $D$.

  3. (3)

In order for $\mathbf{A1}$ to be satisfied, we check the hypotheses of Theorem
3. We need $F^{\prime}(\zeta)\neq 0$ for $\zeta\in(-1,1]$ and
$G^{\prime\prime}(z)$ strictly positive or negative on $(0,\pi)$. Now,
$F(\zeta)=\frac{1}{s}\arccos(\zeta)^{s}$. Therefore,
$F^{\prime}(\zeta)=-\arccos(\zeta)^{s-1}/\sqrt{1-\zeta^{2}}$. Thus, we have
$F^{\prime}(\zeta)\neq 0$ for $\zeta\leq 1$ if $s\leq 2$. Since
$G^{\prime\prime}(z)=(s-1)z^{s-2}$, we have that $G^{\prime\prime}(z)$ is
strictly positive on $(0,\pi)$ for $s>1$. Therefore, condition A1 is satisfied
for $1<s\leq 2$.

  4. (4)

The condition $\mathbf{A2}$, by Equation (77) is satisfied for $1<s\leq 2$.

  5. (5)

In order to verify condition $\mathbf{As}$, we compute the term
$R_{1}(\left\|p\right\|)$ from Equation (21). This is simple, using the fact
that $\beta$ and $G^{\prime}$ are inverses and therefore
$\left\|p\right\|=G^{\prime}(z)=z^{s-1}$. Therefore,
$z=\left\|p\right\|^{1/(s-1)}$. Therefore, we get
$R_{1}(\left\|p\right\|)=\cos\left(\left\|p\right\|^{\frac{1}{s-1}}\right)$.
We must check that $f_{4}(p)$ and $f_{2}(p)+f_{4}(p)$ are strictly convex as a
function of $\left\|p\right\|$, as shown in Section 3.6. We have

(115) $f_{4}(p)=-\zeta
F^{\prime}(\zeta)=\frac{\zeta\arccos(\zeta)^{s-1}}{\sqrt{1-\zeta^{2}}}$

where $\zeta=x\cdot y$. Thus,

(116)
$f_{4}(p)=\frac{\left\|p\right\|\cos\left(\left\|p\right\|^{\frac{1}{s-1}}\right)}{\sqrt{1-\cos^{2}\left(\left\|p\right\|^{\frac{1}{s-1}}\right)}},\
\ \ \left\|p\right\|\in[0,\pi^{s-1}]$

and

(117) $f_{2}(p)=F^{\prime\prime}(\zeta)/(F^{\prime}(\zeta))^{2}=\\\
\left\|p\right\|^{\frac{s}{s-1}}\left((s-1)-\left\|p\right\|^{\frac{1}{s-1}}\frac{\cos\left(\left\|p\right\|^{\frac{1}{s-1}}\right)}{\sqrt{1-\cos^{2}\left(\left\|p\right\|^{\frac{1}{s-1}}\right)}}\right),\
\ \ \left\|p\right\|\in[0,\pi^{s-1}]$

(118) $=(s-1)\arccos(\zeta)^{-s}-\arccos(\zeta)^{2-2s}f_{4}(p)$

(119)
$=\frac{s-1}{\left\|p\right\|^{\frac{s}{s-1}}}-\frac{1}{\left\|p\right\|^{2}}f_{4}(p)$

and therefore,

(120)
$f_{2}(p)+f_{4}(p)=\frac{s-1}{\left\|p\right\|^{\frac{s}{s-1}}}+\frac{\left\|p\right\|^{2}-1}{\left\|p\right\|^{2}}f_{4}(p)$

This can be used to confirm, for example, that condition $\mathbf{As}$ does
not hold for $s>2$, since for small values of $\left\|p\right\|$, we have:

(121)
$f_{2}(p)+f_{4}(p)=(s-2)\left\|p\right\|^{\frac{-s}{s-1}}+o\left(\left\|p\right\|^{\frac{-s}{s-1}}\right)$

so the cancellation of the highest order term only happens when $s=2$. The
fact that the power costs for $s\neq 2$ do not satisfy the $\mathbf{As}$
condition agrees with the result for Euclidean power costs in [15].

###  4.2. Conditions on Source and Target Mass Density Functions that
Constrain the Mapping

Now that we have verified the MTW conditions on a subdomain $D_{\gamma}$, we
need to make sure that our source and target densities $f$ and $g$ are such
that they do not require mass to move too far, as was the case in the example
in Section 3.2. The work in [11] showed that certain technical conditions on
the source and target mass distributions restricted mass to move away from the
cut locus on the sphere. Therefore, in that paper, the MTW conditions were
shown to be satisfied on the subdomain
$D\subset\mathbb{S}^{2}\times\mathbb{S}^{2}$, where the domain $D$ was chosen
to be $D=\mathbb{S}^{2}\times\mathbb{S}^{2}\setminus\text{antidiag}$, where
the set antidiag is defined as $(x,-x)\in\mathbb{S}^{2}\times\mathbb{S}^{2}$.

In this section, we show that more strict conditions can be found on the
source and target mass distributions that restrict the Optimal Transport
mapping to satisfy $d_{\mathbb{S}^{2}}(x,T)\leq\gamma$ for any desired
$\gamma<\pi$, where the Optimal Transport mapping $T$ arises from certain
defective cost functions, such as the cost function arising in the lens
refractor problem. This then allows us to use the regularity framework in [11]
to build a regularity theory for a large class of defective cost functions.

Define $G(\phi)(x)=\text{exp}_{x}(\nabla\phi(x))$, which is the Optimal
Transport mapping arising from the squared geodesic cost function. Also define
$T(\phi)(x)=\text{exp}_{x}\left(\frac{\nabla\phi(x)}{\left\|\nabla\phi(x)\right\|}\beta\left(\nabla\phi(x)\right)\right)$
where $\beta$ is the function in Theorem 3 arising from a defective cost
function. As before, let $z^{*}$ denote the smallest value for which
$G^{\prime\prime}(z)=0$ or $G^{\prime}(z)=\infty$. First, we begin with a
result due to Delanoë and Loeper [3].

######  Theorem 17 (Delanoë and Loeper).

Let $\phi:\mathbb{S}^{2}\rightarrow\mathbb{R}$ be a $C^{3}$ functions such
that $G(\phi)$ is a diffeomorphism and let
$\rho:\mathbb{S}^{2}\rightarrow\mathbb{R}$ be the positive $C^{1}$ functions
defined by:

(122) $G(\phi)_{\\#}\text{dVol}=\rho\text{dVol},$

where dVol is the standard $2$-form on the unit sphere $\mathbb{S}^{2}$
induced by the standard Euclidean metric on $\mathbb{R}^{3}$. Then,

(123) $\max_{\mathbb{S}^{2}}\left|\nabla\phi\right|\leq
C\max_{\mathbb{S}^{2}}\left|\nabla[\rho^{-1}]\right|.$

We continue by defining a $c$-convex function.

######  Definition 18.

A function $\phi:\mathbb{S}^{2}:\rightarrow\mathbb{R}$ is $c$-convex if at
each point $x\in\mathbb{S}^{2}$ there exists a point $y\in\mathbb{S}^{2}$ and
a value $\phi^{c}(y)$ such that:

(124) $\displaystyle-\phi^{c}(y)-c(x,y)$ $\displaystyle=\phi(x)$ (125)
$\displaystyle-\phi^{c}(y)-c(x^{\prime},y)$
$\displaystyle\leq\phi(x^{\prime}),\ \ \ \forall x^{\prime}\in\mathbb{S}^{2}.$

######  Theorem 19.

Let $\phi:\mathbb{S}^{2}\rightarrow\mathbb{R}$ be a $C^{3}$, $c$-convex
function, for the defective cost function $c$, such that $G(\phi)$ and
$T(\phi)$ are diffeomorphisms and the Jacobian $\left|\nabla T\circ
G^{-1}\right|$ has a bounded derivative and let
$\tilde{\rho}:\mathbb{S}^{2}\rightarrow\mathbb{R}$ be the positive $C^{1}$
function defined by:

(126) $T(\phi)_{\\#}\text{dVol}=\tilde{\rho}\text{dVol}.$

Then, there exists a constant $\tilde{C}>0$ such that:

(127)
$\max_{\mathbb{S}^{2}}\left|\nabla\phi\right|\leq\tilde{C}\left(\max\tilde{\rho}+\max\left|\nabla\tilde{\rho}\right|\right).$

###### Proof.

Note: if $\rho=\tilde{\rho}=1$, then $\phi$ is $c$-convex, $G(\phi)$ and
$T(\phi)$ are diffeomorphisms and $\left|\nabla T\circ G^{-1}\right|$ has a
bounded derivative.

From Theorem 17, the $c$-convex function is $C^{3}$ and $G(\phi)$ is a
diffeomorphism. Therefore, $\phi$ satisfies the bound:

(128) $\max_{\mathbb{S}^{2}}\left|\nabla\phi\right|\leq
C\max_{\mathbb{S}^{2}}\left|=\nabla[\rho^{-1}]\right|.$

where $\rho$ is defined by:

(129) $G(\phi)_{\\#}\text{dVol}=\rho\text{dVol},$

Since $T(\phi)$ is a diffeomorphism, it satisfies $0<\left|\nabla
T(\phi)(x)\right|<\infty$. The function $\tilde{\rho}$ then satisfies
$\tilde{\rho}=(T\circ G^{-1})_{\\#}\rho$ and likewise $\rho=(G\circ
T^{-1})_{\\#}\tilde{\rho}$. We denote $S=T\circ
G^{-1}(x)=\text{exp}_{x}\left(\frac{\nabla\phi(x)}{\left\|\nabla\phi(x)\right\|}\left(-\nabla\phi(x)+\beta\left(\nabla\phi(x)\right)\right)\right)$.
We have then, that:

(130) $\tilde{\rho}=\rho(S^{-1}(x))\left|\nabla S^{-1}(x)\right|.$

Since $S$ is invertible and a diffeomorphism, since it is the composition of
diffeomorphisms, then we have $c\leq\left|\nabla S^{-1}(x)\right|\leq C$, and
since $G$ is a diffeomorphism, we have $0<m\leq\rho\leq M$. Therefore,

(131) $\tilde{\rho}\geq cm.$

Then, since

(132) $\rho=\tilde{\rho}(S(x))\left|\nabla S(x)\right|,$

we get:

(133) $\rho\geq C\min\tilde{\rho}\geq cCm,$

and therefore,

(134) $\frac{1}{\rho}\leq\frac{1}{cCm}.$

Now, since $\left|\nabla S\right|$ has bounded derivatives, there exists a
$\mu>0$ such that:

(135) $\left|\nabla\rho\right|=\left|\nabla
S^{-1}_{\\#}\tilde{\rho}\right|=\left|\nabla\left[\tilde{\rho}(S(x))\left|\nabla
S(x)\right|\right]\right|$

(136) $=\left|\nabla\tilde{\rho}(S(x))\left|\nabla
S(x)\right|+\tilde{\rho}(S(x))\nabla\left|\nabla
S(x)\right|\right|\leq\left|c\nabla\tilde{\rho}(S(x))+\mu\tilde{\rho}(S(x))\right|,$

Therefore, we get:

(137)
$\left|\nabla\phi\right|\leq\left|\nabla[\rho^{-1}]\right|\leq\tilde{C}\left(\max\tilde{\rho}+\max\left|\nabla\tilde{\rho}\right|\right).$

∎

This result shows that provided that we control the ratio of $g/f$ and the
ratio of the derivative of $g/f$, then the derivative of the potential
function $\phi$ can be controlled as well. This, in turn, via the fact that
for a defective cost function, as long as $\left\|p\right\|<p^{*}$, means that
the Optimal Transport mapping will not transport a distance greater than
$z^{*}$ from Definition 7. At this point, we can adapt the results of Loeper
in [11] to the case of defective cost functions to get the following result:

######  Theorem 20.

Let $\mu$, $\nu$ be two $C^{\infty}$ probability measures that are strictly
bounded away from zero and $u$ a $c$-convex potential function, where $c$ is a
defective cost function, such that, defining
$T(u)(x)=\text{exp}_{x}\left(\frac{\nabla u(x)}{\left\|\nabla
u(x)\right\|}\beta(\left\|\nabla u(x)\right\|)\right)$, we have
$T(u)_{\\#}\mu=\nu$ and also the bound $\left|\nabla u\right|<p^{*}$, where
$p^{*}$ is defined in Definition 7. Additionally, let $c$ satisfy assumption
$\mathbf{As}$ on $D_{\gamma}$. Then $u\in C^{\infty}(\mathbb{S}^{2})$.

###### Proof.

The proof follows the line of reasoning in Section $5$ of Loeper [11]. First,
we have shown by the estimate in Theorem 19 that there exist potential
functions whose derivatives can be bounded by $g/f$ and the derivative of
$g/f$. Assuming that $\mathbf{As}$ is satisfied, we use a $C^{2}$ estimate
established in [12] on $u$. Then, we use the method of continuity to construct
smooth solutions for any $\mu$, $\nu$ for any smooth positive densities
satisfying the hypotheses of Theorem 19. ∎

By applying the above theorem to the case of refractor problem II, and keeping
in mind that the condition As is only true for the cost function
$c(x,y)=-\log(\kappa x\cdot y-1)$ for all $\left\|p\right\|<p^{*}$, we arrive
at the following result for the refractor problem II:

######  Corollary 21.

Let $\mu$, $\nu$ be two $C^{\infty}$ probability measures that are strictly
bounded away from zero. Then, the potential function $u$ for the refractor
problem II is $C^{\infty}$.

######  Remark 22.

One interesting consequence of Theorem 20 is that it suggests that some lens
refractor problems may be solved using a series of $N$ lenses
$\left\\{L_{i}\right\\}_{i=1,\dots,N}$, in the case that it is impossible to
solve using a single lens.

##  5\. Conclusion

In examining the far-field lens refractor problem and its solvability
condition, we found that we could answer many theoretical and computational
questions by formulating a theory for defective cost functions on the unit
sphere. Using these definitions, we examined when the MTW conditions held for
such cost functions and derived formulas for verifying the cost-sectional
curvature for defective cost functions. We discussed how to extend these ideas
to Euclidean space and some ideas to general Riemannian manifolds. We examined
a sufficient condition on the source and target densities $f$ and $g$,
respectively, that ensured the mapping $T$ satisfied the solvability bound,
which was achieved by bounding the ratio $g/f$ and the derivative of $g/f$ in
the sup-norm. Using this, we could employ the regularity framework of Loeper
in order to establish $C^{\infty}$ smoothness for the potential function for
defective cost functions satisfying the strictly positive cost-sectional
curvature condition.

Acknowledgements: I would like to especially thank Brittany Hamfeldt for
introducing me to the lens refractor problem and working alongside for some of
the initial exploratory computations. I would also like to thank Rene Cabrera
for listening to and discussing the ideas as they developed.

Data Availability Statement: The author declares that the data supporting the
findings of this study are available within the paper.

## References

  * [1] Martijn J. H. Anthonissen, Lotte B. Romijn, Jan H. M. ten Thije Boonkkamp, and Wilbert L. IJzerman.  Unified mathematical framework for a class of fundamental freeform optical systems.  Optics Express, 29(20), September 2021.
  * [2] M. Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transportation distances.  Advances in Neural Information Processing Systems, 26, June 2013\.
  * [3] P. Delanoë and G. Loeper.  Gradient estimates for potentials of invertible gradient-mappings on the sphere.  Calculus of Variations and Partial Differential Equations, 26(3):297–311, 2006.
  * [4] A. Figalli, L. Rifford, and C. Villani.  On the Ma-Trudinger-Wang curvature on surfaces.  Calculus of Variations, 39:307–332, 2010.
  * [5] Wilfrid Gangbo and Robert J McCann.  The geometry of optimal transportation.  Acta Mathematica, 177:113–161, 1996.
  * [6] C. E. Guitiérrez and Q. Huang.  The near field refractor.  Annaless de l’I.H.P. Analyse non linéaire, 31(4):655–684, 2014\.
  * [7] C. E. Gutiérrez and Q. Huang.  The refractor problem in reshaping light beams.  Archive for Rational Mechanics and Analysis, 193(2):423–443, 2009\.
  * [8] Brittany Froese Hamfeldt and Axel G. R. Turnquist.  A convergent finite difference method for optimal transport on the sphere.  Journal of Computational Physics, 445, November 2021.
  * [9] Aram L. Karakhanyan.  On the regularity of weak solutions to refractor problems.  Armenian Journal of Mathematics, 2(1):28–37, 2009.
  * [10] G. Loeper.  On the regularity of solutions of optimal transportation problems.  Acta Mathematica, 202:241–283, 2009.
  * [11] G. Loeper.  Regularity of optimal maps on the sphere: the quadratic cost and the reflector antenna.  Archive for rational mechanics and analysis, 199(1):269–289, 2011\.
  * [12] X.-N. Ma, N. X. Trudinger, and X.-J. Wang.  Regularity of potential functions of the optimal transportation problem.  Archive for Rational Mechanics and Analysis, 177(2):151–183, 2005\.
  * [13] R. J. McCann.  Polar factorization of maps on Riemannian manifolds.  Geometric and Functional Analysis, 11:589–608, 2001.
  * [14] L. Romijn, J. ten Thije Boonkkamp, and W. IJzerman.  Freeform lens design for a point source and far-field target.  Journal of the Optical Society of America A, 36(11):1926–1939, 2019\.
  * [15] Neil S. Trudinger and X.-J. Wang.  On the second boundary value problem for monge-ampère type equations and optimal transportation.  Ann. Scuola Norm. Sup. Pisa Cl. Sci., 8(5):143–174, 2009.
  * [16] C. Villani.  Regularity of optimal transport and cut locus: From nonsmooth analysis to geometry to smooth analysis.  Discrete & Continuous Dynamical Systems, 30(2):559–571, 2011.
  * [17] X.-J. Wang.  On the design of a reflector antenna.  IOP Science, 12:351–375, 1996.
  * [18] X.-J. Wang.  On the design of a reflector antenna II.  Calculus of Variations and Partial Differential Equations, 20(3):329–341, 2004.
  * [19] N. K. Yadav.  Monge-Ampère Problems with Non-Quadratic Cost Function: Application to Freeform Optics.  PhD thesis, Technische Universiteit Eindhoven, Eindhoven, Netherlands, 2018.


# A question about points on an elliptic curve with prime denominator

Simon L Rydin Myerson Mathematics Insitute, University of Warwick; for contact
details see https://maths.fan

(7 July 2023
Problem session of Journées Arithmetiques 2023
Session chair: Michel Waldschmidt)

###### Abstract

Let $E$ be an elliptic curve defined by a Weierstrass equation with integer
coefficients. Any $P\in E(\mathbb{Q})\setminus\\{O\\}$ is of the form
$P=(\frac{x(P)}{z(P)^{2}},\frac{y(P)}{z(P)^{3}})$ where
$x(P),y(P)\in\mathbb{Z}$, $z(P)\in\mathbb{N}$ and in addition both
$\gcd(x(P),z(P))=1$ and $\gcd(y(P),z(P))=1$ hold. The question is: how often
is $z(P)$ a prime number?

##  1 Preamble

This is an open problem submitted by the author and presented by Michel
Waldschmidt to the JA2023 conference. I have done my best to include a number
of different plausible strategies and questions for what I feel is a family of
underexplored questions badly in want of new insights. Interested readers who
may have questions or ideas to make progress are encouraged to contact the
author; I would be delighted to do what I can to help workers advance into
this field.

##  2 What we expect

Let $E(\mathbb{Q})\simeq T\oplus\mathbb{Z}^{R}$ where $T$ is the torsion
subgroup. Let $P_{1},\dotsc,P_{r}$ generate the free part of $E(\mathbb{Q})$
and write $\vec{n}\cdot\vec{P}=n_{1}P_{1}+\dotsb+n_{r}P_{r}$.

A general point of $E(\mathbb{Q})$ is of the form $P=\vec{n}\cdot\vec{P}+S$
where $\vec{n}\in\mathbb{Z}^{r},S\in T$. For some constants $0<c<C$ depending
on $E$, we have $e^{c|\vec{n}|^{2}}<z(P)<e^{C|\vec{n}|^{2}}$. Heuristically
the probability that $z(P)$ is prime might be around $1/\log
z(P)\approx|\vec{n}|^{-2}$.

Meanwhile the number of points $P$ with $N\leq|\vec{n}|\leq 2N$ is around
$N^{r}$.

Let $Z$ be a large parameter. Based on the above, we guess the total number of
points with $z(P)$ prime and $z(P)\leq Z$ is around

$\\#\\{P\in E(\mathbb{Q}):z(P)\text{ prime},z(P)\leq Z\\}\\\
\approx\sum_{\begin{subarray}{c}N=2^{k},\,k\in\mathbb{Z}\\\
e^{cN^{2}}<Z\end{subarray}}N^{r-2}\approx\begin{cases}\text{constant}&r\leq
1,\\\ \log\log Z&r=2,\\\ (\log Z)^{\frac{r}{2}-1}&r>2.\end{cases}$

##  3 What is known

######  Lemma 1.

Suppose $E$ is given by a minimal Weierstrass equation. If $z(nP)$ is prime
and $n$ is sufficiently large, then $n$ is prime.

With a little work the condition that the Weierstrass equation is minimal can
be removed (M. Verzobio, personal communication).

###### Proof sketch..

Suppose that $z(P)$ is prime and $n=rs$ with $r\geq s>1$. Then because
$rP\equiv O\bmod z(rP)$ we must have $nP\equiv O\bmod z(rP)$, which implies
that $z(rP)$ divides the denominator $z(nP)$. If $z(nP)$ is prime then
$z(rP)=z(nP)$ or $z(rP)=1$. But $r\geq\sqrt{n}$ and $n$ is assumed large, so
$rP$ cannot be integral and $z(rP)\neq 1$. Thus $z(rP)=z(nP)$.

But if $n$ is sufficiently large, then since the Weierstrass equation is
minimal, $z(nP)$ has a prime divisor which is coprime to $z(rP)$, by a result
of Silverman [4, pp232-233]. This is impossible, so $n$ must have been prime.
∎

If $E(\mathbb{Q})\simeq\mathbb{Z}$, and $Q$ is a generator of $E(\mathbb{Q})$,
then

$\\{P\in E(\mathbb{Q}):z(P)\leq Z\\}\subseteq\\{nQ:n\ll\sqrt{\log Z}\\}.$

From this and the lemma above we get:

######  Corollary 1.

Suppose $E$ is given by a minimal Weierstrass equation, and that
$E(\mathbb{Q})\simeq\mathbb{Z}$. Then

$\frac{\\#\\{P\in E(\mathbb{Q}):z(P)\text{ prime},z(P)\leq Z\\}}{\\#\\{P\in
E(\mathbb{Q}):z(P)\leq Z\\}}\ll\frac{1}{\log\log Z}.$

Actually there are some cases where much more is known. There are several
results saying that for a suitable isogeny $\phi:E^{\prime}\to E$ and a
suitable field $K$, there are only finitely many points $P\in
E(\mathbb{Q})\cap\phi(E^{\prime}(K))$ such that $z(P)$ is prime.

As I understand it the proof idea is always something like this: Imagine for
simplicity that $K=\mathbb{Q}$. We try to show that the denominator
$z(\phi(Q))$ has to be divisible by more primes than the denominator $z(Q)$,
so if $P=\phi(Q)$ has prime denominator then $Q$ is integral and there are
only finitely many such $Q$. When $K\neq\mathbb{Q}$ one has to look at the
number of places at which $Q$ has positive valuation, instead of primes
dividing $z(Q)$. I do not intend to provide a complete list of work of this
kind, but I am thinking in particular of [1, 2, 3].

##  4 Specific questions

  * •

Show that Corollary 1 still holds if $E$ has torsion, i.e. if
$E(\mathbb{Q})\simeq T\oplus\mathbb{Z}$.

  * •

Do better than Corollary 1. That is take an elliptic curve for which

$\\#\\{P\in E(\mathbb{Q}):z(P)\text{ prime}\\}$

is not already known to be finite, and prove that

$\frac{\\#\\{P\in E(\mathbb{Q}):z(P)\text{ prime},z(P)\leq Z\\}}{\\#\\{P\in
E(\mathbb{Q}):z(P)\leq Z\\}/\log\log Z}\to 0\qquad(Z\to\infty).$

## Acknowledgements

This question was inspired by my collaboration with Subham Bhakta, Dan
Loughran, and Masahiro Nakahara.

My subsequent exploration of these ideas has been carried out in the company
of Subham Bhakta and Matteo Verzobio, to whom the author owes many of the
suggestions and hints in this document.

I thank the organisers of Journées Arithmétiques for inviting me to
participate, for their gracious hospitality and for their hard work to make
the meeting run smoothly; I am grateful in particular to Cécile Dartyge,
Youness Lamzouri and Paola Schneider.

## References

  * [1] Manfred Einsiedler, Graham Everest and Thomas Ward  “Primes in Elliptic Divisibility Sequences”  In _LMS Journal of Computation and Mathematics_ 4 Cambridge University Press, 2001, pp. 1–13  DOI: 10.1112/S1461157000000772
  * [2] Graham Everest and Helen King  “Prime Powers in Elliptic Divisibility Sequences”  In _Mathematics of Computation_ 74.252 American Mathematical Society, 2005, pp. 2061–2071  URL: http://www.jstor.org/stable/4100227
  * [3] Valéry Mahé  “Prime power terms in elliptic divisibility sequences”  In _Math. Comp._ 83.288, 2014, pp. 1951–1991  DOI: 10.1090/S0025-5718-2013-02790-1
  * [4] Joseph H. Silverman  “Wieferich’s criterion and the abc-conjecture”  In _Journal of Number Theory_ 30.2, 1988, pp. 226–237  DOI: https://doi.org/10.1016/0022-314X(88)90019-4