# Convex envelopes for ray-concave functions

Javiera Barrera  Faculty of Engineering and Sciences, Universidad Adolfo
Ibáñez, Santiago, Chile<EMAIL_ADDRESS>,  Eduardo Moreno  Faculty of
Engineering and Sciences, Universidad Adolfo Ibáñez, Santiago, Chile
<EMAIL_ADDRESS>and  Gonzalo Muñoz  Institute of Engineering Sciences,
Universidad de O’Higgins, Rancagua, Chile<EMAIL_ADDRESS>
###### Abstract.

Convexification based on convex envelopes is ubiquitous in the non-linear
optimization literature. Thanks to considerable efforts of the optimization
community for decades, we are able to compute the convex envelopes of a
considerable number of functions that appear in practice, and thus obtain
tight and tractable approximations to challenging problems. We contribute to
this line of work by considering a family of functions that, to the best of
our knowledge, has not been considered before in the literature. We call this
family _ray-concave_ functions. We show sufficient conditions that allow us to
easily compute closed-form expressions for the convex envelope of ray-concave
functions over arbitrary polytopes. With these tools, we are able to provide
new perspectives to previously known convex envelopes and derive a previously
unknown convex envelope for a function that arises in probability contexts.

###### Key words and phrases:

Convex envelopes, Nonlinear programming, Convex optimization

###### 2020 Mathematics Subject Classification:

Primary: 90C26, 90C25

The research leading to these results received funding from grants
ANID/CONICYT-Fondecyt Regular 1200809 (J.B., E.M.) and ANID/CONICYT-Fondecyt
Iniciación 11190515 (G.M.).

##  1\. Introduction

Strong convex relaxations of complex optimization problems is a key component
in the development of tractable computational techniques in the field. In this
regard, a popular approach has been the study of _convex underestimators_ of
functions, that is, given an arbitrary function $f$, find a convex function
$f^{\prime}$ such that $f^{\prime}(x)\leq f(x)$ $\forall x\in P$, where $P$ is
a given convex set. Such function can be used to _relax_ a sub-level set
$\\{x\in P:\,f(x)\leq 0\\}$ with the convex set $\\{x\in
P\,:\,f^{\prime}(x)\leq 0\\}$, and thus obtain a computationally tractable
approximation. The pointwise largest convex underestimator is known as the
_convex envelope of $f$ over $P$_, and the optimization community has
allocated considerable efforts on finding such envelopes for various classes
of functions $f$ and sets $P$.

######  Definition 1.

The convex envelope of a function $f$ on a subset $P$ is given by

$\displaystyle f_{\text{vex}}(x)$ $\displaystyle=\sup\\{g(x)\,:\,g\text{ is
convex and }g(x)\leq f(x)\ \forall x\in P\\}$

In this work, we consider $P$ a polytope and study the convex envelope of a
family of functions that are required to be convex of the facets of $P$, and
what we term as _ray-concavity_.

######  Definition 2.

A function $f:P\to\mathbb{R}$ is _ray-concave_ over $P$ if, for every $x\in
P$, the function $f$ restricted to $\\{\alpha x\,:\,\alpha\geq 0\\}\cap P$ is
concave.

We present sufficient conditions for deriving simple closed-form formulas of
the convex envelopes of ray-concave functions over arbitrary polytopes in any
dimension.

Our result is closely related to known results for general functions over
polytopes. To the best of our knowledge, the vast majority of the work
producing closed-form formulas of convex envelopes in arbitrary dimension
either require a rectangular domain, or require $f$ to be _edge-concave_ , in
which case the convex envelope is polyhedral111A function is _polyhedral_ if
its epigraph is a polyhedron.. With our result, through the concept of ray-
concavity, we are able to explicitly construct convex envelopes which are not
necessarily polyhedral, in any dimension, for a new family of functions that
has not been explicitly exploited before in the literature.

Our result yields a previously unknown convex envelope of a function that
appears in probability contexts.

######  Example 1.

The function $f(x,y)=-\frac{x\cdot y}{x+y-x\cdot y}$ is ray-concave over any
box $[0,u_{x}]\times[0,u_{y}]$ with $u_{x},u_{y}\leq 1$.

This function is one of the the main motivations behind this work.
Additionally, many functions for which their convex envelope formulas are
known exhibit ray-concavity (e.g., $f(x_{1},x_{2})=-x_{1}x_{2}$ or
$f(x_{1},x_{2})=x_{1}/x_{2}$ for $x_{1},x_{2}>0$), and our result provide a
new perspective on these expressions and alternative derivations.

##  2\. Literature review

The literature of convex envelopes is vast. Probably the most well-known and
used convex envelope is that of the bilinear function
$f(x_{1},x_{2})=x_{1}x_{2}$ over a rectangular region, for which its convex
(and concave) envelope is obtained through the McCormick envelopes [18, 1].

To the best of our knowledge, the first method capable of constructing the
convex envelope for a family of functions (as opposed to a particular
function) is provided in [30]. Based on disjunctive programming, they show a
general expression of the convex envelope for functions that are concave on
one variable, convex on the rest, and defined over a rectangular region. Later
on, in [5] the authors show how to compute the evaluation of $f_{\text{vex}}$
when $f$ is an $(n-1)$-convex function (i.e., $f$ is convex whenever one
variable is fixed to any value) over a rectangular domain. The function
evaluation requires the resolution of a convex optimization problem. In [7,
8], the authors formulate the convex envelope of a lower semi-continuous
function over a compact set as a convex optimization problem. They use this to
compute, explicitly, the convex envelope for various functions that are the
product of a convex function and a component-wise concave function, over a
box. We remark that in all the aforementioned cases, the convex envelopes may
be non-polyhedral, and that explicit calculations consider hyper-rectangular
domains.

Considerable efforts have been put in the case of polyhedral convex envelopes.
In [27, 28], it is shown that _edge-concavity_ of a function $f$ (i.e.,
concavity over all edge directions of a polytope $P$) implies that the convex
envelope of $f$ over $P$ is polyhedral. The construction of these convex
envelopes is studied in [20]. In [22], necessary and sufficient conditions for
the convex envelope to be polyhedral are also provided, and they are used to
obtain the convex envelope of a multilinear function over the unit box (see
also [25, 23]). In [19], the authors provide explicit expressions for the
facets of the convex envelope of a trilinear monomial over a box. In [3], the
authors design a cutting plane approach to generate, on-the-fly, the convex
envelope of a bilinear function over a box. The strength of the convex
underestimator of a bilinear function that is obtained from using a term-wise
convex envelopes is analyzed in [17].

Other known results include the convex envelopes of odd-degree monomials over
an interval [11] and the fractional function $f(x_{1},x_{2})=x_{1}/x_{2}$ over
a rectangle [31, 30]. Recently, in [15] the author computed the convex
envelope of cubic functions in two dimensions, over a rectangular region.

While a big portion of these works involve rectangular regions, there exist
important work considering sets beyond boxes in two dimensions. In [26], the
authors derive explicit formulas for the convex envelope of bilinear bivariate
functions over a class special polytopes called $D$-polytopes. The case of the
fractional function $x_{1}/x_{2}$ over a trapezoid is studied in [9]. This was
expanded in [4], where convex envelopes for bilinear and fractional bivariate
functions over quadrilaterals are constructed. The convex envelope of a
bilinear bivariate function over a triangle has been carefully studied in [26,
12, 2]. Such envelopes were tested computationally in [12] within a branching
scheme for QCQPs with positive results. In [16] it is shown how to evaluate
the convex envelope, and obtain a supporting hyperplane, for bivariate
functions over arbitrary polytopes. This approach involves solving a low-
dimensional convex problem. This procedure was refined in [14], by shifting
the calculations to the solution of a KKT system. These last techniques were
extensively tested in [21] to improve general-purpose optimization routines.
In [13], the author characterizes the convex envelope of various bivariate
functions (including the bilinear and fractional functions) over arbitrary
polytopes using a polyhedral sub-division of the polytopes. In some cases, the
convex envelope in each element of the sub-division can be given explicitly.

To the best of our knowledge, there is no construction that can provide a
closed-form formula for the convex envelope of Example 1, and almost no
construction allowing the explicit computation of non-polyhedral convex
envelopes over polytopes beyond boxes in dimension $n\geq 3$. The only
exception that we are aware of is [29]. In this work, the authors derive
explicit convex and concave envelopes of several functions on sub-sets of a
hyper-rectangle, which are obtained through polyhedral subdivisions. In this
case the authors can obtain, in closed form, the convex envelope of
disjunctive functions of the form $xf(y)$, and the concave envelope of
_concave-extendable_ supermodular functions. This may produce non-polyhedral
envelopes. We remark some similarities with their construction below, however,
it is worth noting that the results in [29] cannot directly provide a formula
for $f_{\text{vex}}$ for the function $f$ in Example 1. On one hand, such
function $f$ does not fit the disjunctive framework of [29], so we cannot
apply their convex envelope construction. On the other hand, one could
consider using their concave envelope results with $-f$, thus effectively
constructing $f_{\text{vex}}$. However, we show below that the convex envelope
of such $f$ requires a polyhedral partition that introduces new vertices in
the box, while the construction of the concave envelope for concave-extendable
functions (see [29, Corollary 2.8]) is based solely on the original vertices
of the polytope.

##  3\. Convex envelopes for ray-concave functions

Overall, we consider a polytope $P\subset\mathbb{R}^{n}$ with non-empty
interior.

######  Definition 3.

For any $v\in\mathbb{R}^{n}\setminus\\{0\\}$ such that $\exists\,\alpha\geq
0,\alpha v\in P$ (i.e., the ray defined by $v$ intersects the polytope) we
define

$\displaystyle v^{+}$ $\displaystyle=\alpha^{+}v\text{, where
}\alpha^{+}=\arg\max\\{\alpha\,:\,\alpha\geq 0,\,\alpha v\in P\\}$ and
$\displaystyle v^{-}$ $\displaystyle=\alpha^{-}v\text{, where
}\alpha^{-}=\arg\min\\{\alpha\,:\,\alpha\geq 0,\,\alpha v\in P\\}.$

In simple words, $v^{+}$ and $v^{-}$ are the intersections of the ray given by
$v$ with the boundary of $P$ (see Figure 2). Note that if $0\in P$ then
$v^{-}=0$ for all $v\in P$.

We remark that $v^{\pm}$ are continuous functions of $v$. Below, we emphasize
this functional aspect when taking derivatives.

Using this definition, a function $f:P\rightarrow\mathbb{R}$ is ray-concave
iff $f$ restricted to the segment $[v^{-},v^{+}]$ is concave for all $v$ where
$v^{\pm}$ is well defined. Our main results provides an explicit
characterization for the convex envelope of ray-concave functions that are
convex on the facets of $P$.

######  Theorem 1.

Let $f:P\to\mathbb{R}$ be a continuously differentiable and ray-concave
function over a polytope $P$, such that $f$ is convex over the facets of $P$.
Let $g:P\to\mathbb{R}$ be defined as

(1) $\displaystyle g(v)=\alpha_{v}f(v^{-})+(1-\alpha_{v})f(v^{+}),$

where $\alpha_{v}\in[0,1]$ is such that
$v=\alpha_{v}v^{-}+(1-\alpha_{v})v^{+}$. If $g$ is positively homogeneous,
then $f_{\text{vex}}=g$.

######  Remark 1.

In Section 3.3 we provide more insights on the positively homogeneous
requirement. For example, we show that whenever $0\in P$, $g$ is positively
homogeneous iff $f(0)=0$. The latter is not a restrictive requirement, as we
can compute the convex envelope of $f-f(0)$ instead.

A linear interpolation of a similar type as (1) has been considered in
multiple articles. The general result in [5], for example, shows that to
evaluate $f_{\text{vex}}$ for an edge-convex function $f$ over a box, it
suffices to consider the lines passing through $x$ where the function $f$ is
concave, similarly to our result. Each evaluation involves solving an
optimization problem (see [5, Theorem 3.1]). Another example is given by [29],
who construct envelopes explicitly using secants of a similar type. In [12],
the author also uses such lines in his construction of convex envelopes of the
bilinear function over triangles.

In our case, by considering ray-concavity, we only need to consider secants on
the rays emanating from the origin in the envelope construction.

To prove Theorem 1, we first provide three lemmas about the convexity of the
function $g$ over different regions of the domain. We divide the polytope $P$
into subregions using the rays that pass through the vertices of $P$.

######  Definition 4.

Let $\mathcal{F}$ be the set of facets of $P$. If $0\notin P$, for each pair
of facets $F_{i},F_{j}\in\mathcal{F}$ we define the region

$B_{ij}=\\{v\in P:v^{-}\in F_{i},v^{+}\in F_{j}\\}.$

We refer to the hyperplane containing the facet $F_{i}$ as the _in-hyperplane
of $B_{ij}$_, and to the hyperplane containing the facet $F_{j}$ as the _out-
hyperplane of $B_{ij}$_. Alternatively, if $0\in P$, for each facet
$F_{j}\in\mathcal{F}$ we define the region

$B_{0j}=\\{v\in P:v^{+}\in F_{j}\\}.$

In this case we only define the _out-hyperplane_ of $B_{0j}$. We denote by
$\mathcal{B}$ the set of all full-dimensional regions $B_{ij}$.

In Figure 1 we illustrate the regions we consider in $\mathcal{B}$, which
clearly form a sub-division of $P$. Note that if $B_{ij}\neq\emptyset$ then
$B_{ji}=\emptyset$.

Figure 1. Polyhedral sub-division of $P$ into regions $\mathcal{B}$ according
to intersection of rays with the boundary

Also note that every $B_{ij}\in\mathcal{B}$ is polyhedral: for example, in the
case $0\not\in P$, it is not hard to see that

(2) $B_{ij}=\mbox{cone}(F_{i})\cap\mbox{cone}(F_{j})\cap P.$

Polyhedrality follows since both $F_{i}$ and $F_{j}$ are polyhedra.

######  Remark 2.

For a given region $B_{ij}\in\mathcal{B}$ we can provide an explicit formula
for $v^{\pm}$. In fact, note that we can assume that the out-hyperplane of
$B_{ij}$ has the form ${a^{+}}^{\intercal}\vec{x}=1$. Since $v$ and $v^{+}$
lie on the same ray, we obtain $v^{+}=\frac{1}{{a^{+}}^{\intercal}v}v$ for any
$v\in B_{ij}$. Similarly, for the case $0\notin P$, we may assume that the in-
hyperplane of $B_{ij}$ has the form ${a^{-}}^{\intercal}\vec{x}=1$, and then
$v^{-}=\frac{1}{{a^{-}}^{\intercal}v}v$ for any $v\in B_{ij}$.

Moreover, since $v=\alpha_{v}v^{-}+(1-\alpha_{v})v^{+}$, this implies that

(3)
$\frac{\alpha_{v}}{{a^{-}}^{\intercal}v}+\frac{1-\alpha_{v}}{{a^{+}}^{\intercal}v}=1$

###  3.1. Convexity and differentiability over a single region

To show convexity of $g$, we first prove that under the homogeneity assumption
of Theorem 1, $g$ is convex in each region $B\in\mathcal{B}$.

######  Lemma 1.

Let $B\in\mathcal{B}$ and let $g:B\subset\mathbb{R}^{n}\rightarrow\mathbb{R}$
as defined in (1). If $g$ is positively homogeneous, then $g$ is convex in
$B$.

###### Proof.

Let $v,w\in B$ and let $z=\lambda v+(1-\lambda)w$ for $\lambda\in[0,1]$. By
convexity of the region, $z\in B$ as well. To prove the convexity of $g$ over
$B$, we show that $g(z)\leq\lambda g(v)+(1-\lambda)g(w)$.

Recall that $v^{+}$ and $w^{+}$ belong to the same facet defining $B$, and
$v^{-}$ and $w^{-}$ are either $0$ (if $0\in P$) or belong to the same facet
defining $B$ (if $0\notin P$). Hence, there exist
$\gamma,\varepsilon,\delta\in[0,1]$ such that:

$\displaystyle z$ $\displaystyle=\gamma z^{-}+(1-\gamma)z^{+}$ $\displaystyle
z^{-}$ $\displaystyle=\varepsilon v^{-}+(1-\varepsilon)w^{-}$ $\displaystyle
z^{+}$ $\displaystyle=\rho v^{+}+(1-\rho)w^{+}$

In Figure 2 we illustrate these vectors.

Figure 2. Notation for Lemma 1

Since $g(z)=\gamma f(z^{-})+(1-\gamma)f(z^{+})$ and $f$ is convex on the
facets containing $\\{v^{+},z^{+},w^{+}\\}$ and $\\{v^{-},z^{-},w^{-}\\}$ (if
$0\notin P$), we know that

(4) $\displaystyle g(z)$ $\displaystyle=\gamma f(z^{-})+(1-\gamma)f(z^{+})$
(5) $\displaystyle\leq\gamma(\varepsilon
f(v^{-})+(1-\varepsilon)f(w^{-}))+(1-\gamma)(\rho f(v^{+})+(1-\rho)f(w^{+}))$
(6) $\displaystyle=\gamma\varepsilon f(v^{-})+(1-\gamma)\rho
f(v^{+})+\gamma(1-\varepsilon)f(w^{-})+(1-\gamma)(1-\rho)f(w^{+})$

Let ${a^{+}}^{\intercal}\vec{x}=1$ be the out-hyperplane of $B$, i.e., the
hyperplane that contains $v^{+}$, $w^{+}$, and $z^{+}$. By Remark 2, we know
that

$z^{+}=\frac{1}{{a^{+}}^{\intercal}z}z=\frac{1}{{a^{+}}^{\intercal}z}\left(\lambda
v+(1-\lambda)w\right)=\underbrace{\lambda\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}}_{\rho}v^{+}+\underbrace{(1-\lambda)\frac{{a^{+}}^{\intercal}w}{{a^{+}}^{\intercal}z}}_{1-\rho}w^{+}$

where we deduce
$\rho=\lambda\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}$ because
${a^{+}}^{\intercal}z=\lambda{a^{+}}^{\intercal}v+(1-\lambda){a^{+}}^{\intercal}w$.
In a similar way, when $0\notin P$ we can apply the same for $z^{-}$ we obtain

(7)
$\varepsilon=\lambda\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}\qquad\text{and}\qquad
1-\varepsilon=(1-\lambda)\frac{{a^{-}}^{\intercal}w}{{a^{-}}^{\intercal}z}$

where ${a^{-}}^{\intercal}\vec{x}=1$ is the in-hyperplane of $B$. If $0\in P$,
then $v^{-}=w^{-}=z^{-}=0$ and thus $\varepsilon$ can take any value in
$[0,1]$. To simplify the proof, we abuse notation and consider
$\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}=\frac{{a^{-}}^{\intercal}w}{{a^{-}}^{\intercal}z}=1$
for this case, so (7) still holds.

Replacing the values of $\rho$ and $\varepsilon$ in (6), we obtain

$\displaystyle g(z)$ $\displaystyle\leq\gamma\varepsilon
f(v^{-})+(1-\gamma)\rho
f(v^{+})+\gamma(1-\varepsilon)f(w^{-})+(1-\gamma)(1-\rho)f(w^{+})$ (10)
$\displaystyle=\begin{multlined}\gamma\lambda\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}f(v^{-})+(1-\gamma)\lambda\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}f(v^{+})\\\
+\gamma(1-\lambda)\frac{{a^{-}}^{\intercal}w}{{a^{-}}^{\intercal}z}f(w^{-})+(1-\gamma)(1-\lambda)\frac{{a^{+}}^{\intercal}w}{{a^{+}}^{\intercal}z}f(w^{+})\end{multlined}\gamma\lambda\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}f(v^{-})+(1-\gamma)\lambda\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}f(v^{+})\\\
+\gamma(1-\lambda)\frac{{a^{-}}^{\intercal}w}{{a^{-}}^{\intercal}z}f(w^{-})+(1-\gamma)(1-\lambda)\frac{{a^{+}}^{\intercal}w}{{a^{+}}^{\intercal}z}f(w^{+})$
(13)
$\displaystyle=\begin{multlined}\lambda\left(\gamma\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}f(v^{-})+(1-\gamma)\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}f(v^{+})\right)\\\
+(1-\lambda)\left(\gamma\frac{{a^{-}}^{\intercal}w}{{a^{-}}^{\intercal}z}f(w^{-})+(1-\gamma)\frac{{a^{+}}^{\intercal}w}{{a^{+}}^{\intercal}z}f(w^{+})\right).\end{multlined}\lambda\left(\gamma\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}f(v^{-})+(1-\gamma)\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}f(v^{+})\right)\\\
+(1-\lambda)\left(\gamma\frac{{a^{-}}^{\intercal}w}{{a^{-}}^{\intercal}z}f(w^{-})+(1-\gamma)\frac{{a^{+}}^{\intercal}w}{{a^{+}}^{\intercal}z}f(w^{+})\right).$

What follows uses that $g$ is positively homogeneous in order to rewrite the
last inequality. To do so, note that

(14)
$\displaystyle\gamma\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}\cdot
v^{-}+(1-\gamma)\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}\cdot v^{+}$
$\displaystyle=\gamma\frac{1}{{a^{-}}^{\intercal}z}v+(1-\gamma)\frac{1}{{a^{+}}^{\intercal}z}v=v$

because $z=\gamma
z^{-}+(1-\gamma)z^{+}=\left(\gamma\frac{1}{{a^{-}}^{\intercal}z}+(1-\gamma)\frac{1}{{a^{+}}^{\intercal}z}\right)z$.
Let
$\Omega=\gamma\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}+(1-\gamma)\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}$
—this is simply the sum of the weights in the leftmost linear combination of
(14). By definition of $g$, and because we are assuming it to be positively
homogeneous, we have that

$\displaystyle g(v)$
$\displaystyle=g\left(\gamma\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}\cdot
v^{-}+(1-\gamma)\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}\cdot
v^{+}\right)$ (due to (14)) $\displaystyle=\Omega\cdot
g\left(\frac{\gamma\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}}{\Omega}\cdot
v^{-}+\frac{(1-\gamma)\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}}{\Omega}\cdot
v^{+}\right)$ (pos. homog.)
$\displaystyle=\Omega\cdot\left(\frac{\gamma\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}}{\Omega}\cdot
f(v^{-})+\frac{(1-\gamma)\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}}{\Omega}\cdot
f(v^{+})\right)$ (def. of $g$)
$\displaystyle=\gamma\frac{{a^{-}}^{\intercal}v}{{a^{-}}^{\intercal}z}\cdot
f(v^{-})+(1-\gamma)\frac{{a^{+}}^{\intercal}v}{{a^{+}}^{\intercal}z}\cdot
f(v^{+})$

A similar relation can be deduced for $w$ obtaining

$\gamma\frac{{a^{-}}^{\intercal}w}{{a^{-}}^{\intercal}z}w^{-}+(1-\gamma)\frac{{a^{+}}^{\intercal}w}{{a^{+}}^{\intercal}z}w^{+}=\left(\gamma\frac{1}{{a^{-}}^{\intercal}z}+(1-\gamma)\frac{1}{{a^{+}}^{\intercal}z}\right)w=w$

which implies

$g(w)=g\left(\gamma\frac{{a^{-}}^{\intercal}w}{{a^{-}}^{\intercal}z}w^{-}+(1-\gamma)\frac{{a^{+}}^{\intercal}w}{{a^{+}}^{\intercal}z}w^{+}\right)=\gamma\frac{{a^{-}}^{\intercal}w}{{a^{-}}^{\intercal}z}f(w^{-})+(1-\gamma)\frac{{a^{+}}^{\intercal}w}{{a^{+}}^{\intercal}z}f(w^{+}).$

Using these expressions for $g(v)$ and $g(w)$ in (13) we obtain

$g(z)=g\left(\lambda v+(1-\lambda)w\right)\leq\lambda g(v)+(1-\lambda)g(w).$

This shows that $g$ is convex in $B$. ∎

The previous lemma shows that $g$ is convex in each region $B\in\mathcal{B}$.
Before moving to convexity toward $P$, we show differentiability of $g$ in
each region and compute the corresponding gradient, which we rely on in the
next section.

######  Lemma 2.

Let $B\in\mathcal{B}$ and
$g:B\subseteq\mathbb{R}^{n}_{+}\rightarrow\mathbb{R}$ as defined in (1). Let
${a^{\pm}}^{\intercal}x=1$ be the in-hyperplane and out-hyperplane of $B$.
Then, $g$ is a differentiable function in $\mbox{int}(B)$. Moreover, the
gradient is given by

(15) $\nabla
g(v)=\frac{\alpha_{v}}{{a^{-}}^{\intercal}v}\left(\delta^{-}a^{-}+\nabla
f(v^{-})\right)+\frac{1-\alpha_{v}}{{a^{+}}^{\intercal}v}\left(\delta^{+}a^{+}+\nabla
f(v^{+})\right)$

where

$\displaystyle\delta^{-}$ $\displaystyle=\left(\nabla f(v^{*})-\nabla
f(v^{-})\right)^{\intercal}v^{-}\leq 0$ $\displaystyle\delta^{+}$
$\displaystyle=\left(\nabla f(v^{*})-\nabla
f(v^{+})\right)^{\intercal}v^{+}\geq 0$

for a vector $v^{*}$ contained on the segment $[v^{-},v^{+}]$, and
$\frac{1}{{a^{-}}^{\intercal}v}\coloneqq 0$ in the case $0\in P$.

###### Proof.

Consider $v\in\mbox{int}(B)$ arbitrary. In this proof, to aid the reader, we
emphasize that $v^{\pm}$ and $\alpha_{v}$ are functions of $v$ by referring to
them as $v^{\pm}(v)$ and $\alpha_{v}(v)$, respectively.

Since $g(v)=\alpha_{v}(v)f(v^{-}(v))+(1-\alpha_{v}(v))f(v^{+}(v))$ and $f$ is
differentiable, $g$ is also differentiable in the interior of $B$. The
gradient of $g$ is given by

$\displaystyle\nabla g(v)$ $\displaystyle=\nabla\alpha_{v}(v)\cdot
f(v^{-}(v))+\alpha_{v}(v)\nabla f(v^{-}(v))$
$\displaystyle\qquad+\nabla(1-\alpha_{v}(v))\cdot
f(v^{+}(v))+(1-\alpha_{v}(v))\nabla f(v^{+}(v))$
$\displaystyle=\nabla\alpha_{v}(v)\cdot\left(f(v^{-}(v))-f(v^{+}(v))\right)+\alpha_{v}(v)\nabla
f(v^{-}(v))$ (16) $\displaystyle\qquad+(1-\alpha_{v}(v))\nabla f(v^{+}(v))$

where

(17) $\nabla f(v^{\pm}(v))={\mathcal{D}v^{\pm}(v)}^{\intercal}\nabla
f(v)\rvert_{v=v^{\pm}}$

and $\mathcal{D}v^{\pm}(v)$ is the Jacobian matrix of $v^{\pm}(v)$. Recall
that we are assuming $v^{\pm}(v)$ intersects a facet of $P$ contained in an
hyperplane of equation ${a^{\pm}}^{\intercal}\vec{x}=1$. Hence, by Remark 2
and defining $\frac{1}{{a^{-}}^{\intercal}v}\coloneqq 0$ when $0\in P$,

(18) $v^{\pm}(v)=\frac{1}{{a^{\pm}}^{\intercal}v}v,\ \text{ so }\
\nabla\left(\frac{1}{{a^{\pm}}^{\intercal}v}\right)=\frac{-a^{\pm}}{({a^{\pm}}^{\intercal}v)^{2}}$

and

$\displaystyle\mathcal{D}v^{\pm}(v)$
$\displaystyle=v\nabla\left(\frac{1}{{a^{\pm}}^{\intercal}v}\right)^{\intercal}+\frac{1}{{a^{\pm}}^{\intercal}v}I_{n}$
$\displaystyle=\frac{-v{a^{\pm}}^{\intercal}}{({a^{\pm}}^{\intercal}v)^{2}}+\frac{1}{{a^{\pm}}^{\intercal}v}I_{n}=\frac{-1}{{a^{\pm}}^{\intercal}v}\left(v^{\pm}(v)a^{\pm\intercal}-I_{n}\right).$

Replacing in (17) we obtain

$\displaystyle\nabla f(v^{\pm}(v))$
$\displaystyle=\left(\frac{-1}{{a^{\pm}}^{\intercal}v}\left(v^{\pm}(v){a^{\pm}}^{\intercal}-I_{n}\right)\right)^{\intercal}\nabla
f(v^{\pm}(v))$ $\displaystyle=\frac{-1}{{a^{\pm}}^{\intercal}v}\left((\nabla
f(v^{\pm}(v))^{\intercal}v^{\pm}(v))a^{\pm}-\nabla f(v^{\pm}(v))\right)$

Note that ${\nabla v^{\pm}(v)}^{\intercal}v=0$ and $\nabla
f(v^{\pm}(v))^{\intercal}v=0$. This is expected because varying $v$ over its
ray does not change the position of $v^{\pm}(v)$ nor the value of
$f(v^{\pm}(v))$. On the other hand, applying the gradient to (3) we obtain

(19) $\displaystyle\nabla\alpha_{v}(v)$
$\displaystyle=\left(\frac{1}{{a^{-}}^{\intercal}v}-\frac{1}{{a^{+}}^{\intercal}v}\right)^{-1}\left(\alpha_{v}(v)\frac{a^{-}}{({a^{-}}^{\intercal}v)^{2}}+(1-\alpha_{v}(v))\frac{a^{+}}{({a^{+}}^{\intercal}v)^{2}}\right)$

Finally, the mean value theorem ensures there exist $v^{*}$ on the segment
$[v^{-}(v),v^{+}(v)]$ such that

(20) $f(v^{-}(v))-f(v^{+}(v))=\nabla
f(v^{*})^{\intercal}(v^{-}(v)-v^{+}(v))=\nabla
f(v^{*})^{\intercal}v\left(\frac{1}{{a^{-}}^{\intercal}v}-\frac{1}{{a^{+}}^{\intercal}v}\right)$

Grouping all the terms into (16), we obtain an explicit formula for the
gradient of $g$ at $v$ as

$\displaystyle\nabla g(v)$
$\displaystyle=\nabla\alpha_{v}(v)\cdot\left(f(v^{-}(v))-f(v^{+}(v))\right)+\alpha_{v}(v)\nabla
f(v^{-}(v))+(1-\alpha_{v}(v))\nabla f(v^{+}(v))$
$\displaystyle=\left(\alpha_{v}(v)\frac{a^{-}}{({a^{-}}^{\intercal}v)^{2}}+(1-\alpha_{v}(v))\frac{a^{+}}{({a^{+}}^{\intercal}v)^{2}}\right)\nabla
f(v^{*})^{\intercal}v$
$\displaystyle\qquad+\alpha_{v}(v)\left(\frac{-1}{{a^{-}}^{\intercal}v}\left(\left(\nabla
f(v^{-}(v))^{\intercal}v^{-}(v)\right)a^{-}-\nabla f(v^{-}(v))\right)\right)$
$\displaystyle\qquad+(1-\alpha_{v}(v))\left(\frac{-1}{{a^{+}}^{\intercal}v}\left(\left(\nabla
f(v^{+}(v))^{\intercal}v^{+}(v)\right)a^{+}-\nabla f(v^{+}(v))\right)\right)$
$\displaystyle=\frac{\alpha_{v}(v)}{{a^{-}}^{\intercal}v}\left(\left(\left(\nabla
f(v^{*})^{\intercal}v\right)\frac{1}{{a^{-}}^{\intercal}v}-\left(\nabla
f(v^{-}(v))^{\intercal}v^{-}(v)\right)\right)a^{-}+\nabla f(v^{-}(v))\right)$
$\displaystyle\qquad+\frac{1-\alpha_{v}(v)}{{a^{+}}^{\intercal}v}\left(\left(\left(\nabla
f(v^{*})^{\intercal}v\right)\frac{1}{{a^{+}}^{\intercal}v}-\left(\nabla
f(v^{+}(v))^{\intercal}v^{+}(v)\right)\right)a^{+}+\nabla f(v^{+}(v))\right)$
$\displaystyle=\frac{\alpha_{v}(v)}{{a^{-}}^{\intercal}v}\left(\underbrace{\left(\left(\nabla
f(v^{*})-\nabla
f(v^{-}(v))\right)^{\intercal}v^{-}(v)\right)}_{\delta^{-}}a^{-}+\nabla
f(v^{-}(v))\right)$
$\displaystyle\qquad+\frac{1-\alpha_{v}(v)}{{a^{+}}^{\intercal}v}\left(\underbrace{\left(\left(\nabla
f(v^{*})-\nabla
f(v^{+}(v))\right)^{\intercal}v^{+}(v)\right)}_{\delta^{+}}a^{+}+\nabla
f(v^{+}(v))\right)$
$\displaystyle=\frac{\alpha_{v}(v)}{{a^{-}}^{\intercal}v}\left(\delta^{-}a^{-}+\nabla
f(v^{-}(v))\right)+\frac{1-\alpha_{v}(v)}{{a^{+}}^{\intercal}v}\left(\delta^{+}a^{+}+\nabla
f(v^{+}(v))\right)$

Note that since $f$ is ray-concave, then it is concave on the segment
$[v^{-}(v),v^{+}(v)]$, so for any $v^{\prime}$ on the segment

$\displaystyle\left(\nabla f(v^{*})-\nabla
f(v^{-}(v))\right)^{\intercal}v^{\prime}$ $\displaystyle\leq 0$
$\displaystyle\left(\nabla f(v^{*})-\nabla
f(v^{+}(v))\right)^{\intercal}v^{\prime}$ $\displaystyle\geq 0$

because $\vec{0},v,v^{+}(v)$ and $v^{-}(v)$ are colinear. Therefore,
$\delta^{+}\geq 0$ and $\delta^{-}\leq 0$. ∎

As a side note, in the last lemma we only used differentiability of $f$ to
show the formula (15), therefore such formula is always valid for $g$ defined
as (1). Ray-concavity of $f$ was only used to show the signs of
$\delta^{\pm}$, and facet-convexity of $f$ was not needed.

###  3.2. Convexity over the polytope

We know provide the last step which proves that $g$ is convex in $P$.

######  Lemma 3.

Let $g:P\subset\mathbb{R}^{n}_{+}\rightarrow\mathbb{R}$ as defined in (1). If
$g$ is convex over each region $B\in\mathcal{B}$, then it is convex in $P$.

###### Proof.

Our strategy to show convexity is to show _mid-point local convexity_ , that
is, for each $v\in P$, we show there is a neighborhood of $v$ where $g$ is
mid-point convex. We remind the reader that mid-point convexity reads

$g\left(\frac{1}{2}(v_{1}+v_{2})\right)\leq\frac{1}{2}\left(g(v_{1})+g(v_{2})\right)\qquad\forall
v_{1},v_{2}\in P.$

Mid-point convexity does not always imply convexity, but in this case it
suffices as the function $g$ is continuous. Therefore, establishing local mid-
point convexity implies local convexity [6]. And since local convexity implies
convexity (see e.g. [10]), we conclude that $g$ is convex.

We now proceed to proving local mid-point convexity of $g$. Let us consider
$v\in P$, $d\in\mathbb{R}^{n}$ and $\varepsilon>0$. We would like to show that

(21) $g(v)\leq\frac{1}{2}(g(v-\varepsilon d)+g(v+\varepsilon d)).$

If $v\pm\varepsilon d\in\mbox{int}(B)$ the inequality follows from convexity
of $g$ within a region. Therefore, we may assume $v\in B_{s}\cap B_{t}$,
$v-\varepsilon d\in B_{s}$ and $v+\varepsilon d\in B_{t}$ for some
$B_{s},B_{t}\in\mathcal{B}$.

Let ${a_{s}^{+}}^{\intercal}x=1$ be the out-hyperplane of $B_{s}$, and
${a_{s}^{-}}^{\intercal}x=1$ be its in-hyperplane. Similarly, we define
$a^{\pm}_{t}$. Thus,
$v^{\pm}=\frac{1}{{a_{s}^{\pm}}^{\intercal}v}v=\frac{1}{{a_{t}^{\pm}}^{\intercal}v}v$.
Let $\nabla g_{B_{s}}$ and $\nabla g_{B_{t}}$ be the gradients of $g$ in
$\mbox{int}(B_{s})$ and $\mbox{int}(B_{t})$ respectively (see Lemma 2). Since
these gradients are continuous, we can extend their formula (15) to $B_{s}$
and $B_{t}$. From here, we obtain

(22) $\displaystyle\nabla(g_{B_{s}}(v)-g_{B_{t}}(v))$
$\displaystyle=\frac{\alpha_{v}}{{a_{s}^{-}}^{\intercal}v}\delta^{-}(a_{s}^{-}-a_{t}^{-})+\frac{1-\alpha_{v}}{{a_{s}^{+}}^{\intercal}v}\delta^{+}(a_{s}^{+}-a_{t}^{+})$

Now we focus on showing that $\nabla g_{B_{s}}(v)^{\intercal}d\leq\nabla
g_{B_{t}}(v)^{\intercal}d$. Since $v+\varepsilon d\in B_{t}$,

$(v+\varepsilon d)^{\pm}=\frac{v+\varepsilon
d}{{a_{t}^{\pm}}^{\intercal}(v+\varepsilon d)}.$

We start exploring the facet contained in ${a_{t}^{+}}^{\intercal}x=1$. By
convexity of the polytope $P$, we know that
${a_{s}^{+}}^{\intercal}(v+\varepsilon d)^{+}\leq 1$. Hence,
${a_{s}^{+}}^{\intercal}(v+\varepsilon
d)\leq{a_{t}^{+}}^{\intercal}(v+\varepsilon d)$ and since
${a_{s}^{+}}^{\intercal}v={a_{t}^{+}}^{\intercal}v$ we conclude that

$(a_{s}^{+}-a_{t}^{+})^{\intercal}d\leq 0.$

In a similar way, for the facet contained in ${a_{t}^{-}}^{\intercal}x=1$, by
convexity of the polytope we get that ${a_{s}^{-}}^{\intercal}(v+\varepsilon
d)^{-}\geq 1$. So, ${a_{s}^{-}}^{\intercal}(v+\varepsilon
d)\geq{a_{t}^{-}}^{\intercal}(v+\varepsilon d)$ and we conclude that

$(a_{s}^{-}-a_{t}^{-})^{\intercal}d\geq 0.$

As $\delta^{-}\leq 0$ and $\delta^{+}\geq 0$ (Lemma 2), we obtain that

$\nabla(g_{B_{s}}(v)-g_{B_{t}}(v))^{\intercal}d=\frac{\alpha_{v}}{{a_{s}^{-}}^{\intercal}v}\underbrace{\delta^{-}(a_{s}^{-}-a_{t}^{-})^{\intercal}d}_{\leq
0}+\frac{1-\alpha_{v}}{{a_{s}^{+}}^{\intercal}v}\underbrace{\delta^{+}(a_{s}^{+}-a_{t}^{+})^{\intercal}d}_{\leq
0}\leq 0.$

so we conclude that

$\nabla g_{B_{s}}(v)^{\intercal}d\leq\nabla g_{B_{t}}(v)^{\intercal}d$

Finally, we can use the first order characterization of convexity within each
region and obtain

$\displaystyle g(v+\varepsilon d)=g_{B_{t}}(v+\varepsilon d)\geq
g(v)+\varepsilon\nabla g_{B_{t}}(v)^{\intercal}d\geq g(v)+\varepsilon\nabla
g_{B_{s}}(v)^{\intercal}d$ $\displaystyle g(v-\varepsilon
d)=g_{B_{s}}(v-\varepsilon d)\geq g(v)-\varepsilon\nabla
g_{B_{s}}(v)^{\intercal}d\geq g(v)-\varepsilon\nabla
g_{B_{t}}(v)^{\intercal}d$

These two inequalities imply (21). This completes the proof of local mid-point
convexity of $g$ which, as discussed at the beginning of this proof, implies
convexity of $g$ in $P$. ∎

Note that, similarly to Lemma 2, the latter proof does not explicitly rely on
facet-convexity. The result mainly uses that $g$ is convex on each region and
that $f$ is ray-concave (in order to use the signs of $\delta^{\pm}$ in the
gradient formula).

Knowing that $g$ defines a convex function over the domain, we can prove our
main theorem, showing that it corresponds to the convex envelope of $f$ over
the polytope $P$.

###### Theorem 1.

By previous lemma, we know that $g$ is a convex function over the domain $P$.
We show that $g$ is an underestimator of $f$, that is, $g(v)\leq f(v)$ for all
$v\in P$. For $v=0$ it clearly holds. If $v\neq 0$, $v\in P$ implies that
$\alpha_{v}\in[0,1]$. Additionally, since $f$ is concave over $[v^{-},v^{+}]$
we know that

$f(v)=f(\alpha_{v}v^{-}+(1-\alpha_{v})v^{+})\geq\alpha_{v}f(v^{-})+(1-\alpha_{v})f(v^{+})=g(v).$

Finally, we argue why $g$ is the largest convex function that underestimates
$f$. Let $h$ another convex function that underestimates $f$ and let $v\in P$
such that $h(v)>g(v)$. Restricted to the segment $[v^{-},v^{+}]$, the function
$h$ is also convex. But this is a contradiction, because $f$ is concave on
$[v^{-},v^{+}]$, so the largest convex function underestimating $f$ on this
segment is the line interpolating $f(v^{-})$ and $f(v^{+})$, which is exactly
$g$. ∎

###  3.3. On the positively homogeneous condition

In this section we present characterizations for when the function $g$
constructed in (1) is positively homogeneous.

######  Lemma 4.

If $0\in P$, then $g$ is positively homogeneous if and only if $f(0)=0$. In
this case,

$g(v)={a^{+}}^{\intercal}v\cdot f(v^{+}),$

where ${a^{+}}^{\intercal}x=1$ is the out-hyperplane of the region $B\ni v$
(see Remark 2).

###### Proof.

If $0\in P$ then $v=\alpha_{v}\cdot
0+(1-\alpha_{v})v^{+}=\frac{1-\alpha_{v}}{{a^{+}}^{\intercal}v}v$, so

$g(v)=(1-{a^{+}}^{\intercal}v)f(0)+{a^{+}}^{\intercal}v\cdot f(v^{+})$

If $g$ is positively homogeneous, then $g(0)=0=f(0)$. To prove the other
direction, if $f(0)=0$ then $g(v)={a^{+}}^{\intercal}v\cdot f(v^{+})$, which
is homogeneous because for any $\lambda>0$ such that $\lambda v\in P$,
$(\lambda v)^{+}=v^{+}$ so $g(\lambda v)={a^{+}}^{\intercal}(\lambda v)\cdot
f(v^{+})=\lambda g(v)$. ∎

As mentioned in Remark 1, the condition $f(0)=0$ is not restrictive in the
construction of convex envelopes when $0\in P$. If $f(0)\neq 0$, it suffices
to define $\hat{f}=f-f(0)$ and use our construction to derive
$\hat{f}_{\text{vex}}$. The desired convex envelope simply follows from noting
that $f_{\text{vex}}=\hat{f}_{\text{vex}}+f(0)$. We illustrate the use of this
transformation in the upcoming examples section.

######  Lemma 5.

If $0\notin P$, then $g$ is positively homogeneous iff, for every $v\in P$,
${{a^{-}}^{\intercal}v}\cdot f(v^{-})={{a^{+}}^{\intercal}v}\cdot f(v^{+})$,
where ${a^{\pm}}^{\intercal}x=1$ are the in-hyperplane and out-hyperplane of a
region $B\ni v$ (see Remark 2). In this case,

$g(v)={{a^{-}}^{\intercal}v}\cdot f(v^{-})={{a^{+}}^{\intercal}v}\cdot
f(v^{+}).$

###### Proof.

Since $v^{\pm}=\frac{1}{{a^{\pm}}^{\intercal}v}v$, if $g$ is homogeneous then

$g(v)=g\left(({a^{\pm}}^{\intercal}v)\cdot
v^{\pm}\right)={a^{\pm}}^{\intercal}v\cdot
g(v^{\pm})={{a^{\pm}}^{\intercal}v}\cdot f(v^{\pm}).$

For the other direction, if ${{a^{-}}^{\intercal}v}\cdot
f(v^{-})={{a^{+}}^{\intercal}v}\cdot f(v^{+})$, by (3) we obtain

$\displaystyle g(v)$ $\displaystyle=\alpha_{v}f(v^{-})+(1-\alpha_{v})f(v^{+})$
$\displaystyle=\alpha_{v}f(v^{-})+\left(1-\frac{\alpha_{v}}{{a^{-}}^{\intercal}v}\right)({a^{+}}^{\intercal}v)f(v^{+})$
$\displaystyle=\alpha_{v}f(v^{-})+\left(1-\frac{\alpha_{v}}{{a^{-}}^{\intercal}v}\right)({a^{-}}^{\intercal}v)f(v^{-})$
$\displaystyle=\alpha_{v}f(v^{-})+({a^{-}}^{\intercal}v)f(v^{-})-\alpha_{v}f(v^{-})$
$\displaystyle=({a^{-}}^{\intercal}v)f(v^{-})=({a^{+}}^{\intercal}v)f(v^{+})$

So, $g$ is homogeneous because for any $\lambda>0$ such that $\lambda v\in P$,
$g(\lambda v)={a^{\pm}}^{\intercal}(\lambda v)\cdot f((\lambda
v)^{\pm})=\lambda({a^{\pm}}^{\intercal}v)\cdot f(v^{\pm})=\lambda g(v)$. ∎

We note that our results have an unexpected consequence: when $f$ is a
homogeneous function, convexity of $f$ over the facets of $P$ imply convexity
of $f$ over all $P$.

######  Corollary 1.

Let $f:P\rightarrow\mathbb{R}$ is continuously differentiable and convex
(concave) over the facets of $P$. If $f$ is positively homogeneous, then $f$
is convex (concave) over $P$.

###### Proof.

We show the proof for $f$ convex on the facets; the concave case is almost
identical. If $f$ is positively homogeneous then in particular is ray-linear.
Hence $f(v)=({a^{+}}^{\intercal}v)\cdot f(v^{+})=g(v)$. In addition, since $f$
is convex on the facets of $P$, by Theorem 1 $g=f_{\text{vex}}$, so $f$ is
convex over $P$. ∎

##  4\. Examples of ray-concave functions and their envelopes

In this section, we provide the convex envelopes of various explicit
functions. Some of these are new, and some have been provided in the
literature before. In the latter case, our result provides new perspectives,
and in some cases simpler derivations.

######  Example 2.

Consider the function $f(x,y)=-x\cdot y$, whose convex envelope over
$[l_{x},u_{x}]\times[l_{y},u_{y}]$ is well-known. In order to construct its
convex envelope using Theorem 1, we first shift the domain by considering the
function

$\hat{f}(x,y)=f(x+l_{x},y+l_{y})+l_{x}\cdot
l_{y}=-(x+l_{x})\cdot(y+l_{y})+l_{x}\cdot l_{y}$

over the box $[0,u_{x}-l_{x}]\times[0,u_{y}-l_{y}]$. It is easy to verify that
$\hat{f}$ is ray-concave and linear on the facet of any box
$[0,u_{x}-l_{x}]\times[0,u_{y}-l_{y}]$. Theorem 1 implies that
$\hat{f}_{\text{vex}}(v)={a^{+}}^{\intercal}v\cdot\hat{f}(v^{+})$ and thus

$\displaystyle\hat{f}_{\text{vex}}(x,y)$
$\displaystyle=\begin{cases}\frac{y}{u_{y}-l_{y}}\cdot\hat{f}\left(x\cdot\frac{u_{y}-l_{y}}{y},y\cdot\frac{u_{y}-l_{y}}{y}\right)&\text{if
}y\geq\frac{u_{y}-l_{y}}{u_{x}-l_{x}}x\\\
\frac{x}{u_{x}-l_{x}}\cdot\hat{f}\left(x\cdot\frac{u_{x}-l_{x}}{x},y\cdot\frac{u_{x}-l_{x}}{x}\right)&\text{if
}y\leq\frac{u_{y}-l_{y}}{u_{x}-l_{x}}x\\\ \end{cases}$
$\displaystyle=\begin{cases}-u_{y}x-l_{x}y&\text{if
}y\geq\frac{u_{y}-l_{y}}{u_{x}-l_{x}}x\\\ -l_{y}x-u_{x}y&\text{if
}y\leq\frac{u_{y}-l_{y}}{u_{x}-l_{x}}x\end{cases}$

And since
$f_{\text{vex}}(x,y)=\hat{f}_{\text{vex}}(x-l_{x},y-l_{y})-l_{x}\cdot l_{y}$,
we obtain

$f_{\text{vex}}(x,y)=\begin{cases}-u_{y}x-l_{x}y+l_{x}u_{y}&\text{if
}y-l_{y}\geq\frac{u_{y}-l_{y}}{u_{x}-l_{x}}(x-l_{x})\\\
-l_{y}x-u_{x}y+l_{y}u_{x}&\text{if
}y-l_{y}\leq\frac{u_{y}-l_{y}}{u_{x}-l_{x}}(x-l_{x})\end{cases}$

which corresponds to the McCormick envelopes for this function.

######  Example 3.

Let us consider the following example from [16]. Let $f(x,y)=y/x$ and

$P=\left\\{(x,y)\in\mathbb{R}^{2}:-x+2y\leq 2,1\leq x\leq 2,0\leq y\leq
2\right\\}$

We shift the domain by considering the function as $\hat{f}(x,y)=f(x+1,y)$ and
the polytope $\hat{P}=\\{(x,y)\in\mathbb{R}^{2}:(x+1,y)\in P\\}$. Note that
$\hat{f}(x,y)$ is ray-concave because $\hat{f}(x,\lambda
x)=\lambda\frac{x}{x+1}$ is concave for $x\geq 0$ and $\lambda\geq 0$.
Convexity on the facets can be directly verified. Applying Theorem 1, since
the outer facets of $\hat{P}$ are $x=1$ and $-\tfrac{1}{3}x+\tfrac{2}{3}y=1$,
we obtain

$\displaystyle\hat{f}_{\text{vex}}(x,y)$
$\displaystyle=\begin{cases}(-\tfrac{1}{3}x+\tfrac{2}{3}y)\hat{f}\left(\frac{x}{-\tfrac{1}{3}x+\tfrac{2}{3}y},\frac{y}{-\tfrac{1}{3}x+\tfrac{2}{3}y}\right)&\text{if
}y\geq 2x\\\ x\hat{f}\left(\frac{x}{x},\frac{y}{x}\right)&\text{if }y\leq
2x\\\ \end{cases}$
$\displaystyle=\begin{cases}y\frac{-\tfrac{1}{3}x+\tfrac{2}{3}y}{x+\left(-\tfrac{1}{3}x+\tfrac{2}{3}y\right)}=y\frac{-x+2y}{2x+2y}&\text{if
}y\geq 2x\\\ \frac{1}{2}y&\text{if }y\leq 2x\\\ \end{cases}$

Therefore, as $f_{\text{vex}}(x,y)=\hat{f}_{\text{vex}}(x-1,y)$ we obtain

$f_{\text{vex}}(x,y)=\begin{cases}y\frac{1-x+2y}{2(x+y+1)}&\text{if }y\geq
2(x-1)\\\ \frac{1}{2}y&\text{if }y\leq 2(x-1)\\\ \end{cases}$

######  Example 4.

Let us consider the function

(23) $f(x,y)=\frac{xy}{x+y-xy}$

in a box $[0,u_{x}]\times[0,u_{y}]\subseteq[0,1]^{2}$. This function appears
naturally in the context of network reliability optimization. In fact, if
$X,Y$ are independent Bernoulli random variables indicating the current state
of two serial component, with reliabilities $p_{X}:=\mathbb{P}(X=1)$ and
$p_{Y}:=\mathbb{P}(Y=1)$ then

$f(p_{X},p_{Y})=\mathbb{P}\left(X\cdot Y=1|X+Y\geq 1\right).$

corresponds to the resulting reliability of a _degree-2_ reduction [24].

We compute the concave envelope of (23) via the convex envelope of
$\hat{f}=-f$. The function $\hat{f}$ can be directly verified to be convex on
the facets of $[0,u_{x}]\times[0,u_{y}]$. For instance

$h(x)=\hat{f}(x,u_{y})=-\frac{xu_{y}}{x+u_{y}-xu_{y}},$

and a simple calculation shows

$h^{\prime\prime}(x)=\frac{2(1-u_{y})u_{y}^{2}}{(1-(1-u_{y})(1-x))^{3}}\geq
0.$

As for ray-concavity, we compute

$\hat{f}(x,\lambda x)=-\frac{\lambda x}{1+\lambda-\lambda
x}\quad\Rightarrow\quad\frac{\partial^{2}\hat{f}}{\partial x^{2}}(x,\lambda
x)=-\frac{2\lambda^{2}(1+\lambda)}{(1+\lambda(1-x))^{3}},$

therefore $\frac{\partial^{2}f}{\partial x^{2}}(x,\lambda x)\geq 0$ for $x\leq
1$ and $\lambda\geq 0$. By Theorem 1, the concave envelope of $f(x,y)$,
denoted $f_{\text{ave}}$, is given by

$\displaystyle f_{\text{ave}}(x,y)=-\hat{f}_{\text{vex}}(x,y)$
$\displaystyle=-({a^{+}}^{\intercal}v)\hat{f}(v^{+}_{x},v^{+}_{y})$
$\displaystyle=\begin{cases}\frac{y}{u_{y}}\frac{x\frac{u_{y}}{y}\cdot
u_{y}}{x\frac{u_{y}}{y}+u_{y}-x\frac{u_{y}}{y}\cdot u_{y}}&\text{if
}y\geq\frac{u_{y}}{u_{x}}x\\\ \frac{x}{u_{x}}\frac{u_{x}\cdot
y\frac{u_{x}}{x}}{u_{x}+y\frac{u_{x}}{x}-u_{x}\cdot y\frac{u_{x}}{x}}&\text{if
}y\leq\frac{u_{y}}{u_{x}}x\\\ \end{cases}$
$\displaystyle=\begin{cases}\frac{x\cdot y}{x+y-x\cdot u_{y}}&\text{if
}y\geq\frac{u_{y}}{u_{x}}x\\\ \frac{x\cdot y}{x+y-u_{x}\cdot y}&\text{if
}y\leq\frac{u_{y}}{u_{x}}x\\\ \end{cases}$

Note that this procedure also computes, for free, the concave envelope of $f$
on the non-rectangular polytopes
$\\{(x,y)\in[0,u_{x}]\times[0,u_{y}]\,:\,y\leq\frac{u_{y}}{u_{x}}x\\}$ and
$\\{(x,y)\in[0,u_{x}]\times[0,u_{y}]\,:\,y\geq\frac{u_{y}}{u_{x}}x\\}$.

######  Example 5.

Let us consider the function

$f(x,y)=-\frac{y\left(x^{3}y^{2}+2x^{4}y-3x^{3}y-x^{2}y+x^{5}-3x^{4}+2x^{3}-2xy^{2}-y^{3}\right)}{x(x+y)^{2}}$

over the region $P=\left\\{(x,y)\in\mathbb{R}^{2}_{+}:1\leq x+y\leq
2\right\\}$.

Note that $f(x,0)=0$, $\lim_{\epsilon\to 0}f(\epsilon,y)$ is linear,
$f(x,1-x)=(x-1)^{2}/x$ and $f(x,2-x)=(x-2)^{2}/x$, so $f$ is convex on the
facets. On the other hand, over the ray $y=\lambda x$ for $\lambda\geq 0$ we
obtain

$\frac{\partial^{2}f}{\partial x^{2}}(x,\lambda
x)=-6\lambda\frac{(1+\lambda)x-1}{1+\lambda}$

If $(x,\lambda x)\in P$ then
$x\in[\tfrac{1}{1+\lambda},\tfrac{2}{1+\lambda}]$, so $f$ is ray-concave for
any $\lambda\geq 0$. Applying Theorem 1, we obtain:

$\displaystyle f_{\text{vex}}(x,y)$ $\displaystyle=(2-x-y)\cdot
f\left(\frac{x}{x+y},\frac{y}{x+y}\right)+(x+y-1)\cdot
f\left(\frac{2x}{x+y},\frac{2y}{x+y}\right)$
$\displaystyle=(2-x-y)\cdot\frac{y^{2}}{x(x+y)}+(x+y-1)\cdot\frac{2y^{2}}{x(x+y)}=\frac{y^{2}}{x}$

which corresponds to the convex envelope because ${y^{2}}/{x}$ is positively
homogeneous.

The following example shows how Corollary 1 can be used to prove the convexity
of positively homogeneous functions.

######  Example 6.

Let $f$ be a 3-dimensional Cobb-Douglas function

$f(x_{1},x_{2},x_{3})=Ax_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}x_{3}^{\alpha_{3}}$

where $A,\alpha_{1},\alpha_{2},\alpha_{3}>0$,
$\alpha_{1}+\alpha_{2}+\alpha_{3}=1$ and $\vec{x}\in\mathbb{R}^{3}_{+}$.

It is known that the 2-dimensional Cobb-Douglas function is concave if
$\alpha_{i}+\alpha_{j}<1$, hence $f$ is concave over the facets of the box
$P=[l_{1},u_{1}]\times[l_{2},u_{2}]\times[l_{3},u_{3}]\subset\mathbb{R}^{3}_{+}$.

Since $\sum_{i=1}^{3}\alpha_{i}=1$, $f$ is positively homogeneous, so by
Corollary 1 we conclude that $f$ is concave over $P$.

## References

  * [1] Faiz A Al-Khayyal and James E Falk, _Jointly constrained biconvex programming_ , Mathematics of Operations Research 8 (1983), no. 2, 273–286.
  * [2] Kurt M Anstreicher and Samuel Burer, _Computable representations for convex hulls of low-dimensional quadratic forms_ , Mathematical Programming 124 (2010), no. 1, 33–43.
  * [3] Xiaowei Bao, Nikolaos V Sahinidis, and Mohit Tawarmalani, _Multiterm polyhedral relaxations for nonconvex, quadratically constrained quadratic programs_ , Optimization Methods & Software 24 (2009), no. 4-5, 485–504.
  * [4] Harold P Benson, _On the construction of convex and concave envelope formulas for bilinear and fractional functions on quadrilaterals_ , Computational Optimization and Applications 27 (2004), no. 1, 5–22.
  * [5] Matthias Jach, Dennis Michaels, and Robert Weismantel, _The convex envelope of (n–1)-convex functions_ , SIAM Journal on Optimization 19 (2008), no. 3, 1451–1466.
  * [6] Johan Ludwig William Valdemar Jensen, _Om konvekse funktioner og uligheder imellem middelvaerdier_ , Nyt tidsskrift for matematik 16 (1905), 49–68.
  * [7] Aida Khajavirad and Nikolaos V Sahinidis, _Convex envelopes of products of convex and component-wise concave functions_ , Journal of Global Optimization 52 (2012), no. 3, 391–409.
  * [8] by same author, _Convex envelopes generated from finitely many compact convex sets_ , Mathematical Programming 137 (2013), no. 1, 371–408.
  * [9] Takahito Kuno, _A branch-and-bound algorithm for maximizing the sum of several linear ratios_ , Journal of Global Optimization 22 (2002), no. 1, 155–174.
  * [10] Yuan-Chuan Li and Cheh-Chih Yeh, _Some characterizations of convex functions_ , Computers & Mathematics with applications 59 (2010), no. 1, 327–337.
  * [11] Leo Liberti and Constantinos C Pantelides, _Convex envelopes of monomials of odd degree_ , Journal of Global Optimization 25 (2003), no. 2, 157–168.
  * [12] Jeff Linderoth, _A simplicial branch-and-bound algorithm for solving quadratically constrained quadratic programs_ , Mathematical Programming 103 (2005), no. 2, 251–282.
  * [13] Marco Locatelli, _Polyhedral subdivisions and functional forms for the convex envelopes of bilinear, fractional and other bivariate functions over general polytopes_ , Journal of Global Optimization 66 (2016), no. 4, 629–668.
  * [14] by same author, _Convex envelopes of bivariate functions through the solution of KKT systems_ , Journal of Global Optimization 72 (2018), no. 2, 277–303.
  * [15] by same author, _Convex envelope of bivariate cubic functions over rectangular regions_ , Journal of Global Optimization 76 (2020), no. 1, 1–24.
  * [16] Marco Locatelli and Fabio Schoen, _On convex envelopes for bivariate functions over polytopes_ , Mathematical Programming 144 (2014), no. 1, 65–91.
  * [17] James Luedtke, Mahdi Namazifar, and Jeff Linderoth, _Some results on the strength of relaxations of multilinear functions_ , Mathematical Programming 136 (2012), no. 2, 325–351.
  * [18] Garth P McCormick, _Computability of global solutions to factorable nonconvex programs: Part I – Convex underestimating problems_ , Mathematical Programming 10 (1976), no. 1, 147–175.
  * [19] Clifford A Meyer and Christodoulos A Floudas, _Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes_ , Journal of Global Optimization 29 (2004), no. 2, 125–155.
  * [20] by same author, _Convex envelopes for edge-concave functions_ , Mathematical Programming 103 (2005), no. 2, 207–224.
  * [21] Benjamin Muller, Felipe Serrano, and Ambros Gleixner, _Using two-dimensional projections for stronger separation and propagation of bilinear terms_ , SIAM Journal on Optimization 30 (2020), no. 2, 1339–1365.
  * [22] Anatoliy D Rikun, _A convex envelope formula for multilinear functions_ , Journal of Global Optimization 10 (1997), no. 4, 425–437.
  * [23] Hong Seo Ryoo and Nikolaos V Sahinidis, _Analysis of bounds for multilinear functions_ , Journal of Global Optimization 19 (2001), no. 4, 403–424.
  * [24] Appajosyula Satyanarayana and R Kevin Wood, _A linear-time algorithm for computing k-terminal reliability in series-parallel networks_ , SIAM Journal on Computing 14 (1985), no. 4, 818–832.
  * [25] Hanif D Sherali, _Convex envelopes of multilinear functions over a unit hypercube and over special discrete sets_ , Acta Mathematica Vietnamica 22 (1997), no. 1, 245–270.
  * [26] Hanif D Sherali and Amine Alameddine, _An explicit characterization of the convex envelope of a bivariate bilinear function over special polytopes_ , Annals of Operations Research 25 (1990), no. 1, 197–209.
  * [27] Fabio Tardella, _On the existence of polyhedral convex envelopes_ , Frontiers in Global Optimization, Springer, 2004, pp. 563–573.
  * [28] by same author, _Existence and sum decomposition of vertex polyhedral convex envelopes_ , Optimization Letters 2 (2008), no. 3, 363–375.
  * [29] Mohit Tawarmalani, Jean-Philippe P Richard, and Chuanhui Xiong, _Explicit convex and concave envelopes through polyhedral subdivisions_ , Mathematical Programming 138 (2013), no. 1, 531–577.
  * [30] Mohit Tawarmalani and Nikolaos V Sahinidis, _Semidefinite relaxations of fractional programs via novel convexification techniques_ , Journal of Global Optimization 20 (2001), no. 2, 133–154.
  * [31] Juan M Zamora and Ignacio E Grossmann, _A branch and contract algorithm for problems with concave univariate, bilinear and linear fractional terms_ , Journal of Global Optimization 14 (1999), no. 3, 217–249.


11institutetext:  Purple Mountain Observatory, Chinese Academy of Sciences,
No. 10 Yuanhua Road, Nanjing 210033, China
11email<EMAIL_ADDRESS>22institutetext: School of Astronomy and Space
Science, University of Science and Technology of China, No. 96 Jinzhai Road,
Hefei, Anhui 230026, China 33institutetext: Hunan Astronomical Association,
Changsha, Hunan 410000, China 44institutetext: Department of Space Sciences
and Astronomy, Hebei Normal University, No. 20 Road East. 2nd Ring South,
Shijiazhuang, Hebei 050024, China 55institutetext: Shenzhen Astronomical
Observatory, Tianwen Road, Shenzhen, Guangdong 518040, China 66institutetext:
Shenzhen Astronomical Society, 22c Seascape Square Taizi Road, Shenzhen,
Guangdong 518040, China 77institutetext: Shanghai Astronomical Observatory,
Chinese Academy of Sciences, No. 80 Nandan Road, Shanghai 200030, China
88institutetext: Nanjing Amateur Astronomers Association, Nanjing 210000,
China

# Reconciling results of 2019 and 2020 stellar occultations on Pluto’s
atmosphere

New constraints from both the 5 September 2019 event and consistency analysis

Ye Yuan (袁烨) 11 0000-0002-4686-8548 Fan Li (李凡) 11 0000-0002-2997-3098 Yanning
Fu (傅燕宁) 11 Jian Chen (陈健) 1122 0009-0006-9369-0806 Wei Tan (谭巍) 33 Shuai
Zhang (张帅) 44 0000-0002-4799-0780 Wei Zhang (张伟) 11 Chen Zhang (张晨) 1122
0000-0002-9583-263X Qiang Zhang (张强) 44 Jiahui Ye (叶嘉晖) 55 Delai Li (李德铼) 55
Yijing Zhu (朱一静) 66 Zhensen Fu (傅震森) 7744 Ansheng Zhu (朱安生) 1122
0009-0007-5505-526X Yue Chen (陈悦) 1122 0000-0002-6050-9920 Jun Xu (许军) 88 Yang
Zhang 11

A stellar occultation by Pluto on 5 September 2019 yielded positive detections
at two separate stations. Using an approach consistent with comparable
studies, we derived a surface pressure of $11.478\pm 0.55\leavevmode\nobreak\
\mathrm{\upmu bar}$ for Pluto’s atmosphere from the observations of this
event. In addition, to avoid potential method inconsistancies highlighted by
Sicardy et al. when comparing with historical pressure measurements, we
reanalyzed the data by 15 August 2018 and 17 July 2019 events, respectively.
All the new measurements provide a bridge between the two different
perspectives on the pressure variation since 2015: a rapid pressure drop from
previous studies of the 15 August 2018 and 17 July 2019 events and a plateau
phase from that of the 6 June 2020 event. The pressure measurement from the 5
September 2019 event aligns with those from 2016, 2018, and 2020, supporting
the latter perspective. While the measurements from the 4 June 2011 and 17
July 2019 events suggest probable V-shaped pressure variations unaccounted for
by the volatile transport model (VTM) from Meza et al., the VTM remains
applicable on average. And, the validity of the V-shaped variations is
debatable due to the stellar faintness of the 4 June 2011 event and the
grazing single-chord geometry of the 17 July 2019 event. To reveal and
understand all significant pressure variations of Pluto’s atmosphere, it is
essential to provide constraints on both short-term and long-term evolutions
of the interacting atmosphere and surface by continuous pressure monitoring
through occultation observations, whenever possible, complemented by frequent
spectroscopy and photometry of the surface.

###### Key Words.:

Kuiper belt objects: individual: Pluto – planets and satellites: atmospheres –
planets and satellites: physical evolution – occultations – techniques:
photometric

##  1 Introduction

Pluto’s atmosphere was discovered during the 1985 stellar occultation (Brosch,
1995), and since then, stellar occultations have played a crucial role in
studying its structure, composition, and how it evolves over time (Hubbard et
al., 1988; Elliot et al., 1989, 2003; Yelle & Elliot, 1997; Sicardy et al.,
2003, 2011, 2016, 2021; Pasachoff et al., 2005, 2017; Young et al., 2008,
2021; Rannou & Durry, 2009; Person et al., 2013, 2021; Olkin et al., 2015;
Bosh et al., 2015; Gulbis et al., 2015; Dias-Oliveira et al., 2015; Meza et
al., 2019; Arimatsu et al., 2020).

A compilation of twelve occultations observed between 1988 and 2016 revealed a
three-fold monotonic increase in the atmospheric pressure of Pluto during that
period (Meza et al., 2019). This increase can be explained by the volatile
transport model (VTM) of the Laboratoire de Météorologie Dynamique (LMD)
(Bertrand & Forget, 2016; Forget et al., 2017; Bertrand et al., 2018, 2019)
fine-turned by Meza et al. (2019). This model provides a framework for
simulating the volatile cycles on Pluto’s globe over both seasonal and
astronomical timescales, allowing for an understanding of the long-term
evolution of Pluto’s atmosphere and its response to seasonal variations over
its 248-year heliocentric orbital period (Meza et al., 2019). According to the
LMD VTM in Meza et al. (2019) (VTM19, hereafter), Pluto’s atmospheric pressure
is expected to reach its peak around 2020. The pressure increase is attributed
to the progression of summer over the northern hemisphere of Pluto, exposing
Sputnik Planitia (SP)111Sputnik Planitia is a large plain or basin on Pluto
discovered by the New Horizons spacecraft (Stern et al., 2015). to solar
radiation. The surface of SP, composed of nitrogen ($\mathrm{N_{2}}$), methane
($\mathrm{CH_{4}}$), and carbon monoxide ($\mathrm{CO}$) ices, is believed to
sublimate and release volatile gases into the atmosphere during this period,
leading to a pressure increase. After reaching its peak, the model predicts a
gradual decline in pressure over the next two centuries under the combined
effects of Pluto’s recession from the Sun and the prevalence of the winter
season over SP.

On one hand, the VTM19 remains consistent with Sicardy et al. (2021)’s
analysis of the 6 June 2020 occultation observed at Devasthal, where two co-
located telescopes were used. Their analysis suggests that Pluto’s atmosphere
has been in a plateau phase since mid-2015. This aligns with the model
predictions that the atmospheric pressure reached its peak around 2020.

On the other hand, Arimatsu et al. (2020)’s analysis of the 17 July 2019
occultation observed by a single telescope (TUHO) suggests a rapid pressure
decrease between 2016 and 2019. They detected a significant pressure drop at
the $2.4\sigma$ level. However, it is worth noting that the geometry of this
occultation is grazing. This may have introduced larger correlations between
the pressure and the geocentric closest approach distance to Pluto’s shadow
axis, leading to the lack of the necessary precision to confidently support
the claim of a large pressure decrease followed by a return to the pressure
level close to that of 2015 in 2020 (Sicardy et al., 2021).

These contrasting results highlight the need for occultation observations
between 2019 and 2020 to better understand the behavior and evolution of
Pluto’s atmosphere during this time period. Furthermore, while Young et al.
(2021) support the presence of a pressure drop based on their analysis of the
15 August 2018 occultation, Sicardy et al. (2021) suggest that careful
comparisons between measurements by independent teams should be made before
drawing any conclusions on the pressure evolution.

Observations of the 5 September 2019 occultation which has not been reported
by other teams are presented in Section 2, followed by a description of the
light-curve fitting methods in Section 3. These unique observations allow us
to track the changes in Pluto’s atmosphere during the time period between the
events studied by Arimatsu et al. (2020) and Sicardy et al. (2021). Results
are detailed in Section 4, and the pressure evolution is discussed in Section
5, including comparisons with the reanalyzed 15 August 2018 and 17 July 2019
events. Conclusions and recommendations are provided in Section 6.

##  2 Occultation observations

Two observation campaigns were organized in China for occultations in 2019
(see Appendix A). One occurred on 17 July 2019, which has been studied by
Arimatsu et al. (2020), and the other on 5 September 2019, which is reported
in the present paper for the first time. Due to bad weather conditions in many
areas, no effective light curves were observed by our stations for the first
occultation, and only two for the second.

Table 1 lists the circumstances of the 5 September 2019 event. Figure 1
presents all the observation stations and the reconstructed path of the shadow
of Pluto222The occulted star is Gaia DR3 $6771712487062767488$, of which the
astrometric and photometric parameters are obtained from VizieR (Gaia
Collaboration, 2022). during this event. Table 3 lists the circumstances of
stations with positive detections. Their station codes are DWM and HNU.

Table 1: Circumstances and light-curve fitting results of the 5 September 2019
event.

Occulted star
---
Identification (Gaia DR3aa$a$Gaia Collaboration et al. (2022), used to derive
$(\alpha_{\mathrm{s}},\delta_{\mathrm{s}})$.
) | $6771712487062767488$
Geocentric astrometric position | $\alpha_{\mathrm{s}}=19^{\mathrm{h}}29^{\mathrm{m}}11\aas@@fstack{s}1996$
at observational epoch (ICRFbb$b$International Celestial Reference Frame.
) | $\delta_{\mathrm{s}}=-22\degr 21^{\prime}39\aas@@fstack{\prime\prime}880$
Pluto’s body
Masscc$c$Stern et al. (2015), where $G$ is the gravitational constant.
, $GM_{\mathrm{p}}$ ($\mathrm{km}^{3}\cdot\mathrm{s}^{-2}$) | $869.6$
Radiuscc$c$Stern et al. (2015), where $G$ is the gravitational constant.
, $R_{\mathrm{p}}$ ($\mathrm{km}$) | $1187$
Pluto’s atmosphere
$\mathrm{N}_{2}$ molecular massdd$d$Assumed to be the only constituent in the
light-curve model, following, e.g., Sicardy et al. (2021).
, $\mu$ (kg) | $4.652\times 10^{-26}$
$\mathrm{N}_{2}$ molecular | $1.091\times 10^{-23}$
refractivityee$e$Washburn et al. (1930), where $\lambda_{\mathrm{\upmu m}}$ is
the wavelength expressed in microns. For stations use no filters (clear), a
nominal central wavelength ($\lambda_{\mathrm{\upmu m}}=0.55$) was used.
, $K$ ($\mathrm{cm}^{3}$) | $\leavevmode\nobreak\ \leavevmode\nobreak\ +6.282\times 10^{-26}/\lambda_{\mathrm{\upmu m}}^{2}$
Boltzmann constantff$f$Newell & Tiesinga (2019).
, $k_{\mathrm{B}}$ ($\mathrm{J}\cdot\mathrm{K}^{-1}$) | $1.380649\times 10^{-23}$
Given reference radiusgg$g$Using a ratio $p_{\mathrm{surf}}/p_{0}=1.837$ given
by the template model (Meza et al., 2019; Sicardy et al., 2021).
, $r_{0}$ (km) | $1215$
Results of atmospheric fit (with $1\sigma$ error bars)
Pressure at $r_{0}$, $p_{0}$ ($\mathrm{\upmu bar}$) | $6.248\pm 0.30$
Surface pressuregg$g$Using a ratio $p_{\mathrm{surf}}/p_{0}=1.837$ given by
the template model (Meza et al., 2019; Sicardy et al., 2021).
at $R_{\mathrm{p}}$, $p_{\mathrm{surf}}$ ($\mathrm{\upmu bar}$) | $11.478\pm 0.55$
Geocentric closest approach distance |
to shadow centerhh$h$Positive (resp. negative) values mean that the shadow
center passes north (resp. south) of the geocenter.
, $\rho_{\mathrm{cag}}$ (km) | $+3644\pm 25$
Geocentric closest approach time |
to shadow centerii$i$Timings by QHY174GPS camera are used as time references,
considering their reliability and accuracy.
, $t_{\mathrm{cag}}$ (UTCjj$j$Coordinated Universal Time.) |  $15\mathrm{:}01\mathrm{:}19.1\pm 0.38$ s
333

Figure 1: Reconstructed occultation map of the 5 September 2019 event.

To ensure accurate and precise timing in stellar occultations, some stations
(e.g., DWM as shown in Table 3) were equipped with QHY174GPS cameras. These
cameras, manufactured by QHYCCD444https://www.qhyccd.com, offer precise
recording of observation time and location for each frame using a GPS-based
function. And, they have been used in many stellar occultation studies (e.g.,
Buie et al., 2020a, b; Morgado et al., 2021, 2022; Pereira et al., 2023). In
the light-curve fitting procedures described in Section 3.2, the time
recording offsets of the QHY174GPS cameras are fixed to zero, considering
their reliability and accuracy as time references.

All observational data were captured in the FITS format. These data were
processed using the Tangra occultation photometric
tool555http://www.hristopavlov.net/Tangra3/ (Pavlov, 2020) and our data
reduction code (see Appendix B). It was ensured that the targets and reference
stars in all the images we used were not overexposed. The resulting light
curves from the observations, after being normalized, are presented in Figure
2. Each data point on the light curves is represented by
$f_{i}(t)\pm\sigma_{i}(t)$, where $i$ indicates the quantities associated with
a specific station, $t$ represents the recorded timing of each frame, $f$ the
normalized total observed flux of the occulted star and the Pluto’s system,
and $\sigma$ the measurement error associated with each data point.

(a)

(b)

(c)

Figure 2: Occultation observations and the best-fitting light-curve model of
the 5 September 2019 event. Panel a: observed and simultaneously fitted light
curves. Panels b and c: Reconstructed stellar paths seen by DWM and HNU,
respectively.

##  3 Light-curve fitting methods

###  3.1 Light-curve model

In order to simulate observed light curves, we implemented a light-curve
model, $\phi(t;A,s,\Delta t,\Delta\tau,\Delta\rho,p_{0})$, described in
Appendix C and consistent with the DO15 (Dias-Oliveira et al., 2015; Sicardy
et al., 2016; Meza et al., 2019; Sicardy et al., 2021). As a function of model
parameters, its time-dependent Jacobian matrix was also implemented to
represent the sensitivity of the model with respect to the corresponding
parameters that will be estimated through fitting procedures.

The light-curve model of a given station can be formally written as:

$\phi_{i}(t)=A_{i}\cdot\left(s_{i}\cdot\psi_{i}(t;\Delta
t_{i},\Delta\tau,\Delta\rho,p_{0})+(1-s_{i})\right),$ (1)

where $i$ indicates the quantities associated with the station, and for
details the reader is referred to Appendix C. Here, the reference ephemerides
we use are the NIMAv9666https://lesia.obspm.fr/lucky-star/obj.php?p=818
asteroidal ephemeris (Desmars et al., 2015, 2019) for the orbit of the Pluto
system barycenter with respect to the Sun, the
PLU058777https://ssd.jpl.nasa.gov/ftp/eph/satellites/bsp/plu058.bsp satellite
ephemerides (Brozović et al., 2015; Jacobson et al., 2019) for the orbit of
Pluto with respect to the Pluto system barycenter, and the
DE440888https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440.bsp planetary
ephemerides (Park et al., 2021) for the orbits of the Earth and the Sun with
respect to the Solar system barycenter. The reference star catalog where the
data of the occulted star are obtained is the Gaia DR3.

###  3.2 Fitting procedure

The light-curve model was fitted to the normalized observed light curves
simultaneously by nonlinear least squares, returning a $\chi^{2}$ type value
of goodness-of-fit. The goal was to minimize the objective function given by:

$\chi^{2}_{\mathrm{obs}}=\sum_{i,j}\frac{(\phi_{i}(t_{ij})-f_{i}(t_{ij}))^{2}}{\sigma_{i}^{2}(t_{ij})},$
(2)

where $t_{ij}$ represents the mid-exposure time of the $j$-th observation of
the station $i$.

In addition, with the used reference ephemerides and star catalog, a piece of
a priori information on $\Delta\rho$ can be obtained:

$\Delta\rho=0\leavevmode\nobreak\ \mathrm{km}\pm\sigma_{\rho},$ (3)

where the uncertainty $\sigma_{\rho}$ is set to $72$ km using the positional
uncertainties listed in NIMAv9’s “orbit quality” table and in Gaia DR3. This
$\sigma_{\rho}$ value corresponds to about $3$ mas on the sky at the
geocentric distance of Pluto. The a priori information can be treated as an
independent observational data and used in the model fitting, with the
objective function modified as:

$\chi^{2}_{\mathrm{apr}}=\chi^{2}_{\mathrm{obs}}+\frac{\Delta\rho^{2}}{\sigma_{\rho}^{2}}.$
(4)

The fitting steps are as following:

  * •

In order to find all local minima, at some of which a nonlinear least-squares
fitting could get stuck, we explored the two parameter space
$(\Delta\rho,p_{0})$ by generating the variation of $\chi^{2}$ as a function
of them. Figure 3 presents such two $\chi^{2}$ maps, labeled (a) and (b),
which will be analyzed in Section 4. The maps are generated by minimizing
$\chi^{2}_{\mathrm{obs}}$ or $\chi^{2}_{\mathrm{apr}}$ at each fixed
$(\Delta\rho,p_{0})$ point on a regularly spaced grid. The Levenberg-Marquardt
(LM) method, implemented in the LMFIT package999https://lmfit.github.io/lmfit-
py/, was used in each fitting procedure. The free parameters to be adjusted
are $\Delta\tau$101010As a technical note, if no equipments for reliable time
references are used, $\Delta\tau$ would be fixed as zero, avoiding fitting
problems caused by its linear relationship with $t_{i}$., $\Delta t_{i}$ of
any station having no reliable time reference system like QHY174GPS, and
$s_{i}$ and $A_{i}$ of each station.

  * •

For a more accurate best-fitting solution for $(\Delta\rho,p_{0})$, the LM
method is used again, with $\Delta\rho$ and $p_{0}$ adjusted with initial
guesses located at all known local minima of each $\chi^{2}$ map.

  * •

Each $\chi^{2}$ map, which provides information about the quality of the fit,
is used to define confidence limits based on constant $\chi^{2}$ boundaries
(Press et al., 2007).

(a)

(b)

Figure 3: The $\chi^{2}$ maps of the 5 September 2019 event. Panels a and b:
the $\chi^{2}_{\mathrm{obs}}$ and $\chi^{2}_{\mathrm{apr}}$ maps,
respectively. The $\chi^{2}_{\mathrm{obs}}$ denotes the goodness-of-fit only
using observational data, while the $\chi^{2}_{\mathrm{apr}}$ with the
additional a priori information represented as the $\chi^{2}$ type value,
$(\Delta\rho/\sigma_{\rho})^{2}$. These maps are used to derive the best-
fitting atmospheric pressure $p_{0}$ at the reference radius $r_{0}$ of 1215
km and the cross track correction $\Delta\rho$ to the ephemerides, of which
the a priori uncertainty $\sigma_{\rho}$ is $72$ km. The surface pressure
$p_{\mathrm{surf}}$ and the geocentric closest approach distance
$\rho_{\mathrm{cag}}$ to the shadow center are obtained by linear
transformations of $\Delta\rho$ and $p_{0}$, respectively. The best-fitting
$\chi^{2}_{\mathrm{obs}}$ and $\chi^{2}_{\mathrm{apr}}$ values per degree of
freedom are $1.160$ and $1.155$, respectively.

##  4 Results

Figure 3(a) shows the $\chi^{2}_{\mathrm{obs}}$ map for the 5 September 2019
occultation. Two local minima are observed. However, considering the
significant $\chi^{2}$ difference of $9$ between the two local minima, the
global minimum is more likely the correct solution. In addition, the
$\Delta\rho$ value at the global minimum is more consistent with the NIMAv9
solution, $\Delta\rho=0$ km, at the $0.16\leavevmode\nobreak\ \sigma_{\rho}$
level, compared with the other local one at the $2.44\leavevmode\nobreak\
\sigma_{\rho}$ level.

In an effort to mitigate or, at least, further weaken the presence of multiple
local minima, we calculated the $\chi^{2}_{\mathrm{apr}}$ map by adding the
$\chi^{2}$ type value of the a priori informaiton,
${(\Delta\rho/\sigma_{\rho})^{2}}$, into the $\chi^{2}_{\mathrm{obs}}$ map.
Figure 3(b) presents the results, showing two local minima still being
present, but with a $\chi^{2}$ difference of about 14.5, which is larger than
that of the $\chi^{2}_{\text{obs}}$ map. So, the global minimum is confidently
accepted as the solution for $(\rho_{\mathrm{cag}},p_{\mathrm{surf}})$, as
provided in Table 1.

Moreover, Figure 3 presents the consistency of our derived $p_{\mathrm{surf}}$
across the two different local minima. It demonstrates that the specific
choice of local minima does not significantly affect the value of
$p_{\mathrm{surf}}$, further supporting the reliability of our solution for
$p_{\mathrm{surf}}$.

##  5 Pressure evolution

###  5.1 Comparisons and necessary reanalyses of historical events

In Figure 4, the red plot represents our $p_{\mathrm{surf}}$ measurement from
the 5 September 2019 occultation. We also include other published measurements
(Hinson et al., 2017; Meza et al., 2019; Arimatsu et al., 2020; Young et al.,
2021; Sicardy et al., 2021) and the pressure evolution predicted by the VTM19
to provide a comprehensive view of the pressure variations on Pluto.

Figure 4: The pressure evolution of Pluto over a 248-year heliocentric orbital
period predicted by the VTM19, along with the measured pressures with
$1\sigma$ error bars. The $3\sigma$ error bars of our new measurements and of
the previous measurements of the 19 July 2016, 17 July 2019, and 6 June 2020
events are also presented, using the published $\chi^{2}$ maps.
$A_{\mathrm{N_{2}}}$ denotes the albedo of nitrogen ice.

To avoid potential inconsistencies arising from different analysis methods, as
discussed by Sicardy et al. (2021), we reanalyzed the 15 August 2018 event
studied by Young et al. (2021) using the IXON observational data of Silva-
Cabrera et al. (2022). The derived pressure measurement, presented in Appendix
D.1, is ${12.027}_{-0.08}^{+0.09}\leavevmode\nobreak\ \mathrm{\upmu bar}$. In
addition, we also reanalyzed the 17 July 2019 event in Appendix D.2, deriving
a pressure of $p_{\mathrm{surf}}=9.421_{-0.75}^{+0.68}\leavevmode\nobreak\
\mathrm{\upmu bar}$, similar to that of Arimatsu et al. (2020),
${9.56}_{-0.34}^{+0.52}\leavevmode\nobreak\ \mathrm{\upmu bar}$. This
similarity is expected because the same DO15 method is used. Since the method
is used by Meza et al. (2019) and Sicardy et al. (2021), their pressure
measurements, along with that of Arimatsu et al. (2020), can be fully compared
with our new ones. Both the remeasurements are plotted in black in Figure 4.

The pressure measurement from the 5 September 2019 event shows alignments with
those from the 19 July 2016, 15 August 2018, and 6 June 2020 events within
their combined $1\sigma$ levels. Our new measurement from the 15 August 2018
event shows no significant pressure drop previously reported by Young et al.
(2021). The previously reported pressure drop between the 19 July 2016 and 17
July 2019 events is still detected at the same level as in Arimatsu et al.
(2020).

###  5.2 Discussions on pressure variations

While the VTM19, on average, remains applicable and capable of predicting the
main atmospheric behavior during the observed years, there are also two
probable V-shaped pressure variations observed from 2010 to 2015 and from 2015
to 2020, especially when considering the measurements from the 4 June 2011 and
17 July 2019 events. These V-shaped variations suggest the presence of
additional and unaccounted factors. Specifically, short-term changes in
Pluto’s surface ices and their interaction with the atmosphere are likely
contributing to the variation. In fact, spectral monitoring of the surface
composition has revealed some short-term changes in the ices over several
Earth years (e.g., Grundy et al., 2014; Lellouch et al., 2022; Holler et al.,
2022).

However, the validity of the V-shaped variations is debatable due to the
stellar faintness of the 4 June 2011 event and the grazing single-chord
geometry of 17 July 2019 event. If the debatable measurement from the 17 July
2019 event is discarded, no significant changes would be observed between 2016
and 2020. This more likely supports the plateau phase since 2015 predicted by
the VTM19. In order to better understand the relationship between these
factors, further observations using multiple observational techniques
(occultation, spectroscopy, and photometry), as well as simulations with a
refined VTM, are required.

##  6 Conclusions

The unique observations of the 5 September 2019 occultation have provided a
surface pressure of $p_{\mathrm{surf}}=11.478\pm 0.55\leavevmode\nobreak\
\mathrm{\upmu bar}$. In order to avoid potential method inconsistencies in
comparing with historical pressure measurements (Sicardy et al., 2021), we
also reanalyzed the 15 August 2018 and 17 July 2019 events based on publicly
available data (Silva-Cabrera et al., 2022; Arimatsu et al., 2020). All
measurements are presented in Figure 4.

The VTM19 remains applicable on average. Besides, we also observed unaccounted
V-shaped pressure variations with the previously reported pressure drop being
a part of these variations, which, however, are debatable. To better
understand all significant pressure variations of Pluto, continuous pressure
monitoring through occultation observations, whenever possible, is essential.
Also, simultaneous and frequent spectroscopic and photometric monitoring of
its surface ice changes are important, as such comprehensive monitoring will
provide more short-term and long-term evolution constraints of Pluto’s
interacting atmosphere and surface.

###### Acknowledgements.

We acknowledge Bruno Sicardy for his useful comments that helped improving
this manuscript. This work has been supported by the National Natural Science
Foundation of China (Grant Nos. 12203105 and 12103091). We acknowledge the
science research grants from the China Manned Space Project with NO.CMS-
CSST-2021-A12 and NO.CMS-CSST-2021-B10. This research has made use of data
from the 40 cm DOB telescope at the Dawei Mountain Observatory of Hunan
Astronomical Association and from the 50 cm RC telescope at the Observatory of
Hebei Normal University. We ackonwledge the support of Chinese amateur
astronomers from Hunan Astronomical Association, Nanjing Amateur Astronomers
Association, and Shenzhen Astronomical Observatory Team.

## References

  * Arimatsu et al. (2020) Arimatsu, K., Hashimoto, G. L., Kagitani, M., et al. 2020, A&A, 638, L5
  * Bertrand & Forget (2016) Bertrand, T. & Forget, F. 2016, Nature, 540, 86
  * Bertrand et al. (2018) Bertrand, T., Forget, F., Umurhan, O. M., et al. 2018, Icarus, 309, 277
  * Bertrand et al. (2019) Bertrand, T., Forget, F., Umurhan, O. M., et al. 2019, Icarus, 329, 148
  * Bosh et al. (2015) Bosh, A. S., Person, M. J., Levine, S. E., et al. 2015, Icarus, 246, 237
  * Brosch (1995) Brosch, N. 1995, MNRAS, 276, 571
  * Brozović et al. (2015) Brozović, M., Showalter, M. R., Jacobson, R. A., & Buie, M. W. 2015, Icarus, 246, 317
  * Buie et al. (2020a) Buie, M. W., Leiva, R., Keller, J. M., et al. 2020a, AJ, 159, 230
  * Buie et al. (2020b) Buie, M. W., Porter, S. B., Tamblyn, P., et al. 2020b, AJ, 159, 130
  * Desmars et al. (2015) Desmars, J., Camargo, J. I. B., Braga-Ribas, F., et al. 2015, A&A, 584, A96
  * Desmars et al. (2019) Desmars, J., Meza, E., Sicardy, B., et al. 2019, A&A, 625, A43
  * Dias-Oliveira et al. (2015) Dias-Oliveira, A., Sicardy, B., Lellouch, E., et al. 2015, ApJ, 811, 53
  * Elliot et al. (2003) Elliot, J. L., Ates, A., Babcock, B. A., et al. 2003, Nature, 424, 165
  * Elliot et al. (1989) Elliot, J. L., Dunham, E. W., Bosh, A. S., et al. 1989, Icarus, 77, 148
  * Forget et al. (2017) Forget, F., Bertrand, T., Vangvichith, M., et al. 2017, Icarus, 287, 54
  * Gaia Collaboration (2022) Gaia Collaboration. 2022, VizieR Online Data Catalog, I/355
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211
  * Grundy et al. (2014) Grundy, W. M., Olkin, C. B., Young, L. A., & Holler, B. J. 2014, Icarus, 235, 220
  * Gulbis et al. (2015) Gulbis, A. A. S., Emery, J. P., Person, M. J., et al. 2015, Icarus, 246, 226
  * Hinson et al. (2017) Hinson, D. P., Linscott, I. R., Young, L. A., et al. 2017, Icarus, 290, 96
  * Holler et al. (2022) Holler, B. J., Yanez, M. D., Protopapa, S., et al. 2022, Icarus, 373, 114729
  * Hubbard et al. (1988) Hubbard, W. B., Hunten, D. M., Dieters, S. W., Hill, K. M., & Watson, R. D. 1988, Nature, 336, 452
  * Jacobson et al. (2019) Jacobson, R. A., Brozovic, M., Showalter, M., et al. 2019, in Pluto System After New Horizons, Vol. 2133, 7031
  * Lellouch et al. (2022) Lellouch, E., Butler, B., Moreno, R., et al. 2022, Icarus, 372, 114722
  * Meza et al. (2019) Meza, E., Sicardy, B., Assafin, M., et al. 2019, A&A, 625, A42
  * Morgado et al. (2022) Morgado, B. E., Gomes-Júnior, A. R., Braga-Ribas, F., et al. 2022, AJ, 163, 240
  * Morgado et al. (2021) Morgado, B. E., Sicardy, B., Braga-Ribas, F., et al. 2021, A&A, 652, A141
  * Newell & Tiesinga (2019) Newell, D. & Tiesinga, E. 2019, in The International System of Units (SI)
  * Olkin et al. (2015) Olkin, C. B., Young, L. A., Borncamp, D., et al. 2015, Icarus, 246, 220
  * Park et al. (2021) Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. 2021, AJ, 161, 105
  * Pasachoff et al. (2017) Pasachoff, J. M., Babcock, B. A., Durst, R. F., et al. 2017, Icarus, 296, 305
  * Pasachoff et al. (2005) Pasachoff, J. M., Souza, S. P., Babcock, B. A., et al. 2005, AJ, 129, 1718
  * Pavlov (2020) Pavlov, H. 2020, Tangra: Software for video photometry and astrometry, Astrophysics Source Code Library, record ascl:2004.002
  * Pereira et al. (2023) Pereira, C. L., Sicardy, B., Morgado, B. E., et al. 2023, A&A, 673, L4
  * Person et al. (2021) Person, M. J., Bosh, A. S., Zuluaga, C. A., et al. 2021, Icarus, 356, 113572
  * Person et al. (2013) Person, M. J., Dunham, E. W., Bosh, A. S., et al. 2013, AJ, 146, 83
  * Press et al. (2007) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical recipes in C++ : the art of scientific computing (3rd Edition), 3rd edn. (Cambridge University Press)
  * Rannou & Durry (2009) Rannou, P. & Durry, G. 2009, Journal of Geophysical Research (Planets), 114, E11013
  * Sicardy (2022) Sicardy, B. 2022, arXiv e-prints, arXiv:2206.06236
  * Sicardy et al. (2021) Sicardy, B., Ashok, N. M., Tej, A., et al. 2021, ApJ, 923, L31
  * Sicardy et al. (2011) Sicardy, B., Bolt, G., Broughton, J., et al. 2011, AJ, 141, 67
  * Sicardy et al. (2016) Sicardy, B., Talbot, J., Meza, E., et al. 2016, ApJ, 819, L38
  * Sicardy et al. (2003) Sicardy, B., Widemann, T., Lellouch, E., et al. 2003, Nature, 424, 168
  * Silva-Cabrera et al. (2022) Silva-Cabrera, J. S., Castro-Chacón, J. H., Reyes-Ruiz, M., et al. 2022, MNRAS, 511, 5550
  * Stern et al. (2015) Stern, S. A., Bagenal, F., Ennico, K., et al. 2015, Science, 350, aad1815
  * Washburn et al. (1930) Washburn, E. W., West, C. J., Dorsey, N. E., et al. 1930, International Critical Tables of Numerical Data, Physics, Chemistry and Technology (Washington, DC: The National Academies Press)
  * Yelle & Elliot (1997) Yelle, R. V. & Elliot, J. L. 1997, in Pluto and Charon, ed. S. A. Stern & D. J. Tholen, 347
  * Young et al. (2021) Young, E., Young, L. A., Johnson, P. E., & PHOT Team. 2021, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 53, AAS/Division for Planetary Sciences Meeting Abstracts, 307.06
  * Young et al. (2008) Young, E. F., French, R. G., Young, L. A., et al. 2008, AJ, 136, 1757

##  Appendix A The two occultation campaigns in 2019

The 17 July 2019111111https://lesia.obspm.fr/lucky-star/occ.php?p=13163 and
the 5 September 2019121212https://lesia.obspm.fr/lucky-star/occ.php?p=13166
stellar occultations by Pluto were originally predicted by the ERC Lucky Star
project131313https://lesia.obspm.fr/lucky-star.

###  A.1 The 17 July 2019 occultation campaign

Figure 5(a) presents the reconstructed path of the shadow of Pluto141414The
occulted star is Gaia DR3 $6772059623498733952$, of which the astrometric and
photometric parameters are obtained from VizieR. during the 17 July 2019
event, which has been studied by Arimatsu et al. (2020). We also organized an
observation campaign, involving stations detailed in Table 2. All the above-
mentioned stations are also presented in Figure 5(a). Unfortunately, all our
stations encountered weather problems.

(a)

(b)

(c)

(d)

Figure 5: Reanalysis of the 17 July 2019 event based on the TUHO observational data from Figure 1 of Arimatsu et al. (2020) (credit: Arimatsu et al, A&A, 638, L5, 2020, reproduced with permission © ESO). Panel a: reconstructed occultation map. Panel b: observed and simultaneously fitted light curves. Panel c: reconstructed stellar paths seen by TUHO. Panel d: the $\chi^{2}$ map, where $\chi^{2}$ denotes the goodness-of-fit value using these data. The best-fitting $\chi^{2}$ value per degree of freedom is $0.860$. Table 2: Stations that encountered weather problems on 17 July 2019 Station | Longitude (N) | Latitude (E) | Altitude (m) | Observers
---|---|---|---|---
PMOXL | $118\degr 57^{\prime}30\aas@@fstack{\prime\prime}0$ | $32\degr 07^{\prime}33\aas@@fstack{\prime\prime}0$ | $74$ | Jian Chen, Ansheng Zhu, Yue Chen
XYOS | $118\degr 27^{\prime}50\aas@@fstack{\prime\prime}0$ | $32\degr 44^{\prime}03\aas@@fstack{\prime\prime}7$ | $227$ | Wei Zhang
YAOS | $101\degr 10^{\prime}52\aas@@fstack{\prime\prime}0$ | $25\degr 31^{\prime}43\aas@@fstack{\prime\prime}0$ | $1943$ | Ye Yuan, Fan Li
GMG | $100\degr 01^{\prime}51\aas@@fstack{\prime\prime}6$ | $26\degr 42^{\prime}33\aas@@fstack{\prime\prime}1$ | $3185$ | Xiaoli Wang, Jianguo Wang
SZAO | $114\degr 33^{\prime}21\aas@@fstack{\prime\prime}0$ | $22\degr 28^{\prime}56\aas@@fstack{\prime\prime}0$ | $220$ | Jiahui Ye, Yijing Zhu, Delai Li, Lin Mei

###  A.2 The 5 September 2019 occultation campaign

The occultation map of the 5 September 2019 event is given by Figure 1. Table
3 lists the circumstances of stations with positive detections. Table 4 lists
the circumstances of stations that encountered weather problems.

Table 3: Stations with positive detections on 5 September 2019 Station | Longitude (N) | Telescope | Exposure (s) | Observers
---|---|---|---|---
| Latitude (E) | Camera | Circle (s) |
| Altitude (m) | Filter |  |
DWM | $114\degr 06^{\prime}44\aas@@fstack{\prime\prime}5$ | DOB 0.4 m | 0.2 | Ye Yuan, Jian Chen,
| $28\degr 25^{\prime}28\aas@@fstack{\prime\prime}0$ | QHY174GPS | 0.2 | Wei Tan, Ying Wang,
| $1388$ | clear |  | Hui Zhang, Shun Yao
HNU | $114\degr 30^{\prime}58\aas@@fstack{\prime\prime}0$ | RC 0.5 m | 9 | Qiang Zhang,
| $37\degr 59^{\prime}55\aas@@fstack{\prime\prime}0$ | CCD | 14 | Fan Li,
| $69$ | clear |  | Zhensen Fu
Table 4: Stations that encountered weather problems on 5 September 2019 Station | Longitude (N) | Latitude (E) | Altitude (m) | Observers
---|---|---|---|---
PMOXL | $118\degr 57^{\prime}30\aas@@fstack{\prime\prime}0$ | $32\degr 07^{\prime}33\aas@@fstack{\prime\prime}0$ | $74$ | Jian Chen, Ansheng Zhu, Yue Chen
NJXD | $118\degr 25^{\prime}17\aas@@fstack{\prime\prime}5$ | $32\degr 04^{\prime}40\aas@@fstack{\prime\prime}8$ | $74$ | Jun Xu
XYOS | $118\degr 27^{\prime}50\aas@@fstack{\prime\prime}0$ | $32\degr 44^{\prime}03\aas@@fstack{\prime\prime}7$ | $227$ | Wei Zhang
YAOS | $101\degr 10^{\prime}52\aas@@fstack{\prime\prime}0$ | $25\degr 31^{\prime}43\aas@@fstack{\prime\prime}0$ | $1943$ | Chen Zhang, Fan Li, Ye Yuan
GMG | $100\degr 01^{\prime}51\aas@@fstack{\prime\prime}6$ | $26\degr 42^{\prime}33\aas@@fstack{\prime\prime}1$ | $3185$ | Xiliang Zhang, Xiaoli Wang, Jianguo Wang
SZAO | $114\degr 33^{\prime}21\aas@@fstack{\prime\prime}0$ | $22\degr 28^{\prime}56\aas@@fstack{\prime\prime}0$ | $220$ | Jiahui Ye, Delai Li, Lin Mei
QHOS | $97\degr 33^{\prime}36\aas@@fstack{\prime\prime}0$ | $37\degr 22^{\prime}24\aas@@fstack{\prime\prime}0$ | $3200$ | Wei Zhang
QDAS | $119\degr 56^{\prime}38\aas@@fstack{\prime\prime}0$ | $36\degr 07^{\prime}20\aas@@fstack{\prime\prime}0$ | $107$ | Kun Zhou

##  Appendix B Data processing

The observational data were initially obtained in the FITS format, a commonly
used format for astronomical data. To process and analyze these data, we used
the Tangra occultation photometric tool developed by Pavlov (2020). This tool
offers various functionalities and algorithms specifically designed for
analyzing occultation events. When calibration images including bias, dark,
and flat-field frames are available, they are used for image correction. The
default measurement type and tracking method provided by Tangra were applied.
The reduction method was aperture photometry with median background
subtraction. The combined image of the occulted star and Pluto’s system is
circled as the target.

For each station, a selection of suitable guiding and nearby reference stars
(typically three) was made. All images of target and reference stars are
confirmed not overexposed. The signal flux counts, denoted as
$I_{\mathrm{signal}}$, and the background flux counts, denoted as
$I_{\mathrm{bkg}}$, were measured for the target and each reference star in
all FITS files. The timing information and measurement results were exported
to a CSV file specific to each station, from which the observed light curves
would be generated.

Further data processing and analysis were performed using our data reduction
code implemented in Python. Each flux measurement for the target or a
reference star was calculated as $I=I_{\mathrm{signal}}-I_{\mathrm{bkg}}$. The
measurement error, $\sigma_{I}$, was modeled using the classical form given
by:

$\sigma_{I}=\sqrt{\sigma_{\mathrm{bkg}}^{2}+I\cdot g_{\mathrm{eff}}^{-1}},$
(5)

where $\sigma_{\mathrm{bkg}}$ is the standard deviation of the background
($I_{\mathrm{bkg}}$); and $g_{\mathrm{eff}}$ is the effective gain that
represents the number of electrons per flux count.

The relative photometric result of the target is defined as

$F=\frac{I_{T}}{I_{RS}},$ (6)

where $I_{T}$ is the $I$ of target and $I_{RS}$ the sum of those of all
reference stars. This step helps eliminate low-frequency variations caused by
atmospheric and instrumental effects. The corresponding measurement error,
$\sigma_{F}$, was derived using the usual formula for propagation of errors
(Press et al. 2007).

If the effective gain $g_{\mathrm{eff}}$ was not precisely known,
$g_{\mathrm{eff}}^{-1}$ was adjusted with the non-negative constraint to fit
the derived $\sigma_{F}$ model ($\sigma_{F}(g_{\mathrm{eff}}^{-1})$ where
$g_{\mathrm{eff}}^{-1}\geq 0$) to the standard deviation of $F$ outside the
occultation part.

Finally, the normalized total observed flux, denoted as $f$, and its error,
$\sigma$, were derived by dividing $F$ and $\sigma_{F}$, respectively, by the
median of $F$ outside the occultation part. As a result, $f$ outside the
occultation part approximates unity.

##  Appendix C Synthetic light curve model for a stellar occultation by Pluto

For a given station, the light-curve model is represented as:

$\phi(t)=A\cdot\left(s\cdot\psi(t)+(1-s)\right),$ (7)

where $t$ is the recorded timing of the observation and $\phi$ is the received
flux, as a function of $t$; $A$ is the total flux of the star and Pluto’s
system when the star is not occulted; $s$ is the flux ratio of the occulted
star to $A$; and $\psi$ is the normalized (between zero and unity) flux of the
occulted star, which represents the total flux of the primary (sometimes
called the near-limb) and secondary (or far-limb, when available) images
produced by Pluto’s spherical planetary atmosphere (Sicardy 2022) as shown in
Figure 2.

To calculate $\psi$, the following equation is used:

$\psi(t)=\left|\frac{r_{+}}{z}\right|\left|\frac{dr_{+}}{dz}\right|+\left|\frac{r_{-}}{z}\right|\left|\frac{dr_{-}}{dz}\right|,\quad(z>0,r_{\pm}>R_{\mathrm{p}})$
(8)

where $z$ is the distance of the station to the shadow axis, that is, the line
passing through the center of Pluto and parallel to the light ray from the
star; $r_{+}$ and $r_{-}$ represent the closest approach distances of the
primary and secondary images, respectively to the shadow axis; and
$R_{\mathrm{p}}$ is the radius of Pluto’s body. The dividers, $|r_{\pm}/z|$,
and derivatives, $|dr_{\pm}/dz|$, account for the compression and stretching
effects produced by the spherical planetary atmosphere, respectively, caused
by limb curvature and differential refraction.

Given reference ephemerides and a star catalog, the distance $z$ is modeled in
the International Celestial Reference Frame (ICRF) by using the following
equation:

$z(t)=\left|\left(\bm{\rho}_{\mathrm{p}}(t+\Delta
t_{i}+\Delta\tau)+\Delta\rho\cdot\hat{\bm{n}}_{\mathrm{cag}}-\bm{\rho}_{i}(t+\Delta
t_{i})\right)\bm{\times}\hat{\bm{\ell}}_{\mathrm{s}}\right|,$ (9)

where $t+\Delta t_{i}$ represents the real observation epoch corrected for the
camera time recording offset $\Delta t_{i}$; $\hat{\bm{\ell}}_{\mathrm{s}}$ is
the unit vector in the direction of the geocentric astrometric position of the
occulted star; $\bm{\rho}_{i}$ is the geocentric geometric position of $i$;
$\bm{\rho}_{\mathrm{p}}$ is the geocentric astrometric position of Pluto
obtained at $t+\Delta t_{i}+\Delta\tau$, with $\Delta\tau$ accounting for the
ephemeris offset along the geocentric motion of Pluto;
$\hat{\bm{n}}_{\mathrm{cag}}$ is the unit vector in the northern direction
along the predicted geocentric closest approach to the shadow axis, multiplied
by $\Delta\rho$ to account for the ephemeris offset across the geocentric
motion of Pluto. The two ephemeris offset parameters, $\Delta\tau$ and
$\Delta\rho$, represent the corrections to the epoch $t_{\mathrm{cag}}$ and
distance $\rho_{\mathrm{cag}}$ of the geoncetric closest approach predicted
with the given reference ephemerides and star catalog (see Section 3.1).

Then, the values of $r_{\pm}$ can be numerically solved by using the
relationship between $z$ and $r_{\pm}$:

$z=\pm(r_{\pm}+D\cdot\omega(r_{\pm})),$ (10)

where $D$ represents the light travel distance, which can be approximated by
the geocentric distance of Pluto, and $\omega(r)$ is the total deviation angle
at $r$.

As detailed in Dias-Oliveira et al. (2015) and Sicardy (2022), $\omega(r)$ can
be obtained by a ray-tracing code that follows these steps:

  * •

set the temperature profile $T(r)$ of Pluto using Equation (4) from Dias-
Oliveira et al. (2015), with parameters obtained from Sicardy et al. (2016);

  * •

set the gas molecular refractivity $K$ corresponding to a stellar wavelength
of $\lambda_{\mathrm{\upmu m}}$;

  * •

set the atmospheric pressure $p_{0}$ at the reference radius $r_{0}$ of 1215
km either directly, or determined from a surface pressure $p_{\mathrm{surf}}$
(at $R_{\mathrm{p}}=1187$ km) with the ratio $p_{\mathrm{surf}}/p_{0}=1.837$
used by previous studies (Meza et al. 2019; Sicardy et al. 2021);

  * •

set the needed boundary condition for the gas molecular density,
$n_{0}=p_{0}/(k_{\mathrm{B}}\cdot T(r_{0}))$, where $k_{\mathrm{B}}$ is
Boltzmann’s constant;

  * •

integrate the first-order differential equation (Equations (5) and (6) of
Dias-Oliveira et al. (2015)) to obtain the gas molecular density profile
$n(r)$;

  * •

transform $n(r)$ to the gas refractivity $\nu(r)$ using $\nu(r)=K\cdot n(r)$
(Equation (7) from Dias-Oliveira et al. (2015));

  * •

calculate the total deviation based on the straight line approximation
(Equation (16) of Sicardy (2022), with the upper limit of $r$ set as 2600 km).

Finally, an interpolator of $\omega(r)$ for specific $p_{0}$ and $K$ values is
built in advance. This allows for efficient computation of $\omega(r)$ during
the analysis of the occultation observations.

##  Appendix D Results of the reanalyzed events necessary for comparisons

###  D.1 The 15 August 2018 event

Based on the IXON observational data from the online supplemental data
file151515“stac401_Supplement_File - zip file”, via the paper link
https://doi.org/10.1093/mnras/stac401. of Silva-Cabrera et al. (2022), we
reanalyzed the 15 August 2018 event161616The occulted star is Gaia DR3
$6772629170525258240$, of which the astrometric and photometric parameters are
obtained from VizieR. using our fitting procedure. As the raw data had been
normalized, with the flux contribution of the Pluto system carefully removed
by Silva-Cabrera et al. (2022), all $A_{i}$ and $s_{i}$ were fixed as unity in
our fitting procedure. Moreover, as the uncertainties used for weighting are
not given, we used the standard deviation of the data outside the occultation
part for each station. Results are presented in Figure 6. Our surface pressure
remeasurement is ${12.027}_{-0.08}^{+0.09}\leavevmode\nobreak\ \mathrm{\upmu
bar}$.

(a)

(b)

(c)

(d)

Figure 6: Reanalysis of the 15 August 2018 event based on the IXON
observational data from the online supplemental data file of Silva-Cabrera et
al. (2022). The raw data had been normalized, with the flux contribution of
the Pluto system carefully removed by Silva-Cabrera et al. (2022). The
uncertainties used for weighting are the standard deviation of the data
outside the occultation part for each station. Panel a: reconstructed
occultation map. Panel b: observed and simultaneously fitted light curves.
Panel c: reconstructed stellar paths seen by IXON. Panel d: the $\chi^{2}$
map, where $\chi^{2}$ denotes the goodness-of-fit value using these data. The
best-fitting $\chi^{2}$ value per degree of freedom is ${1.046}$.

###  D.2 The 17 July 2019 event

Based on the TUHO observational data and uncertainties from Figure 1 of
Arimatsu et al. (2020) (credit: Arimatsu et al, A&A, 638, L5, 2020, reproduced
with permission © ESO), we reanalyzed the 17 July 2019 event using our fitting
procedure. Results are presented in Figure 5. Our surface pressure
remeasurement is $9.421_{-0.75}^{+0.68}\leavevmode\nobreak\ \mathrm{\upmu
bar}$, similar to the result obtained by Arimatsu et al. (2020),
${9.56}_{-0.34}^{+0.52}\leavevmode\nobreak\ \mathrm{\upmu bar}$.