\section{Introduction}
Classical models in financial mathematics work similar in principle. They fix a certain underlying probability space and assume that the random future behaviour of the underlying asset price process is specified somehow, for example as the solution of a stochastic differential equation. Further assuming no-arbitrage and completeness of the considered financial market, a unique equivalent martingale measure, i.e. a probability measure such that the discounted asset price process is a martingale, exists by the fundamental theorem of asset pricing. Then using the law of one price, it is possible to derive the uniquely determined price of an exotic option written on the underlying by calculating either the expected discounted payoff of the exotic option with respect to the equivalent martingale measure or the price of a self-financing, replicating hedging strategy. Most of the models are of parametric form, where the parameters are determined by calibrating the model to observable market prices of certain options.

The assumptions which are made are simplifying and often quite  unrealistic thus leading to model prices which are likely to be inaccurate and unreliable. Indeed, various studies, see e.g. \cite{schoutens2003perfect}, observe a great range of option prices when calibrating several different models to the same underlying market.

On the other hand, over the years, trading of call options became so liquid that  \cite{dupire1993model, dupire1994pricing} and others argued that they should rather be considered as contingent claims with exogenously fixed prices. Thus, the prices of European call options became available as information for pricing other, more complicated exotic options. In model-independent finance, this is used under the idealizing assumption that for the maturities of interest, the call option prices are observable for a continuum of strike prices.

The target of model-independent finance is now to price exotic options such that the prices satisfy no-arbitrage and are consistent with  observable call option prices. Therefore, no specific martingale measure but a set of different consistent martingale measures emerges from the analysis. Thus, no unique option price but a range of possible option prices may be derived. In return, the risk of model misspecification is eliminated. This intuition was first formalized by Hobson in his famous paper \cite{hobson1998robust}, where he used the no-arbitrage assumption and the knowledge about the call option prices to derive upper and lower price bounds for the lookback option in continuous time.

Such price bounds can be derived in various ways. There is one stream of literature which uses Skorokhod-type stopping problems to derive upper bounds (see e.g. \cite{hobson1998robust}) and another one using methods from optimal transport, see e.g. \cite{beiglbock2013model}. We will pursue the latter approach in this paper.

 In optimal transport, the problem is to minimize the cost that transportation of mass from one point to another generates in the sense that a cost-minimal transport allocation is aimed for. Mathematically, we may specify the mass at the origins and the destinations by measures. Then, minimizing the transport cost is equivalent to minimizing the integral over a usually bivariate function representing the cost of transporting a unit of mass from one point to another with respect to the set of all couplings or so-called transport plans which have the specified measures as marginals. The problem was originally introduced by Monge \cite{monge1781memoire} in 1781 and then refined by Kantorovich \cite{kantorovich1948problem, kantorovich2006problem} in 1948. A great variety of researchers considered the optimal transport problem and in the course of their research many important results on optimal transport were established, see for example  \cite{rachev1998mass, rachev1998application} or  \cite{villani2008optimal} for excellent monographes on the topic.

  Observing that there is an analogy between model-independent finance and optimal transport, as in both areas the marginals of the distribution over which some function is optimized are specified, \cite{beiglbock2013model} introduced a new research field that we refer to as martingale optimal transport. Reinterpreting the transport cost function as the payoff function of an exotic option and implementing the usual martingale condition of mathematical finance, the minimization problem of optimal transport cost evolves to the lower price bound problem of model-independent finance. Also, properly implementing the martingale condition in the dual problem of optimal transport, a pricing-hedging duality is shown using only the usual assumptions of model-independent finance, no-arbitrage and consistency with call options prices.

In this paper we consider the martingale optimal transport with discrete margins. This is in contrast to previous papers, see e.g. \cite{henry2016explicit} where continuous margins are  assumed. Moreover, in the discrete case it is possible to solve the problem which reduces to a linear program  with simple techniques. We do not rely on any previous result in this direction and present a stand-alone work. We are also able to generalize some notions like the martingale Spence Mirrlees condition which appears in \cite{henry2016explicit}. For payoff functions with this property we prove the optimality of a left-monotone transport plan, show its uniqueness and present an algorithm to compute the optimal solution. Under the martingale Spence Mirrlees condition the optimality of a left-monotone transport plan has for continuous marginals been shown in \cite{henry2016explicit} and for general marginals in \cite{beiglbock2017monotone}. Also note that the recent paper \cite{hobno18} considers the construction of optimal left-monotone transport plans for general margins. However, their results are far less explicit. The duality with the superhedging problem is trivial in the discrete case.

Our paper is organized as follows: After a mathematical rigorous introduction of the problem and the presentation of some well-known facts about the convex order, we present the linear programming formulation in Section \ref{sec:model}. In Section \ref{sec:spence} we define the generalized Spence Mirrlees condition and prove the optimality of the left-monotone transport plan. In the next section we present an algorithm which solves the primal problem. In Section \ref{sec:dual} we consider an algorithm which solves the dual problem. The last section contains an example.

\section{Preliminary Results and Original Problem}
Consider a financial market with one risk-free asset and one risky asset.
We allow only for two trading times which we denote by $0<t<T$.  The risk-free asset has no interest and its price is given by by $B=(B_t,B_T)=(1,1)$ with $B_0=1$ and for the risky asset we write $S=(S_t,S_T)=(X,Y)$ with $S_0 =s_0 \in \mathbb{R}$. Throughout we will assume that $X$ and $Y$ are discrete random variables and we consider a canonical construction, i.e. for the probability space we may choose $\Omega:=\Omega_1\times \Omega_2:= \{x_1,\ldots,x_N\}\times \{y_1,\ldots ,y_M\}$ and $X(\omega)=X(x,y)=x$ and $Y(\omega)=Y(x,y)=y.$ The $\sigma$-algebra is given by the power set. In what follows we denote by $\mathcal{P}(\Omega)$ the set of all probability measures on $\Omega$. Finally, we denote by $C_t(k)$ the price of a call-option with strike $k$ and maturity $t$, i.e. the payoff is given by $(X-k)^{+}$ and by $C_T(k)$ the price of a call-option with strike $k$ and maturity $T$, i.e. the payoff is given by $(Y-k)^{+}$. Calls are frequently traded and we assume that these prices can be observed at the market. The question in model-independent finance now is to find bounds on consistent prices for other derivatives on this market. A general payoff function of such a derivative in our framework is given by
  \[
       c:\Omega\to \mathbb{R},\quad (x,y) \mapsto c(x,y),
  \]
  i.e.\ the payoff depends only on $S_t$ and $S_T$.
The no-arbitrage assumption generally implies the existence of a probability measure $\mathbb{Q}$ such that the price process $(X,Y)$ is a martingale under $\mathbb{Q}$ and the risk-neutral pricing formula holds, i.e. the price of a derivative with payoff $c$ is
\begin{equation}\label{eq:priceformula}
\mathbb{E}_\mathbb{Q}[c(X,Y)].
\end{equation}  An interesting result in this direction is given in \cite{bl} where it is shown that the observation of call prices $C_t(k)$  for all strikes $k$ implies the knowledge of the marginal at time $t$ of the consistent pricing measure $\mathbb{Q}$. More precisely the following is true.

  \begin{lemma}[{\cite[Sec. 2]{bl}}]\label{Lem Breeden Litzenbeger}
    Let $\mathbb{Q} \in \mathcal{P}(\Omega)$ be consistent with the price functions of call options, i.e. for $t,T$ and all $k \in \mathbb{R}$, we have
    \begin{eqnarray*}\label{Lem Breeden Litzenberger Consitency}
      C_t(k)       = \int_{\mathbb{R}}(x-k)_+ \mathbb{Q}(d (x,y)),\\
       C_T(k)       = \int_{\mathbb{R}}(y-k)_+ \mathbb{Q}(d (x,y)).\
    \end{eqnarray*}
    Then we have
    \begin{eqnarray*}\label{Lem Breeden Litzenberger Distribution}
	  \mathbb{Q}(X\leq k)=1+C_t'(k+),\\
	  \mathbb{Q}(Y\leq k)=1+C_T'(k+)
	\end{eqnarray*}
	for the distribution function of $X$ and $Y$ under $\mathbb{Q}$, where $C_t'(\cdot +)$ denotes the right side derivative of $C_t$ and likewise for $C_T$.
  \end{lemma}

Thus for pricing further derivatives we assume that the marginal distributions of the risky asset price process are known under the consistent pricing measure. But of course this does not imply the joint distribution $\mathbb{Q}$ which is necessary to compute general prices in \eqref{eq:priceformula}. However we have the further information that $(X,Y)$ has to be a martinagle under $\mathbb{Q}$ which restricts the distributions.

  \begin{definition}\label{Def MTP}
    Let $\mu,\nu$ be probability measures on $\Omega_1$  and  $\Omega_2$ respectively. The elements of the set
    \begin{eqnarray*
      \mathcal{M}_2(\mu,\nu):=
       \left\{ \mathbb{Q} \in \mathcal{P}(\Omega)\mid \mu(A) = \mathbb{Q}(A\times \Omega_2), \nu(B) = \mathbb{Q}(\Omega_1\times B), \mathbb{E}_{\mathbb{Q}}[Y\mid X] = X  \ \mathbb{Q}-\text{a.s.}  \right\}
    \end{eqnarray*}
     are called \textit{martingale transport plans} or \textit{potential pricing measures}.
  \end{definition}

Hence $\mathcal{M}_2(\mu,\nu)$  is the set of all possible consistent pricing measures when we assume that $\mu$ is the distribution of $X$ and $\nu$ is the distribution of $Y$.
 Then we define for a given derivative with payoff $c$ the upper price bound problem
  \begin{equation}\label{Eq Upper Price Bound Standard Case}
    P(\mu,\nu):= \sup_{\mathbb{Q} \in \mathcal{M}_2(\mu,\nu)}\mathbb{E}_\mathbb{Q}[c(X,Y)].
    \end{equation}
    and the lower price bound problem
    \begin{equation}\label{Eq Lower Price Bound Standard Case}
   \underline P(\mu,\nu):= \inf_{\mathbb{Q} \in \mathcal{M}_2(\mu,\nu)}\mathbb{E}_\mathbb{Q}[c(X,Y)].
 \end{equation}
 In what follows we will concentrate on problem \eqref{Eq Upper Price Bound Standard Case}. Problem \eqref{Eq Lower Price Bound Standard Case} is in some sense symmetric (see Remark \ref{rem:lower}).
We will later see that in our setting this problem reduces to a simple linear program which can be solved efficiently. However in case $c$ satisfies a certain property, the solution of this linear program is given by a special structure and we will provide an algorithm for the solution.

In what follows when we consider measures $\mu$  on $\mathbb{R}$, we always assume that they are finite, i.e.   $\mu(\mathbb{R})<\infty$ and that the integral exists, i.e.  $\int |x|\mu(dx)<\infty$. Also the next definition of the convex order is introduced for measures (see e.g. \cite{s51}) and not only probability measures as it is often done.

 \begin{definition}\label{Def Convex Order}
    Two measures $\mu,\nu $ on $\mathbb{R}$ are said to be in \textit{convex order}, denoted by $\mu\leq_c\nu$, if for any convex function $f:\mathbb{R} \to \mathbb{R}$ such that the integrals exist,
    \[
      \int_{\mathbb{R}} f(x) \mu( d x) \leq \int_{\mathbb{R}} f(x) \nu( d x).
    \]
      \end{definition}
Since both $f(x)=x$ and $f(x)=-x$ are convex as well as $f(x)=1$ and $f(x)=-1$, the property $\mu\leq_c\nu$ implies that $\int x\mu(dx)=\int x\nu(dx)$ and $\mu(\mathbb{R})=\nu(\mathbb{R})$.

  \begin{definition}\label{Def Call Option Price Func}
    For a measure  $\mu $ on $\mathbb{R}$ the corresponding \textit{call option price function}  is defined by
    \[
          C_\mu:\mathbb{R} \to \mathbb{R}_+, \quad  k \mapsto \int_{\mathbb{R}}(x-k)_{+} \mu( d x).
        \]
  \end{definition}

In arbitrage-free markets it is well-known that call option prices increase with the maturity, i.e. we have $C_\mu(k) \le C_\nu(k)$. Recall here that  $\mu$ is the distribution of $X$ and $\nu$ is the distribution of $Y$ for all $k$, otherwise arbitrage opportunities exist. Now we obtain the following relation which can be found in \cite{chong74}:

  \begin{lemma}\label{Prop Char Conv Ord MultiMarginal}
    Let $\mu, \nu$ be measures on $\mathbb{R}$ with $\frac{1}{\mu(\mathbb{R})}\int x\mu(dx)= \frac{1}{\nu(\mathbb{R})}\int x\nu(dx)$. Then the following are equivalent.
    \begin{enumerate}
          \item[1.] $\mu \leq_c  \nu$.
          \item[2.] $C_{\mu} \leq  C_{\nu}$.
    \end{enumerate}
      \end{lemma}

Thus the condition $\mu \leq_c  \nu$ is a natural assumption in our setting, since $\mu$ and $\nu$ are probability measures and hence $\mu(\mathbb{R})=\nu(\mathbb{R})=1$, $(X,Y)$ is a martingale and thus has constant expectation $\int x\mu(dx)= \int x\nu(dx)$ and due to no arbitrage requirements we have that $C_{\mu} \leq  C_{\nu}$.  Hence Lemma \ref{Prop Char Conv Ord MultiMarginal} implies the convex order between the marginal distributions $\mu$ and $\nu$.

Problem \eqref{Eq Upper Price Bound Standard Case} has a dual problem which is given by (see e.g. \cite{beiglbock2013model})
  \begin{align}\label{Eq Super Hedging Standard Case}
 & \inf_{(\varphi,\psi,h) \in \mathcal{D}_{2}}\left\{ \int_{\mathbb{R}}\varphi(x)\mu(d x) + \int_{\mathbb{R}}\psi(y)\nu(d y) \right\}
    \\ \notag
    &=\inf_{(\varphi,\psi,h) \in \mathcal{D}_{2}}\left\{ \mathbb{E}_{\mu}[\varphi(X)]+\mathbb{E}_{\nu}[\psi(Y)] \right\},
  \end{align}
  where
\begin{align*}
    \mathcal{D}_{2}:= \{(\varphi,\psi,h) \;|\; &\varphi,h :\Omega_1 \to \mathbb{R}, \psi : \Omega_2\to \mathbb{R},
    \varphi(x) + \psi(y)+ h(x)(y-x) \geq c(x,y), \ (x,y) \in \Omega \}.
  \end{align*}
 This problem may be interpreted as finding the cheapest super-replication strategy for the payoff $c$.    Hedging strategies of this form are called semi-static hedging strategies, as $\varphi$ and $\psi$ may be interpreted as static investments in European options with maturity $t$ and $T$ respectively, while $h$ may be understood as a dynamic investment in the underlying asset.   Again we will see that this problem reduces to a linear program and it is indeed the dual linear program to \eqref{Eq Upper Price Bound Standard Case}.

\section{The Problem with Discrete Marginals}\label{sec:model}
In what follows we assume that $\mu$ and $\nu $ are discrete probability distributions, i.e. there are $N,M \in \mathbb{N} $ such that
    \[
        \mu=\sum_{j=1}^{N}\omega_{j}\delta_{x_j} \quad \text{ and }\quad  \nu=\sum_{i=1}^{M}\vartheta_i\delta_{y_i},
    \]
    where $\omega_j, \vartheta_i \geq 0 $, $x_j,y_i \in \mathbb{R}$ for all $j=1,\ldots,N $ and all $i=1,\ldots , M$ and $\sum_{j=1}^{N}\omega_j=\sum_{i=1}^{M}\vartheta_i=1$. $\delta_x$ is the Dirac measure in point $x$. Moreover we assume that $\mu \leq_c \nu$.

 Under this assumptions, martingale transport plans $\mathbb{Q} \in \mathcal{M}_{2}(\mu,\nu)$ are of the form
  \[
      \mathbb{Q}=\sum_{j=1}^{N}\sum_{i=1}^{M} q_{j,i}\delta_{(x_j,y_i)},
  \]
  where the following additional constraints have to be satisfied.
  \begin{enumerate}
       \item The masses of $\mathbb{Q}$ are non-negative, i.e. we have $q_{j,i} \geq 0$, for all $j=1,\ldots, N $ and all $i=1,\ldots,M$.
      \item The marginal distributions of $\mathbb{Q}$ are $\mu$ and $\nu$, i.e. we have
           \begin{align*}
                       \sum_{i=1}^{M}q_{j,i}=\omega_j, \quad   j=1,\ldots,N,
                       \\
                       \sum_{j=1}^{N} q_{j,i}=\vartheta_i, \quad    i=1,\ldots, M.
                     \end{align*}
          This implies $\sum_{j=1}^{N}\sum_{i=1}^{M}q_{j,i}=1$ such that $\mathbb{Q}$ is indeed a probability measure.
      \item The measure $\mathbb{Q}$ satisfies the martingale condition. Transferring the classic condition $\mathbb{E}_{\mathbb{Q}}\left[Y\mid X \right]=X$ to the discrete situation, we have \[
            \sum_{i=1}^{M}\frac{q_{j,i}}{\omega_j}y_i=x_j, \quad  j = 1 ,\ldots, N,
          \]
          as $\frac{q_{j,i}}{\omega_j} = \frac{\mathbb{Q}((X,Y)=(x_j,y_i))}{\mathbb{Q}(X=x_j)}$ is the  conditional distribution $\mathbb{Q}(Y=y_i|X=x_j)$. Rewriting this as
           \[
                       \sum_{i=1}^{M}q_{j,i}y_i=\omega_jx_j\iff \sum_{i=1}^{M}q_{j,i}y_i - \sum_{i=1}^{M}q_{j,i}x_j = 0 \iff \sum_{i=1}^{M} q_{j,i}(y_i-x_j)=0,
                     \]
                     we find an alternative condition.
       \end{enumerate}

   Altogether, the upper price bound problem in \eqref{Eq Upper Price Bound Standard Case} reduces in the discrete case to the following  linear program
  \begin{align}\label{Prob Upper Price Bound Discrete}\tag{P}
    \max \quad   &\sum_{j=1}^{N} \sum_{i=1}^{M}q_{j,i}c(x_j,y_i)  := \sum_{j=1}^{N} \sum_{i=1}^{M}q_{j,i}c_{j,i}
    \\ \notag
    \text{s.t.} \quad  &\sum_{i=1}^{M}q_{j,i}=\omega_j, \quad  j=1,\ldots, N ,
    \\ \notag
   &\sum_{j=1}^{N}q_{j,i}=\vartheta_i, \quad  i=1,\ldots,M,
    \\ \notag
    & \sum_{i=1}^{M} q_{j,i}(y_i-x_j)=0, \quad j =1 ,\ldots, N,
    \\ \notag
    &   q_{j,i} \geq 0, \quad j=1,\ldots,N, i=1,\ldots, M.
 \end{align}
 The first two equations guarantee that we have the correct  marginal distributions, the last equation is the martingale condition.
The dual problem is given by
  \begin{align}\label{Prob Super Hedging Discrete}\tag{D}
    \min \quad &\sum_{j=1}^{N}\omega_j \varphi_j + \sum_{i=1}^{M} \vartheta_i \psi_i
    \\
    \notag \text{s.t.}\quad  &\varphi_j + \psi_i + h_j (y_i-x_j)\geq c_{j,i},\quad  j=1,\ldots,N, i=1,\ldots,M,
    \\
    \notag& \varphi_j,h_j,\psi_i \in \mathbb{R}, \quad  j=1,\ldots,N, i=1,\ldots,M,
  \end{align}
and can be interpreted as  a super hedging problem.

The next result is due to \cite{s51} and has later been generalized  in various ways (see also \cite{str}).

\begin{lemma}\label{lem:sherman}
It holds that $\mu \leq_c \nu$ if and only if there exists $q_{j,i} \geq 0$ such that
  \begin{align}
    &\sum_{i=1}^{M}q_{j,i}=\omega_j, \quad  j=1,\ldots, N ,
    \\ \notag
   &\sum_{j=1}^{N}q_{j,i}=\vartheta_i, \quad  i=1,\ldots,M,
    \\ \notag
    & \sum_{i=1}^{M} q_{j,i}(y_i-x_j)=0, \quad j =1 ,\ldots, N.
 \end{align}

\end{lemma}

Thus, we immediately obtain the following existence result.

\begin{theorem}\label{theo:existence}
The linear programs (P) and (D) have optimal solutions $q^*$ and $(\varphi^*,h^*,\psi^*)$ and $\sum_{i=1}^{M}q^*_{j,i}c_{j,i} = \sum_{j=1}^{N}\omega_j \varphi^*_j + \sum_{i=1}^{M} \vartheta_i \psi^*_i$, i.e. the optimal value of the objective functions coincide.
\end{theorem}

\begin{proof}
Lemma \ref{lem:sherman} implies that (P) has feasible points. Since $q_{j,i} \in [0,1]$ the target function is bounded and thus by the classical existence theorem for linear programs an optimal solution for (P) exists. But by the strong duality theorem of linear programming this implies that (D) has a solution and the optimal values of the objective functions coincide.
\end{proof}

\section{Left-Monotonicity and a kind of Martingale Spence Mirrlees Condition}\label{sec:spence}
We will consider the linear program (P) with functions $c$ having a special property.
In this case it is possible to prove that a solution $\mathbb{Q}^*=q^*$ of (P) has a certain property.
The following definition is due to \cite{beiglboeck2016problem}.

 \begin{definition}\label{def:leftmo}
  A martingale transport plan ${\mathbb{Q} \in \mathcal{M}_2(\mu,\nu)}$ is called \textit{left monotone}, if there is a set $\Gamma \subseteq \mathrm{supp}(\mu)\times \mathrm{supp}(\nu)$ with ${\mathbb{Q}( \Gamma)=1}$ and such that for $(x,y^-),(x,y^+),(x',y')\in \Gamma$ with $x<x'$, we have $y' \notin (y^-,y^+)$.
  \end{definition}

    \begin{figure}[!htbp]\begin{center}

       \begin{tikzpicture}[scale=1]

          \node  at (0,0.4) {$x$};
          \node at (2,0.4) {$x'$};
         \fill (0,0) circle (0.05);
         \fill (2,0) circle (0.05);

          \node at (-2,-3.4) {$y^{-}$};
      \node at (0,-3.4) {$y'$};
          \node at (2,-3.4) {$y^{+}$};
           \fill (-2,-3) circle (0.05);
                  \fill (0,-3) circle (0.05);
                  \fill (2,-3) circle (0.05);

                  \draw(0,0)--(-2,-3);
                     \draw(0,0)--(2,-3);
                     \draw(2,0)--(0,-3);

        \end{tikzpicture}
        \caption{Forbidden configuration for left monotonicity.}
  \end{center}
    \end{figure}

\begin{theorem}\label{theo:lmopt}
Let $c:\mathbb{R}^2 \to \mathbb{R}$ be a payoff function such that for all $x'>x$ and $y^{+}>y'>y^{-}$ with $x,x' \in \mathrm{supp}(\mu)$ and $ y^{-}, y',y^{+} \in \mathrm{supp}(\nu)$ we have
   \begin{equation}\label{eq:optcondition}
      \lambda\left[ c(x',y^{+})-c(x,y^{+}) \right] + (1-\lambda)\left[ c(x',y^{-})-c(x,y^{-}) \right] - \left[ c(x',y')-c(x,y') \right] > 0,
    \end{equation}
    where $\lambda = \frac{y'-y^{-}}{y^{+}-y^{-}}\in[0,1]$. Then an optimal $q^*=\mathbb{Q}^*$ of (P) is left-monotone. The (partial) converse is also true: If  there exists a unique solution of (P) which is left-monotone, then $c$ has to satisfy \eqref{eq:optcondition}.

\end{theorem}

\begin{proof}
Suppose $x,x' \in \mathrm{supp}(\mu)$ and $ y^{-}, y',y^{+} \in \mathrm{supp}(\nu)$ are such that $x'>x$ and $y^{+}>y'>y^{-}$ and suppose that ${\mathbb{Q} \in \mathcal{M}_2(\mu,\nu)}$ satisfies
\begin{eqnarray*}
\mathbb{Q}(x,y^-) := \theta^->0\\
\mathbb{Q}(x,y^+) := \theta^+>0\\
\mathbb{Q}(x',y') := \theta'>0,
\end{eqnarray*}
i.e. the condition in Definition \ref{def:leftmo} is not satisfied.
Let $\lambda = \frac{y'-y^{-}}{y^{+}-y^{-}} $ which implies that $y'=\lambda y^++(1-\lambda)y^-$.  We will now define a new $\tilde{\mathbb{Q}}\in \mathcal{M}_2(\mu,\nu)$ such that the objective function of (P) attains a higher value which implies that $\mathbb{Q}$ cannot be optimal.
We consider the following new arrangement $\tilde{\mathbb{Q}}$ of the probability mass
\begin{eqnarray*}
\tilde{\mathbb{Q}}(x',y^-) &:=& \mathbb{Q}(x',y^-)+\theta'(1-\lambda)\\
\tilde{\mathbb{Q}}(x',y^+) &:=& \mathbb{Q}(x',y^+)+\theta'\lambda\\
\tilde{\mathbb{Q}}(x',y') &:=& 0\\
\tilde{\mathbb{Q}}(x,y^-) &:=& \theta^--\theta'(1-\lambda)\\
\tilde{\mathbb{Q}}(x,y^+) &:=& \theta^+-\theta'\lambda\\
\tilde{\mathbb{Q}}(x,y') &:=& \mathbb{Q}(x,y')+\theta'
\end{eqnarray*}
All other assignments are left unchanged.
In a first step we assume here that $\theta^--\theta'(1-\lambda)\ge0$ and $\theta^+-\theta'\lambda\ge0$ and claim that $\tilde{\mathbb{Q}}$ is again admissible for (P). That the marginal distributions are preserved is easy to see.  In order to show that the martingale condition holds we have to show that the values $\sum_i \mathbb{Q}(x,y_i)y_i$ and $\sum_i \mathbb{Q}(x',y_i)y_i$ do not change under the new probability assignment. First consider $x$. Here we obtain
\begin{eqnarray*}
\sum_i \mathbb{Q}(x,y_i)y_i-\sum_i \tilde{\mathbb{Q}}(x,y_i)y_i &=&
y^- \theta^- + y^+ \theta^+ -\big( y^- (\theta^--\theta'(1-\lambda)) + y^+(\theta^+-\theta'\lambda)+y'\theta'\big)\\
&=&  y^- \theta^- + y^+ \theta^+ -\big(y^- \theta^-+ y^+\theta^++\theta'(y'-y^-(1-\lambda)-y^+\lambda) \big)\\
&=&0
\end{eqnarray*}
which is true since the term in the inner brackets vanishes due to the definition of $\lambda$. Next consider $x'$. Here we obtain:
\begin{eqnarray*}
\sum_i \mathbb{Q}(x',y_i)y_i-\sum_i \tilde{\mathbb{Q}}(x',y_i)y_i &=&
y' \theta'-\big( y^- \theta' (1-\lambda)  + y^+\theta'\lambda\big) \\
&=& \theta' \big[ y'-(y^-(1-\lambda)+y^+\lambda)\big]=0
\end{eqnarray*}
which is again true by the definition of $\lambda$. Thus we have that $\tilde{\mathbb{Q}}\in \mathcal{M}_2(\mu,\nu)$. Now we finally prove that the value of the objective function under $\tilde{\mathbb{Q}}$ is larger than under $\mathbb{Q}$. We obtain for the difference of the objective functions:
\begin{eqnarray*}
\sum_{j=1}^{N} \sum_{i=1}^{M}\big(\tilde{\mathbb{Q}}(x_j,y_i)-\mathbb{Q}(x_j,y_i)\big)c(x_j,y_i) &=& c(x',y^-) \theta'(1-\lambda)+ c(x',y^+) \theta'\lambda-c(x',y')\theta'\\[-0.4cm]
&& -c(x,y^-) \theta'(1-\lambda)-c(x,y^+) \theta'\lambda+c(x,y')\theta'\\
&=& \lambda\theta' \left[ c(x',y^{+})-c(x,y^{+}) \right] \\&&+ (1-\lambda)\theta'\left[ c(x',y^{-})-c(x,y^{-}) \right] \\
&&- \theta'\left[ c(x',y')-c(x,y') \right]> 0
\end{eqnarray*}
by our assumption. Also note that the new assignment dissolves the forbidden configuration. In case the conditions $\theta^--\theta'(1-\lambda)\ge0$ and $\theta^+-\theta'\lambda\ge0$ are not satisfied, we keep part of the probability mass $\theta'$ on $(x',y')$ and shift only the amount $\tilde{\theta}$ such that $\min\{\theta^--\tilde{\theta}(1-\lambda), \theta^+-\tilde{\theta}\lambda \}=0$ in the same way as before. In this case we dissolve  again the forbidden configuration and obtain a higher value for the objective function using the same arguments as before. From the proof so far, we see that the converse statement is also true: Suppose there are $x,x' \in \mathrm{supp}(\mu)$ and $ y^{-}, y',y^{+} \in \mathrm{supp}(\nu)$  such that $x'>x$ and $y^{+}>y'>y^{-}$ and $c$ does not satisfy \eqref{eq:optcondition} on these points. Then inverting the transport that we have constructed in the first part we see that $\mathbb{Q}$ is not worse than $\tilde{\mathbb{Q}}$ which implies that $\tilde{\mathbb{Q}}$ cannot be the unique solution.
\end{proof}

The following notion has been introduced in \cite{henry2016explicit} and it has been shown there that this property of $c$ implies that the optimal transport plan is left-monotone in the case of continuous margins. The same has been shown in \cite{beiglbock2017monotone} for general margins. In the next lemma we show that this condition implies \eqref{eq:optcondition}. Of course \eqref{eq:optcondition} does not need a differentiable function $c$ and is thus more general.

  \begin{definition}\label{Def Martingale Spence Mirrlees}
    A function $c:\mathbb{R}^2 \to \mathbb{R}$ satisfies the \textit{martingale Spence Mirrlees condition}, if the partial derivative $c_{xyy}$ exists and satisfies $c_{xyy}>0$.
  \end{definition}

 \begin{lemma}
    Let $c:\mathbb{R}^2 \to \mathbb{R}$ be a function satisfying the martingale Spence Mirrlees condition. Then \eqref{eq:optcondition} holds for all $x'>x$, $y^{+}>y'>y^{-}$.
  \end{lemma}

  \begin{proof}
 First note that
    \[
      \lambda\left[ c(x',y^{+})-c(x,y^{+}) \right] + (1-\lambda)\left[ c(x',y^{-})-c(x,y^{-}) \right] - \left[ c(x',y')-c(x,y') \right] > 0,
    \]
         is  equivalent to
    \begin{align*}
      &\lambda\left[ c(x',y^{+})-c(x',y')-c(x,y^{+})+c(x,y') \right]
      \\
      &-  (1-\lambda)\left[ c(x',y') - c(x',y^{-})-c(x,y')+c(x,y^{-}) \right]>0.
    \end{align*}

    If we plug in $\lambda = \frac{y'-y^{-}}{y^{+}-y^{-}} $ and multiply by $y^{+}-y^{-}$, we obtain
    \begin{align}\nonumber
    &    {[c(x',y^{+})-c(x',y')-c(x,y^{+})+c(x,y')]}(y'-y^{-}) \\ \label{eq:msmc}
       -&  {[c(x',y') - c(x',y^{-})-c(x,y')+c(x,y^{-})]}(y^{+}-y')>0.
     \end{align}
Since $c_{xyy}>0$ and $s \ge y' \ge u$ for all $s \in [y',y^+]$ and $u \in [y^-,y']$ we obtain:
\begin{align*}
  0 &< \int_{x}^{x'}\int_{y'}^{y^{+}}\int_{y^{-}}^{y'}\int_{u}^{s} c_{xyy}(t,v) \mathrm{d}v\mathrm{d}u\mathrm{d}s\mathrm{d}t
  \\&= \int_{x}^{x'}\int_{y'}^{y^{+}}\int_{y^{-}}^{y'} c_{xy}(t,s)-c_{xy}(t,u) \mathrm{d}u\mathrm{d}s\mathrm{d}t
  \\&= \int_{x}^{x'}\int_{y'}^{y^{+}} c_{xy}(t,s)(y'-y^{-})-(c_{x}(t,y')-c_{x}(t,y^{-})) \mathrm{d}s\mathrm{d}t
  \\&= \int_{x}^{x'} (c_{x}(t,y^{+})-c_{x}(t,y'))(y'-y^{-})-(c_{x}(t,y')-c_{x}(t,y^{-}))(y^{+}-y') \mathrm{d}t
  \\&=  [c(x',y^{+})-c(x,y^{+})-(c(x',y')-c(x,y'))](y'-y^{-})
    \\&\quad    - [c(x',y')-c(x,y') - (c(x',y^{-})-c(x,y^{-}))](y^{+}-y').
\end{align*}
Comparing it with \eqref{eq:msmc} yields the assertion.
  \end{proof}

\begin{corollary}
If $c$ satifies the martingale Spence Mirrlees condition, then the optimal $q^*=\mathbb{Q}^*$ of (P) is left-monotone.
\end{corollary}

\begin{remark}\label{rem:lower}
As far as the lower price bound in problem \eqref{Eq Lower Price Bound Standard Case} is concerned, we obtain a similar result: If the martingale Spence Mirrlees condition is satisfied, then a right-monotone transport plan is optimal.  A $\mathbb{Q} \in \mathcal{M}_2(\mu,\nu)$ is called \textit{right monotone}, if there is a Borel set $\Gamma  \subset
\mathrm{supp}(\mu)\times \mathrm{supp}(\nu) $ with $\mathbb{Q}( \Gamma)=1$ and such that for all $(x,y_1),(x,y_2),(x',y')\in \Gamma$ with $x>x'$, we have $y' \notin (y_1,y_2)$. In case $c_{xyy}<0$, the optimality properties reverse, i.e. the right monotone transport plan is optimal for \eqref{Eq Upper Price Bound Standard Case} and the left monotone transport plan for \eqref{Eq Lower Price Bound Standard Case}.
\end{remark}

\section{An Algorithm for the Construction of the Optimal Solution of (P)}\label{sec:algo}
In this section we derive an algorithm which computes an optimal solution of (P). We assume that $\mu,\nu$ satisfy the conditions $\mu \leq_c \nu$, $\mu= \sum_{j=1}^{N}\omega_j \delta_{x_j}$ and $\nu=\sum_{i=1}^{M}\vartheta_i\delta_{y_i}$ and the function $c$ has property \eqref{eq:optcondition}. Theorem \ref{theo:lmopt} and Theorem \ref{theo:existence} imply that a left-monotone transport plan exists. We denote by $\mathbb{Q}_l(\mu,\nu)$ any such left-monotone transport plan.

Before we formalize the algorithm to determine $\mathbb{Q}_l(\mu,\nu)$, let us develop an intuition how it could work. As $\mathbb{Q}_l(\mu,\nu)$ is in particular a coupling of its margins $\mu$ and $\nu$, we have to specify how much of the mass of each $x_j$ is transported to each $y_i$. Of course, this transport specification has to satisfy some prerequisites. Clearly, the coupling can be done for one atom of $\mu$ at a time. For this purpose, by (a)-(c) in Section \ref{sec:model}, it is obvious how an algorithm should proceed in order to satisfy the marginal and martingale conditions. Thus, the central idea stems from the desired left monotonicity. In order not to introduce a contradiction to left monotonicity when coupling an atom of $\mu$ with atoms of $\nu$, we should couple the smallest atom of $\mu$ martingale-consistently with atoms of $\nu$ such that no mass of $\nu$ is left over in the convex hull of those atoms. This does indeed guarantee that in a next step, coupling the second smallest atom of $\mu$, no contradictions to the left monotonicity can be introduced.

In order to get a non-trivial problem we assume now that $\mu$ and $\nu$ satisfy the following conditions.
    \begin{enumerate}
       \item $x_1<\ldots <x_N$, $y_1< \ldots < y_M$.
      \item $\mu \neq \nu$.
      \item $\mu$ has at least two different atoms of positive mass.
    \end{enumerate}

First we would like to note the following obvious fact.

\begin{lemma}\label{lem:cximpl}
The assumption $\mu\le_c\nu$ implies that for all $x_j\in \mathrm{supp}(\mu)$ there either exists $\ell \in \{1,\ldots,M\}$ such that $x_j=y_\ell$ or there exists an $\ell \in \{1,\ldots,M-1\}$ such that $y_\ell<x_j<y_{\ell+1}$.
\end{lemma}

\begin{proof}
This property follows directly from the martingale condition:
$$\sum_{i=1}^M q_{j,i}(y_i-x_j)=0$$
because the atoms of $\nu$ cannot all lay on one side of $x_j$.
\end{proof}

The next two properties have to be satisfied by any left-monotone transport plan.

  \begin{lemma}\label{lem:c1}
    Suppose $x_1 \in \mathrm{supp}(\nu)$, say $x_1=y_\ell$ for some $\ell \in \{1,\ldots,M\}.$
    \begin{enumerate}
      \item Then $(x_1,y_\ell)$ is an atom of $\mathbb{Q}_{l}(\mu,\nu)$.
      \item Further suppose $\ell =1$. Then $\omega_1 \leq \vartheta_1$.
    \end{enumerate}
  \end{lemma}

  \begin{proof}
    \begin{enumerate}
      \item In order to get a contradiction, assume that $x_1$ is not coupled with $y_\ell$ under $\mathbb{Q}_{l}(\mu,\nu)$. Then there are $y^-,y^+ \in \mathrm{supp}(\nu)$ with $y^-<y_\ell < y^+$ and such that $x_1$ is coupled with $y^-$ and $y^+$. Also, there is some $x' \in \mathrm{supp}(\mu)$ with $x'>x_1$ and such that $x'$ is coupled with $y_{\ell}$. This contradicts the left monotonicity.
      \item In order to get a contradiction, assume that $\omega_1 >\vartheta_1$. After coupling $x_1$ with $y_1$ there is mass of at least $\omega_1 - \vartheta_1$ left in $x_1$. This has then to be coupled with some atoms from $\mathrm{supp}(\nu)\setminus \{y_1\}$. As $y_i > x_1$ for all $i=2,\ldots,M$, this contradicts the martingale property.
    \end{enumerate}
      \end{proof}

 \begin{lemma}\label{lem:c2}
 Suppose $y_{\ell}<x_1 < y_{\ell+1}$ for some $\ell \in \{1,\ldots,M-1\}$. Then $(x_1,y_\ell)$ and $(x_1,y_{\ell +1})$ are atoms of $\mathbb{Q}_{l}(\mu,\nu)$.
 \end{lemma}

 \begin{proof}
   In order to get a contradiction, assume without loss of generality that $x_1$ is not coupled with $y_{\ell+1}$. Then there is a $y^+ \in \mathrm{supp}(\nu)$ with $y^+ > y_{\ell +1}$ and such that $x_1$ is coupled with at least $y_{\ell}$ and $y^+$. Also, there is an $x' \in \mathrm{supp}(\mu)$ with $x' >x_1$ and such that $x'$ is coupled with $y_{\ell+1}$. This contradicts the left monotonicity.
 \end{proof}

Before we start with the algorithm to construct a left-monotone transport plan we show that it is unique. This has already been shown in \cite{beiglboeck2016problem} for general margins using  rather complex methods. In contrast, our proof is very easy and does not rely on sophisticated concepts.

\begin{theorem}
There exists one unique left-monotone transport plan.
\end{theorem}

\begin{proof}
Suppose  there are at least two different left-monotone transport plans $q, \tilde{q} \in \mathcal{M}_2(\mu,\nu)$. Let $x_{j^*}$ be smallest atom of $\mu$ where there is an $i^{*}\in \{1,\ldots, M\}$ such that $q_{j^*,i^{*}}\neq \tilde{q}_{j^*,i^{*}}.$ We eliminate the variables $q_{j,i}, \tilde{q}_{j,i}$ for $j<j^*$ and $i=1,\ldots,M$ from the transport problem and consider the remaining constraints which for $q$ read:
  \begin{align*}
  &\sum_{i=1}^{M}q_{j,i}=\omega_j, \quad  j=j^*,\ldots, N ,
    \\ \notag
   &\sum_{j=j^*}^{N}q_{j,i}=\vartheta_i - \sum_{j=1}^{j^*-1}q_{j,i} , \quad  i=1,\ldots,M,
    \\ \notag
    & \sum_{i=1}^{M} q_{j,i}(y_i-x_j)=0, \quad j =j^* ,\ldots, N,
    \\ \notag
    &   q_{j,i} \geq 0, \quad j=j^*,\ldots,N, i=1,\ldots, M.
 \end{align*}
The remaining transport plans $q$ and $\tilde{q}$ are again left-monotone. Now consider the atom $x_{j^*}$  with the probability mass $\omega_{j^*}$. By the left-monotonicity $x_{j^*}$  has to be coupled with atoms $y_\ell <y_{\ell+1}<\ldots <y_{\ell+m}$ of $\nu$ such that  $q_{j^*,k}= \vartheta_k$ for $k=\ell+1,\ldots ,\ell+m-1$, $ \sum_{k=0}^m q_{j^*,\ell+k}=\omega_{j^*}$ and $$ x_{j^*} \omega_{j^*} = \sum_{k=0}^m q_{j^*,\ell+k} y_{\ell+k}.$$
However these conditions imply a unique solution. Thus $q$ and $\tilde{q}$ cannot be different.
\end{proof}

Next we write down an algorithm which produces the left-monotone transport plan. Note that Lemma \ref{lem:cximpl} implies that once we fix an $x_j$ then only two situations can occur: either  there exists an $\ell \in \{1,\ldots,M\}$ such that $x_j=y_\ell$ or there exists an $\ell \in \{1,\ldots,M-1\}$ such that $y_\ell<x_j<y_{\ell+1}$. Note that the algorithm works in a recursive way: the mass from $\mu$ to $\nu$ is transported step by step. After the transport we continue with the remaining masses. In particular the index $1$ always refers to the smallest remaining atom of  $\mu$ and thus may be different from the real index.\\

{\bf Algorithm (primal problem):} Start with $x_1$.
\begin{description}
\item[Case I]  There exists $\ell \in \{1,\ldots,M\}$ such that $x_1=y_\ell$
\begin{description}
\item[Case 1] $\omega_1\le \vartheta_\ell$.

Then set $q_{1,\ell}:= \omega_1, q_{1,i}:= 0, i\neq \ell$ and define $ \vartheta_\ell := \vartheta_\ell-\omega_1, \omega_1 := 0$. Cross $x_1$ (and possibly $y_{\ell}$) from the list. All other probability masses are kept unchanged. Continue with the smallest remaining atom of $\mu$.

   \begin{figure}[!htbp]\begin{center}

       \begin{tikzpicture}[scale=1]

               \draw(0,0)--(7,0);
               \draw(0,2)--(7,2);
              \draw(1,2.1)--(1,1.9);
         \node  at (1,1.7) {$x_1$};
          \node at (1,2.3) {$0$};
            \draw(1,0.1)--(1,-0.1);
           \node  at (1,0.3) {$y_\ell$};
            \node at (1,-0.3) {$\vartheta_\ell-\omega_1$};
         { \color{red}{   \draw(1,0)--(1,2);}
            \node  at (1.3,1) {$\omega_1$};}
            \draw(2,2.1)--(2,1.9);
           \node  at (2,1.7) {$x_2$};
          \draw(5,2.1)--(5,1.9);
           \node  at (5,1.7) {$x_N$};
                  \draw(6,0.1)--(6,-0.1);
           \node  at (6,0.3) {$y_M$};

        \end{tikzpicture}
        \caption{Case I, Case 1. The figure shows the new probability masses of the involved atoms.}
  \end{center}
    \end{figure}

 \item[Case 2] $\omega_1> \vartheta_\ell$.

Then set $q_{1,\ell}:= \vartheta_\ell, q_{j,\ell}:= 0, j\neq 1$ and define $\omega_1 := \omega_1-\vartheta_\ell, \vartheta_\ell := 0$. Cross $y_\ell$ from the list. All other probability masses are kept unchanged. Continue with the smallest remaining atom of $\mu$.

      \begin{figure}[!htbp]\begin{center}

       \begin{tikzpicture}[scale=1]

               \draw(0,0)--(7,0);
               \draw(0,2)--(7,2);
              \draw(1,2.1)--(1,1.9);
         \node  at (1,1.7) {$x_1$};
          \node at (1,2.3) {$\omega_1-\vartheta_\ell$};
            \draw(1,0.1)--(1,-0.1);
           \node  at (1,0.3) {$y_\ell$};
            \node at (1,-0.3) {$0$};
         { \color{red}{   \draw(1,0)--(1,2);}
            \node  at (1.3,1) {$\vartheta_\ell$};}
            \draw(2,2.1)--(2,1.9);
           \node  at (2,1.7) {$x_2$};
          \draw(5,2.1)--(5,1.9);
           \node  at (5,1.7) {$x_N$};
                  \draw(6,0.1)--(6,-0.1);
           \node  at (6,0.3) {$y_M$};

        \end{tikzpicture}
        \caption{Case I, Case 2. The figure shows the new probability masses of the involved atoms}
  \end{center}
    \end{figure}

\end{description}

\item[Case II]  There exists $\ell \in \{1,\ldots,M-1\}$ such that $y_\ell<x_1<y_{\ell+1}$.

Then there exist $\vartheta'_\ell, \vartheta'_{\ell+1} \ge 0$ such that
\begin{eqnarray}
\label{eq:leq} \vartheta'_\ell + \vartheta'_{\ell+1} &=& \omega_1\\ \nonumber
 \vartheta'_\ell  y_\ell +\vartheta'_{\ell+1} y_{\ell+1} &=& \omega_1x_1.
\end{eqnarray}
Note that these linear equations have a unique solution.
\begin{description}
 \item[Case 1] $\vartheta'_\ell \le \vartheta_\ell, \vartheta'_{\ell+1}\le \vartheta_{\ell+1}.$

 Then set $q_{1,\ell} := \vartheta'_\ell, q_{1,\ell+1}:= \vartheta'_{\ell+1} , q_{1,i}:=0, i\neq \ell,\ell+1$ and define $\omega_1 := 0, \vartheta_\ell := \vartheta_\ell-\vartheta'_\ell, \vartheta_{\ell+1}:= \vartheta_{\ell+1}-\vartheta'_{\ell+1}.$ Cross $x_1$ from the list (and possibly $y_{\ell}, y_{\ell+1}$). All other probability masses are kept unchanged. Continue with the smallest remaining atom of $\mu$.

      \begin{figure}[!htbp]\begin{center}

       \begin{tikzpicture}[scale=1]

               \draw(-0.5,0)--(7,0);
               \draw(-0.5,2)--(7,2);
              \draw(1,2.1)--(1,1.9);
         \node  at (1,1.7) {$x_1$};
          \node at (1,2.3) {$0$};
            \draw(-0.1,0.1)--(-0.1,-0.1);
           \node  at (-0.1,0.3) {$y_\ell$};
            \node at (-0.1,-0.3) {$\vartheta_\ell-\vartheta'_\ell$};
                  \draw(2.1,0.1)--(2.1,-0.1);
           \node  at (2.1,0.3) {$y_{\ell+1}$};
            \node at (2.1,-0.3) {$\vartheta_{\ell+1}-\vartheta'_{\ell+1}$};

         { \color{red}{   \draw(-0.1,0)--(1,2);}
            \node  at (0,1.1) {$\vartheta'_\ell$};
            \draw(2.1,0)--(1,2);
             \node  at (2,1.1) {$\vartheta'_{\ell+1}$};}
           \node  at (2,1.7) {$x_2$};
          \draw(5,2.1)--(5,1.9);
           \node  at (5,1.7) {$x_N$};
                  \draw(6,0.1)--(6,-0.1);
           \node  at (6,0.3) {$y_M$};

        \end{tikzpicture}
        \caption{Case II, Case 1. The figure shows the new probability masses of the involved atoms}
  \end{center}
    \end{figure}
 \item[Case 2] $\vartheta'_\ell > \vartheta_\ell, \vartheta'_{\ell+1}\le \vartheta_{\ell+1}.$

 In this case we multiply the linear equations of  \eqref{eq:leq} with $\frac{\vartheta_\ell}{\vartheta'_\ell}$ to obtain:
 \begin{eqnarray*}
 \vartheta_\ell + \frac{\vartheta_\ell}{\vartheta'_\ell} \vartheta'_{\ell+1} &=& \omega_1 \frac{\vartheta_\ell}{\vartheta'_\ell}\\
 \vartheta_\ell  y_\ell +\frac{\vartheta_\ell}{\vartheta'_\ell} \vartheta'_{\ell+1} y_{\ell+1} &=& \omega_1 \frac{\vartheta_\ell}{\vartheta'_\ell}x_1.
\end{eqnarray*}

 Then set $q_{1,\ell} := \vartheta_\ell, q_{1,\ell+1}:= \frac{\vartheta_\ell}{\vartheta'_\ell}\vartheta'_{\ell+1} , q_{i,\ell}:=0, i\neq 1$ and define $\omega_1 := \omega_1(1-\frac{\vartheta_\ell}{\vartheta'_\ell}), \vartheta_{\ell+1}:= \vartheta_{\ell+1}-\frac{\vartheta_\ell}{\vartheta'_\ell}\vartheta'_{\ell+1}$, $\vartheta_\ell := 0$. Cross $y_\ell$ from the list. All other probability masses are kept unchanged. Continue with the smallest remaining atom of $\mu$. Note that after this step there is still probability mass on $x_1$ and $y_{\ell+1}$ which could result in the next step in another transport from $x_1$ to $y_{\ell+1}$. At the end these transported masses have to be added.

  \begin{figure}[!htbp]\begin{center}

       \begin{tikzpicture}[scale=1]

               \draw(-0.5,0)--(7,0);
               \draw(-0.5,2)--(7,2);
              \draw(1,2.1)--(1,1.9);
         \node  at (1,1.7) {$x_1$};
          \node at (1,2.3) {$\omega_1(1-\frac{\vartheta_\ell}{\vartheta'_\ell})$};
            \draw(-0.1,0.1)--(-0.1,-0.1);
           \node  at (-0.1,0.3) {$y_\ell$};
            \node at (-0.1,-0.3) {$0$};
                  \draw(2.1,0.1)--(2.1,-0.1);
           \node  at (2.1,0.3) {$y_{\ell+1}$};
            \node at (2.1,-0.5) {$\vartheta_{\ell+1}-\frac{\vartheta_\ell}{\vartheta'_\ell}\vartheta'_{\ell+1}$};

         { \color{red}{   \draw(-0.1,0)--(1,2);}
            \node  at (0,1.1) {$\vartheta_\ell$};
            \draw(2.1,0)--(1,2);
             \node  at (2.2,1.1) {$\frac{\vartheta_\ell}{\vartheta'_\ell}\vartheta'_{\ell+1}$};}
           \node  at (2,1.7) {$x_2$};
          \draw(5,2.1)--(5,1.9);
           \node  at (5,1.7) {$x_N$};
                  \draw(6,0.1)--(6,-0.1);
           \node  at (6,0.3) {$y_M$};

        \end{tikzpicture}
        \caption{Case II, Case 2. The figure shows the new probability masses of the involved atoms}
  \end{center}
    \end{figure}

  \item[Case 3] $\vartheta'_\ell \le \vartheta_\ell, \vartheta'_{\ell+1}> \vartheta_{\ell+1}.$

 In this case we multiply the linear equation \eqref{eq:leq} with $\frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}$ to obtain:
 \begin{eqnarray*}
 \vartheta'_\ell \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}} + \vartheta_{\ell+1} &=& \omega_1 \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}\\
 \vartheta'_\ell \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}  y_\ell + \vartheta_{\ell+1} y_{\ell+1} &=& \omega_1 \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}x_1.
\end{eqnarray*}

 Then set $q_{1,\ell} := \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}} \vartheta'_\ell, q_{1,\ell+1}:=\vartheta_{\ell+1} , q_{i,\ell+1}:=0, i\neq 1$ and define $\omega_1 := \omega_1(1-\frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}), \vartheta_\ell := \vartheta_\ell- \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}\vartheta'_\ell, \vartheta_{\ell+1}:= 0.$ Cross $y_{\ell+1}$ from the list. All other probability masses are kept unchanged. Continue with the smallest remaining atom of $\mu$. Note that after this step there is still probability mass on $x_1$ and $y_{\ell}$ which could result in the next step in another transport from $x_1$ to $y_{\ell}$. At the end these transported masses have to be added.
  \begin{figure}[!htbp]\begin{center}

       \begin{tikzpicture}[scale=1]

               \draw(-0.5,0)--(7,0);
               \draw(-0.5,2)--(7,2);
              \draw(1,2.1)--(1,1.9);
         \node  at (1,1.7) {$x_1$};
          \node at (1,2.3) {$\omega_1(1-\frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}})$};
            \draw(-0.1,0.1)--(-0.1,-0.1);
           \node  at (-0.1,0.3) {$y_\ell$};
            \node at (-0.1,-0.5) {$\vartheta_\ell- \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}\vartheta'_\ell$};
                  \draw(2.1,0.1)--(2.1,-0.1);
           \node  at (2.1,0.3) {$y_{\ell+1}$};
            \node at (2.1,-0.3) {$0$};

         { \color{red}{   \draw(-0.1,0)--(1,2);}
            \node  at (-0.2,1.1) {$\frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}\vartheta'_\ell$};
            \draw(2.1,0)--(1,2);
             \node  at (2.2,1.1) {$\vartheta_{\ell+1}$};}
           \node  at (2,1.7) {$x_2$};
          \draw(5,2.1)--(5,1.9);
           \node  at (5,1.7) {$x_N$};
                  \draw(6,0.1)--(6,-0.1);
           \node  at (6,0.3) {$y_M$};

        \end{tikzpicture}
        \caption{Case II, Case 3. The figure shows the new probability masses of the involved atoms}
  \end{center}
    \end{figure}

   \item[Case 4] $\vartheta'_\ell > \vartheta_\ell, \vartheta'_{\ell+1}> \vartheta_{\ell+1}.$

Define here $k := \min\big\{\frac{\vartheta_{\ell}}{\vartheta'_{\ell}}, \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}} \big\}.$ If $k= \frac{\vartheta_{\ell}}{\vartheta'_{\ell}}$ proceed as in case 2. If $k=\frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}$ proceed as in case 3.
\end{description}
\end{description}

\begin{theorem}
The algorithm produces the left-monotone transport plan in a finite number of steps.
\end{theorem}

\begin{proof}
In each step of the algorithm at least one atom is crossed from the list. In the final step at least two atoms vanish. Thus, the algorithm stops in at most $N+M-1$ steps. Moreover, the algorithms produces a $\mathbb{Q}\in  \mathcal{M}_{2}(\mu,\nu)$  since
\begin{enumerate}
\item the whole probability mass $\omega_j$ is transported from $x_j$ and the whole mass $\vartheta_i$ is received by $y_i$.
\item the martingale property is satisfied. This is obvious in Case I and guaranteed by the second linear equations  in Case II.
\end{enumerate}
What is left to show is that the algorithm continues after each step, i.e. that the remaining probability masses still satisfy $\mu\le_c \nu.$ For this property indeed it is crucial that the algorithm starts with the smallest atom $x_1$. We will show this statement for each case separately by using the characterization of Lemma \ref{Prop Char Conv Ord MultiMarginal} of the convex order. Note that the condition $\frac{1}{\mu(\mathbb{R})}\int x\mu(dx)= \frac{1}{\nu(\mathbb{R})}\int x\nu(dx)$ is satisfied in each step, since the same probability masses are subtracted on both sides and the integrals remain the same due to the second equation.

\begin{description}
\item[Case I]  There exists $\ell \in \{1,\ldots,M\}$ such that $x_1=y_\ell$.
Initially we have by assumption that \begin{equation}\label{eq:cnew}
\mathbb{E}[(X-t)_+]= \sum_{j=1}^N(x_j-t)_+\omega_j \le\sum_{i=1}^M(y_i-t)_+\vartheta_i= \mathbb{E}[(Y-t)_+].\end{equation} In Case I the new probability mass is constructed by subtracting on both sides the same mass which belongs to the same value. Thus, obviously the inequality still holds true and the remaining masses still satisfy $\mu\le_c\nu.$
\item[Case II]  There exists $\ell \in \{1,\ldots,M-1\}$ such that $y_\ell<x_1<y_{\ell+1}$.
\begin{description}
 \item[Case 1] $\vartheta'_\ell \le \vartheta_\ell, \vartheta'_{\ell+1}\le \vartheta_{\ell+1}.$

 First assume that $t\le y_\ell$. In what follows we will always set $m_t := \min\{i: y_i > t\}$. By assumption we have that
 $$\mathbb{E}[(X-t)_+]= \sum_{j=1}^N(x_j-t)\omega_j \le\sum_{i=m_t}^M(y_i-t)\vartheta_i= \mathbb{E}[(Y-t)_+].$$ This implies
 \begin{eqnarray*}
 \sum_{j=2}^N(x_j-t)\omega_j  &\le & \sum_{i=m_t}^M(y_i-t)\vartheta_i-(x_1-t) \omega_1\\
 &=& \sum_{i=m_t}^M(y_i-t)\vartheta_i-(y_\ell-t)\vartheta'_\ell -(y_{\ell+1}-t)\vartheta'_{\ell+1}
 \end{eqnarray*}
  which implies that the new masses still satisfy the inequality.
  Now assume that $t\ge y_{\ell+1}$. In this case all atoms where masses are changed do not enter the expectation. Hence obviously the inequality \eqref{eq:cnew} is satisfied for the new masses. Finally assume that $t\in (y_\ell,y_{\ell+1})$ and let $n_t := \min\{j: x_j > t\}$. Here we have to show that
  $$  \sum_{j=n_t}^N(x_j-t)\omega_j  \le \sum_{i=\ell+1}^M(y_i-t)\vartheta_i-(y_{\ell+1}-t)\vartheta'_{\ell+1}$$
which is equivalent to showing that $f(t)\ge 0$ for $ t\in (y_\ell,y_{\ell+1})$ where
$$ f(t) := \sum_{i=\ell+1}^My_i\vartheta_i-y_{\ell+1}\vartheta'_{\ell+1}-
\sum_{j=n_t}^Nx_j\omega_j + t\Big( \sum_{j=n_t}^N \omega_j - \sum_{i=\ell+1}^M \vartheta_i+\vartheta'_{\ell+1}\Big).
$$
Note that
\begin{enumerate}
\item $f$ is continuous and piecewise linear (it changes its slope at $x_j$).
\item $f(y_\ell)=f(y_{\ell+1})\ge 0.$
\item the slope is decreasing in $t$, hence $f$ is concave.
\end{enumerate}
Thus we can conclude that $f(t)\ge 0$ for $ t\in (y_\ell,y_{\ell+1})$.
  \item[Case 2] $\vartheta'_\ell > \vartheta_\ell, \vartheta'_{\ell+1}\le \vartheta_{\ell+1}.$

   First assume that $t\le y_\ell$. By assumption we have that
 $$\mathbb{E}[(X-t)_+]= \sum_{j=1}^N(x_j-t)\omega_j \le\sum_{i=m_t}^M(y_i-t)\vartheta_i= \mathbb{E}[(Y-t)_+].$$ This implies
 \begin{eqnarray*}
 \sum_{j=1}^N(x_j-t)\omega_j  -(x_1-t)\omega_1 \frac{\vartheta_\ell}{\vartheta'_{\ell}}&\le & \sum_{i=m_t}^M(y_i-t)\vartheta_i-(x_1-t)\omega_1 \frac{\vartheta_\ell}{\vartheta'_{\ell}}\\ &=&\sum_{i=m_t}^M(y_i-t)\vartheta_i-\\ &&-(y_\ell-t)\vartheta_\ell-(y_{\ell+1}-t)  \frac{\vartheta_\ell}{\vartheta'_{\ell}} \vartheta'_{\ell+1}
 \end{eqnarray*}
  which implies that the new masses still satisfy the inequality.   Now assume that $t\ge y_{\ell+1}$. In this case all atoms where masses are changed do not enter the expectation. Hence obviously the inequality \eqref{eq:cnew} is satisfied for the new masses. Finally assume that $t\in (y_\ell,y_{\ell+1})$. Define
 \begin{eqnarray*}f(t) &:=& \sum_{i=\ell+1}^My_i\vartheta_i-y_{\ell+1}\frac{\vartheta_\ell}{\vartheta'_\ell}\vartheta'_{\ell+1}-
\sum_{j=n_t}^Nx_j\omega_j -1_{[n_t=1]}x_1\omega_1\frac{\vartheta_\ell}{\vartheta'_\ell}\\&&+ t\Big( \sum_{j=n_t}^N \omega_j - 1_{[n_t=1]}\omega_1\frac{\vartheta_\ell}{\vartheta'_\ell}-\sum_{i=\ell+1}^M \vartheta_i+\frac{\vartheta_\ell}{\vartheta'_\ell}\vartheta'_{\ell+1}\Big).
\end{eqnarray*}
Note that
\begin{enumerate}
\item $f$ is continuous and piecewise linear (it changes its slope at $x_j$).
\item $f(y_\ell)=f(y_{\ell+1})\ge 0.$
\item the slope is decreasing in $t$, hence $f$ is concave.
\end{enumerate}
  Thus we can conclude that $f(t)\ge 0$ for $ t\in (y_\ell,y_{\ell+1})$.
    \item[Case 3] $\vartheta'_\ell \le \vartheta_\ell, \vartheta'_{\ell+1}> \vartheta_{\ell+1}.$

      First assume that $t\le y_\ell$. By assumption we have that
 $$\mathbb{E}[(X-t)_+]= \sum_{j=1}^N(x_j-t)\omega_j \le\sum_{i=m_t}^M(y_i-t)\vartheta_i= \mathbb{E}[(Y-t)_+].$$ This implies
 \begin{eqnarray*}
 \sum_{j=1}^N(x_j-t)\omega_j  -(x_1-t)\omega_1 \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}&\le & \sum_{i=m_t}^M(y_i-t)\vartheta_i-(x_1-t)\omega_1 \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}\\ &=&\sum_{i=m_t}^M(y_i-t)\vartheta_i-\\ &&-(y_\ell-t)\vartheta'_\ell \frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}     -(y_{\ell+1}-t)  \vartheta_{\ell+1}
 \end{eqnarray*}
  which implies that the new masses still satisfy the inequality.  Now assume that $t\ge y_{\ell+1}$. In this case all atoms where masses are changed do not enter the expectation. Hence obviously the inequality \eqref{eq:cnew} is satisfied for the new masses. Finally assume that $t\in (y_\ell,y_{\ell+1})$. Define
 \begin{eqnarray*}f(t) &:=& \sum_{i=\ell+2}^My_i\vartheta_i-
\sum_{j=n_t}^Nx_j\omega_j -1_{[n_t=1]}x_1\omega_1\frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}\\&&+ t\Big( \sum_{j=n_t}^N \omega_j - 1_{[n_t=1]}\omega_1\frac{\vartheta_{\ell+1}}{\vartheta'_{\ell+1}}-\sum_{i=\ell+2}^M \vartheta_i\Big).
\end{eqnarray*}
  Note that
\begin{enumerate}
\item $f$ is continuous and piecewise linear (it changes its slope at $x_j$).
\item $f(y_\ell)=f(y_{\ell+1})\ge 0.$
\item the slope is decreasing in $t$, hence $f$ is concave.
\end{enumerate}
  Thus we can conclude that $f(t)\ge 0$ for $ t\in (y_\ell,y_{\ell+1})$.

       \item[Case 4] $\vartheta'_\ell > \vartheta_\ell, \vartheta'_{\ell+1}> \vartheta_{\ell+1}.$

       This case either leads to case 2 or case 3.
       \end{description}
\end{description}
In total we have shown that the algorithm can be continued the same way in the next step.
\end{proof}

\section{An Algorithm for the Construction of the Optimal Solution of (D)}\label{sec:dual}
In order to construct an optimal solution for the dual problem it is clear due to the complementary slackness condition that whenever $q_{j,i}>0$ we must have
$$ \varphi_j + \psi_i + h_j (y_i-x_j)= c_{j,i}.$$
Thus, when performing the algorithm for the primal problem we can at the same time fix some variables for the optimal solution of the hedging problem. This algorithm is again recursive. Note that we have to replace $1$ and $\ell$ by the true indices of these atoms.

{\bf Algorithm (dual problem):}

\begin{description}
\item[Case I]  There exists $\ell \in \{1,\ldots,M\}$ such that $x_1=y_\ell$.
\begin{description}
\item[Case 1] $\omega_1\le \vartheta_\ell$.

In this case $x_1$ is crossed from the list and we fix $\varphi_1 := c_{1,\ell}-\psi_\ell$ since it will not appear in the remaining equations.
 \item[Case 2] $\omega_1> \vartheta_\ell$.
 In this case $y_\ell$ is crossed from the list and we fix $\psi_\ell := c_{1,\ell}-\varphi_1$ since it will not appear in the remaining equations.
 \end{description}

\item[Case II]  There exists $\ell \in \{1,\ldots,M-1\}$ such that $y_\ell<x_1<y_{\ell+1}$.
\begin{description}
 \item[Case 1] $\vartheta'_\ell \le \vartheta_\ell, \vartheta'_{\ell+1}\le \vartheta_{\ell+1}.$

 In this case we cross $x_1$ from the list and we fix $\varphi_1, h_1$ as the solution of
 \begin{eqnarray*}
 && \varphi_1 + \psi_\ell + h_1 (y_\ell-x_1)= c_{1,\ell}\\
 && \varphi_1 + \psi_{\ell+1} + h_1 (y_{\ell+1}-x_1)= c_{1,\ell+1}
 \end{eqnarray*}
 since this atom will not appear in the remaining equations.

 \item[Case 2] $\vartheta'_\ell > \vartheta_\ell, \vartheta'_{\ell+1}\le \vartheta_{\ell+1}.$

 In this case we cross $y_\ell$ from the list and we fix $\psi_\ell$ as
 \begin{eqnarray*}
 &&  \psi_\ell = c_{1,\ell} - \varphi_1 - h_1 (y_\ell-x_1)
 \end{eqnarray*}
 since this atom will not appear in the remaining equations.

    \item[Case 3] $\vartheta'_\ell \le \vartheta_\ell, \vartheta'_{\ell+1}> \vartheta_{\ell+1}.$

  In this case we cross $y_{\ell+1}$ from the list and we fix $\psi_{\ell+1}$ as
 \begin{eqnarray*}
 &&  \psi_{\ell+1} = c_{1,\ell+1} - \varphi_1 - h_1 (y_{\ell+1}-x_1)
 \end{eqnarray*}
 since this atom will not appear in the remaining equations.

       \item[Case 4] $\vartheta'_\ell > \vartheta_\ell, \vartheta'_{\ell+1}> \vartheta_{\ell+1}.$

       This case either leads to case 2 or case 3.
       \end{description}
\end{description}

Note that in each step we can choose the free variable since we only fix variables when the corresponding atom is crossed from the list. Thus, during the algorithm all $\varphi_j$ are determined and some of the $h_j$ and $\psi_i$. The remaining $h_j$ and $\psi_i$ have to be chosen such that the inequalities
 $$ \varphi_j + \psi_i + h_j (y_i-x_j)\geq c_{j,i}$$
are satisfied and the objective function is minimized. This is possible due to complementary slackness.

\section{An Example}\label{sec:example}

In this section we illustrate our algorithms to solve the problems \ref{Prob Upper Price Bound Discrete} and \eqref{Prob Super Hedging Discrete} for a payoff function $c:\mathbb{R}^2 \to \mathbb{R}$ such that \eqref{eq:optcondition} is satisfied. We consider the margins
\[
  \mu=\frac{1}{2}(\delta_{1}+\delta_3) \leq_c \nu= \frac{1}{2}\delta_0 + \frac{1}{6}\delta_2 + \frac{1}{3}\delta_5
\]
and begin illustrating the algorithm that determines $\mathbb{Q}_l(\mu,\nu)$ and thus solves \eqref{Prob Upper Price Bound Discrete}. First, we couple the mass of $\delta_1$. Therefore we solve
\begin{eqnarray}
\nonumber \vartheta'_1 + \vartheta'_{2} &=& \omega_1=\frac12\\ \nonumber
 \vartheta'_1 \cdot 0 +\vartheta'_{2} \cdot 2 &=& \omega_1 \cdot 1 = \frac12,
\end{eqnarray}
which yields $\vartheta'_1=\vartheta'_2=\frac{1}{4}.$ As $\vartheta'_1 \le \vartheta_1, \vartheta'_2> \vartheta_2$ we are in \textbf{case II, 3}. Thus we define $q_{1,1}=q_{1,2}=\frac{1}{6}$, $\omega_1=\frac{1}{6}$, $\vartheta_1=\frac{1}{3}$ and $\vartheta_2=0$.

Then we proceed to couple the remaining mass of $\delta_1$ and thus solve
\begin{eqnarray}
\nonumber \vartheta'_1 + \vartheta'_{3} &=& \omega_1=\frac16\\ \nonumber
 \vartheta'_1 \cdot 0 +\vartheta'_{3} \cdot 5 &=& \omega_1 \cdot 1=\frac16,
\end{eqnarray}
which yields $\vartheta'_1= \frac{2}{15},\vartheta'_3=\frac{1}{30}.$ As $\vartheta'_1 \le \vartheta_1, \vartheta'_3\le \vartheta_3$ we are in \textbf{case  II, 1}. Thus we define $q_{1,1}=\frac{2}{15}, q_{1,3}=\frac{1}{30}$, $\omega_1=0$, $\vartheta_1=\frac{1}{5}$ and $\vartheta_3=\frac{3}{10}$.

Finally, we couple the mass of $\delta_3$. Therefore we solve
\begin{eqnarray}
\nonumber \vartheta'_1 + \vartheta'_{3} &=& \omega_2=\frac12\\ \nonumber
 \vartheta'_1 \cdot 0 +\vartheta'_{3} \cdot 5 &=& \omega_2 \cdot 3=\frac32,
\end{eqnarray}
which clearly yields $\vartheta'_1= \frac{1}{5},\vartheta'_3=\frac{3}{10}$ and as $\vartheta'_1 \le \vartheta_1, \vartheta'_3\le \vartheta_3$ we are in \textbf{case II, 1}. Thus we define $q_{2,1}=\frac{1}{5}, q_{2,3}=\frac{3}{10}$, $\omega_2=0$, $\vartheta_1=0$ and $\vartheta_3=0$. Note that in this example mass between the same atoms is transported in two different steps. These masses have to be added in the end.

In total this yields the following unique left-monotone transport plan \[
  \mathbb{Q}_l(\mu,\nu)=\frac{3}{10}\delta_{(1,0)} + \frac{1}{6}\delta_{(1,2)} + \frac{1}{30}\delta_{(1,5)} + \frac{1}{5}\delta_{(3,0)}+ \frac{3}{10}\delta_{(3,5)}.
\]

    \begin{figure}[!htbp]\begin{center}

       \begin{tikzpicture}[scale=1]

          \node  at (0,0.4) {$1$};
          \node at (2,0.4) {$3$};
         \fill (0,0) circle (0.05);
         \fill (2,0) circle (0.05);

          \node at (-2,-3.4) {$0$};
      \node at (1,-3.4) {$2$};
          \node at (3,-3.4) {$5$};
           \fill (-2,-3) circle (0.05);
                  \fill (1,-3) circle (0.05);
                  \fill (3,-3) circle (0.05);

                  \draw(0,0)--(-2,-3);
                     \draw(0,0)--(1,-3);
                      \draw(0,0)--(3,-3);
                     \draw(2,0)--(-2,-3);
                      \draw(2,0)--(3,-3);

                        \node  at (-1.4,-1.5) {$\frac{3}{10}$};
                         \node  at (2.8,-1.5) {$\frac{3}{10}$};
                          \node  at (1.2,-2.5) {$\frac{1}{6}$};
                            \node  at (-0.9,-2.5) {$\frac{1}{5}$};
                              \node  at (1.2,-1.5) {$\frac{1}{30}$};

        \end{tikzpicture}
        \caption{Left-monotone transport plan in the example.}
  \end{center}
    \end{figure}

Parallel to determining $\mathbb{Q}_{l}(\mu,\nu)$ we may determine the hedging functions $\varphi_1,\varphi_2$, $\psi_2$, and $h_1,h_2$ as follows. In the first step, where we are in \textbf{case II, 3}, we define \[
  \psi_2=c_{1,2}-\varphi_1-h_1(y_2-x_1)=c_{1,2}-\varphi_1-h_1.
\]
In the second step, where we are in \textbf{case II,1}, we solve the given linear equation system which yields\[
  h_1=\frac{c_{1,3}-c_{1,1}-(\psi_3-\psi_1)}{y_3-y_1}=\frac{c_{1,3}-c_{1,1}-(\psi_3-\psi_1)}{5}
\]
and
\[
  \varphi_1 = \frac{y_3-x_1}{y_3-y_1}(c_{1,1}-\psi_1)+\frac{x_1-y_1}{y_3-y_1}(c_{1,3}-\psi_3) = \frac{4}{5}(c_{1,1}-\psi_1)+\frac{1}{5}(c_{1,3}-\psi_3).
\]
Finally, in the third step we analogously determine
\[
  h_2=\frac{c_{2,3}-c_{2,1}-(\psi_3-\psi_1)}{y_3-y_1}=\frac{c_{2,3}-c_{2,1}-(\psi_3-\psi_1)}{5}
\]
and
\[
  \varphi_2 = \frac{y_3-x_2}{y_3-y_1}(c_{2,1}-\psi_1)+\frac{x_2-y_1}{y_3-y_1}(c_{2,3}-\psi_3) = \frac{2}{5}(c_{2,1}-\psi_1)+\frac{3}{5}(c_{2,3}-\psi_3).
\]
Then by construction we have \[
  \varphi_j+\psi_i+h_j(y_i-x_j)=c_{j,i}
\]
for all $j,i$ such that $(x_j,y_i) \in \mathrm{supp}(\mathbb{Q}_{l}(\mu,\nu)).$ The only remaining inequality to be satisfied is
\[
  \varphi_2+\psi_2+h_2(y_2-x_2)\geq c_{2,2}.
\]
Plugging in all defined notions, after a short calculation we obtain that the above is equivalent to\[
  \frac{3}{5}(c_{2,1}-c_{1,1})+\frac{2}{5}(c_{2,3}-c_{1,3})-(c_{2,2}-c_{1,2}) \geq 0,
\] which is satisfied by assumption. That is, we may define $\psi_1=\psi_3=0$ and thus obtain
$
  \varphi_1=\frac{4}{5}c_{1,1}+\frac{1}{5}c_{1,3},\ h_1=\frac{c_{1,3}-c_{1,1}}{5},\ \varphi_2=\frac{2}{5}c_{2,1}+\frac{3}{5}c_{2,3},\ h_2=\frac{c_{2,3}-c_{2,1}}{5},$ and $\psi_{2}=c_{1,2}-\frac{3}{5}c_{1,1}-\frac{2}{5}c_{1,3}.
$
Note that the solution of the dual problem is not unique.
For further examples we refer the reader to \cite{schm18}.


\section{Introduction and formalism}

Transverse momentum dependent parton distributions (TMDs) are nowadays believed to represent fundamental pieces in our
comprehension of the internal structure of the nucleons.
They are also crucial to explain very interesting spin-polarization phenomena like
single spin asymmetries (SSAs). The Sivers function~\cite{Sivers:1989cc},
related to the asymmetry in the azimuthal distribution of unpolarized partons
inside a high-energy transversely polarized proton, is of great interest.
Its existence is well established and many phenomenological extractions of its quark component are available.
In contrast, the gluon Sivers function (GSF) is much more difficult to access and still poorly known.

In Ref.~\cite{DAlesio:2015fwo} a first extraction of the GSF from mid-rapidity data on
SSAs in $p^\uparrow p \to \pi^0\,X $~\cite{Adare:2013ekj}
has been obtained within a TMD approach, the so-called generalized parton model (GPM).
In this scheme, applied to single-scale processes, TMDs are taken to be universal.
Although not formally proven, it has been shown to be phenomenologically very successful
(see Refs.~\cite{DAlesio:2007bjf, Aschenauer:2015ndk}).

A QCD extension of this model, the Color Gauge Invariant GPM (CGI-GPM),
was developed in Refs.~\cite{Gamberg:2010tj,DAlesio:2011kkm},
focusing on the process dependence of quark Sivers function.
Here, the effects of initial- (ISI) and final- (FSI) state interactions
are taken into account, within a one-gluon exchange approximation.
In particular, the CGI-GPM can reproduce the opposite relative sign of the quark Sivers
functions in SIDIS and in the Drell-Yan processes~\cite{Collins:2002kn,Brodsky:2002cx}.

In Ref.~\cite{DAlesio:2017rzj} the CGI-GPM has been extended, for the first time, to the
gluon Sivers effect in SSAs for inclusive $J/\psi$ and $D$-meson production
in $pp$ collisions, while inclusive pion and photon production
has been considered in Ref.~\cite{DAlesio:2018rnv}.
In this scheme, two different classes of modified partonic cross sections enter,
each convoluted with a different GSF. These universal and independent distributions
are named the $f$-type and $d$-type gluon Sivers functions.

We present here a phenomenological analysis, giving only few basic expressions
and referring to Refs.~\cite{DAlesio:2017rzj,DAlesio:2018rnv} for all details.
The SSA for the inclusive process $p^\uparrow p \to h \, X$, defined as
\begin{equation}
A_N \equiv \frac{\d\sigma^\uparrow-\d\sigma^\downarrow}{\d\sigma^\uparrow+\d\sigma^\downarrow}
\equiv\, \frac{ \d\Delta\sigma}{ 2 \d\sigma}\,,
\end{equation}
within a TMD scheme and in suitable kinematical regions, is mainly sensitive to the Sivers function:
\begin{equation}
\Delta \hat f_{a/p^\uparrow}\,(x_a, \bm{k}_{\perp a})
\label{defsiv}
= \Delta^N f_{a/p^\uparrow}\,(x_a, k_{\perp a}) \> \cos\phi_a
=  -2 \, \frac{k_{\perp a}}{M_p} \, f_{1T}^{\perp a} (x_a, k_{\perp a}) \> \cos\phi_a \, .
\end{equation}
In particular, for $D$-meson, $J/\psi$ or mid-rapidity $\pi\;(\gamma)$ production
in $pp$ collisions only the gluon Sivers effect can be sizeable,
being all these processes dominated by gluon initiated sub-processes.

Within the CGI-GPM the $f$- and $d$-type contributions to $\d\Delta\sigma$, \emph{to be added up}, read
\begin{align}
  \left ( -\frac{k_{\perp g}}{M_p} \right ) f^{\perp g\, (f,d)}_{1 T}(x_g, k_{\perp g})\cos\phi_g \otimes
\> f_{b/p}(x_b, \bm{k}_{\perp b}) \otimes H^{{\rm Inc}\, (f,d)}_{gb \to cd}
\otimes \>  D_{h/c}(z, \bm{k}_{\perp h}) \>,
\label{sivgen}
\end{align}
where $H^{\rm Inc}$ include ISIs and FSIs.
For $\gamma$ production $D_{h/c}$ is replaced by a Dirac delta-function on $z$ and $\bm{k}_{\perp h}$,
and for $J/\psi$ production by a proper non perturbative factor.
In the GPM approach only one GSF appears and $H^{\rm inc}$ is replaced
by the ordinary unpolarized partonic contribution $H^U$.

The first $k_\perp$-moment of the Sivers function is also of relevance:
\begin{equation}
\Delta^N \! f_{g/p^\uparrow}^{(1)}(x) =
\int \d^2 \bm{k}_\perp \frac{k_\perp}{4 M_p} \Delta^N \! f_{g/p^\uparrow}(x,k_\perp) \equiv - f_{1T}^{\perp (1) g}(x) \, .
\label{siversm1}
\end{equation}
For the gluon Sivers function we adopt a factorized Gaussian-like form:
\begin{equation}
\Delta^N\! f_{g/p^\uparrow}(x,k_\perp) =
\left (-2\frac{k_\perp}{M_p}  \right )f_{1T}^{\perp\,g} (x,k_\perp)  = 2 \, {\cal N}_g(x)\,f_{g/p}(x)\,
h(k_\perp) \,\frac{e^{-k_\perp^2/\langle k_\perp^2 \rangle_g}}
{\pi \langle k_\perp^2 \rangle_g}\,,
\label{eq:siv-par-1}
\end{equation}
where $f_{g/p}(x)$ is the unpolarized collinear gluon distribution and
\begin{equation}
{\cal N}_g(x) = N_g x^{\alpha}(1-x)^{\beta}\,
\frac{(\alpha+\beta)^{(\alpha+\beta)}}
{\alpha^{\alpha}\beta^{\beta}}\;\;\;\;\;\; h(k_\perp) = \sqrt{2e}\,\frac{k_\perp}{M'}\,e^{-k_\perp^2/M'^2}\,,
\label{eq:nq-coll}
\end{equation}
with $|N_g|\leq 1$. We also define the parameter $ \rho = M'^2/(\langle k_\perp^2 \rangle_g +M'^2)$, with $0 < \rho < 1$.
Notice that the collinear component of all TMDs defined above evolve through a DGLAP evolution.

In the sequel, by a combined analysis of SSAs for mid-rapidity pion and $D$-meson production, we will show how
one can constrain the two GSFs in the CGI-GPM.
Predictions for $p^\uparrow p\to J/\psi\,X$
and $p^\uparrow p \to \gamma\, X$ are also given, as well as the corresponding results within the GPM.

\section{Results and comparison with data}

In the GPM approach, via
the analysis of the SSA in $pp\to\pi^0 X$ at mid-rapidity~\cite{DAlesio:2015fwo}, one can directly extract the GSF
and then give predictions to be compared against other sets of data.
The parameters obtained, adopting for $\langle k_\perp^2\rangle_g = 1$ GeV$^2$, are~\cite{DAlesio:2018rnv}:
\begin{equation}
\label{eq:par_gsf_GPM}
N_g=0.25\,, \hspace*{1cm} \alpha = 0.6\,, \hspace*{1cm}\beta = 0.6\,, \hspace*{1cm}\rho = 0.1  \,.
\end{equation}
The analysis within the CGI-GPM, for which we have two independent GSFs, appears more difficult.
The first issue we address is the effective role of the $f$- and $d$-type contributions.
To do this, we maximize the effects of the two functions, by saturating the positivity bound
for their $x$-dependent parts (i.e.~${\cal N}_g(x)=\pm 1$) and adopting the value $\rho = 2/3$.

In Fig.~\ref{fig:AN-pi} (left panel) we present the maximized (${\cal N}_g(x) = +1$) gluon Sivers
contributions to $A_N$ for the process $p^\uparrow p\to \pi^0\, X$ at $\sqrt s=200$ GeV
and mid-rapidity as a function of $p_T$, together with the PHENIX data~\cite{Adare:2013ekj}:
$f$-type (red solid line), $d$-type (blue dot-dashed line) and GPM (green dashed line).
The quark Sivers contribution (red dotted line) is totally negligible.
More importantly, the $d$-type term appears strongly suppressed.
This is due to the cancellation between the $gq$ and $g\bar q$ channels, entering with a relative opposite sign,
and the absence of the $gg\to gg$ contribution,
relevant at moderate values of $p_T$. In contrast, the GPM approach gives the largest maximized effect,
since its partonic contributions are all positive and unsuppressed.

\begin{figure}[t]
\begin{center}
\includegraphics[trim = 1.cm 0cm 1cm 0cm, width=7cm]{AN_pi0_PHENIX_xF0_GPM_CGI.pdf}
\includegraphics[trim = 1.cm 0cm 1cm 0cm, width=7cm]{AN_pi0_PHENIX_xF0_CGI_full.pdf}
\caption{Left panel: maximized gluon Sivers contributions (${\cal N}_g(x) =+1$) to $A_N$ for
$p^\uparrow p\to \pi^0\, X$ as a function of $p_T$. The quark Sivers effect, as extracted from SIDIS fits,
is also shown (red dotted line). Right panel: $A_N$ estimates obtained adopting
a reduced $f$-type GSF ((${\cal N}_g^{(f)}(x) =0.1$)) and a negative saturated $d$-type GSF (${\cal N}_g^{(d)}(x) =-1$).
Shaded area represents a $\pm 20$\% uncertainty on ${\cal N}_g^{(f)}$. Data are from Ref.~\cite{Adare:2013ekj}.}
\label{fig:AN-pi}
\end{center}
\end{figure}

To obtain SSAs comparable with the very tiny and \emph{positive} data values,
one would conservatively expect a relative cancellation between the $f$- and $d$-type pieces,
with a strongly reduced and \emph{positive} $f$-type GSF.
The results, for ${\cal N}_g^{(f)}(x)=+0.1$ and ${\cal N}_g^{(d)}(x)=-1$, are
shown in the right panel of Fig.~\ref{fig:AN-pi},
together with an overall uncertainty band of about $\pm$20\% on ${\cal N}_g^{(f)}$.

Let us now consider $A_N$ for $D^0$ production at $\sqrt s= 200$ GeV
and the corresponding analysis for its muon decays,
for which data are available~\cite{Aidala:2017pum}.
In Fig.~\ref{fig:AN-D0-sat} we show the results for $A_N$ as a function of $x_F$ for different $p_T$ values,
obtained by separately maximizing the $d$- (left panel) and $f$-type (right panel) contributions.
In this case, in the forward region the $d$-type term is sizeable, while the $f$-type is suppressed.
The reason is that, for $D^0$ production, at leading order only the charm-quark
fragmentation function enters and no cancellation, as in the previous case, occurs.
Moreover, the hard partonic parts favor the $d$-type w.r.t.~the $f$-type term at $x_F>0$, while
for $x_F<0$ both effects are almost washed out by the integration over the Sivers azimuthal phase.

Adopting now the results of the analysis of $\pi^0$ SSA data, the $f$-type piece should be accordingly reduced by
0.1, thus becoming practically negligible.
This means that for $D^0$ production we can safely consider
only the $d$-type GSF. Similar considerations apply also to $\bar D^0$ production.

\begin{figure}[t]
\begin{center}
\includegraphics[trim = 1.cm 0cm 1cm 0cm, width=7cm]{AN_D0_PHENIX_d-type-sat.pdf}
\includegraphics[trim = 1.cm 0cm 1cm 0cm,width=7cm]{AN_D0_PHENIX_f-type-sat.pdf}
\caption{Maximized (${\cal N}_g(x) =+1$) $A_N$ for $p^\uparrow p\to D^0\, X$ at different $p_T$
values (1 -- 6 GeV) as a function of $x_F$, within the CGI-GPM approach:
$d$-type (left panel) and $f$-type (right panel) contributions.}
\label{fig:AN-D0-sat}
\end{center}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[trim = 1.cm 0cm 1cm 0cm, width=7cm]{AN_muplus_d-type_vs_xF-v2.pdf}
\includegraphics[trim = 1.cm 0cm 1cm 0cm, width=7cm]{AN_muminus_d-type_vs_xF-v2.pdf}
\caption{$d$-type gluon contributions to $A_N$ for $p^\uparrow p\to \mu^+\, X$ (left panel) and
$p^\uparrow p\to \mu^-\, X$ (right panel) from $D$-meson decay as a function of $x_F$. Maximized effect:
${\cal N}_g^{(d)}(x) =+1$ (thin red solid lines), ${\cal N}_g^{(d)}(x) =-1$ (thin blue dot-dashed lines).
GPM predictions (green dashed lines) are also shown. Data are from Ref.~\cite{Aidala:2017pum}.}
\label{fig:AN-muon}
\end{center}
\end{figure}

For a direct comparison with data, the estimates for $D$-meson production have
to be converted to the SSAs for their muon decay products~\cite{Aidala:2017pum}.
Although the amount and precision of muon data do not allow for a true fit,
important information can be gathered even adopting a very simple scenario with ${\cal N}_g^{(f,d)}(x)\equiv N_g^{(f,d)}$.

As shown in Fig.~\ref{fig:AN-muon} the maximized $d$-type contributions
(thin red solid line: ${\cal N}_g^{(d)} =+1$, thin blue dot-dashed line: ${\cal N}_g^{(d)} = -1$)
largely overestimate the positive $x_F$ data in size. Notice that ${\cal N}_g^{(d)} = -1$,
together with ${\cal N}_g^{(f)} = + 0.1$, was adopted to reproduce the $\pi^0$ SSA data (see Fig.~\ref{fig:AN-pi}, right panel).
To get a fair account of the muon SSA data, one can take  $|{\cal N}_g^{(d)}| \leq 0.15$,
with a mild preference for positive values, because of the positive $\mu^+$ data point.
In Fig.~\ref{fig:AN-muon}, the results with ${\cal N}_g^{(d)} = +0.15(-0.15)$ are shown
as thick red solid lines (thick blue dot-dashed lines)
both for $\mu^+$ (left panel) and $\mu^-$ (right panel) production.

Going back to the pion SSA, we find that by varying ${\cal N}_g^{(d)}$
in the range $-0.15\,\div\,+0.15$, while keeping $\rho=2/3$, a very good description of both
the $\mu^\pm$ and $\pi^0$ data can be obtained if ${\cal N}_g^{(f)}$
varies in the \emph{corresponding} range $+0.05\,\div\,-0.01$:
\begin{equation}
\label{eq:par_gsf_CGI}
{\cal N}_g^{(d)} = - 0.15\,  \rightarrow  {\cal N}_g^{(f)} = +0.05 \hspace*{1cm}
{\cal N}_g^{(d)} = + 0.15\,  \rightarrow  {\cal N}_g^{(f)} = -0.01\,.
\end{equation}
A stronger reduction of the $f$-type GSF is then
required by the combined analysis of $\pi$ and $\mu$ data.
In contrast, in the GPM approach the parametrization of the GSF
extracted from the $\pi^0$ SSA data, see Eq.~(\ref{eq:par_gsf_GPM}),
allows us to describe the SSA for $\mu^\pm$ fairly well (Fig.~\ref{fig:AN-muon}, green dashed lines).

\noindent
Note that a comparable good agreement with data is reached within the twist-three approach~\cite{Koike:2011mb}.

\begin{figure}[b]
\begin{center}
\includegraphics[width=8cm]{gsf-1stmom_CGI-GPM.pdf}
\caption{Upper values for the first $k_\perp$-moments of the gluon Sivers functions at $Q^2 = 2$ GeV$^2$.
}
\label{fig:1stmom}
\end{center}
\end{figure}

The absolute value of the first $k_\perp$-moment of
the GSFs is shown in Fig.~\ref{fig:1stmom} for the GPM (green dashed line)
and the CGI-GPM approaches, $d$-type (blue dot-dashed line) and $f$-type ($N_g^{(f)}=0.05$, red solid line),
together with the positivity bound (black dotted line).  Notice that
the GSF in the GPM approach, for which the hard partonic parts are all positive, is the smallest.

\section{Predictions}

$A_N$ for $J/\psi$ production is directly sensitive to the gluon Sivers function,
and, within the CGI-GPM approach and the Color Singlet model,
only to the $f$-type distribution~\cite{DAlesio:2017rzj}.

In Fig.~\ref{fig:Jpsi} (left panel) we show our estimates against PHENIX data~\cite{Aidala:2018gmp}
for $A_N$ in $p^\uparrow p\to J/\psi\, X$:
GPM (green dashed line), CGI-GPM (red band: $-0.01\le {\cal N}_g^{(f)} \le 0.05$).
With the exception of one point at $x_F<0$,
data are compatible with zero and well described by our estimates.

The corresponding analysis of $A_N$ at LHC in the fixed polarized target
mode ($\sqrt s=115$ GeV) could be also extremely useful.
Here one could probe even larger $x$ in the polarized proton,
accessing the gluon TMDs in a complementary region.
Notice that in such a configuration the backward rapidity region refers to the forward region
for the polarized proton target.
Predictions within the GPM (thick green dashed lines) and the CGI-GPM (red bands) approaches,
together with the corresponding upper/lower positivity bounds (thin lines)
are shown in the right panel of Fig.~\ref{fig:Jpsi}.

SSAs for direct photon production is another ideal tool to access the GSF.
In Fig.~\ref{fig:AN-gamma} we show our predictions obtained adopting Eqs.~(\ref{eq:par_gsf_GPM}), (\ref{eq:par_gsf_CGI}).
Despite of their tiny values, a measure of $A_N$ for direct photon production
would be extremely important to test the consistency of this picture.

\begin{figure}[t]
\begin{center}
\includegraphics[trim = 1.cm 0cm 1cm 0cm,width=7cm]{AN_Jpsi_vs_xF_pt165-CGI-GPM-v2.pdf}
\includegraphics[trim = 1.cm 0cm 1cm 0cm,width=7cm]{AN_Jpsi_vs_pT-y2-CGI-GPM-LHC-v2.pdf}
\caption{Left panel: Predictions of $A_N$ for $p^\uparrow p\to J/\psi\, X$
at $\sqrt s=200$ GeV as a function of $x_F$ (RHIC set-up). Data are from Ref.~\cite{Aidala:2018gmp}.
Right panel: $A_N$ for $pp^\uparrow \to J/\psi\, X$ at $\sqrt{s} = 115$ GeV and $y=-2$
as a function of $p_T$ (LHC fixed target set-up): predictions (thick lines) and maximized contributions (thin lines).
}
\label{fig:Jpsi}
\end{center}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[trim = 1.cm 0cm 1cm 0cm,width=7cm]{AN_gamma_PHENIX_xF0-v2.pdf}
\includegraphics[trim = 1.cm 0cm 1cm 0cm,width=7cm]{AN_gamma_PHENIX_xF-01-v2.pdf}
\caption{Predictions of $A_N$ for $p^\uparrow p\to \gamma\, X$
as a function of $p_T$ (see Eqs.~(\ref{eq:par_gsf_GPM}), (\ref{eq:par_gsf_CGI})).
}
\label{fig:AN-gamma}
\end{center}
\end{figure}

Summarizing, the present analysis can be considered the first detailed attempt
towards a quantitative extraction of the process dependent GSFs in the CGI-GPM scheme.
Even if it is not yet possible to clearly discriminate between the GPM and the CGI-GPM approaches, the results are encouraging.
More data, and more precise, will certainly help in shedding light on this important TMD,
as well as on the validity of the whole approach.

\providecommand{\href}[2]{#2}\begingroup\raggedright