# The tropical polytope is the set of all weighted tropical Fermat-Weber
points

Shelby Cox and Mark Curiel

###### Abstract.

Let $v_{1},\ldots,v_{m}$ be points in $\mathbb{R}^{n}$, and let
$w_{1},\ldots,w_{m}$ be positive real weights. The weighted Fermat-Weber
points are those points $z$ which minimize $\sum w_{i}d(v_{i},z)$. Following
Comăneci and Joswig, we study the weighted Fermat-Weber points with respect to
an asymmetric tropical metric. In the unweighted case (when
$w_{1}=\cdots=w_{m}=1$), Comăneci and Joswig showed that the set of Fermat-
Weber points is the “central” cell of the tropical convex hull of
$v_{1},\ldots,v_{m}$. We show that for any fixed data points
$v_{1},\ldots,v_{m}$, as the weights $w_{i}$ vary, the set of all Fermat-Weber
points is the entire tropical convex hull of the $v_{i}$.

##  1\. Introduction

Phylogenetics is a field of computational biology concerned with converting
molecular data into trees that capture evolutionary relationships. While there
are several, often combinatorial, methods for computing these phylogenetic
trees, typically there is much discrepancy among the various resulting trees.
Thus, a key hurdle to overcome is to find a tree which best represents the
true evolutionary data. A common strategy is to gather information from many
trees and compile this information into a single tree. The resulting tree
obtained from such a strategy is called a consensus tree and the algorithm
used to compute a consensus tree is called a consensus method. Still there are
many consensus methods for computing a consensus tree. For a survey of
consensus methods see [1]. While many consensus methods are combinatorial, in
this paper our approach is through tropical geometry.

A phylogenetic tree is a metric tree with $N$ leaves. We embed the set of
phylogenetic trees on N leaves into $\mathbb{R}^{\binom{N}{2}}$, by sending a
tree $T$ to the list of $\binom{N}{2}$ distances between each of the pairs of
leaves of $T$. This space is known as the space of phylogenetic trees [8]. In
this paper, we only consider trees up to adding a constant to all pendant
(leaf-adjacent) edge lengths. In this setting, the space of phylogenetic trees
can be realized as a tropical linear space, the space of ultrametrics, in the
tropical projective torus $\mathbb{R}^{\binom{N}{2}}/\mathbb{R}\mathbbm{1}$
[8]. In particular, this tree space is tropically convex – for more background
see [8, §4.3].

A problem of computing a consensus tree can now be stated as the following
optimization problem: for a fixed set of points in phylogenetic tree space,
locate a point which minimizes the sum of distances to that point. This type
of problem is more broadly known as a Fermat-Weber problem. A Fermat-Weber
problem is a geometric problem seeking the median of a collection of data
points $v_{1},\ldots,v_{m}$ taken from a metric space $X$, specifically the
problem asks to find a point $x$ such that sum $\sum_{i=1}^{m}d(x,v_{i})$ is
minimized. Such a median $x$ is called a Fermat-Weber point. This optimization
problem arises in many different contexts such as linear programming via
transportation problems [2], economics via location theory [4], and of course
phylogenetics via tropical geometry [2, 7].

In the Euclidean case, there is a single Fermat-Weber point, and it is
contained in the convex hull of the data points. In [2], Comăneci and Joswig
show that with an asymmetrical tropical metric on phylogenetic tree space, the
Fermat-Weber points form a cell in the tropical convex hull of the data
points. Since phylogenetic tree space is tropically convex, this means that
each Fermat-Weber point can be interpreted as a tree itself.

Classically, the Fermat-Weber problem was first posed by Fermat before 1640 to
compute the Fermat-Weber point when $P$ is a triangle in Euclidean space. It
was solved geometrically by Evangelista Torricelli in 1645. We note that the
sum $\sum_{i=1}^{m}d(x,v_{i})$ assumes that the distances are weighted
equally. Motivated by unequal attracting forces on particles, a generalization
of the problem was introduced and solved by Thomas Simpson in 1750, later
popularized by Alfred Weber in 1909, by considering weighted distances. In
that spirit, this paper is concerned with generalizing the result by Comăneci
and Joswig, namely, we are interested in locating the specific cells for which
the Fermat-Weber points live in $P$ by minimizing the sum
$\sum_{i=1}^{m}w_{i}d(x,v_{i})$ for some choice of positive real weights
$w_{i}$. Our main theorem states that the Fermat-Weber set is a cell of the
tropical convex hull $P$ and that, by choosing weights appropriately, it can
be any cell of $P$.

######  Theorem 1.1.

Given data points $v_{1},\ldots,v_{m}$, the collection of asymmetric tropical
weighted Fermat-Weber points over all possible positive real weights $w_{i}$
is $P=\mathrm{tconv}\\{v_{1},\ldots,v_{m}\\}$.

This paper is organized as follows: in Section 2 we provide background in
tropical geometry and polyhedral geometry necessary to understand our
approach. Of particular interest we recall the Cayley trick which gives a
correspondence between mixed subdivisions of the Minkowski sum of polytopes
and subdivisions of the corresponding Cayley polytope. Further, at the end of
Section 2 we formulate the Fermat-Weber problem in the language of tropical
convexity. In Section 3, we prove our main result Theorem 1.1 as a corollary
to Theorem 3.2.

##  2\. Background

In this section we set up the framework of the Fermat-Weber problem and
introduce the algebraic and geometric tools that we use to study it.

###  2.1. Classical Polytopes

In this subsection we recall objects from convex geometry that play keys roles
in computing Fermat-Weber sets. A polyhedron $P\subset\mathbb{R}^{n}$ is an
intersection of finitely many half spaces. We call $P$ a polytope if this
intersection is bounded. In this case, we write $P=\mbox{conv}(A)$ to mean $P$
is the convex hull of some finite set $A\subset\mathbb{R}^{n}$. For instance
the standard simplex is the convex hull of the $n$ standard basis vectors in
$\mathbb{R}^{n}$ and we denote it by $\Delta^{n-1}$. A face of a polytope $P$
is the collection of points in $P$ that minimizes the dot product with a fixed
vector $u$, specifically

(1) $\mbox{face}_{\mathbf{u}}(P):=\\{\mathbf{x}\in
P\mid\mathbf{u}\cdot\mathbf{x}\leq\mathbf{u}\cdot\mathbf{y},\forall\mathbf{y}\in
P\\}.$

A polyhedral complex is a collection of polyhedra $S=\\{C_{i}\\}$ with the
following properties: (1) $S$ is closed under taking faces, and (2) $C_{i}\cap
C_{j}$ is a face of both $C_{i}$ and $C_{j}$, or is empty. The normal cone to
$F$ in $P$, denoted $\sigma_{F}$, is the set of all vectors
$u\in\left(\mathbb{R}^{m}\right)^{\ast}$ such that
$\mbox{face}_{\mathbf{u}}(P)\supseteq F$. Alternatively, $\sigma_{F}$ is the
closure of
$\\{u\in\left(\mathbb{R}^{m}\right)^{\ast}\mid\mbox{face}_{\mathbf{u}}(P)=F\\}$.
The normal fan of a polytope $P$ is the collection of cones
$\\{\sigma_{F}\\}_{F\text{ a face of }P}$.

For the remainder of this section we are concerned with subdivisions of
polytopes. Consider $m$ polytopes $P_{1},\ldots,P_{m}\subset\mathbb{R}^{n}$
where $P_{i}=\mathrm{conv}(A_{i})$ for some finite sets
$A_{i}\subset\mathbb{R}^{n}$. A subdivision of a polytope $P$ is a polyhedral
complex $S=\\{C_{i}\\}$ such that $P=\bigcup_{i}C_{i}$.

######  Definition 2.1.

Given a polytope
$P=\mbox{conv}\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{k}\\}\subset\mathbb{R}^{n}$,
and weights $w_{i}\in\mathbb{R}$ on $\mathbf{v}_{i}$, the lift of $P$ with
respect to $\mathbf{w}$ is

(2) $\tilde{P}:=\mbox{conv}\\{(\mathbf{v}_{i},w_{i})\mid
i=1,\ldots,k\\}\subset\mathbb{R}^{n+1}$

A face $F=\mbox{face}_{\mathbf{u}}(\tilde{P})$ of $\tilde{P}$ is an upper face
if $u_{n+1}<0$. The regular subdivision of $P$ with respect to $\mathbf{w}$ is
the projection of the upper faces of $\tilde{P}$ onto $P$ (forgetting the last
coordinate). An example is given in Example 2.24.

######  Definition 2.2.

The normal complex of a regular subdivision is the projection of the cones of
the normal fan of $\tilde{P}$ that are normal to upper faces of $\tilde{P}$.

Notation: A subdivision of $P$ will be denoted $\underline{P}$. Regular
subdivisions induced by a piecewise-linear convex function $\lambda$ will be
denoted $\underline{P}_{\lambda}$.

######  Definition 2.3.

Let $P_{1},\ldots,P_{m}$ be polytopes in $\mathbb{R}^{n}$. The Cayley
polytope, $\mbox{Cayley}(P_{1},\ldots,P_{m})$, is the convex hull of
$\bigcup_{i=1}^{m}e_{i}\times P_{i}$ in $\mathbb{R}^{m}\times\mathbb{R}^{n}$.
In the special case where $P_{1}=P_{2}=\cdots=P_{m}$, the Cayley polytope is
$\Delta^{r-1}\times P$. For an example see Figure 5,

######  Definition 2.4.

The Minkowski sum of $P_{1},\ldots,P_{m}$ is the set

$P=P_{1}+\cdots+P_{m}=\\{\mathbf{p}_{1}+\cdots\mathbf{p}_{m}\mid\mathbf{p}_{i}\in
P_{i}\\}\subset\mathbb{R}^{n}.$

The Minkowski sum of polytopes is indeed a polytope since its vertices are
necessarily sums of vertices of the summands, hence the Minkowski sum is a
convex hull of a finite set.

######  Example 2.5 (Minkowski Sum).

Let $A_{1}=\\{a,b,c,d\\}$ and $A_{2}=\\{e,f,g\\}$ be the vertex sets in
$\mathbb{R}^{2}$ of $P_{1}$ and $P_{2}$ respectively. The minkowski sum
$P_{1}+P_{2}$ is the convex hull of $\\{a+e,b+f,c+f,c+g,d+g\\}$ and is shown
in Figure 1.

Figure 1. A square (left), a triangle (middle), and their Minkowski sum
(right).

######  Definition 2.6.

A cell of the minkowski sum $P=P_{1}+\ldots+P_{m}$ is a tuple
$C=(C_{1},\ldots,C_{m})$ where $C_{i}\subseteq A_{i}$ for all $i$.

######  Remark 2.7.

Each cell $(C_{1},C_{2},\ldots,C_{m})$ gives a polytope $\sum
C_{i}\subseteq\sum P_{i}$, and we will often abuse notation by identifying
$(C_{1},C_{2},\ldots,C_{m})$ with this sum. If $C$ and $C^{\prime}$ are two
such cells, then $C\cap C^{\prime}$ refers to the intersection
$\left(\sum_{i}\mathrm{conv}(C_{i})\right)\cap\left(\sum_{i}\mathrm{conv}(C^{\prime}_{i})\right)$.
Additionally, a cell of $P$ may be the Minkowski sum of two (or more)
different ordered $m$-tuples, and we would like to consider these as different
cells.

The Cayley trick is a correspondence between mixed subdivisions of the
Minkowski sum $P_{1}+\cdots+P_{m}$ and subdivisions of
$\mbox{Cayley}(P_{1},\ldots,P_{m})$, illustrated in Figure 5. The top right
polytope in Figure 5 is
$\mbox{Cayley}(\Delta^{2},\Delta^{2})\cong\Delta^{1}\times\Delta^{2}$.
Explicitly, a subdivision of the Cayley polytope gives rise to a subdivision
of the Minkowski sum after forgetting the first $m$ coordinates. For more
details, see [11, Section 5] for coherent/regular subdivisions, and [5,
Theorem 3.1] for all subdivisions. Although often stated as a theorem, we will
use the Cayley trick to define mixed subdivisions.

######  Definition 2.8.

A mixed subdivision of $P=(P_{1},\ldots,P_{m})$ is a collection of cells
$C^{1},\ldots,C^{k}$ so that
$\\{\mbox{Cayley}(\mbox{conv}(C^{j}_{1}),\ldots,\mbox{conv}(C^{j}_{m}))\mid
j=1,\ldots,k\\}$ is a subdivision of $\mbox{Cayley}(P_{1},\ldots,P_{m})$.

######  Example 2.9.

A subdivision of the polytope $P_{1}+P_{2}$ of Example 2.5 consists of a
collection of cells $\\{C^{1},C^{2},C^{3},C^{4}\\}$ with

$\begin{array}[]{lll}C^{1}=\mathrm{conv}\\{a,b,c,d\\}+\mathrm{conv}\\{e\\}&&C^{3}=\mathrm{conv}\\{b,c\\}+\mathrm{conv}\\{e,f\\}\\\
C^{2}=\mathrm{conv}\\{c,d\\}+\mathrm{conv}\\{e,g\\}&&C^{4}=\mathrm{conv}\\{c\\}+\mathrm{conv}\\{e,f,g\\}\\\
\end{array}$

It is mixed since the cells $\mbox{Cayley}(\\{a,b,c,d\\},\\{e\\})$,
$\mbox{Cayley}(\\{b,c\\},\\{e,f\\})$, $\mbox{Cayley}(\\{c,d\\},\\{e,g\\})$,
$\mbox{Cayley}(\\{c\\},\\{e,f,g\\})$ give a subdivision of
$\mbox{Cayley}(P_{1},P_{2})$. Another subdivision of $P_{1}+P_{2}$ can be
achieved with the cells

$\displaystyle C^{1}$
$\displaystyle=\mathrm{conv}\\{a,b,c,d\\}+\mathrm{conv}\\{e\\}$ $\displaystyle
C^{2}$ $\displaystyle=\mathrm{conv}\\{c,d\\}+\mathrm{conv}\\{e,g\\}$
$\displaystyle C^{3}$
$\displaystyle=\mathrm{conv}\\{c\\}+\mathrm{conv}\\{e,f,g\\}$ $\displaystyle
C^{4}$ $\displaystyle=\mathrm{conv}\\{a,b,c\\}+\mathrm{conv}\\{f\\}$
$\displaystyle C^{5}$
$\displaystyle=\mathrm{conv}\\{a,c,d\\}+\mathrm{conv}\\{f\\}$

However, it is not mixed since for example the cells
$\mbox{Cayley}(\\{a,b,c,d\\},\\{e\\})$ and
$\mbox{Cayley}(\\{a,b,c\\},\\{f\\})$ intersect on their interior.

Figure 2. A mixed subdivision of $P_{1}+P_{2}$ (left) and a subdivision of
$P_{1}+P_{2}$ that is not mixed (right).

We will refer to Definition 2.8 above as the combinatorial Cayley trick. In
addition to the combinatorial correspondence above, there is also an explicit
geometric correspondence between a subdivision of the Cayley polytope and a
mixed subdivision, sometimes called the geometric Cayley trick.

######  Theorem 2.10 ([5, Theorem 3.1]).

Let $\underline{C}$ be a subdivision of $\mbox{Cayley}(P_{1},\ldots,P_{m})$.
Then the corresponding mixed subdivision of $P=\sum_{i=1}^{m}P_{i}$ is
$\underline{P}=n\cdot\underline{C}\cap\\{\frac{1}{m}\mathbbm{1}_{m}\\}\times\mathbb{R}^{n}$.

Proposition 2.11 is essentially stated in [10, §1.3]. It allows us to think
about mixed subdivisions of a weighted Minkowski sum,
$P^{\mathbf{w}}=\sum_{i=1}^{m}w_{i}P_{i}$, in terms of the Cayley polytope
without weights, $\mbox{Cayley}(P_{1},\ldots,P_{m})$, by slicing at
$\\{\mathbf{w}\\}\times\mathbb{R}^{n}$. We provide a proof that explicitly
states the map that induces the bijection, which we will use later in the
paper.

######  Proposition 2.11.

Let $C=\mbox{Cayley}(P_{1},\ldots,P_{m})$, and let $D=\mbox{cone}(C)$ denote
the cone over it. Let $C^{\mathbf{w}}$, $D^{\mathbf{w}}$, and $P^{\mathbf{w}}$
denote the corresponding weighted versions, for $\mathbf{w}\in\mathbb{R}^{m}$
with $|\mathbf{w}|=1$. Let $\lambda$ be a piecewise-linear convex function on
$\mbox{Cayley}(w_{1}P_{1},\ldots,w_{m}P_{m})$, and let
$g(\mathbf{x},\mathbf{y})=(w_{1}x_{1},\ldots,w_{m}x_{m},\mathbf{y})$. If
$\lambda^{\prime}=\lambda(g^{-1}(x,y))$, then the following diagram commutes.

${\underline{D}_{\lambda}^{\mathbf{w}}}$${\underline{D}_{\lambda^{\prime}}}$${\underline{C}_{\lambda}^{\mathbf{w}}}$${\underline{C}_{\lambda^{\prime}}}$${\underline{P}_{\lambda}^{\mathbf{w}}}$$\scriptstyle{g}$$\scriptstyle{\cap\\{\sum_{i=1}^{m}x_{i}=1\\}}$$\scriptstyle{n\cdot\left(\left(\\{\frac{1}{m}\\}\times\mathbb{R}^{n}\right)\cap\bullet\right)}$$\scriptstyle{\left(\bullet\cap\left(\\{\mathbf{w}\\}\times\mathbb{R}^{n}\right)\right)}$

###### Proof.

The subdivision $\underline{C}_{\lambda}^{\mathbf{w}}$ induces a subdivsion
$\underline{D}_{\lambda}^{\mathbf{w}}$. The function $g$ is an invertible
linear function. In particular, this means $g$ preserves convexity, dimension,
and the containment relations within a polyhedral subdivision. Thus,
$\underline{D}_{\lambda}^{\mathbf{w}}$ induces a subdivision
$\underline{D}_{\lambda^{\prime}}$ via
$\lambda^{\prime}(z)=\lambda(g^{-1}(z))$. This in turn induces a subdivision
$\underline{C}_{\lambda^{\prime}}$ with the weightings
$\lambda^{\prime}(e_{i},v_{ij})=\lambda(e_{i},w_{i}v_{ij})$, which is
combinatorially equivalent to $\underline{C}_{\lambda^{\prime}}^{\mathbf{w}}$.
Moreover,
$n\cdot\underline{C}_{\lambda^{\prime}}^{\mathbf{w}}\cap\left(\\{\frac{1}{m}\\}\times\mathbb{R}^{n}\right)$,
and
$\underline{C}_{\lambda^{\prime}}\cap\left(\\{\mathbf{w}\\}\times\mathbb{R}^{n}\right)$
give the same subdivision of $\underline{P}_{\lambda}^{\mathbf{w}}$. ∎

###  2.2. Tropical Polynomials and Regular Subdivisions

####  2.2.1. Tropical Arithmetic

In the tropical max-plus semi-ring
$(\mathbb{R}\cup\\{-\infty\\},\oplus,\odot)$, tropical addition $\oplus$ and
tropical multiplication $\odot$ are defined by

$a\oplus b=\max\\{a,b\\},\quad a\odot b=a+b\quad\text{where}\
a,b\in\mathbb{R}.$

The multiplicative identity is 0, and the additive identity is $-\infty$.
These operations can be extended component-wise to the tropical projective
torus $\mathbbm{TP}^{n-1}\cong\mathbb{R}^{n}/\mathbb{R}\mathbbm{1}$. For
vectors $\mathbf{v},\mathbf{u}\in\mathbbm{TP}^{n-1}$, the notation
$\mathbf{v}\oplus\mathbf{u}$ and $\mathbf{v}\odot\mathbf{u}$ denotes
component-wise $\max$ and addition, respectively. Tropical scalar
multiplication of a vector amounts to adding a (classical) multiple of the all
ones vector $\mathbbm{1}_{n}$, namely
$\lambda\odot\mathbf{v}=\lambda\mathbbm{1}+\mathbf{v}=(\lambda+v_{1},\ldots,\lambda+v_{n})$
for any $\lambda\in\mathbb{R}$.

######  Example 2.12.

If $\mathbf{v}_{1}=(1,2,-3)$ and $\mathbf{v}_{2}=(-5,3,2)$ are points in
$\mathbbm{TP}^{2}$, then

$\displaystyle\mathbf{v}_{1}\oplus\mathbf{v}_{2}$
$\displaystyle=(\max\\{1,-5\\},\max\\{2,3\\},\max\\{-3,2\\})$
$\displaystyle=(1,3,2)=(-1,1,0)$
$\displaystyle\mathbf{v}_{1}\odot\mathbf{v}_{2}$
$\displaystyle=(1-5,2+3,-3+2)$ $\displaystyle=(-4,5,-1)$ $\displaystyle 3\odot
v_{2}$ $\displaystyle=(3-5,3+3,3+2)$ $\displaystyle=(-2,6,5)$

Later it will be convenient to fix a specific representative of
$v\in\mathbbm{TP}^{n-1}$, namely the one where the sum of the coordinates is
zero. We denote by $H_{0}$ the hyperplane where these points are located:

$H_{0}:=\left\\{\mathbf{z}\in\mathbb{R}^{n}\left|\phantom{.}\mathbf{z}\cdot\mathbbm{1}=\sum_{i}z_{i}=0\right.\right\\}\subset\mathbb{R}^{n}.$

Each point in $\mathbbm{TP}^{n-1}$ has a unique representative in $H_{0}$.

When we draw pictures in $\mathbbm{TP}^{n-1}$, we will tropically scale points
to have last coordinate zero, then project away the last coordinate and draw
the point in $\mathbb{R}^{n-1}$. For example, $(3,1,2)\equiv(1,-1,0)$ will be
drawn in the plane at the location $(1,-1)$.

####  2.2.2. Tropical Polynomials

Let $\mathbf{x},\mathbf{a}\in\mathbb{R}^{n}$, and define the tropical
monomial: $\mathbf{x}^{\mathbf{a}}:=\sum_{i=1}^{n}a_{i}x_{i}$. Note that when
the entries of $\mathbf{a}$ are non-negative integers:

$x_{1}^{a_{1}}\odot\cdots\odot x_{n}^{a_{n}}=\underbrace{x_{1}\odot\cdots\odot
x_{1}}_{a_{1}\ \text{times}}\odot\cdots\odot\underbrace{x_{n}\odot\cdots\odot
x_{n}}_{a_{n}\ \text{times}},$

which explains the notation. Note that for $\mathbf{a}\in\mathbb{R}^{n}$,
$\mathbf{x}^{\mathbf{a}}$ is still a well-defined tropical function, but not a
tropical polynomial.

######  Definition 2.13.

Let $\mathbf{x}\in\mathbbm{TP}^{n-1}$. A tropical signomial in $\mathbf{x}$ is
a finite linear combination of tropical monomials, i.e.

$f(\mathbf{x})=\bigoplus_{\mathbf{a}\in
A}\lambda_{\mathbf{a}}\odot\mathbf{x}^{\mathbf{a}}$

where $A\subset\mathbb{R}^{n}_{\geq 0}$ is finite and
$\lambda_{\mathbf{a}}\in\mathbb{R}$ for all $\mathbf{a}\in A$. If
$A\subset\mathbb{Z}^{n}_{\geq 0}$, then $f(\mathbf{x})$ is a tropical
polynomial.

######  Example 2.14.

Let $f(x)=1\oplus 3\odot x\oplus-1\odot x^{\sqrt{2}}$. In terms of classical
arithmetic operations:

$f(x)=\max\\{1,3+x,-1+\sqrt{2}x\\}.$

The graph of $f(x)$ is depicted in fig. 3; it has three linear pieces:

$f(x)=\begin{cases}1&\text{if }x\leq-2\\\ 3+x&\text{if }-2\leq
x\leq\frac{4}{\sqrt{2}-1}\\\ -1+\sqrt{2}x&\text{if }\frac{4}{\sqrt{2}-1}\leq
x.\end{cases}$

Figure 3. The graph of $f(x)=1\oplus 3\odot x\oplus-1\odot x^{\sqrt{2}}$. The
connection to the Newton polytope of $f$ is explained in section 2.2.4.

####  2.2.3. Optimization

A tropical max-plus signomial is a piecewise-linear, continuous, convex
function on $\mathbb{R}^{n}$. For a convex function, any local minimum is a
global minimum. This minimum can be identified by locating the tangent plane
with zero slope, which we formalize using subgradients.

######  Definition 2.15 ([9, §3.1.5]).

Given a convex function $f:\mathbb{R}^{n}\to\mathbb{R}$, the subdifferential
of $f$ at $x$ is:

(3) $\partial_{f}(x):=\\{u\in\mathbb{R}^{n}\mid\forall
z\in\mbox{dom}(f),f(z)\geq f(x)+u^{\top}\cdot(z-x)\\}.$

A subgradient of $f$ at $x$ is any element of $\partial_{f}(x)$.

The subdifferential of any function is a closed convex set. If $f$ is convex
and differentiable at $x$, then the subdifferential of $f$ at $x$ is a
singleton. In particular, if $f$ is linear then the subdifferential contains
only the slope of $f$. And if $f$ is piecewise-linear, then $\partial_{f}(x)$
is constant on the linear pieces of $f$.

######  Lemma 2.16 ([9, Theorem 3.1.15]).

For any function $f$, the subdifferential at $x$ contains
$\overline{\mathbf{0}}$ if and only if $x$ is a global minimizer for $f$.

###### Proof.

By definition, $\overline{\mathbf{0}}\in\partial_{f}(x)$ if and only if

$\forall z\in\mbox{dom}(f),f(z)\geq
f(x)+\overline{\mathbf{0}}^{\top}\cdot(z-x)\iff\forall
z\in\mbox{dom}(f),f(z)\geq f(x),$

which is if and only if $x$ is a global minimizer for $f$. ∎

######  Example 2.17.

Let $f(x)=1\oplus 3\odot x\oplus-1\odot x^{\sqrt{2}}$. The subdifferential of
$f(x)$ on each linear piece is the slope of that piece. The subdifferential at
$x=-2$ is $\partial_{f}(-2)=[0,1]$, and the subdifferential at
$x=4(1+\sqrt{2})$ is $[1,\sqrt{2}]$. Note that $0$ is in the subdifferential
of the constant (left-most) linear piece, and this is where the global minimum
of $f(x)$ is achieved. See fig. 3.

####  2.2.4. Tropical Hypersurfaces

The results we stated for subgradients hold for any convex function. In this
section we recall further results for tropical signomials. For a tropical
polynomial $f:\mathbb{R}^{n}\to\mathbb{R}$, subdifferentials of linear pieces
of $f$ are encoded by a subdivision of $\mbox{Newt}(f)$. It is this connection
that will allow us to convert the problem of optimizing $f$ into a polyhedral
geometry problem.

######  Definition 2.18.

The tropical vanishing set or tropical hypersurface of a tropical signomial
$f=\bigoplus_{i}c_{i}\odot x^{\alpha_{i}}$, denoted $\mbox{tropV}(f)$, is the
set of $x\in\mathbb{R}^{n}$ for which the $\max$ in $f(x)$ is achieved at
least twice.

(4) $\mbox{tropV}(f)=\\{x\in\mathbb{R}^{n}\mid\max\text{ in }f(x)\text{ is
achieved at least twice}\\}.$

######  Lemma 2.19.

The tropical vanishing set of a product of real positive powers of
polynomials, $f_{i}^{w_{i}}$, is (as a set) the union of tropical vanishing
sets of the $f_{i}$. That is,

$\mbox{tropV}\left(\bigodot_{i=1}^{m}f_{i}^{w_{i}}\right)=\bigcup_{i=1}^{m}\mbox{tropV}(f_{i}),\text{
for }w_{i}>0.$

###### Proof.

Let $f_{i}=\bigoplus_{\alpha\in A_{i}}c_{\alpha}\odot\mathbf{x}^{\alpha}$,
$A_{i}\subset\mathbb{R}^{n}_{>0}$. If the maximum in $f_{j}$ is achieved twice
by $\mathbf{x}$, then
$f_{j}^{w_{j}}(\mathbf{x})=w_{j}(c_{\alpha_{1}}+\alpha_{1}\cdot\mathbf{x})=w_{j}(c_{\alpha_{2}}+\alpha_{2}\cdot\mathbf{x})$
for some $\alpha_{1}\neq\alpha_{2}\in A_{j}$. It follows that the maximum in
$\bigodot_{i=1}^{m}f_{i}^{w_{i}}$ is also achieved at least twice:

$\bigodot_{i=1}^{m}f_{i}^{w_{i}}=w_{j}(c_{\alpha_{1}}+\alpha_{1}\cdot\mathbf{x})+\sum_{i\neq
j}f_{i}^{w_{i}}(\mathbf{x})=w_{j}(c_{\alpha_{2}}+\alpha_{2}\cdot\mathbf{x})+\sum_{i\neq
j}f_{i}^{w_{i}}(\mathbf{x}).$

On the other hand, if the maximum is achieved twice in
$\bigodot_{i=1}^{m}f_{i}^{w_{i}}$ at $\mathbf{x}$, then we must be able to
write the maximum as two distinct sums:
$\sum_{i=1}^{m}w_{i}(c_{\alpha_{i}^{1}}+\alpha_{i}^{1}\cdot\mathbf{x})=\sum_{i=1}^{m}w_{i}(c_{\alpha_{i}^{2}}+\alpha_{i}^{2}\cdot\mathbf{x})$,
where $\alpha_{i}^{1},\alpha_{i}^{2}\in A_{i}$. It follows that for some $j$,
$\alpha_{j}^{1}\neq\alpha_{j}^{2}$, so the maximum in $f_{j}^{w_{j}}$ is
achieved at least twice. ∎

A tropical hypersurface is a polyhedral complex, and it can be understood
combinatorially in terms of the Newton polytope of $f$, defined below.

######  Definition 2.20 (Newton polytope).

Suppose $f(\mathbf{x})=\bigoplus_{\mathbf{a}\in
A}c_{\mathbf{a}}\odot\mathbf{x}^{\mathbf{a}}$ is a multivariate polynomial for
some finite $A\subset\mathbb{R}^{n}$ and $c_{\mathbf{a}}\in\mathbb{R}$. The
support of $f(\mathbf{x})$ is the set, denoted $\mathrm{supp}(f)$, containing
all $\mathbf{a}\in A$ such that $c_{\mathbf{a}}\neq-\infty$. The Newton
polytope $\mathrm{Newt}(f)$ of the polynomial $f(\mathbf{x})$ is the convex
hull of its support, i.e. $\mathrm{Newt}(f)=\mathrm{conv}(\mathrm{supp}(f))$.
If $f,g$ are polynomials, then
$\mathrm{Newt}(fg)=\mathrm{Newt}(f)+\mathrm{Newt}(g)$.

######  Proposition 2.21 ([8, Proposition 3.1.6]).

Given a tropical signomial $f=\sum_{i}c_{i}x^{\alpha_{i}}$, let
$\underline{N}_{f}$ denote the regular subdivision of $\mbox{Newt}(f)$ induced
by the weighting $w(a_{i})=c_{i}$. Then $\mbox{tropV}(f)$ is the codimension-1
skeleton of the normal complex of $\underline{N}_{f}$.

######  Remark 2.22.

First, although Proposition 3.1.6 in [8] is originally stated for tropical
polynomials, the arguments clearly hold for tropical signomials as well.
Moreover, the proof in [8] shows that the subdifferential of a linear piece of
$f$ consists of the points in the corresponding cell of $\underline{N}$.

######  Definition 2.23.

Given a tropical polynomial $f$, let $\underline{N}$ be the regular
subdivision of $\mathrm{Newt}(f)$ induced by the coefficients of $f$. The
normal complex of $f$ is normal complex of $\underline{N}$; it is a
subdivision of $\mathbb{R}^{n}$.

######  Example 2.24.

The following is an example of Proposition 2.21. Let

(5) $f=x^{2}\oplus 4\odot xy\oplus 3\odot xz\oplus y^{2}\oplus 4\odot yz\oplus
3\odot z^{2}.$

Figure 4 depicts $\mbox{tropV}(f)$, the subdivision of the Newton polytope
dual to it, and the lift of the Newton polytope that induces that subdivision.

Figure 4. Left: A lift of $2\Delta^{2}=\mbox{Newt}(f)$ with weights given by
the coefficients of $f$, overlaid with $\mbox{tropV}(f)$ in black; right:
$\mbox{tropV}(f)$.

######  Lemma 2.25.

The minimum of a max-plus tropical polynomial is achieved on the cell dual to
the cell of the Newton polytope containing $\overline{\mathbf{0}}$.

###### Proof.

Let $f=\sum_{i}c_{i}x^{\alpha_{i}}$ be a max-plus tropical polynomial (so $f$
is a piecewise-linear convex function). Let $\underline{N}$ be the regular
subdivision of $\mbox{Newt}(f)$ induced by the weighting
$w(\alpha_{i})=c_{i}$. Let $L$ be a linear piece of $f$, and let $N_{L}$ be
the cell dual to it in $\underline{N}$. If $\overline{\mathbf{0}}\in N_{L}$,
then by Proposition 2.21 $\overline{\mathbf{0}}$ is in the subdifferential of
$f$ at any point in $L$. It then follows from Lemma 2.16 that the minimum of
$f$ is achieved on $L$. ∎

###  2.3. Fermat-Weber Problems

A Fermat-Weber problem is a geometric problem seeking the median of a
collection of data points $V=\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{m}\\}\subset
X$, where $X$ is a metric space with distance $d(\mathbf{x},\mathbf{y})$. We
are particularly interested in the Fermat-Weber points for a collection of
data points in $\mathbbm{TP}^{n-1}$, where the points could represent
phylogenetic trees. The goal of this section is to introduce the Fermat-Weber
problem, and reframe a tropical version as a problem on Newton polytopes.

In general, the median of a collection of points is not unique and hence we
seek the set of all such medians, called the called Fermat-Weber set. The
medians belonging to the Fermat-Weber set are called Fermat-Weber points.
Formally, the Fermat-Weber points are the points $\mathbf{x}\in X$ minimizing
the sum in (6).

######  Definition 2.26.

The Fermat-Weber points on the data
$V=\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{m}\\}\subset X$ are the points
$\mathbf{x}\in X$ minimizing the following sum

(6) $\mbox{FW}(V):=\frac{1}{m}\sum_{i=1}^{m}d(\mathbf{x},\mathbf{v}_{i}).$

In this paper, we are interested in a variant of the Fermat-Weber problem,
called the weighted Fermat-Weber problem. This new problem seeks the points
$\mathbf{x}$ minimizing the sum in Equation 7, where the weights $w_{i}$ are
positive real numbers.

(7)
$\mbox{FW}(V,\mathbf{w}):=\frac{1}{m}\sum_{i=1}^{m}w_{i}d(\mathbf{x},\mathbf{v}_{i}).$

We will use the asymmetric tropical distance first defined by Comăneci and
Joswig in [2].

######  Definition 2.27.

The asymmetric tropical distance, $d_{\Delta}(\mathbf{x},\mathbf{y})$ is:

(8)
$d_{\Delta}(\mathbf{x},\mathbf{y}):=n\max_{i\in[n]}(x_{i}-y_{i})+\sum_{i\in[n]}(y_{i}-x_{i}).$

When the points $\mathbf{x},\mathbf{y}\in\mathbb{R}^{n}/\mathbb{R}\mathbbm{1}$
are given by their unique representative in $H_{0}$ (the subspace where the
coordinates sum to zero), the metric $d_{\Delta}(\mathbf{x},\mathbf{y})$ can
be simplified to the following

(9)
$d_{\Delta}(\mathbf{x},\mathbf{y}):=n\max_{i\in[n]}(x_{i}-y_{i}),\phantom{.}\mathbf{x},\mathbf{y}\in
H_{0}.$

Note that $d_{\Delta}(x,y)$ is invariant under independent scalar
multiplication of the input vectors, so it is well-defined on
$\mathbbm{TP}^{n-1}$. From now on, we will assume that all points in
$\mathbbm{TP}^{n-1}$ are given by their representative in $H_{0}$. With this
assumption, the distance to a point $v_{i}$ can be reinterpreted as a power of
a tropical linear equations, and the sum in equation 7 can be realized as a
tropical product of tropical linear functions (possibly with real exponents).

The distance to a point $v_{i}$, denote by $f_{v_{i}}(\mathbf{x})$ is

(10)
$f_{v_{i}}^{n}(\mathbf{x}):=d_{\Delta}(x,v_{i})=n\max_{j}(x_{j}-v_{ij})=\left(\bigoplus_{j=1}^{n}-v_{ij}\odot
x_{j}\right)^{n}.$

It follows that the sum in (eq. 7) for $d=d_{\Delta}$ is

(11)
$\frac{1}{m}\sum_{i=1}^{m}w_{i}d_{\Delta}(\mathbf{x},\mathbf{v}_{i})=\frac{1}{m}\sum_{i=1}^{m}w_{i}f_{v_{i}}^{nw_{i}}(\mathbf{x})=\frac{n}{m}\bigodot_{i=1}^{m}f_{v_{i}}^{w_{i}}.$

######  Definition 2.28.

We define the tropical signomial associated to data $V$ with weights
$\mathbf{w}$, $f_{V,\mathbf{w}}$, to be the following tropical function:

$f_{V,\mathbf{w}}(\mathbf{x}):=\bigodot_{i=1}^{m}f_{v_{i}}^{w_{i}}=\bigodot_{i=1}^{m}\left(\bigoplus_{j=1}^{n}-v_{ij}\odot
x_{j}\right)^{w_{i}}.$

The tropical hypersurface $\mbox{tropV}(f_{v_{i}})$ is a tropical hyperplane
centered at $v_{i}$; it is the codimension-1 skeleton of the normal fan of the
standard simplex $\Delta^{n-1}$. By lemma 2.19, the hypersurface
$\mbox{tropV}(f_{V,\mathbf{w}})$ is the union of tropical hyperplanes centered
at the data points $v_{i}$. The Newton polytope of
$f_{V,\mathbf{w}}=\bigodot_{i=1}^{m}f_{v_{i}}^{w_{i}}$ is
$\sum_{i=1}^{m}w_{i}\cdot\Delta^{n-1}$.

######  Example 2.29.

The polynomial $f_{V,\mathbf{w}}(\mathbf{x})$ with
$\mathbf{x}\in\mathbbm{TP}^{2}$ and $\mathbf{w}\in\mathbb{R}^{2}$ has nine
terms for generic $V$ and $\mathbf{w}$.

$\displaystyle f_{V,\mathbf{w}}$ $\displaystyle=\left(-v_{11}\odot
x_{1}\oplus-v_{12}\odot x_{2}\oplus-v_{13}\odot
x_{3}\right)^{w_{1}}\odot\left(-v_{21}\odot x_{1}\oplus-v_{22}\odot
x_{2}\oplus-v_{23}\odot x_{3}\right)^{w_{2}}$ $\displaystyle=(-v_{11}\odot-
v_{21})\odot x_{1}^{w_{1}+w_{2}}\oplus(-v_{11}\odot-v_{22})\odot
x_{1}^{w_{1}}x_{2}^{w_{2}}\oplus(-v_{11}\odot-v_{23})\odot
x_{1}^{w_{1}}x_{3}^{w_{2}}$ $\displaystyle\oplus(-v_{12}\odot-v_{21})\odot
x_{1}^{w_{2}}x_{2}^{w_{1}}\oplus(-v_{12}\odot-v_{22})\odot
x_{2}^{w_{1}+w_{2}}\oplus(-v_{12}\odot-v_{23})\odot
x_{2}^{w_{1}}x_{3}^{w_{2}}$ $\displaystyle\oplus(-v_{13}\odot-v_{21})\odot
x_{1}^{w_{2}}x_{3}^{w_{1}}\oplus(-v_{13}\odot-v_{22})\odot
x_{2}^{w_{2}}x_{3}^{w_{1}}\oplus(-v_{13}\odot-v_{23})\odot
x_{3}^{w_{1}+w_{2}}.$

Its Newton polytope is $(w_{1}+w_{2})\cdot\Delta^{2}$, and it is depicted in
the lower right of Figure 5, with $w_{2}>w_{1}$. The image on the lower left
of Figure 5 is the Newton polytope in the special case where $w_{1}=w_{2}=1$.

We now apply the results of the previous subsection to translate the problem
of optimizing $f_{V,\mathbf{w}}$ into a problem on
$\mbox{Newt}(f_{V,\mathbf{w}})$. Since $f_{V,\mathbf{w}}$ is a function from
$\mathbbm{1}^{\perp}$ rather than $\mathbb{R}^{n}$, we need the following
result to apply the results of the previous subsection.

######  Proposition 2.30.

Given a max-plus tropical polynomial $f:\mathbbm{1}^{\perp}\to\mathbb{R}$, let
$\underline{N}$ be the subdivision of $\mathrm{Newt}(f)$ induced by the
coefficients of $f$. The minimum of $f$ is achieved on the cell dual to the
cell of the Newton polytope containing $\lambda\mathbbm{1}$ for any
$\lambda\in\mathbb{R}\setminus\\{0\\}$.

###### Proof.

Consider the following isomorphism
$\mathbb{R}^{n}/\mathbb{R}\mathbbm{1}\cong\mathbbm{1}^{\perp}$.

(12)
$\psi(x_{1},\ldots,x_{n-1})=\left(x_{1},\ldots,x_{n-1},-\sum_{i=1}^{n-1}x_{i}\right)$

The Newton polytope of $f$ lives in $\left(\mathbbm{1}^{\perp}\right)^{\ast}$.
The dual function of $\psi$ is below.

(13)
$\psi^{\ast}(x_{1},\ldots,x_{n})=\left(x_{1}-x_{n},\ldots,x_{n-1}-x_{n}\right)$

Then
$\psi^{\ast}\left(\lambda\mathbbm{1}\right)=\lambda(1-1,\ldots,1-1)=\overline{\mathbf{0}}$.
Lemma 2.16 says that for a tropical function $g:\mathbb{R}^{n}\to\mathbb{R}$,
the linear piece of $g$ whose dual cell contains $\overline{\textbf{0}}$ is
the linear piece achieving the minimum. Combining this result with the map
$\psi$, it follows that $f:\mathbbm{1}^{\perp}\to\mathbb{R}$ achieves its
minimum on the linear piece dual to the cell containing
$\lambda\mathbbm{1}_{n}$. ∎

###  2.4. Tropical Convexity

####  2.4.1. Tropical Factorization

Let $f=\bigodot_{i=1}^{m}f_{i}$ be a tropical signomial that factors into a
product of tropical signomials. The following theorem tells us how to compute
the tropical hypersurface $\mbox{tropV}(f)$ in terms of a subdivision of they
Cayley polytope. Let $P_{i}=\mbox{Newt}(f_{i})$, $P=\sum_{i}P_{i}$, and write
$f_{i}=\sum c_{i,\alpha}\mathbf{x}^{\alpha}$.

######  Theorem 2.31 (Corollary 4.9 in [6]).

Let $\underline{C}$ be the regular subdivision of
$\mbox{Cayley}(P_{1},\ldots,P_{m})$ induced by the weights
$w(e_{i},\alpha)=c_{i,\alpha}$. Then the mixed subdivision of $P$
corresponding to $\underline{C}$ coincides with the regular subdivsion of $P$
induced by the coefficients of $f$.

Recall that in the weighted tropical Fermat-Weber problem, $f_{V,\mathbf{w}}$
factors into linear pieces, and so theorem 2.31 applies.

######  Corollary 2.32.

Let $\underline{N}$ be the subdivision of $\mbox{Newt}(f_{V,\mathbf{w}})$
induced by the coefficients of $f_{V,\mathbf{w}}$. Then
$\underline{N}=\underline{C}\cap\\{\frac{1}{m}\mathbbm{1}_{m}\\}\times\mathbb{R}^{n}$.
In particular, the cell of $\underline{N}$ containing $\lambda\mathbbm{1}_{n}$
corresponds to the cell of $\underline{C}$ containing
$(\frac{1}{m}\mathbbm{1}_{m},\frac{|\mathbf{w}|}{mn}\mathbbm{1}_{n})$.

###### Proof.

The point
$(\frac{1}{m}\mathbbm{1}_{m},\frac{|\mathbf{w}|}{mn}\mathbbm{1}_{n})$ is the
barycenter of $\mbox{Cayley}(P_{1},\ldots,P_{m})$, so in particular, it lies
in
$\mbox{Cayley}(P_{1},\ldots,P_{m})\cap\\{\frac{1}{m}\mathbbm{1}_{m}\\}\times\mathbb{R}^{n}$.
Apply theorem 2.31 and then theorem 2.10. ∎

######  Proposition 2.33.

Let $P_{i}=\mbox{Newt}(f_{v_{i}})$, and let $\underline{C}^{\prime}$ be the
subdivision of $\mbox{Cayley}(P_{1},\ldots,P_{m})$ induced by
$w(\alpha_{j},e_{i})=-v_{ij}$. Then the subdivision of
$P_{\mathbf{w}}=\sum_{i=1}^{m}w_{i}P_{i}$ induced by the coefficients of
$f_{V,\mathbf{w}}$ is
$\underline{C}^{\prime}\cap\\{\mathbf{w}\\}\times\mathbb{R}^{n}$. In
particular, the cell of $\underline{P}$ containing $\mathbbm{1}_{n}$
corresponds to the cell of $\underline{C}^{\prime}$ containing
$(\mathbf{w},\frac{1}{n}\mathbbm{1}_{n})$.

###### Proof.

Apply proposition 2.11 to corollary 2.32. ∎

Figure 5. Subdivisions with weightings. Clockwise starting on top left: Two
tropical hyperplanes in $\mathbb{TR}^{2}$, the corresponding regular
subdivision of $\Delta^{1}\times\Delta^{2}$, the corresponding mixed
subdivision of $(w_{1}+w_{2})\Delta^{2}$ (weighted FW problem), and the
corresponding mixed subdivision of $2\Delta^{2}$ (unweighted FW problem).

####  2.4.2. Tropical Convex Hull

######  Definition 2.34.

The min-tropical convex hull of a set of points $A\subset\mathbbm{TP}^{n-1}$,
denoted $\mathrm{tconv}^{\min}(A)$ or just $\mathrm{tconv}(A)$, is the set of
all tropical linear combinations of points in $A$, that is,

(14)
$\mathrm{tconv}(A):=\\{\lambda_{1}\odot\mathbf{a}_{1}\oplus_{\min}\cdots\oplus_{\min}\lambda_{k}\odot\mathbf{a}_{k}\mid\lambda_{i}\in\mathbb{R},\mathbf{a}_{i}\in
A,k\in\mathbb{N}\\}.$

If $A$ is a finite set, then $\mathrm{tconv}(A)$ is called a tropical
polytope.

Note that the tropical convex hull is independent of the representatives in
$\mathbb{R}^{n}/\mathbb{R}\mathbbm{1}$ we choose for the points in $A$. That
is, if $\mathbf{v}_{i}^{\prime}=\alpha_{i}\odot\mathbf{v}_{i}$,
$\lambda_{i}\in\mathbb{R}$, then with
$\lambda_{i}^{\prime}=\lambda_{i}-\alpha_{i}$

$\lambda_{1}\odot v_{1}\oplus\cdots\oplus\lambda_{m}\odot
v_{m}=(\lambda_{1}-\alpha_{1})\odot
v_{1}^{\prime}\oplus\cdots\oplus(\lambda_{m}-\alpha_{m})\odot
v_{m}=\lambda_{1}^{\prime}\odot
v_{1}^{\prime}\oplus\cdots\oplus\lambda_{m}^{\prime}\odot v_{m}^{\prime}.$

######  Example 2.35.

Let $v_{1}=(0,0,0)$, and $v_{2}=(1,-1,0)$. The tropical polytope
$\mathrm{tconv}(v_{1},v_{2})$ consists of points of the form $\lambda_{1}\odot
v_{1}\oplus\lambda_{2}\odot v_{2}$, and is illustrated in Figure 6. The
tropical polytope consists of three points, connected by two classical line
segments.

Figure 6. The tropical line segment between $v_{1}=(0,0,0)$ and
$v_{2}=(1,-1,0)$.

It turns out that the tropical convex hull of the data points coincides with
the bounded part of the tropical hypersurface $\mbox{tropV}(f_{V})$.

######  Theorem 2.36 (Theorem 5.2.11 in [8]).

The bounded part of the tropical hypersurface $f_{V}$ is
$\mbox{tconv}^{\min}(v_{1},\ldots,v_{m})$.

We now see that $f_{V}$ and $f_{V,\mathbf{w}}$ define the same tropical
hypersurface. Thus, the bounded part $\mbox{tropV}(f_{V,\mathbf{w}})$ is the
tropical convex hull of the $v_{i}$’s.

######  Lemma 2.37.

For any $\mathbf{w}\in\mathbb{R}^{m}$,
$\mbox{tropV}(f_{V})=\mbox{tropV}(f_{V,\mathbf{w}})$.

###### Proof.

By lemma 2.19, $\mbox{tropV}(f_{v}^{w})=\mbox{tropV}(f_{v})$ for any
$v\in\mathbbm{1}^{\perp}$, and any $w>0$. Applying lemma 2.19 to $f_{V}$ and
$f_{V,\mathbf{w}}$, it follows that

$\mbox{tropV}(f_{V})=\bigcup_{i=1}^{m}\mbox{tropV}(f_{v_{i}})=\mbox{tropV}(f_{V,\mathbf{w}}).\qed$

######  Corollary 2.38.

The bounded part of the tropical hypersurface $f_{V,\mathbf{w}}$ is
$\mbox{tconv}^{\min}(v_{1},\ldots,v_{m})$.

##  3\. Solving the Weighted Tropical Fermat-Weber Problem

In this section, we use combinatorics and tropical geometry to solve the
weighted Fermat-Weber problem for $\mathbbm{TP}^{n-1}$ equipped with the
tropical asymmetric distance. We begin by discussing the extrema of tropical
polynomials.

###  3.1. Containment

######  Theorem 3.1.

Given data points $v_{1},\ldots,v_{m}\in\mathbb{R}^{n}/\mathbb{R}\mathbbm{1}$,
and weights $w_{1},\ldots,w_{m}>0$, the weighted Fermat-Weber points under the
tropical asymmetric metric are a cell of $\mbox{tconv}(v_{1},\ldots,v_{m})$.

###### Proof.

According to corollary 2.38, $\mbox{tconv}(v_{1},\ldots,v_{m})$ is the bounded
part of $\mbox{tropV}(f_{V,\mathbf{w}})$. The bounded cells of
$\mbox{tropV}(f_{V,\mathbf{w}})$ are exactly those cells dual to interior
cells of the Newton polytope, so by proposition 2.30, it suffices to show that
$\lambda\mathbbm{1}$ is in the interior of $\mbox{Newt}(f_{V,\mathbf{w}})$.
The vertices of $\mbox{Newt}(f_{V,\mathbf{w}})$ are $|\mathbf{w}|e_{i}$, and
their average, $\frac{|\mathbf{w}|}{n}\mathbbm{1}$, is in the interior of the
Newton polytope. This proves the Fermat-Weber points are achieved on a bounded
cell of $\mbox{tropV}(f_{V,\mathbf{w}})$, and therefore form a cell of
$\mbox{tconv}(v_{1},\ldots,v_{m})$. ∎

###  3.2. Any cell can be the weighted Fermat-Weber cell

The following result shows that we can pick weights $w_{i}$ so that the
weighted barycenter lies in any interior cell of the subdivision of the Cayley
polytope. This finishes the proof of the main theorem.

######  Theorem 3.2.

Given some data points
$\mathbf{v}_{1},\ldots,\mathbf{v}_{m}\in\mathbb{R}^{n}/\mathbb{R}\mathbbm{1}$,
and any simplex $S$ in $\Delta^{r-1}\times\Delta^{n-1}$, which intersects the
relative interior of $\Delta^{r-1}\times\Delta^{n-1}$, there is a choice of
weights $w_{1},\ldots,w_{m}\in[0,1]$ with $\sum w_{i}=1$, so that $S$ contains
the point $(w_{1},\ldots,w_{m},\frac{1}{n}\mathbbm{1}_{n})$.

The proof uses a well-known correspondence between subsets of the vertices of
$\Delta^{m-1}\times\Delta^{n-1}$ and subgraphs of $K_{m,n}$ (the complete
bipartite graph with $m$ left vertices, and $n$ right vertices), which we now
briefly recall (see [3, §6.2.2] for more details). The vertex $(e_{i},e_{j})$
in a simplex $S\subseteq\Delta^{r-1}\times\Delta^{n-1}$ corresponds to the
edge between left vertex $i$ and right vertex $j$ in the bipartite graph.
Thus, a subset of vertices $A\subseteq\Delta^{r-1}\times\Delta^{n-1}$
corresponds to the subgraph of $K_{m,n}$.

######  Example 3.3 (Simplex-Forest Correspondence for $m=n=2$).

The product of two $2$-simplices (i.e. line segments) is a square; the
corresponding bipartite graph has two left vertices and two right vertices.
Both are illustrated in fig. 7. The vertex $(e_{i},e_{j})$ in the simplex
corresponds to the edge $(l_{i},r_{j})$ in the bipartite graph. For example,
the top left vertex of the shaded gray simplex, $(e_{1},e_{2})$, corresponds
to the edge $(l_{1},r_{2})$. The shaded simplex is full-dimensional, so it
corresponds to a spanning tree of $K_{2,2}$ (see lemma 3.4).

Figure 7. Vertices in the product of simplices (left) correspond to the color-
coded edges of the bipartite graph (right).

######  Lemma 3.4 (Lemma 6.2.8 in [3]).

Let $A$ be a subset of the vertices of $\Delta^{m-1}\times\Delta^{n-1}$. Then,

  1. (a)

$\mathrm{conv}(A)$ is a simplex if and only if the corresponding subgraph of
$K_{m,n}$ is a forest.

  2. (b)

$\mathrm{conv}(A)$ is full dimensional if and only if the corresponding
subgraph of $K_{m,n}$ is spanning and connected.

###### Proof of Theorem 3.2.

Let $S$ be a simplex in $\Delta^{r-1}\times\Delta^{n-1}$, and let $F$ be the
corresponding forest in $K_{r,n}$. Assume that
$S\cap\mbox{int}\left(\Delta^{r-1}\times\Delta^{n-1}\right)\neq\emptyset$ (so
$F$ is a spanning forest).

A point $(\mathbf{p},\mathbf{q})\in\mathbb{R}^{m}\times\mathbb{R}^{n}$ lies in
$S$ if it can be written as a convex combination of the vertices of $S$. In
terms of the forest $F$, $(\mathbf{p},\mathbf{q})$ lies in $S$ if there exist
$\lambda(e)>0$ for each edge $e\in E(F)$ such that the sum of edge weights on
any left vertex adds up to the corresponding $\mathbf{p}$ coordinate, and the
sum of edge weights on any right vertex adds up to the corresponding
$\mathbf{q}$ coordinate. Let $r(e)$ be the node on the right side connected to
$e$, and let $\ell(e)$ be the node on the left side connected to $e$. The
choice of $\lambda$’s in (15) leads to a valid choice of weights
$w_{1},\ldots,w_{m}$ (given in (16)) so that $S$ contains the weighted
barycenter.

(15) $\lambda(e):=\frac{1}{n\cdot\deg r(e)}.$ (16) $w_{i}:=\sum_{e\text{ s.t.
}\ell(e)=i}\lambda(e).$

The equations in (17) show that the weights on any right node sum to
$\frac{1}{n}$ (since $F$ is spanning, every vertex has at least one edge); by
definition, the weights on the $i$th left node sum to $w_{i}$. It follows that
$\mathbf{b}=(w_{1},\ldots,w_{m},\frac{1}{n}\mathbbm{1})$ lies in the relative
interior of $S$.

(17)
$\sum_{r(e)=j}\frac{1}{n\cdot\deg(j)}=\frac{1}{n}\sum_{r(e)=j}\frac{1}{\deg(j)}=\frac{1}{n}\deg(j)\frac{1}{\deg(j)}=\frac{1}{n}.$

Moreover, $w_{1},\ldots,w_{m}$ is a valid choice of weights for the Fermat-
Weber problem. The weights $w_{i}$ are positive because $F$ is spanning, so
the sum in (16) is never empty; The equations in (18) show that the $w_{i}$
sum to one.

(18)
$\sum_{i=1}^{m}w_{i}=\sum_{e}\lambda(e)=\sum_{j=1}^{n}\sum_{r(e)=j}\frac{1}{n\cdot\deg(j)}=\sum_{j=1}^{n}\frac{1}{n}\sum_{r(e)=j}\frac{1}{\deg(j)}=\sum_{j=1}^{n}\frac{1}{n}=n\frac{1}{n}=1.\qed$

######  Corollary 3.5.

Given a cell $T$ in the tropical polytope
$\mathrm{tconv}(\mathbf{v}_{1},\ldots,\mathbf{v}_{m})$, there is a choice of
weights $w_{1},\ldots,w_{m}$ so that $T$ is the set of weighted tropical
Fermat-Weber points for $\mathbf{v}_{1},\ldots,\mathbf{v}_{m}$ with weights
$w_{1},\ldots,w_{m}$.

### Acknowledgements

We are grateful to David Speyer and Michael Joswig for helpful conversations.
This work was started at the “Algebra of phylogenetic networks” workshop held
at the University of Hawai‘i at Mānoa from May 23 - 27, 2022 which was
supported by the National Science Foundation under grant DMS-1945584. The
first author was supported by National Science Foundation Graduate Research
Fellowship under Grant No. DGE-1841052, and by the National Science Foundation
under Grant No. 1855135.

## References

  * [1] David Bryant.  A classifcation of consensus methods for phylogenetics.  BioConsensus, DIMACS, 61, 05 2002.
  * [2] Andrei Comăneci and Michael Joswig.  Tropical medians by transportation.  https://arxiv.org/abs/2205.00036, 2022.
  * [3] Jesús A. De Loera, Jörg Rambau, and Francisco Santos.  Triangulations, volume 25 of Algorithms and Computation in Mathematics.  Springer-Verlag, Berlin, 2010.  Structures for algorithms and applications.
  * [4] Frank A. Fetter.  Alfred weber’s theory of the location of industries. alfred weber.  Journal of Political Economy, 38(2):232–234, 1930.
  * [5] Birkett Huber, Jörg Rambau, and Francisco Santos.  The Cayley trick, lifting subdivisions and the Bohne-Dress theorem on zonotopal tilings.  J. Eur. Math. Soc. (JEMS), 2(2):179–198, 2000.
  * [6] Michael Joswig.  Essentials of tropical combinatorics, volume 219 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, [2021] ©2021.
  * [7] Bo Lin and Ruriko Yoshida.  Tropical Fermat-Weber points.  SIAM J. Discrete Math., 32(2):1229–1245, 2018.
  * [8] Diane Maclagan and Bernd Sturmfels.  Introduction to tropical geometry, volume 161 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2015.
  * [9] Yurii Nesterov.  Introductory lectures on convex optimization, volume 87 of Applied Optimization.  Kluwer Academic Publishers, Boston, MA, 2004.  A basic course.
  * [10] Francisco Santos.  The Cayley trick and triangulations of products of simplices.  In Integer points in polyhedra—geometry, number theory, algebra, optimization, volume 374 of Contemp. Math., pages 151–177. Amer. Math. Soc., Providence, RI, 2005.
  * [11] Bernd Sturmfels.  On the Newton polytope of the resultant.  J. Algebraic Combin., 3(2):207–236, 1994.


# Impact of edge morphology and chemistry on nanoribbons’ gapwidth

Elisa Serrano Richaud  Université Paris-Saclay, ONERA, CNRS, Laboratoire
d’etude des microstructures, 92322, Châtillon, France  Sylvain Latil
Université Paris-Saclay, CEA, CNRS, SPEC, 91191 Gif-sur-Yvette, France
Lorenzo Sponza  Université Paris-Saclay, ONERA, CNRS, Laboratoire d’etude des
microstructures, 92322, Châtillon, France  European Theoretical Spectroscopy
Facility (ETSF), B-4000 Sart Tilman, Liège, Belgium

(today)

###### Abstract

In this work, we scrutinise theoretically how the gap of C and BN armchair
nanoribbons changes upon variations of the bond length between edge atoms and
their distance from passivating species. Our DFT calculations indicate that
the gap of C-based nanoribbons is more sensitive to the relaxation of the
bonding length between edge atoms (morphology) whereas in BN-nanorribons it is
more sensitive to the distance between edge atoms and passivating hydrogens
(chemical environment). To understand the origin of these two different
behaviours, we solved a tight-binding ladder model numerically and at the
first-order perturbation theory, demonstrating that the different dependence
is due to the interference of the wavefunctions of the top valence and the
bottom conduction states.

##  I Introduction

In recent decades, graphene and hexagonal boron nitride (BN) have attracted a
great deal of interest because of their remarkable transport and optical
properties Geim and Novoselov (2007); Watanabe  _et al._ (2009); Castro Neto
_et al._ (2009); Weng  _et al._ (2016); Zhang  _et al._ (2017). A much
explored way to modulate them is by adding extra confinement (as in 2D quantum
dots, nanoribbons or nanotubes). The presence of confining edges endows them
with novel size-dependent features dominated by the characteristics of the
edge itself. This is why graphene and BN nanoribbons are often classified
according to their edge shape, which can be zig-zag, armchair, fall in an
intermediate chiral angle, or present structures that require a more general
nomenclature Ezawa (2006). In zig-zag nanoribbons, well localised edge-state
are formed which confer antiferromagnetic properties to C-based zig-zag
nanoribbons Ezawa (2006); Fujita  _et al._ (1996); Nakada  _et al._ (1996);
Behzad and Chegel (2018); Wakabayashi  _et al._ (1999); Yang  _et al._ (2007);
Son  _et al._ (2006). Instead, BN-based zig-zag nanoribbons have an indirect
gap and display an intrinsic dipole moment Behzad and Chegel (2018); Nakamura
_et al._ (2005); Park and Louie (2008); Wang  _et al._ (2011); Shyu (2014);
Topsakal  _et al._ (2009); Zhang and Guo (2008); Ding  _et al._ (2009). At
variance, both graphene Ezawa (2006); Nakada  _et al._ (1996); Behzad and
Chegel (2018); Wakabayashi  _et al._ (1999); Yang  _et al._ (2007); Son  _et
al._ (2006); Nematian  _et al._ (2012); Barone  _et al._ (2006); Jha  _et al._
(2019); Nishad  _et al._ (2020); Prabhakar and Melnik (2019); Ren  _et al._
(2007); Lu  _et al._ (2009); Raza and Kan (2008) and BN Behzad and Chegel
(2018); Park and Louie (2008); Wang  _et al._ (2011); Shyu (2014); Zhang and
Guo (2008); Topsakal  _et al._ (2009) armchair nanoribbons (AGNR and ABNNR),
have no magnetic states and display a direct size-dependent gapwidth To take
full advantage of this richness of properties, several methods have been
explored including the application of external electromagnetic fields Behzad
and Chegel (2018); Wakabayashi  _et al._ (1999); Park and Louie (2008); Zhang
and Guo (2008); Raza and Kan (2008), strain Topsakal  _et al._ (2009);
Prabhakar and Melnik (2019); Li  _et al._ (2010) and edge engineering Topsakal
_et al._ (2009); Barone  _et al._ (2006); Ding  _et al._ (2009); Jha  _et al._
(2019); Nishad  _et al._ (2020); Lu  _et al._ (2009); Prabhakar and Melnik
(2019); Zheng  _et al._ (2013); Ren  _et al._ (2007).

As a matter of fact, the edge characteristics are crucial for the performances
of nanoribbon-based devices such as transistors, interconnects and logical
devicesMurali  _et al._ (2009); Zheng  _et al._ (2013); Das  _et al._ (2021);
Marmolejo-Tejada and Velasco-Medina (2016); Nishad  _et al._ (2020); Saraswat
_et al._ (2021), photovoltaic applications Osella  _et al._ (2012); Saraswat
_et al._ (2021), or chemical sensingMehdi Pour  _et al._ (2017); Saraswat  _et
al._ (2021). Experimentally, edge engineeringSchwab  _et al._ (2015); Wang
_et al._ (2016); Osella  _et al._ (2012), chemical treatmentDauber  _et al._
(2014) or selective passivationZheng  _et al._ (2013) have been demonstrated
to have a significant impact on the device quality, precisely because of their
action on the edges.

Figure 1: Scheme of relaxed AGNR5 and ABNNR5 structures. Unitary cells are
reported as black dashed rectangles and some relevant structural parameters
are labelled in red. In blue, the row index $j$ and the edge or inner
character of atoms are also reported.

Alterations of the electronic structure due to edge modifications can be
divided into morphology effects (variation of the bondlengths) and chemistry
effects (variation of the passivating species and their distance from the
edges)Ezawa (2006); Lu  _et al._ (2009). The sensitivity of AGNR and ABNNR gap
to the passivation has been investigated by many authors Ezawa (2006); Barone
_et al._ (2006); Ding  _et al._ (2009); Jha  _et al._ (2019); Nishad  _et al._
(2020); Topsakal  _et al._ (2009); Lu  _et al._ (2009); Prabhakar and Melnik
(2019); Zheng  _et al._ (2013); Ren  _et al._ (2007) who showed that its
effect depends on the type of atoms involved, and/or on the number and
position of the passivated sites. Most of these first-principle studies Barone
_et al._ (2006); Ding  _et al._ (2009); Jha  _et al._ (2019); Nishad  _et al._
(2020); Topsakal  _et al._ (2009); Prabhakar and Melnik (2019); Ren  _et al._
(2007) discuss the role of passivation on fully relaxed structures, so
morphology and chemistry effects are actually treated on the same footing. At
best of our knowledge, only two studies Ezawa (2006); Lu  _et al._ (2009) have
been conducted in frameworks that separate the two effects, but in both
publications the focus is put on other aspects than the dependence of the
gapwidth on morphology and chemistry modifications.

On the other hand, rare are the studies on genuine morphological effects Ezawa
(2006); Lu  _et al._ (2009); Son  _et al._ (2006) and they are done only on
AGNRs. Actually, in Huang  _et al._ (2009, 2012), a thorough study of edge
reconstruction and edge stress in graphene and BN nanoribbons is actually
carried out, but the investigation stops at a stability level and the relation
to the gapwidth is not explored. However, both effects seem to be decisive in
determining the gap of nanoribbons and we deemed that the subject deserved a
more focused study.

In this article, we employ density functional theory (DFT) to study the
evolution of the gap, the top valence (TV) and the bottom conduction (BC)
states of AGNRs and ABNRs as a function of the nanoribbon size upon variations
of the distance between edge atoms and between these and the passivating
species. Our objective is to compare the effect of morphological and chemical
variations on the gapwidth and understand which of them is dominant and in
which situation. We demonstrate that the response of the gapwidth to changes
of the distance between edge atoms (morphology) or between edge atoms and
passivating atoms (chemical environment) is opposite in the two materials and
we rationalise this different behaviour by means of a tight-binding model
which we solved both numerically and perturbatively.

##  II Structural and computational details

All nanoribbons studied in this article have armchair edges passivated with H
atoms. They form an infinite periodic structure in the $y$ direction and are
confined along $x$. The extension of the periodic cell along $y$ is the cell
parameter $a$, while the width is expressed by the number $N_{a}$ which
indicates the number of dimers aligned along $y$ inside the unitary cell
(number of rows). To indicate a specific structure we will attach the index
$N_{a}$ after the label of the material, as in Figure 1, so for instance AGNR5
designates an armchair graphene nanoribbon of size $N_{a}=5$.

Density functional theory calculations were carried out within the generalized
gradient approximation using the PBE Perdew  _et al._ (1996) exchange
correlation potential as implemented in the Quantum ESPRESSO Giannozzi  _et
al._ (2009) simulation package. Long-range van der Waals corrections were
included via the DFT-D2 method Grimme (2006). To avoid interactions between
consecutive cells, we included 15 Å and 20 Å of empty space in the $z$ and $x$
directions respectively. In electron density calculations and relaxation runs,
the periodic axis was sampled with 20 k-points centered in $\Gamma$
(corresponding to 11 irreducible k-points). This mesh was dense enough to
converge total energies in the smallest nanoribbons. For density of states
(DOS) calculations, a five times denser sampling was adopted for all systems
and the resulting spectra have been broadened with a Gaussian distribution
with a width of $0.02$ eV.

We used norm-conserving pseudopotentials Hamann (2013) and set the kinetic
energy cutoff at 80 Ry in both materials. It is worth stressing that using a
large vertical empty space and a high energy cutoff is essential even in the
relaxation runs in order to prevent nearly free-electron states from hanging
below the $p_{z}$ states hence jeopardizing the gap description. In fact, as
already well known for free-standing layers Posternak  _et al._ (1983, 1984);
Blase  _et al._ (1995); Paleari (2019); Latil  _et al._ (2022) and nanotubes
Blase  _et al._ (1994a, b); Hu  _et al._ (2010) in BN nanomaterials there is a
competition at the bottom conduction between $2p_{z}$ and $3s$ states, whose
right alignment requires a dedicated convergence study. If sometimes one can
overlook this issue in BN layers, because the two competing states originate
direct and indirect band gaps, this is not the case in ABNNRs where both
states give rise to a direct gap at $\Gamma$.

In non-relaxed structures, all atoms occupy the sites of a regular honeycomb
lattice with an inter-atomic distance of 1.42 Å. Structural relaxation runs
have been performed with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm
for all systems with the stopping criterion of all forces being lower than
$5\times 10^{-5}$ eV/Å. We allowed variations of the cell parameter $a$ and
all atomic positions. As clarified in the following, we also run some
calculations letting only specific atoms to move. In Figure 1 we report the
relaxed structures of AGNR and ABNNR at $N_{a}=5$ for sake of example, and we
introduce some notable structural parameters. In the AGNRs, the main
modifications with respect to non-relaxed structures are a contraction of the
distance between edge atoms $d_{E}$ and between C and H $d_{HC}$. In ABNNR, we
observe a similar contraction of the B-N distance on the edges $d_{E}$, and
different contractions of the distances between H-B and H-N ($d_{HB}\neq
d_{HN}$). We observed also that these modifications are basically independent
on the size of the nanoribbon both qualitatively and quantitatively, so the
structural parameters undergo minimal variations when comparing nanoribbons of
different size.

##  III Gap edge states

###  III.1 AGNRs

The electronic structure of AGNRs has been already studied in the past Ezawa
(2006); Nakada  _et al._ (1996); Behzad and Chegel (2018); Wakabayashi  _et
al._ (1999); Yang  _et al._ (2007); Osella  _et al._ (2012); Saraswat  _et
al._ (2021); Nematian  _et al._ (2012); Barone  _et al._ (2006); Jha  _et al._
(2019); Nishad  _et al._ (2020); Lu  _et al._ (2009); Prabhakar and Melnik
(2019); Ren  _et al._ (2007). Both non-relaxed and relaxed ribbons display a
band gap at $\Gamma$ of gapwidth $\Delta_{N_{a}}$. Because of the 1D
confinement, the gapwidth falls in one of the three families $N_{a}=3m-1$,
$3m$ or $3m+1$ (with $m\in\mathbb{N}^{*}$). Each family follows a different
trend which asymptotically tends to zero for growing nanoribbon sizes and
follows the general rule $\Delta_{3m-1}<\Delta_{3m}<\Delta_{3m+1}$. This is
depicted in Figure 2 where we plot the gapwidth of AGNRs versus $N_{a}$ for
both non-relaxed and relaxed structures (red dashed and solid blue curves).
The effect of relaxation is to open the gap by about 0.1 eV in families
$N_{a}=3m+1$ and $3m-1$, while in the $N_{a}=3m$ the opening is observed only
in small nanoribbons, while the gap closes in larger ones. Our results are in
quantitative agreement with previous works both for relaxed Son  _et al._
(2006); Yang  _et al._ (2007); Ren  _et al._ (2007); Barone  _et al._ (2006),
and unrelaxed simulationsLu  _et al._ (2009).

To characterise better the gap states, we analyzed in more detail the nature
of the TV and the BC states at $\Gamma$ in the relaxed structures. In panels
a) and b) of Figure 3, we report the band structure and the density of states
(DOS) of the AGNR8, chosen as a representative example. For sake of
comparison, in panel b) we also report the orbital-projected DOS and the DOS
of an infinite graphene sheet with the same inter-atomic distance. The DOS
around the gap (from -1.5 eV to 1.5 eV) displays neat van Hove singularities
arranged more or less symmetrically with respect to the middle of the gap. As
the inset of panel b) shows clearly, the states composing the gap are entirely
of $p_{z}$ character. They form a $\pi$ bonding with nodes on the $xy$ plane,
as expected. Instead, the first empty $\sigma$ state is found at 3 eV above
the BC. To go deeper in the analysis of the gap-edge states, we look at the
site-projected DOS. We integrated the bare data inside an interval of 0.1 eV
encompassing the TV and BC (shaded bands in the instet of Figure 3b). The
outcome of this analysis is summarised in Figure 3c), where the site-projected
DOS of gap-edge states is reported as a function of the row index (note that
the curves are plotted on the same $y$ axis). At variance from what observed
in zigzag nanoribbons Fujita  _et al._ (1996), the gap states are not
concentrated on the edge atoms, but rather delocalized throughout the full
nanoribbon and present a modulation that nicely displays the characteristics
of a static wave. This observation is confirmed by the wave-like modulation of
the charge probability $|\psi(\mathbf{r})|^{2}$ associated with the TV and BC
states, reported aside panel c). The wavefunction plot shows also that there
is no spill-out on the passivating hydrogens and that, with respect to the
edge-bbondings $d_{E}$, TV and BC states display respectively a bonding and an
antibbonding character.

Figure 2: Energy gap of graphene nanoribbons as a function of the width
$N_{a}$. Relaxed calculations (blue solid line), unrelaxed (red dashed line)
and tight-binding numerical solution (black dotted) with parameters indicated
in Table 1. The three families are reported with different symbols. A blue
arrow at $N_{a}=8$ indicate the nanoribbon chosen for the analysis presented
in Figure 3. Figure 3: Electronic structure of the relaxed AGNR8. a) Band
structure. b) and Inset: Total density of states (thick black) and projected
on $p_{z}$ orbital character (red bullets) compared with the DOS of the
graphene sheet (dashed blue). c) Row-projected DOS from the integration of the
total DOS around the band-edge states (shaded areas of panel b) and charge
density associated with the TV and BC states at $\Gamma$.

###  III.2 ABNNRs

The gapwidth of ABNNRs fall in the same three families with the same
hierarchyZhang and Guo (2008); Li  _et al._ (2010); Topsakal  _et al._ (2009).
This similarity with the graphene ribbons is actually quite general and can be
understood from a simple tight-binding model (see section IV.2). The evolution
of the ABNNRs gapwidth for sizes going from $N_{a}$=5 to 19 in the relaxed and
non-relaxed configurations is presented in Figure 4 by the solid blue and the
red dashed lines. The non-relaxed structures present a gap that monotonically
tends to the limit $N_{a}\rightarrow\infty$ in a way that is similar to non-
passivated calculationsTopsakal  _et al._ (2009). We estimate
$N_{a}\rightarrow\infty=3.885$ eV from the weighted average of the curves
extrapolated at $1/N_{a}=0$ (cfr. inset of the Figure). This value is about
0.8 eV lower than the gapwidth of the isolated BN sheet (4.69 eV in PBE). All
these aspects are consistent because, as it will become clearer later, in non-
relaxed calculation, H atoms are too far to saturate efficiently the dangling
bonds located at the edges of the ribbbon. As a consequence, these form edge
states inside the gap that lower the gapwidth similarly to what happens in
non-passivated (bare) ribbons.

As a result of the structural optimisation, the gapwidth of all families opens
and tends to an asymptotic limit that is still about 0.1 eV lower than in the
isolated monolayer, in agreement with similar calculations Park and Louie
(2008); Topsakal  _et al._ (2009). This discrepancy is ascribed to a non-
negligible edge contribution to the BC state, obviously absent in the isolated
monolayer (cfr. the row-projected DOS analysis here below, andPark and Louie
(2008)). Finally, we note that the first empty $\sigma$ state, i.e. the near
free-electron state, is only 0.5 eV above the BC.

Figure 4: Energy gap of BN nanoribbons as a function of the size $N_{a}$.
Relaxed DFT (blue solid line), unrelaxed (red dashed line) and the numerical
tight-binding solution (Table 1). The three families are reported with
different symbols. Horizontal dashed lines indicate the gapwith of the DFT hBN
sheet (4.69 eV) and the asymptotic $N_{a}=\infty$ limit ($\sim$3.885 eV). The
blue arrow pointing at the ABNNR8 indicates the system analysed in Figure 4.
Inset: extrapolation of non-relaxed calculations at $1/N_{a}=0$. The red arrow
in the inset indicates the $N_{a}=\infty$ limit as the weighted average of the
extrapolation of the three families.

Similarly to what done before, in Figure 5 we report the band structure, the
projected DOS and the row-resolved DOS of the TV and BC states of the
representative ABNNR8 system. We verify that the TV and the BC states are
formed essentially of N-centered and B-centered $p_{z}$ orbitals respectively.
The row-projected DOS of both TV and BC, reported in panel c), shows again a
very nice static-wave-like modulation with nodes in rows 3 and 6, but at
vairance with the AGNR8 case, here the TV and BC states localize differently:
while the TV states are delocalised on the entire nanoribbon as in the
previous case, the BC states are clearly peaked at the edges. The
visualization of the associated charge density confirms that the TV state is
characterised by a wavefunction equally delocalised on all the N atoms except
those on rows 3 and 6. Instead, the BC state presents a wavefunction more
concentrated on the edge B atoms with non negligible tails touching
passivating H and edge nitrogens, in contrast to the isolated monolayer.

The compared study of the TV and BC states of AGNRs and ABNNRs suggests that
the gap of the two materials responds differently to modifications of the
morphology and the passivation of the edges. To test this intuition, we have
performed a detailed analysis by separating the two effects.

Figure 5: Electronic structure of the relaxed ABNNR8. a) band structure; b)
total density of states (thick black) and projected on $p_{z}$ orbital
character (red and green dotted for B and N states) compared to the hBN sheet
DOS (dashed blue). c) Row-projected DOS integrated around the band-edge states
(shaded areas of panel b). Insets: charge density of the TV and BC states at
$\Gamma$.

##  IV Morphology VS Chemistry of the edges

###  IV.1 Distinguishing the effects through selective relaxation in DFT

Several investigations can be found in literature on the effects of edge
reconstruction on the gapwidth of AGNR and ABNNREzawa (2006); Son  _et al._
(2006); Barone  _et al._ (2006); Ding  _et al._ (2009); Jha  _et al._ (2019);
Nishad  _et al._ (2020); Topsakal  _et al._ (2009); Prabhakar and Melnik
(2019); Lu  _et al._ (2009); Ren  _et al._ (2007). However, a study that
systematically compares the effects of passivation and edge morphology is
absent. Here we monitor the gapwidth in the family $N_{a}=3m-1$ by relaxing
separately the H-X distances $d_{HX}$ ($X$ = C, B or N) and the C-C or B-N
distance on the edges $d_{E}$. We did calculate data from the other two
families but we do not report them because they have qualitatively the same
behaviour.

In Figure 6, a variation of $d_{HX}$ is represented by a change in the line’s
type (color and dash), while a variation of $d_{E}$ is represented by a change
in the symbols (colour filled or empty). Let us examine first the case of
AGNRs in panel a). We can start from a non-relaxed configuration where all
atoms are equidistant $d_{HC}$=$d_{E}$=1.42 Å (empty bullets, red dashed
line), then we reduce $d_{HC}$ to its relaxed value 1.08 Å (empty bullets,
blue solid line). We observe that there is basically no variation on the
AGNRs’ gapwidth. Instead, contracting the edge bonds from $d_{E}$=1.42 Å to
$d_{E}$=1.36 Å opens the gap by around 0.15 eV irrespective of the value of
$d_{HC}$. Consequently, we conclude that in AGNRs, the variations of the
gapwidth induced by the relaxation and reported in Figure 2 are essentially
due to changes of bond length $d_{E}$ between C atoms at the edge.
Interestingly, this gap opening is approximately independent on the width of
the ribbon.

Passing now to the study of ABNNRs (bottom panel), we observe an opposite
behaviour. The gapwidth undergoes very small changes upon relaxation of
$d_{E}$, whereas the passage from the unrelaxed H-B and H-N distance (1.42 Å)
to the relaxed values clearly opens the gap by about 0.8 eV. To be more
precise, by changing separately the two distances $d_{HB}$ and $d_{HN}$ (not
shown), we found that it is the bonding between H and B that plays a major
role in the opening of the gapwidth, indicating a dominant contribution from
conduction states consistently with the observations we drew from Figure 5.
According to this analysis, the gapwidth of ABNNRs is more sensitive to the
passivation than to the very morphology of the edge. Once again we notice that
the gap opening is basically independent on $N_{a}$. This clarifies why our
non-relaxed DFT gapwidth look very similar to the non-passivated results of
Topsakal and coworkers Topsakal  _et al._ (2009).

Figure 6: Gapwidth of the $N_{a}=3m-1$ family of a) AGNRs and b) ABNNRs. Full
(empty) symbols stand for relaxed (non-relaxed) edge-atom bondings. Blue solid
(red dashed) lines for relaxed (non-relaxed) passivating-to-edge-atoms
bondings.

###  IV.2 Unperturbed tight-binding model

To investigate further the reasons of this different behaviour, we generalise
a ladder tight-binding model introduced initially for AGNRs to the
heteroatomic case. Changes in the edge passivation and morphology will be
successively introduced through variations of the on-site and hopping
parameters of the model, as suggested in Ezawa (2006); Son  _et al._ (2006),
and the modified Hamiltonian solved both numerically and perturbativelySon
_et al._ (2006).

Figure 7: Scheme of the ladder model of width $N_{a}=8$. The first neighbours
distance is $a$, the index $j$ defines the position of a dimer. Atoms $\mu=1$
are placed above $\mu=2$ if $j$ is even, below if $j$ is odd.

Following references Ezawa (2006); Fujita  _et al._ (1996); Nakada  _et al._
(1996); Wakabayashi  _et al._ (1999); Son  _et al._ (2006); Shyu (2014);
Nematian  _et al._ (2012), the gap of an armchair nanoribbon whose TV and BC
states are formed of $p_{z}$ orbitals can be described with a ladder tight-
binding model as the one reported in Figure 7. The Hamiltonian of the model
reads:

$H^{0}=\sum_{j,\mu}\left(\epsilon_{\mu j}\ket{\Phi_{\mu
j}}+\sum_{j^{\prime},\mu^{\prime}}t_{\mu\mu^{\prime}jj^{\prime}}\ket{\Phi_{\mu^{\prime}j^{\prime}}}\right)\bra{\Phi_{\mu
j}}\,.$ (1)

The index $j\in\left[1,N_{a}\right]$ labels the position of a dimer in the $x$
coordinate (row coordinate), while $\mu=1,2$ indicates the atomic site within
the dimer ($C_{1}$ or $C_{2}$ in AGNRs and $B$ or $N$ in ABNNRs). The basis
function $\braket{r}{\Phi_{\mu j}}=\Phi_{\mu}({\mathbf{r}}-{\mathbf{r}}_{j})$
is the $p_{z}$ orbital of the atom $\mu$ of the dimer placed at
$\mathbf{r}_{j}=\hat{x}(j-1)a$. For $\mu=1$,
$\Phi_{\mu}({\mathbf{r}}-{\mathbf{r}}_{j})$ is centered on the bottom rung if
$j$ is odd and in the upper rung if $j$ is even, and the opposite for $\mu=2$.

At the unperturbed level, $\epsilon_{\mu j}$ does not depend on the row-index
$j$ and is equal to $\epsilon$ for $\mu=1$ and $-\epsilon$ for $\mu=2$, with
$\epsilon\geq 0$. In the first-neighbour approximation, the hopping term
$t_{\mu\mu^{\prime}jj^{\prime}}=t\in\mathbb{R}$ if $\mu\neq\mu^{\prime}$ and
$j-1\leq j^{\prime}\leq j+1$ and vanishes otherwise. The unperturbed solutions
of this model are:

$E^{0}_{n\pm}=\pm\sqrt{\epsilon^{2}+\tau^{2}_{n}}=\pm\mathcal{E}_{n}\,,$ (2)

where $\tau_{n}=t\left[1+2\cos\left(\theta_{n}\right)\right]$, the discrete
index $n$ comes from the confinement in the $x$ direction and
$\theta_{n}=n\pi/(N_{a}+1)$. The eigenfunction associated to these states read

$\Psi_{n\pm}=\sum_{j=1}^{N_{a}}\sum_{\mu=1,2}\sin\left(j\theta_{n}\right)D^{n\pm}_{\mu}\Phi_{\mu}(\mathbf{r}-\mathbf{r}_{j})$
(3)

with

$\begin{split}D_{1}^{n\pm}&=\sqrt{\frac{\mathcal{E}_{n}\pm\epsilon}{(N_{a}+1)\mathcal{E}_{n}}}\\\
D_{2}^{n\pm}&=\pm\text{sgn}\left(\tau_{n}\right)\sqrt{\frac{\mathcal{E}_{n}\mp\epsilon}{(N_{a}+1)\mathcal{E}_{n}}}\end{split}$
(4)

where the function $\text{sgn}\left(x\right)=1$ if $x\geq 0$ and $-1$ if
$x<0$. At this point, it is worth stressing two aspects. First, if one poses
$\tau_{n}=0$, then the Hamiltonian becomes diagonal and equivalent to that of
a non-interacting system. Consistently, the coefficients $D^{n\pm}_{\mu}$
become those of a pure system: $D^{n+}_{1}=-D^{n-}_{2}=\sqrt{2/(N_{a}+1)}$ and
$D^{n-}_{1}=D^{n+}_{2}=0$. If instead one takes the homoatomic limit, i.e.
$\epsilon\rightarrow 0$, then the coefficients become a bonding and
antibonding pair, with $D^{n\pm}_{1}=1/\sqrt{N_{a}+1}$ and
$D^{n\pm}_{n}=\pm\text{sgn}\left(\tau_{n}\right)/\sqrt{N_{a}+1}$.

The last occupied state (TV) $\ket{\tilde{n},-}$ and the first empty state
(BC) $\ket{\tilde{n},+}$ are found at the integer quantum number $\tilde{n}$
that minimizes the quantity $\mathcal{E}_{n}$, i.e. that minimize
$|\tau_{n}|$. If $N_{a}=3m$ or $3m+1$ with $m\in\mathbb{N}^{*}$, then
$\tilde{n}=2m+1$. Note that the interacting term $\tau_{2m+1}$ changes sign in
passing from one family to the other. Instead if $N_{a}=3m-1$, then the
integer $\tilde{n}=2m$ and $\tau_{n}=0$. These considerations leads to the
unperturbed gap of a heteroatomic system ($\epsilon>0$):

$\displaystyle\Delta^{0}_{N_{a}}=\left\\{\begin{array}[]{ll}2\epsilon&\text{for
}N_{a}=3m-1\\\ 2\mathcal{E}_{2m+1}&\text{for the other values of
$N_{a}$}\end{array}\right.$ (7)

and the eigenstates of the TV and BC of the $N_{a}=3m-1$ family are pure
states. The gap of a homoatomic system ($\epsilon=0$) reads:

$\displaystyle\Delta^{0}_{N_{a}}=\left\\{\begin{array}[]{ll}0&\text{for
}N_{a}=3m-1\\\ 2|\tau_{2m+1}|&\text{for the other values of
$N_{a}$}\end{array}\right.$ (10)

and the eigenstates of the TV and BC of the $N_{a}=3m-1$ family are the
bonding and antibonding combinations of $C_{1}$ and $C_{2}$.

###  IV.3 Distinguishing the effects through perturbation theory

As inEzawa (2006); Son  _et al._ (2006), we now add to $H^{0}$ a perturbation
Hamiltonian $\delta H$ which consists in adding $\delta t$ to the hopping term
connecting the atoms of the edge rows ($j=1,N_{a}$) and in changing their on-
side energy by $\delta\epsilon_{\mu}$. The hopping perturbation $\delta t$
accounts for changes in $d_{E}$, so it is more strongly related to the edge
morphology, while the on-site one $\delta\epsilon$ takes into account
variations of $d_{HX}$ and of the passivating species. The perturbative
correction to the energy of the generic state $\ket{n\pm}$ reads

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\braket{n,\pm}{\delta
H}{n,\pm}=2\sin^{2}(\theta_{n})\times$\mbox{}\hfil\\\
&\quad\times\Big{[}(D_{1}^{n\pm})^{2}\delta\epsilon_{1}+(D_{2}^{n\pm})^{2}\delta\epsilon_{2}+2D^{n\pm}_{1}D^{n\pm}_{2}\delta
t\Big{]}\end{split}$ (11)

Figure 8: Deviation from $\Delta^{0}_{N_{a}}$ as a function of $\delta t$
(squares) or $\delta\epsilon$ (circles) in AGNRs (top panels) and ABNNRs
(bottom panels). From left to right: $N_{a}=11,12,13$. Colored symbols result
from equations (12 \- 14); empty symbols from the numerical solution.

In the heteroatomic case $\epsilon>0$, the perturbative correction to the gap
is always $\delta\Delta=\braket{\tilde{n},+}{\delta
H}{\tilde{n},+}-\braket{\tilde{n},-}{\delta H}{\tilde{n},-}$. Using (11), the
coefficients (4) or their appropriate limit, and remembering that
$\Delta^{0}_{N_{a}}=2\mathcal{E}_{\tilde{n}}$, then the gap correction for the
case $\epsilon>0$ reads,

$\delta\Delta=\left(\delta\epsilon_{1}-\delta\epsilon_{2}\right)/m$ (12)

for $N_{a}=3m-1$; and

$\begin{split}\delta\Delta&=\frac{8\sin^{2}\left(\theta_{2m+1}\right)}{(N_{a}+1)\Delta^{0}}\times\\\
&\qquad\times\Big{[}\epsilon\left(\delta\epsilon_{1}-\delta\epsilon_{2}\right)+2\tau_{2m+1}\delta
t\Big{]}\end{split}$ (13)

for $N_{a}=3m$ and $N_{a}=3m+1$. Notice that, by construction, $\tau_{2m+1}$
is the closest to zero among the accessible values, so the term
$2\tau_{2m+1}\delta t$ is always negligible. The result shows that in ABNNRs
the variations of the gap are mostly due to the chemical environment of the
edge atoms. This dependence comes ultimately from an interference between the
TV and the BC wavefunctions. These two states are very close to pure states,
so the mixed products $D_{1}^{+}D_{2}^{+}$ and $D_{1}^{-}D_{2}^{-}$ of
equation (11) are systematically negligible, and they do actually vanish in
the family $N_{a}=3m-1$ where the two states are perfectly pure.

In the homoatomic case ($\epsilon=0$) the corrected gap can be obtained
following the same approach as before, and taking the appropriate limits of
the coefficients (4). However, more attention must be paid in studying the
family $N_{a}=3m-1$. In fact this case corresponds to the double limit
$\epsilon\rightarrow 0$ and $\tau_{n}\rightarrow 0$. Even though the final
eigenvalues do not depend on the order with which the two limits are taken,
the eigenstates do, therefore also the perturbative corrections depend on this
choice. In DFT calculation and experiments, the system itself is well defined
at the very first place, because one works either with ABNNRs or with AGNRs.
So, for comparisons with DFT to make sense, the right order with wich the
limits must be taken is: first $\epsilon\rightarrow 0$, followed by
$\tau_{n}\rightarrow 0$. Finally, one has to pay attention to another point:
in the $N_{a}=3m-1$ family, the TV and the BC states are degenerate and the
unperturbed gap is 0. So there is no reason to define
$\delta\Delta=\braket{\tilde{n},+}{\delta
H}{\tilde{n},+}-\braket{\tilde{n},-}{\delta H}{\tilde{n},-}$ rather than its
opposite. However, the correction must be positive, so the correction must be
defined as the modulus of the difference above. Putting all these things
together, one gets for the homoatomic ($\epsilon=0$) case

$\delta\Delta=\left\\{\begin{array}[]{ll}\frac{2}{m}|\delta t|&\text{for
$N_{a}=3m-1$}\\\
\text{sgn}\left(\tau_{2m+1}\right)\frac{8\sin^{2}\left(\theta_{2m+1}\right)}{(N_{a}+1)}\delta
t&\text{otherwise}\end{array}\right.$ (14)

This result shows that in AGNRs most of the variations of the gap is accounted
by $\delta t$, so by morphological changes of the bonding between edge atoms,
and not by changes of their chemical environment. Once again this result can
be understood from the symmetries of the TV and BC wavefunctions. In fact,
when $\epsilon=0$, the TV and BC states are perfect bonding and antibonding
combinations at any $N_{a}$, so their difference causes the terms in
$(D_{\mu}^{\tilde{n}\pm})^{2}$ of equation (11) to always cancel out. This
result, although in perfect agreement withSon  _et al._ (2006), seems to be in
blatant contradiction with results from 2H-passivated AGNRsLu  _et al._
(2009), where the gap is found independend on the C-C edge distance. Actually,
these systems present a hybridisation of the $sp3$ type and their gapwidth can
not be described by this model.

AGNR
---
$\epsilon_{C}$ | $t_{CC}$ | $\delta t$
0.0 | -2.6 | -0.4

ABNNR
---
$\epsilon_{B}$ | $\epsilon_{N}$ | $t_{BN}$ | $\delta\epsilon$
2.32 | -2.32 | -2.46 | -0.3

Table 1: Parameters of $H=H^{0}+\delta H$ in eV used to plot curves in Figures
2 and 4.

###  IV.4 Validitation of the perturbative approach

Besides the perturbative approach, we also solved the perturbed Hamiltonian
$H=H^{0}+\delta H$ numerically. For the unperturbed problem, we parametrized
the model with values that fit the band structure of the isolated graphene and
hBN monolayers. Instead, the perturbation parameters $\delta\epsilon$ and
$\delta t$ have been adjusted to recover as best as possible the DFT curves
reported in Figures 2 and 4. The best parameters are reported in Table 1.
Successively we explored how the gap changes upon variations of the
perturbative parameters $\delta t$ and $\delta\epsilon_{\mu}$ in the range -1
eV, +1 eV in the nanoribbons of width $N_{a}$=11, 12 and 13, i.e. one
representative nanoribbon per family. Guided by physical intuitions we took
$\delta\epsilon_{1}=\delta\epsilon_{2}=\delta\epsilon$ in the case of AGNRs,
and $\delta\epsilon_{1}=-\delta\epsilon_{2}=\delta\epsilon$ in the case of
ABNNRs.

Globally, the numerical and the perturbative gapwidth are in very good
agreement for both ABNNRs and AGNRs in the range explored, confirming our
conclusions. In all cases, the numerical solution displays a quadratic trend
with respect to $\delta\epsilon$ which adds on top of the invariance (AGNR) or
the linear (ABNNR) dependence predicted by the perturbative approach. The
deviations between the two approaches are larger for this parameter than for
$\delta t$, with the larger deviations of the order of 0.2 eV in the
$N_{a}=3m$ and $N_{a}=3m+1$ families of ABNNRs. Instead, the deviations for
the parameter $\delta t$ are in general very small and never larger than 0.1
eV. Note however that for extreme values of $\delta t$, the numerical solution
may undergo a band crossing in the top valence and the bottom conduction which
would lead to a sudden closing of the gap, as it is the case at $\delta
t=-0.9$ in AGNR13 and $\delta t=0.9$ in AGNR12. This physics is not accessible
in our first order expansion and clearly sets the limit of applicability of
the perturbative approach.

##  V Conclusion

We have calculated with DFT the gapwidth of graphene and boron nitride
armchair nanoribbons (AGNRs and ABNNRs) for ribbon sizes going from $N_{a}=5$
rows to $N_{a}=19$ rows both for relaxed and unrelaxed structures. We have
relaxed selectively specific interatomic distances and reported how the
gapwidth changes upon variations of the bondlength with passivating atoms
(chemistry-driven changes) and between edge atoms (morphology-driven changes).
Thanks to this selective relaxation, we showed that the variations of the
gapwidth in AGNRs are morphology-driven, while in ABNNRs are chemistry-driven.
To understand why, we adopted and extended the tight-binding approach
introduced by Son and coworkers Son  _et al._ (2006) and we demonstrated that
the interference between the wavefunctions of the top valence and the bottom
conduction are at the origin of these two distinct responses.

In the AGNR case, these states are basically a bonding and antibonding pair.
As the two states are equally distributed on the atoms, the difference between
BC and TV leads to a mutual cancellation of on-site changes, and only hopping
terms survive. This explains the stronger dependence of the gapwidth on
interatomic distances and hence on the morphology of the edges rather than the
chemical environment. At variance, in ABNNR case, the TV and the BC states are
basically pure states and the effective Hamiltonian is quasi non-interacting.
As a result, the two states are mostly insensitives to variations in the
hopping term and are instead strongly affected by on-site variations (chemical
environment).

Our results can help pushing further the research on nanoribbon-based devices,
as they clarify the role played by edge-engineering, and selective passivation
and provide the tools to investigate more complex scenarios.

## References

  * Geim and Novoselov (2007) A. K. Geim and K. S. Novoselov, Nature Materials 6, 183 (2007).
  * Watanabe  _et al._ (2009) K. Watanabe, T. Taniguchi, T. Niiyama, K. Miya, and M. Taniguchi, Nature Photonics 3, 591 (2009).
  * Castro Neto  _et al._ (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Reviews of Modern Physics 81, 109 (2009), arXiv:0709.1163 .
  * Weng  _et al._ (2016) Q. Weng, X. Wang, X. Wang, Y. Bando, and D. Golberg, Chemical Society Reviews 45, 3989 (2016).
  * Zhang  _et al._ (2017) K. Zhang, Y. Feng, F. Wang, Z. Yang, and J. Wang, Journal of Materials Chemistry C 5, 11992 (2017).
  * Ezawa (2006) M. Ezawa, Phys. Rev. B 73, 045432 (2006).
  * Fujita  _et al._ (1996) M. Fujita, K. Wakabayashi, K. Nakada, and K. Kusakabe, Journal of the Physical Society of Japan 65, 1920 (1996).
  * Nakada  _et al._ (1996) K. Nakada, M. Fujita, G. Dresselhaus, and M. S. Dresselhaus, Phys. Rev. B 54, 17954 (1996).
  * Behzad and Chegel (2018) S. Behzad and R. Chegel, Diamond and Related Materials 88, 101 (2018).
  * Wakabayashi  _et al._ (1999) K. Wakabayashi, M. Fujita, H. Ajiki, and M. Sigrist, Physical Review B 59, 8271–8282 (1999).
  * Yang  _et al._ (2007) L. Yang, C.-H. Park, Y.-W. Son, M. L. Cohen, and S. G. Louie, Physical Review Letters 99, 186801 (2007), arXiv:0706.1589 .
  * Son  _et al._ (2006) Y.-W. Son, M. L. Cohen, and S. G. Louie, Physical Review Letters 97, 216803 (2006).
  * Nakamura  _et al._ (2005) J. Nakamura, T. Nitta, and A. Natori, Phys. Rev. B 72, 205429 (2005).
  * Park and Louie (2008) C.-H. Park and S. G. Louie, Nano Letters 8, 2200 (2008), pMID: 18593205, https://doi.org/10.1021/nl080695i .
  * Wang  _et al._ (2011) S. Wang, Q. Chen, and J. Wang, Applied Physics Letters 99, 063114 (2011).
  * Shyu (2014) F.-L. Shyu, Physica B: Condensed Matter 452, 7 (2014).
  * Topsakal  _et al._ (2009) M. Topsakal, E. Aktürk, and S. Ciraci, Physical Review B 79, 115442 (2009).
  * Zhang and Guo (2008) Z. Zhang and W. Guo, Physical Review B 77, 075403 (2008).
  * Ding  _et al._ (2009) Y. Ding, Y. Wang, and J. Ni, Applied Physics Letters 94, 233107 (2009).
  * Nematian  _et al._ (2012) H. Nematian, M. Moradinasab, M. Pourfath, M. Fathipour, and H. Kosina, Journal of Applied Physics 111, 093512 (2012).
  * Barone  _et al._ (2006) V. Barone, O. Hod, and G. E. Scuseria, Nano Letters 6, 2748 (2006).
  * Jha  _et al._ (2019) K. K. Jha, N. Tyagi, N. K. Jaiswal, and P. Srivastava, Physics Letters A 383, 125949 (2019).
  * Nishad  _et al._ (2020) V. K. Nishad, A. K. Nishad, S. Roy, B. K. Kaushik, and R. Sharma, Proceedings of the IEEE Conference on Nanotechnology 2020-July, 155 (2020).
  * Prabhakar and Melnik (2019) S. Prabhakar and R. Melnik, Physica E: Low-dimensional Systems and Nanostructures 114, 113648 (2019).
  * Ren  _et al._ (2007) H. Ren, Q. Li, H. Su, Q. W. Shi, J. Chen, and J. Yang, arXiv , 0711.1700 (2007), arXiv:0711.1700 .
  * Lu  _et al._ (2009) Y. Lu, R. Wu, L. Shen, M. Yang, Z. Sha, Y. Cai, and Y. Feng, Applied Physics Letters 94, 122111 (2009).
  * Raza and Kan (2008) H. Raza and E. C. Kan, Phys. Rev. B 77, 245434 (2008).
  * Li  _et al._ (2010) J. Li, L. Z. Sun, and J. X. Zhong, Chinese Physics Letters 27, 077101 (2010).
  * Zheng  _et al._ (2013) P. Zheng, S. E. Bryan, Y. Yang, R. Murali, A. Naeemi, and J. D. Meindl, IEEE Electron Device Letters 34, 707 (2013).
  * Murali  _et al._ (2009) R. Murali, Y. Yang, K. Brenner, T. Beck, and J. D. Meindl, Applied Physics Letters 94, 243114 (2009).
  * Das  _et al._ (2021) S. Das, S. Bhattacharya, D. Das, and H. Rahaman, AIMS Materials Science 8, 247 (2021).
  * Marmolejo-Tejada and Velasco-Medina (2016) J. M. Marmolejo-Tejada and J. Velasco-Medina, Microelectronics Journal 48, 18 (2016).
  * Saraswat  _et al._ (2021) V. Saraswat, R. M. Jacobberger, and M. S. Arnold, ACS Nano 15, 3674 (2021).
  * Osella  _et al._ (2012) S. Osella, A. Narita, M. G. Schwab, Y. Hernandez, X. Feng, K. Mullen, and D. Beljonne, ACS nano 6, 5539 (2012).
  * Mehdi Pour  _et al._ (2017) M. Mehdi Pour, A. Lashkov, A. Radocea, X. Liu, T. Sun, A. Lipatov, R. A. Korlacki, M. Shekhirev, N. R. Aluru, J. W. Lyding, V. Sysoev, and A. Sinitskii, Nature Communications 8, 820 (2017).
  * Schwab  _et al._ (2015) M. G. Schwab, A. Narita, S. Osella, Y. Hu, A. Maghsoumi, A. Mavrinsky, W. Pisula, C. Castiglioni, M. Tommasini, D. Beljonne, X. Feng, and K. Müllen, Chemistry – An Asian Journal 10, 2134 (2015), https://onlinelibrary.wiley.com/doi/pdf/10.1002/asia.201500450 .
  * Wang  _et al._ (2016) G. Wang, S. Wu, T. Zhang, P. Chen, X. Lu, S. Wang, D. Wang, K. Watanabe, T. Taniguchi, D. Shi, R. Yang, and G. Zhang, Applied Physics Letters 109, 053101 (2016), https://doi.org/10.1063/1.4959963 .
  * Dauber  _et al._ (2014) J. Dauber, B. Terrés, C. Volk, S. Trellenkamp, and C. Stampfer, Applied Physics Letters 104, 083105 (2014), https://doi.org/10.1063/1.4866289 .
  * Huang  _et al._ (2009) B. Huang, M. Liu, N. Su, J. Wu, W. Duan, B.-l. Gu, and F. Liu, Phys. Rev. Lett. 102, 166404 (2009).
  * Huang  _et al._ (2012) B. Huang, H. Lee, B. L. Gu, F. Liu, and W. Duan, Nano Research 5, 62 (2012).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, _et al._ , Journal of physics: Condensed matter 21, 395502 (2009).
  * Grimme (2006) S. Grimme, Journal of Computational Chemistry 27, 1787 (2006), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.20495 .
  * Hamann (2013) D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * Posternak  _et al._ (1983) M. Posternak, A. Baldereschi, A. J. Freeman, E. Wimmer, and M. Weinert, Phys. Rev. Lett. 50, 761 (1983).
  * Posternak  _et al._ (1984) M. Posternak, A. Baldereschi, A. J. Freeman, and E. Wimmer, Phys. Rev. Lett. 52, 863 (1984).
  * Blase  _et al._ (1995) X. Blase, A. Rubio, S. G. Louie, and M. L. Cohen, Physical Review B 51, 6868 (1995).
  * Paleari (2019) F. Paleari,  _First-principles approaches to the description of indirect absorption and luminescence spectroscopy: exciton-phonon coupling in hexagonal boron nitride_ , Ph.D. thesis, University of Luxembourg (2019).
  * Latil  _et al._ (2022) S. Latil, H. Amara, and L. Sponza, “Structural classification of boron nitride twisted bilayers and ab initio investigation of their stacking-dependent electronic structure,”  (2022), accepted in SciPost Physics.
  * Blase  _et al._ (1994a) X. Blase, L. X. Benedict, E. L. Shirley, and S. G. Louie, Phys. Rev. Lett. 72, 1878 (1994a).
  * Blase  _et al._ (1994b) X. Blase, A. Rubio, S. G. Louie, and M. L. Cohen, Europhysics Letters (EPL) 28, 335 (1994b).
  * Hu  _et al._ (2010) S. Hu, J. Zhao, Y. Jin, J. Yang, H. Petek, and J. G. Hou, Nano Letters 10, 4830 (2010), pMID: 21049977, https://doi.org/10.1021/nl1023854 .