# Witten deformation and divergence-free symmetric Killing 2-tensors

Kwangho Choi and Junho Lee

###### Abstract

Using a Morse function and a Witten deformation argument, we obtain an upper
bound for the dimension of the space of divergence-free symmetric Killing
$p$-tensors on a closed Riemannian manifold, and calculate it explicitly for
$p=2$.

##  1 Introduction

Let $(M,g)$ be a connected closed Riemannian $n$-manifold and let ${\rm
Sym}^{p}(T^{*}M)$ denote the $p$-th symmetric product of the cotangent bundle
$T^{*}M$. Denote the space of smooth sections of ${\rm Sym}^{p}(T^{*}M)$ by

$S^{p}(M)=\Gamma\big{(}{\rm Sym}^{p}(T^{*}M)\big{)}.$

Here the symmetric product is given by

$\alpha_{1}\cdots\alpha_{p}=\sum_{\sigma\in\mathfrak{S}_{p}}\alpha_{\sigma(1)}\otimes\cdots\otimes\alpha_{\sigma(p)},$

where $\alpha_{1},\cdots,\alpha_{p}\in S^{1}(M)=\Omega^{1}(M)$ and
$\mathfrak{S}_{p}$ denotes the symmetric group on $p$ letters. For simplicity,
we will often write $\alpha\beta$ instead of $\alpha\cdot\beta$.

Let $\nabla$ be the Levi-Civita connection on $TM$, and denote its induced
connection on $S^{p}(T^{*}M)$ by the same $\nabla$. Then the symmetric
covariant derivative is the first order operator

$d^{s}:S^{p}(M)\to S^{p+1}(M):\alpha\mapsto\sum_{k}e^{k}\nabla_{e_{k}}\alpha,$
(1.1)

where $\\{e_{k}\\}$ is an orthonormal frame on $TM$ and its dual coframe
$\\{e^{k}\\}$.

The Riemannian metric on $M$ induces a pointwise inner product on $S^{p}(M)$
by the formula

$\langle\alpha_{1}\cdots\alpha_{p},\beta_{1}\cdots\beta_{p}\rangle=\sum_{\sigma\in\mathfrak{S}_{p}}\langle\alpha_{1},\beta_{\sigma(1)}\rangle\cdots\langle\alpha_{p},\beta_{\sigma(p)}\rangle,$

where $\alpha_{i},\beta_{j}\in S^{1}(M)$. Denote its induced $L^{2}$-inner
product on $S^{p}(M)$ by

$(\alpha,\beta)=\int_{M}\langle\alpha,\beta\rangle.$

Then the formal adjoint of $d^{s}$ is the divergence operator

$d^{s*}:S^{p+1}(M)\to
S^{p}(M):\alpha\mapsto-\sum_{k}e_{k}\mathbin{\lrcorner}\nabla_{e_{k}}\alpha$
(1.2)

(cf. Theorem 3.3 of [8]).

A symmetric tensor $\alpha\in S^{p}(M)$ is called Killing tensor if
$d^{s}\alpha=0$, and is divergence-free if $d^{s*}\alpha=0$. Consider the
following self-adjoint second order operator

$D=d^{s*}d^{s}+d^{s}d^{s*}:S^{p}(M)\to S^{p}(M).$ (1.3)

It is an elliptic operator on $S^{p}(M)$ (see Lemma 2.3) such that

$\mathscr{D}^{p}(M,g):=\ker D=\ker d^{s}\cap\ker d^{s*},$

hence $\mathscr{D}^{p}(M,g)$ is the space of divergence-free symmetric Killing
$p$-tensors.

There is a well-known dimension formula for symmetric Killing $p$-tensors:

$\dim\ker
d^{s}\leq\frac{1}{n}\binom{n+p}{p+1}\mbox{\small$\displaystyle\binom{n+p-1}{p}$}$
(1.4)

(see [1, 5, 11, 19, 21]). The equality holds when $M$ is the $n$-sphere
$S^{n}$ with the standard metric. When $p=1$, Killing 1-forms are traceless,
hence divergence-free (cf. Section 2 of [7]). These are dual to Killing vector
fields, and (1.4) specializes to

$\dim\mathscr{D}^{1}(M,g)\leq\binom{n+1}{2}.$ (1.5)

This is the well-known bound on the dimension of the Lie algebra of the space
of infinitesimal isometries of a Riemannian $n$-manifold.

The aim of this paper is to find a dimension bound, analogous to (1.4), for
the space of divergence-free symmetric Killing $p$-tensors. Our approach is to
use a Morse function $f$ to find a deformation $D_{t}$ of the operator $D$
satisfying the following properties:

  * (P1)

$\dim\ker D_{t}=\dim\ker D$ for all $t$.

  * (P2)

$\dim\ker D_{t}\leq\sum_{i}\dim\ker\bar{K}(x_{i})$ for sufficiently large $t$,
where the sum is over the critical points $\\{x_{i}\\}$ of $f$ and
$\bar{K}(x_{i})$ is a model operator (see Section 3.2.2).

  * (P3)

$\dim\ker\bar{K}(x_{i})$ can be calculated.

For (P1), we follow Witten’s method in [22]; for (P2), we use semiclassical
analysis of model operators in [17] (see also [16]); for (P3), we make the
technical assumption that the metric $g$ is adapted to $f$, that is,
$g_{ij}=\delta_{ij}$ in a Morse chart around each critical point.

Witten deformations and semiclassical analysis are extensively studied and
well-understood (see, for example, [4], [18], [23], and references therein).
However, it is generally difficult to explicitly calculate the kernel or
spectrum of the model operators $\bar{K}(x_{i})$. Thus for a deformation
$D_{t}$ satisfying (P1) and (P2), the main task is to calculate
$\dim\ker\bar{K}(x_{i})$.

######  Theorem 1.1.

Let $(M,g)$ be a connected closed $n$-manifold and $f$ be a Morse function
such that the metric $g$ is adapted to $f$. Then,

$\dim\mathscr{D}^{p}(M,g)\leq\sum\dim\ker\bar{K}(x_{i}),$ (1.6)

where the sum is over all critical points $x_{i}$ of $f$. Furthermore, for the
critical point $x$ of $f$ with index $m$, the integers

$K^{n,p}_{m}=\dim\ker\bar{K}(x)$

depend only on $n$, $p$, and $m$, and satisfy

$\begin{array}[]{l}(a)\ K^{n,p}_{m}=0\ \text{if $m\geq p$ or if
$m=n$}.\\\\[10.0pt] (b)\
K^{n,1}_{0}=\mbox{\small$\displaystyle\binom{n+1}{2}$}\ \ \ (c)\
K^{n,2}_{0}={\small\left\\{\begin{array}[]{cl}1&\text{if}\ n=1\\\ 3&\text{if}\
n=2\\\
\mbox{\small$\displaystyle\frac{1}{n}$}\mbox{\small$\displaystyle\binom{n+2}{3}$}\mbox{\small$\displaystyle\binom{n+1}{2}$}-\mbox{\small$\displaystyle\frac{n(n+3)}{2}$}&\text{if}\
n\geq 3\end{array}\right.}\end{array}$

The inequality (1.6) is a consequence of (P1) and (P2), which are proved in
Proposition 3.2 (a) and Proposition 3.3, respectively, and the integers
$K^{n,p}_{m}$ depend only on $n$, $p$, and $m$ by (3.5). When $p=1$, this
theorem implies the inequality (1.5) under the adapted metric assumption
because every connected closed manifold admits a Morse function with a unique
critical point of index zero. When $p=2$, Theorem 1.1 shows:

######  Corollary 1.2.

If $M$ admits a Morse function with a unique minimum and no critical point of
index $1$ then, for an adapted metric,

$\dim\mathscr{D}^{2}(M,g)\leq K^{n,2}_{0},$

where $K^{n,2}_{0}$ is given above.

A simply connected closed manifold admitting a perfect Morse function
satisfies the condition of Corollary 1.2. Examples of such a manifold include
$S^{n}$ ($n\geq 2$), $\mathbb{C}\mathbb{P}^{n}$, complex Grassmannian
manifolds and their products (cf. [12]).

At the moment, the authors do not know whether the adapted metric assumption
can be eliminated and whether $K^{n,p}_{m}=0$ when $m\neq 0$

The rest of the paper is organized as follows. Section 2 reviews basic
operators on ${\rm Sym}^{p}(TM)$ and calculates the symbol of $D$. In Section
3, we define a deformation $D_{t}$ satisfying (P1) and find its model
operators $\bar{K}(x)$ satisfying (P2). We then employ Hermite polynomials and
operators satisfying the canonical commutation relations to prove Theorem 1.1
(a) in Section 4, and to compute $K^{n,1}_{0}$ and $K^{n,2}_{0}$ in Section 5;
this completes the proof of Theorem 1.1, hence shows (P3).

Acknowledgement: The authors sincerely thank Thomas H. Parker for his
insightful discussion and helpful comments on this paper. The second author
also thanks Seoul National University and the Korea Institute for Advanced
Study for their support and hospitality during his visit when work on this
paper was undertaken.

##  2 Preliminaries

This section reviews the canonical commutation relations (2.1) below, fixes
notations, and calculates the symbol of the operator $D$ defined in (1.3).

###  2.1 Canonical commutation relations

######  Definition 2.1.

For each 1-form $\gamma\in S^{1}(M)=\Omega^{1}(M)$, define operators

$\displaystyle c_{\gamma}:S^{p}(M)$ $\displaystyle\to S^{p+1}(M)\ \ \
\text{by}\ \ \ c_{\gamma}(\alpha)=\gamma\alpha,$ $\displaystyle
a_{\gamma}:S^{p}(M)$ $\displaystyle\to S^{p-1}(M)\ \ \ \text{by}\ \ \
a_{\gamma}(\alpha)=\gamma^{*}\\!\mathbin{\lrcorner}\alpha,$

where $\gamma^{*}$ is the metric dual of $\gamma$, i.e.,
$\gamma(v)=g(\gamma^{*},v)$.

Let $\\{e^{k}\\}$ be an orthonormal frame on $T^{*}M$. If we regard $S^{p}(M)$
as the space of polynomials of degree $p$ in the variables
$e^{1},\cdots,e^{n}$, then $c_{\gamma}$ is the multiplication by $\gamma$ and
$a_{\gamma}$ is the directional derivative in the direction of $\gamma$. For
instance, if $\gamma=\sum_{k}\gamma_{k}e^{k}$, $\alpha=(e^{1})^{3}e^{2}$, and
$\beta=(e^{1})^{2}e^{2}$ then

$\langle a_{\gamma}\alpha,\beta\rangle=\langle
3\gamma_{1}(e^{1})^{2}e^{2}+\gamma_{2}(e^{1})^{3},\beta\rangle=6\gamma_{1}=\langle\alpha,c_{\gamma}\beta\rangle.$

The following well-known fact plays an essential role in our discussions.

######  Lemma 2.2.

For any $\gamma,\mu\in S^{1}(M)$, we have:

  * (a)

The adjoint of $a_{\gamma}$ is $a_{\gamma}^{*}=c_{\gamma}$.

  * (b)

The operators $c_{\gamma}$ and $a_{\gamma}$ satisfy the canonical commutation
relations

$c_{\gamma}c_{\mu}=c_{\mu}c_{\gamma},\ \ \ \
a_{\gamma}a_{\mu}=a_{\mu}a_{\gamma},\ \ \ \
c_{\gamma}a_{\mu}-a_{\mu}c_{\gamma}=-\langle\gamma,\mu\rangle Id.$ (2.1)

Let $\nabla$ denote the connection on ${\rm Sym}^{p}(T^{*}M)$ induced from the
Levi-Civita connection. Given an orthonormal frame $\\{e_{k}\\}$ on $TM$ and
its dual coframe $\\{e^{k}\\}$, we set

$\nabla_{k}=\nabla_{e_{k}},\ \ \ \ c_{k}=c_{e^{k}},\ \ \ \ a_{k}=a_{e^{k}}.$

Then, we can write the symmetric covariant derivative $d^{s}$ in (1.1) and the
divergence operator $d^{s*}$ in (1.2) as

$d^{s}=\sum_{k}c_{k}\nabla_{k}\ \ \ \ \text{and}\ \ \ \
d^{s*}=-\sum_{k}a_{k}\nabla_{k}.$

###  2.2 The symbol map

Let $E$ and $F$ be vector bundles over $M$. Recall that the (principal) symbol
of a $k^{th}$-order liner differential operator $D:\Gamma(E)\to\Gamma(F)$ is
the section

$\sigma_{D}\in\Gamma\big{(}{\rm Hom}(\pi^{*}E,\pi^{*}F)\big{)}\ \ \ \
\text{(where $\pi:T^{*}M\to M$)}$

such that for $(x,w)\in T^{*}M$, the symbol map $\sigma_{D}(x,w):E_{x}\to
F_{x}$ is given by

$\sigma_{D}(x,w)(\xi)=\mbox{\small$\displaystyle\frac{1}{k!}$}D(f^{k}\xi)(x),$

where $f$ is a function satisfying $f(x)=0$ and $df(x)=w$. A linear
differential operator $D$ is called elliptic if the symbol map
$\sigma_{D}(x,w)$ is an isomorphism for every $(x,w)\in
T^{*}M\setminus\\{0\\}$.

######  Lemma 2.3.

  * (a)

$d^{*s}d^{s}$ is elliptic, hence $\ker d^{s}=\ker d^{s*}d^{s}$ has a finite
dimension.

  * (b)

$d^{s*}d^{s}\pm d^{s}d^{s*}$ are elliptic.

Proof. The symbols of $d^{s}$, $d^{s*}$, $d^{s*}d^{s}$, and $d^{s}d^{s*}$ are
given as

$\begin{array}[]{ll}\sigma_{d^{s}}(x,w)=c_{w},&\sigma_{d^{s^{*}}}(x,w)=-a_{w},\\\
\sigma_{d^{s*}d^{s}}(x,w)=-a_{w}c_{w},&\sigma_{d^{s}d^{s*}}(x,w)=-c_{w}a_{w},\end{array}$
(2.2)

where $(x,w)\in T^{*}M\setminus\\{0\\}$. Since $c_{w}$ is injective and
$a_{w}^{*}=c_{w}$, the symbol $-a_{w}c_{w}$ is also injective, hence an
isomorphism. This shows (a). For (b), first note that

$c_{w}a_{w}-a_{w}c_{w}=-|w|^{2}Id$

by (2.1). Also, note that $c_{w}a_{w}+a_{w}c_{w}$ is an isomorphism because

$\big{\langle}(a_{w}c_{w}+c_{w}a_{w})\alpha,\alpha\big{\rangle}=\langle
c_{w}\alpha,c_{w}\alpha\rangle+\langle a_{w}\alpha,a_{w}\alpha\rangle$

and $c_{w}$ is injective. Therefore, (b) follows. ∎

######  Remark 2.4.

In general, the divergence operator $d^{s*}$ has an infinite-dimensional
kernel since its symbol is given by contraction $a_{w}$. One can find such an
example in Section 3.3 of [3] in the context of symplectic Dirac operators.

######  Remark 2.5.

The elliptic operator $d^{s*}d^{s}-d^{s}d^{s*}$ on symmetric tensors, first
studied by Sampson [15], is a zeroth-order perturbation of the trace
Laplacian; hence it is a symmetric analog of the Hodge Laplacian on
differential forms.

##  3 Witten deformation

This section applies Witten’s method in [22] to perturb the operator
$D=d^{s*}d^{s}+d^{s}d^{s*}$ in (1.3) via a Morse function $f$ and then uses
Corollary 1.3 of [17] to obtain an upper bound for $\dim\ker D$ given by local
data (or model operators) at the critical points of $f$.

###  3.1 Perturbation

We will choose a deformation $D_{t}$ of $D$ satisfying (P1) and (P2) in the
Introduction as follows. Fix a smooth function $f$ on $M$ and for
$t\in\mathbb{R}$, define

$d_{t}:S^{p}(M)\to S^{p+1}(M)\ \ \ \text{by}\ \ \ d_{t}=e^{-tf}d^{s}e^{tf}.$

######  Definition 3.1.

For $t\in\mathbb{R}$, let $D_{t}=d_{t}^{*}d_{t}+d_{-t}d_{-t}^{*}$.

Since $d_{t}=d^{s}+tc_{df}$ and
$d_{t}^{*}=e^{tf}d^{s*}e^{-tf}=d^{s*}+ta_{df}$, we have

$\displaystyle D_{t}=D+tB+t^{2}V,\ \ \ \text{where}$ (3.1) $\displaystyle
B=d^{s*}c_{df}-c_{df}d^{s*}+a_{df}d^{s}-d^{s}a_{df},\ \ \
V=a_{df}c_{df}+c_{df}a_{df}.$

Note that $D_{t}$, $D$, $B$, and $V$ are all self-adjoint. The proposition
below shows (P1) and is crucial for obtaining (P2).

######  Proposition 3.2.

  * (a)

$\dim\ker D_{t}=\dim\ker D$ for all $t\in\mathbb{R}$.

  * (b)

$B$ is a bundle map, i.e., its symbol $\sigma_{B}=0$.

  * (c)

$V=2a_{df}^{*}a_{df}+|df|^{2}Id$, hence $\langle
V\alpha,\alpha\rangle\geq|df|^{2}|\alpha|^{2}$, $\forall\alpha\in S^{p}(M)$.

Proof. Define an isomorphism $L_{t}:S^{p}(M)\to S^{p}(M)$ by

$L_{t}(\alpha)=e^{-tf}\alpha.$

Then, for $d^{s},d_{t}:S^{p}(M)\to S^{p+1}(M)$ and
$d^{s*},d_{-t}^{*}:S^{p}(M)\to S^{p-1}(M)$,

$L_{t}(\ker d^{s})=\ker d_{t}\ \ \ \ \text{and}\ \ \ \ L_{t}(\ker d^{s*})=\ker
d^{*}_{-t}.$

Now, (a) follows from

$\ker D_{t}=\ker d_{t}\cap\ker d_{-t}^{*}=L_{t}\big{(}\ker d^{s}\cap\ker
d^{s*}\big{)}=L_{t}\big{(}\ker D\big{)}.$

(b) follows from (2.2) and (2.1):

$\sigma_{B}(x,w)=-a_{w}c_{df}+c_{df}a_{w}+a_{df}c_{w}-c_{w}a_{df}=-\langle
df,w\rangle Id+\langle w,df\rangle Id=0.$

Lastly, (c) follows from Lemma 2.2. ∎

Because $B$ is a bundle map, $D_{t}$ is a zeroth-order perturbation of $D$
such that, as $\langle D\alpha,\alpha\rangle\geq 0$ and $\langle
V\alpha,\alpha\rangle\geq|df|^{2}|\alpha|^{2}$, elements of $\ker D_{t}$
concentrate at critical points of $f$ as $t\to\infty$ (cf. Lemma 9.2 of [9] or
Proposition 2.4 of [10]). For further examples of concentration by perturbing
elliptic operators, see Section 7 of [20], Section 9 of [9], [14], and [10].

###  3.2 Local data

In this subsection, we compute model operators under the adapted metric
assumption below and apply Corollary 1.3 of [17] to show (P2).

####  3.2.1 Adapted metric

Suppose $f$ is a Morse function and $x$ is a critical point of index $m$.
Then, by Morse Lemma, there exists local coordinates $\\{y_{k}\\}$ around $x$
with $y_{k}(x)=0$ such that

$f(y)=f(0)-\mbox{\small$\displaystyle\frac{1}{2}$}y_{1}^{2}-\cdots-\mbox{\small$\displaystyle\frac{1}{2}$}y_{m}^{2}+\mbox{\small$\displaystyle\frac{1}{2}$}y_{m+1}^{2}+\cdots+\mbox{\small$\displaystyle\frac{1}{2}$}y_{n}^{2}.$

Now assume that the Riemannian metric $g$ on $M$ is adapted to $f$, that is,
the metric $g$ is given by

$g=(dy_{1})^{2}+\cdots+(dy_{n})^{2}.$

(cf. [12]). Such a metric always exists and the 1-forms $e^{k}=dy_{k}$ form an
orthonormal frame of $T^{*}M$.

Let $\alpha=\sum_{J}\alpha_{J}e^{J}$, where $J=(j_{1},\cdots,j_{n})$ is a
multi-index, that is, $e^{J}=(e^{1})^{j_{1}}\cdots(e^{n})^{j_{n}}.$ Then, in
the local coordinates $\\{y_{k}\\}$,

$\displaystyle
d^{s}\alpha=\sum_{J}\sum_{k}\mbox{\small$\displaystyle\frac{\partial\alpha_{J}}{\partial
y_{k}}$}c_{k}e^{J},\ \ \ \
d^{s*}\alpha=-\sum_{J}\sum_{k}\mbox{\small$\displaystyle\frac{\partial\alpha_{J}}{\partial
y_{k}}$}a_{k}e^{J},$ (3.2) $\displaystyle
df=-y^{1}e^{1}\cdots-y^{m}e^{m}+y^{m+1}e^{m+1}+\cdots+y^{n}e^{n},$
$\displaystyle\nabla_{k}df=\sum_{k}s_{k}e^{k}\ \ \ \text{where}\ \ \
s_{k}=\left\\{\begin{array}[]{rl}-1&k\leq m,\\\ 1&k>m\end{array}\right.$
$\displaystyle
V(y)=\sum_{k,\ell}(s_{k}s_{\ell})y_{k}y_{\ell}(a_{k}c_{\ell}+c_{k}a_{\ell}).$

Observe that by Proposition 3.2 (c), we have

$V(0)=0\ \ \ \ \text{and}\ \ \ \ \langle
V(y)\alpha,\alpha\rangle\geq|y|^{2}|\alpha|^{2}.$ (3.3)

####  3.2.2 Model operators

Recall that the rank of ${\rm Sym}^{p}(T^{*}M)$ is the dimension of the space
of homogeneous polynomials of degree $p$ on $\mathbb{R}^{n}$, which is

$n_{p}:=\binom{n+p-1}{p}.$ (3.4)

In the local coordinates $\\{y_{k}\\}$ and a trivialization of ${\rm
Sym}^{p}(T^{*}M)$ given by the local frame $\\{e^{k}=dy_{k}\\}$, the operator
$D_{t}=D+tB+t^{2}V$ becomes an $n_{p}\times n_{p}$ matrix differential
operator in $\mathbb{R}^{n}$ and $D$ becomes a matrix differential operator of
the form

$D=\sum_{k,\ell}\frac{\partial^{2}}{\partial y_{k}\partial
y_{\ell}}D_{k\ell}(y)+\text{lower order terms}.$

The model operator $\bar{K}(x)$ at the critical point $x$ is the sum of
operators

$\bar{K}(x)=\bar{D}+\bar{B}+\bar{V},$

where the operators $\bar{D}$, $\bar{B}$, and $\bar{V}$ are defined as:

  * •

${\displaystyle\bar{D}=\sum_{k,\ell}\mbox{\small$\displaystyle\frac{\partial^{2}}{\partial
y_{k}\partial y_{\ell}}$}D_{k\ell}(0)}$ is the second order homogeneous matrix
differential operator obtained by taking the second order terms of $D$ with
the coefficients frozen at the point $x$.

  * •

$\bar{B}=B(0)$ is an endomorphism of the fiber ${\rm Sym}^{p}(T^{*}_{x}M)$.

  * •

${\displaystyle\bar{V}=\mbox{\small$\displaystyle\frac{1}{2}$}\sum_{k,\ell}\mbox{\small$\displaystyle\frac{\partial^{2}V}{\partial
y_{k}\partial y_{\ell}}$}(0)y_{k}y_{\ell}}$ is the quadratic part of $V$ near
$x$.

(See [17] and also [16].) In our case, by (2.1) and (3.2),

$\displaystyle\bar{D}$
$\displaystyle=-\sum_{k,\ell}\partial_{k}\partial_{\ell}(a_{k}c_{\ell}+c_{k}a_{\ell}),\
\ \ \ \bar{B}=-\sum_{k}s_{k}(a_{k}c_{k}+c_{k}a_{k}),$ (3.5)
$\displaystyle\bar{V}$
$\displaystyle=\sum_{k,\ell}(s_{k}s_{\ell})y_{k}y_{\ell}(a_{k}c_{\ell}+c_{k}a_{\ell}).$

The model operator $\bar{K}(x)$ is a self-adjoint second order elliptic
operator acting in $L^{2}(\mathbb{R}^{n})\otimes{\rm Sym}^{p}(T^{*}_{x}M)$, as
a closure from the domain $C_{0}^{\infty}(\mathbb{R}^{n})\otimes{\rm
Sym}^{p}(T^{*}_{x}M)$. Moreover, $\bar{K}(x)$ is positive semi-definite (see
Remark 4.7 below). Observe that (3.5) implies the dimension
$K^{n,p}_{m}=\dim\ker\bar{K}(x)$ depends only on $n$, $p$, and $m$,
independent of the critical point $x$.

####  3.2.3 Semiclassical analysis

Let $x_{1},\cdots,x_{N}$ be the critical points of the Morse function $f$.
Then the model operator for $D_{t}$ on $M$ is the elliptic operator

$K=\bigoplus_{i=1}^{N}\bar{K}(x_{i})$

acting in the orthogonal direct sum
$\bigoplus_{i}\big{(}L^{2}(\mathbb{R}^{n})\otimes
S^{p}(T^{*}_{x_{i}}M)\big{)}$. Since $\langle D\alpha,\alpha\rangle\geq 0$ and
$B$ is a bundle map, (3.3) implies that as $t\to\infty$, the eigenvalues
(multiplicities counted) of the operator $t^{-1}D_{t}$ concentrate near the
eigenvalues of the model operator $K$ (see Proposition 1.2 of [17]). As a
corollary, we obtain (P2) in the Introduction as follows.

######  Proposition 3.3 (Corollary 1.3 of [17]).

There exists $t_{0}>0$ such that for any $t>t_{0}$, we have

$\dim\ker D_{t}=\dim\ker t^{-1}D_{t}\leq\dim\ker
K=\sum_{i}\dim\ker\bar{K}(x_{i}).$

##  4 Hermite polynomials and operator formalism

Throughout the remaining sections, let $x$ denote a critical point of the
Morse function $f$ with Morse index $m$ as in Section 3.2.1. For notational
simplicity, we set

$\bar{K}=\bar{K}(x)\ \ \ \ \text{and}\ \ \ \ S^{p}={\rm Sym}^{p}(T^{*}_{x}M).$

As the model operator $\bar{K}$ acts in $L^{2}(\mathbb{R}^{n})\otimes S^{p}$,
Hermite functions on $\mathbb{R}^{n}$ are useful for calculating $\ker\bar{K}$
because they constitute an orthogonal basis for $L^{2}(\mathbb{R}^{n})$ (cf.
[6]). However, we will use Hermite polynomials on $\mathbb{R}^{n}$ and
relevant weighted operators for computational convenience. In this section, we
write weighted operators in terms of operators convenient for handling Hermite
polynomials and use this operator formalism to prove Theorem 1.1 (a).

Recall that for a multi-index $I=(i_{1},\cdots,i_{n})$,

$I!=i_{1}!\cdots i_{n}!\ \ \ \ \text{and}\ \ \ \ |I|={\textstyle\sum}i_{k}.$

Also recall that $\\{e^{J}\,|\,e^{k}=dy_{k},|J|=p\\}$ is an orthogonal basis
for $S^{p}$. Let $1_{k}$ denote the multi-index with $i$-th component
$\delta_{ik}$. Note that $I=\sum i_{k}1_{k}$.

###  4.1 Hermite polynomials on $\mathbb{R}^{n}$

The classical Hermite polynomials on $\mathbb{R}$ are given by

$H_{k}(t)=(-1)^{k}e^{t^{2}}\mbox{\small$\displaystyle\frac{d^{k}}{dt^{k}}$}e^{-t^{2}}.$

Here the degree of $H_{k}(t)$ is $k$. These Hermite polynomials on
$\mathbb{R}$ satisfy

$H_{k}^{\prime}(t)=2kH_{k-1}\ \ \ \ \text{and}\ \ \ \
tH_{k}(t)=kH_{k-1}(t)+\mbox{\small$\displaystyle\frac{1}{2}$}H_{k+1}.$ (4.1)

######  Definition 4.1.

  * (a)

The Hermite polynomials on $\mathbb{R}^{n}$ are products of Hermite
polynomials on $\mathbb{R}$, that is,

$H_{I}(y)=H_{i_{1}}(y_{1})\cdots H_{i_{n}}(y_{n}),$

where $I\in\mathbb{Z}^{n}$. If some $i_{k}<0$, we set $H_{I}=0$.

  * (b)

Denote the space generated by degree $q$ Hermite polynomials by

${\cal H}^{q}={\rm span}\\{H_{I}:|I|=q\\}.$

Below are well-known properties of the Hermite polynomials $H_{I}$ that are
frequently used for our subsequent discussions.

  * •

For $d\mu=e^{-|y|^{2}}dy$, $\\{H_{I}\\}$ is an orthogonal basis for
$L^{2}(\mathbb{R}^{n},d\mu)$; hence the vectors $H_{I}e^{J}$ form an
orthogonal basis for $L^{2}(\mathbb{R}^{n},d\mu)\otimes S^{p}$ such that

$\big{(}H_{I}e^{J},H_{I^{\prime}}e^{J^{\prime}}\big{)}=\int_{\mathbb{R}^{n}}H_{I}H_{I^{\prime}}d\mu\,\langle
e^{J},e^{J^{\prime}}\rangle=\delta_{I,I^{\prime}}\delta_{J,J^{\prime}}\big{(}\sqrt{\pi}\big{)}^{n}2^{|I|}I!J!.$
(4.2)

  * •

The dimension of the space ${\cal H}^{q}$ is the number of multi-indices $I$
with $|I|=q$, so $\dim{\cal H}^{q}=n_{q}$, where $n_{q}$ is given in (3.4).

  * •

By (4.1), the Hermite polynomials $H_{I}$ satisfy

$\partial_{k}H_{I}=2i_{k}H_{I-1_{k}}\ \ \ \ \text{and}\ \ \ \
y_{k}H_{I}=i_{k}H_{I-1_{k}}+\mbox{\small$\displaystyle\frac{1}{2}$}H_{I+1_{k}}.$
(4.3)

###  4.2 Weighted operators

Noting that $\phi\in L^{2}(\mathbb{R}^{n},d\mu)$ if and only if
$e^{-\frac{1}{2}|y|^{2}}\phi\in L^{2}(\mathbb{R}^{n})$, we define:

######  Definition 4.2.

For the exponential weight function $w=e^{-\frac{1}{2}|y|^{2}}$, let

$\bar{K}_{w}=w^{-1}\bar{K}w,\ \text{that is},\
K_{w}(\beta)=w^{-1}\bar{K}(w\beta).$

The weighted operator $\bar{K}_{w}$ acts in $L^{2}(\mathbb{R}^{n},d\mu)\otimes
S^{p}$ such that

$\dim\ker\bar{K}_{w}=\dim\ker\bar{K}\ \ \ \ \text{and}\ \ \ \
\bar{K}_{w}=w^{-1}\bar{D}w+\bar{B}+\bar{V},$

where by a simple calculation

$w^{-1}\bar{D}w=-\sum_{k,\ell}(-\delta_{k\ell}+y_{k}y_{\ell}-y_{k}\partial_{\ell}-y_{\ell}\partial_{k}+\partial_{k}\partial_{\ell})(a_{k}c_{\ell}+c_{k}a_{\ell}).$
(4.4)

###  4.3 Canonical commutation relations

Without further confusion, we will often regard $y_{k}$ as an operator

$y_{k}:g(y)\mapsto y_{k}g(y).$

Noting (4.3), we obtain handy operators as follows.

######  Definition 4.3.

For $k=1,\cdots,n$, set

$C_{k}=2y_{k}-\partial_{k}\ \ \ \ \text{and}\ \ \ \
A_{k}=\mbox{\small$\displaystyle\frac{1}{2}$}\partial_{k}.$

$C_{k}$ and $A_{k}$ are analogous to $c_{k}$ and $a_{k}$ in the following
sense. Recall that

$\displaystyle c_{k}:S^{p}$ $\displaystyle\to S^{p+1}:e^{J}\mapsto
e^{J+1_{k}},$ (4.5) $\displaystyle a_{k}:S^{p}$ $\displaystyle\to
S^{p-1}:e^{J}\mapsto j_{k}e^{J-1_{k}},$

Analogously, when restricted to ${\cal H}^{q}$, by (4.3) $C_{k}$ and $A_{k}$
are given as

$\displaystyle C_{k}:{\cal H}^{q}\to{\cal H}^{q+1}:H_{I}\mapsto H_{I+1_{k}},$
(4.6) $\displaystyle A_{k}:{\cal H}^{q}\to{\cal H}^{q-1}:H_{I}\mapsto
i_{k}H_{I-1_{k}}.$

Furthermore, one has:

######  Lemma 4.4.

  * (a)

$A_{k}^{*}=\frac{1}{2}C_{k}$ with respect to the inner product of
$L^{2}(\mathbb{R}^{n},d\mu)$.

  * (b)

The operators $C_{k}$ and $A_{k}$ satisfy the canonical commutation relations:

$C_{k}C_{\ell}=C_{\ell}C_{k},\ \ \ \ A_{k}A_{\ell}=A_{\ell}A_{k},\ \ \ \
C_{k}A_{\ell}-A_{\ell}C_{k}=-\delta_{k\ell}I.$ (4.7)

The proof is immediate from definitions.

###  4.4 Operator formalism

######  Definition 4.5.

  * (a)

Recalling $s_{k}$ in (3.2) is defined by the Morse index $m$, we set

$(Cc)^{-}=\sum_{k\leq m}C_{k}c_{k},\ \ \ \
(Cc)^{+}=\sum_{m<\ell}C_{\ell}c_{\ell},\ \ \ \ (Ca)^{-}=\sum_{k\leq
m}C_{k}a_{k},$

and likewise for $(Ca)^{+},(Ac)^{-},(Ac)^{+},(Aa)^{-},(Aa)^{+}$.

  * (b)

Let $P=(Ca)^{-}-2(Aa)^{+}$ and $Q=(Cc)^{-}-2(Ac)^{+}$.

The proposition below is our key fact for analyzing $\ker\bar{K}_{w}$. In
particular, it shows $\bar{K}_{w}$ is self-adjoint and positive semi-definite.
Observe that as $A_{k}^{*}=\frac{1}{2}C_{k}$ and $a_{k}^{*}=c_{k}$, we have

$P^{*}=2(Ac)^{-}-(Cc)^{+}\ \ \ \ \text{and}\ \ \ \ Q^{*}=2(Aa)^{-}-(Ca)^{+}.$

######  Proposition 4.6.

$\bar{K}_{w}=P^{*}P+Q^{*}Q$ and

$Q^{*}Q-QQ^{*}=2\sum_{k}C_{k}A_{k}+2\sum_{k\leq
m}c_{k}a_{k}-2\sum_{m<\ell}c_{\ell}a_{\ell}+2m.$ (4.8)

Proof. Since $\partial_{k}=2A_{k}$ and $y_{k}=A_{k}+\frac{1}{2}C_{k}$, by
(4.7) and (4.4) one has

$\displaystyle\bar{K}_{w}$ $\displaystyle=\sum_{k,\ell\leq
m}(C_{k}A_{\ell}+C_{\ell}A_{k})(a_{k}c_{\ell}+c_{k}a_{\ell})+\sum_{m<k,\ell}(C_{k}A_{\ell}+C_{\ell}A_{k})(a_{k}c_{\ell}+c_{k}a_{\ell})$
$\displaystyle+2\sum_{k\leq m}(a_{k}c_{k}+c_{k}a_{k})-2\sum_{k\leq
m<\ell}(4A_{k}A_{\ell}+C_{k}C_{\ell})(a_{k}c_{\ell}+c_{k}a_{\ell}).$ (4.9)

Since
$C_{\ell}A_{k}(a_{k}c_{\ell}+c_{k}a_{\ell})=C_{\ell}A_{k}(c_{\ell}a_{k}+a_{\ell}c_{k})$
by (2.1), the equation (4.4) becomes

$\displaystyle\bar{K}_{w}$
$\displaystyle=2(Ca)^{-}(Ac)^{-}+2(Cc)^{-}(Aa)^{-}+2(Ca)^{+}(Ac)^{+}+2(Cc)^{+}(Aa)^{+}$
$\displaystyle-4(Aa)^{-}(Ac)^{+}-(Cc)^{-}(Ca)^{+}-4(Ac)^{-}(Aa)^{+}-(Ca)^{-}(Cc)^{+}$
$\displaystyle+2\sum_{k\leq m}(a_{k}c_{k}+c_{k}a_{k}).$ (4.10)

On the other hand, by (2.1) and (4.7), one obtain

$\displaystyle(Ca)^{-}(Ac)^{-}$ $\displaystyle=(Ac)^{-}(Ca)^{-}+\sum_{k\leq
m}A_{k}C_{k}-\sum_{k\leq m}c_{k}a_{k}-m,$ (4.11)
$\displaystyle(Cc)^{-}(Aa)^{-}$ $\displaystyle=(Aa)^{-}(Cc)^{-}-\sum_{k\leq
m}A_{k}C_{k}-\sum_{k\leq m}a_{k}c_{k}+m,$ $\displaystyle(Ca)^{-}(Cc)^{+}$
$\displaystyle=(Cc)^{+}(Ca)^{-},\ \ \ \ (Cc)^{-}(Ca)^{+}=(Ca)^{+}(Cc)^{-}.$

Now, plugging (4.11) to the equation (4.4) and then factoring gives the first
assertion. For the second assertion, apply (2.1) and (4.7) to have

$2(Ca)^{+}(Ac)^{+}=2(Ac)^{+}(Ca)^{+}+2\sum_{m<\ell}A_{\ell}C_{\ell}-2\sum_{m<\ell}c_{\ell}a_{\ell}.$

This and the second equation of (4.11) then yields the second assertion. ∎

######  Remark 4.7.

The weighted operator $\bar{K}_{w}=w^{-1}\bar{K}w$ is positive semi-definite
with respect to the inner product on $L^{2}(\mathbb{R}^{n},d\mu)\otimes
S^{p}$, so the model operator $\bar{K}$ is positive semi-definite with respect
to the inner product on $L^{2}(\mathbb{R}^{n})\otimes S^{p}$.

As a corollary to Proposition 4.6, we have:

######  Corollary 4.8.

If $m=n$, then $\ker\bar{K}_{w}=0$.

Proof. Let $m=n$. Then $(Aa)^{+}=(Ac)^{+}=0$, hence we have

$\bar{K}_{w}=(Ca)^{*}Ca+(Cc)^{*}Cc,$

where $Ca=\sum_{k}C_{k}a_{k}$ and $Cc=\sum_{k}C_{k}c_{k}$. As $Cc$ is
injective, we conclude $\bar{K}_{w}$ is positive definite, hence
$\ker\bar{K}_{w}=0$. ∎

###  4.5 Proof of Theorem 1.1 (a)

By Corollary 4.8 and $\dim\ker\bar{K}=\dim\ker\bar{K}_{w}$, the proposition
below completes the proof of Theorem 1.1 (a).

######  Proposition 4.9.

If $m\geq p$, then $\ker\bar{K}_{w}=0$.

Proof. Suppose $\beta\in\ker\bar{K}_{w}$. Then by Proposition 4.6,
$Q^{*}Q(\beta)=0$. Noting that $QQ^{*}$ is positive semi-definite, by (4.8) we
have

$\displaystyle 0$
$\displaystyle\geq\Big{(}\scalerel*[5.5pt]{\Big{(}}{\ensurestackMath{\addstackgap[1.5pt]{\big{(}}}}\sum_{k}C_{k}A_{k}+\sum_{k\leq
m}c_{k}a_{k}-\sum_{m<\ell}c_{\ell}a_{\ell}+m\scalerel*[5.5pt]{\Big{)}}{\ensurestackMath{\addstackgap[1.5pt]{\big{)}}}}\beta,\beta$
$\displaystyle=2\sum_{k}\big{(}A_{k}\beta,A_{k}\beta)+\sum_{k\leq
m}(a_{k}\beta,a_{k}\beta)-\sum_{m<\ell}\big{(}a_{\ell}\beta,a_{\ell}\beta)+m||\beta||^{2}.$
(4.12)

Since $\beta$ can be written as $\beta=\sum_{q,J}\beta^{q}_{J}e^{J}$, where
$\beta^{q}_{J}e^{J}\in{\cal H}^{q}\otimes S^{p}$, by (4.2), (4.5) and (4.6),
the inequality (4.5) implies that

$0\geq\sum_{q,J}\Big{(}q+\sum_{k\leq
m}j_{k}-\sum_{m<\ell}j_{\ell}+m\Big{)}||\beta^{q}_{J}e^{J}||^{2}.$

As $\sum j_{\ell}\leq|J|=p\leq m$, it follows that

$0\geq\sum_{q,J}\Big{(}q+\sum_{k\leq
m}j_{k}\Big{)}||\beta^{q}_{J}e^{J}||^{2}.$

This inequality implies that $\beta^{q}_{J}=0$ for all $q\neq 0$, hence
$\beta=\sum_{J}\beta^{0}_{J}e^{J}\in{\cal H}^{0}\otimes S^{p}$ such that

$0=Q(\beta)=(Cc)^{-}\beta-2(Ac)^{+}\beta=(Cc)^{-}\beta.$

Since $(Cc)^{-}$ is injective on ${\cal H}^{0}\otimes S^{p}$ (as $m\neq 0$),
we conclude that $\beta=0$. ∎

##  5 Dimension count

This section proves Theorem 1.1 (b) and (c) in the Introduction.

###  5.1 Weighted operator with Morse index zero

Let $\bar{K}_{w}$ be the weighted operator as in Section 4. From now on, we
assume that the Morse index $m=0$. In this case, as $(Ca)^{-}=(Cc)^{-}=0$, by
Proposition 4.6,

$\bar{K}_{w}=(2Aa)^{*}2Aa+(2Ac)^{*}2Ac,$

where $Aa=\sum_{k}A_{k}a_{k}$ and $Ac=\sum_{k}A_{k}c_{k}$. Consider

---
$\textstyle{{\cal H}^{q+1}\otimes
S^{p-1}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{Ac}$$\textstyle{{\cal
H}^{q}\otimes
S^{p}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{Ac}$$\scriptstyle{Ca}$$\scriptstyle{Aa}$$\textstyle{{\cal
H}^{q-1}\otimes
S^{p+1}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{Ca}$$\textstyle{{\cal
H}^{q-1}\otimes S^{p-1}}$

and observe that:

  * (a)

$\bar{K}_{w}({\cal H}^{q}\otimes S^{p})\subset{\cal H}^{q}\otimes S^{p}$.

  * (b)

If $\beta=\sum_{q}\beta_{q}\in\ker\mathbb{K}_{w}$, where $\beta_{q}\in{\cal
H}^{q}\otimes S^{p}$, then by (a),

$\beta_{q}\in\ker\bar{K}_{w},\ \ \ \forall q\geq 0.$

  * (c)

$\ker\bar{K}_{w}=\ker(Aa)\cap\ker(Ac)$. In particular, ${\cal H}^{0}\otimes
S^{p}\subset\ker\bar{K}_{w}$.

  * (d)

By (4.5) and (4.6), it follows from (4.8) that

$(CaAc-AcCa)|_{{\cal H}^{q}\otimes S^{p}}=(q-p)Id.$ (5.1)

  * (e)

Since both $CaAc$ and $AcCa$ are positive semi-definite, (5.1) shows

$q>p\Rightarrow\ker Ac|_{{\cal H}^{q}\otimes S^{p}}=0\ \ \ \text{and}\ \ \
q<p\Rightarrow\ker Ca|_{{\cal H}^{q}\otimes S^{p}}=0.$ (5.2)

With this understood, we define:

######  Definition 5.1.

Let $\bar{K}_{w}$ be as above with $m=0$ and define

$\mathbb{K}^{q,p}=\ker(Aa|_{{\cal H}^{q}\otimes S^{p}})\cap\ker(Ac|_{{\cal
H}^{q}\otimes S^{p}}).$

One may then use the above facts (a) – (e) to conclude:

######  Lemma 5.2.

Let $\bar{K}_{w}$ be as above with $m=0$. Then

$\ker\bar{K}_{w}=\bigoplus_{q\leq p}\mathbb{K}^{q,p}\ \ \ \ \text{and}\ \ \ \
\mathbb{K}^{0,p}={\cal H}^{0}\otimes S^{p}.$

###  5.2 Proof of Theorem 1.1 (b)

Let $p=1$ and $m=0$. In this case, by Lemma 5.2,

$K^{n,1}_{0}=\dim\ker\bar{K}=\dim\ker\bar{K}_{w}=\dim{\cal H}^{0}\otimes
S^{1}+\dim\mathbb{K}^{1,1}.$

Since $\dim{\cal H}^{0}\otimes S^{1}=n$, Theorem 1.1 (b) follows from the
following lemma.

######  Lemma 5.3.

${\displaystyle\dim\mathbb{K}^{1,1}=\binom{n}{2}}$.

Proof. Write $\beta\in{\cal H}^{1}\otimes S^{1}$ as
$\beta=\sum_{k,\ell}\beta^{k}_{\ell}H_{k}e^{\ell}$. Then

$Ac(\beta)=\sum_{k,\ell}\beta^{k}_{\ell}H_{0}e^{k}e^{\ell}=0\ \
\Longleftrightarrow\ \ \beta^{k}_{\ell}+\beta^{\ell}_{k}=0,\ \forall k,\ell.$

This implies that $\ker Ac|_{{\cal H}^{1}\otimes S^{1}}$ is spanned by the
vectors $H_{k}e^{\ell}-H_{\ell}e^{k}$. The lemma now follows since
$Aa(H_{k}e^{\ell}-H_{\ell}e^{k})=0$ for all $k\neq\ell$. ∎.

###  5.3 Proof of Theorem 1.1 (c)

For $p=2$ and $m=0$, by Lemma 5.2,

$K^{n,2}_{0}=\dim\ker\bar{K}=\dim\ker\bar{K}_{w}=\dim{\cal H}^{0}\otimes
S^{2}+\dim\mathbb{K}^{1,2}+\dim\mathbb{K}^{2,2}.$

Since ${\displaystyle\dim{\cal H}^{0}\otimes S^{2}=\binom{n+1}{2}}$, Theorem
1.1 (c) follows from Lemma 5.4 and Lemma 5.5 below.

######  Lemma 5.4.

$\mathbb{K}^{1,2}=0$ for $n=1$, and if $n\geq 2$, we have

$\dim\mathbb{K}^{1,2}=n{n+1\choose 2}-{n+2\choose 3}-n.$

Proof. Since the map $Ca:{\cal H}^{0}\otimes S^{3}\to{\cal H}^{1}\otimes
S^{2}$ is injective by (5.2), its adjoint map is surjective. Hence,

$\dim\ker Ac|_{{\cal H}^{1}\otimes S^{2}}=n{n+1\choose 2}-{n+2\choose 3}.$
(5.3)

If $n=1$, $\ker Ac|_{{\cal H}^{1}\otimes S^{2}}$ is trivial by (5.3), hence so
is $\mathbb{K}^{1,2}$. For $n\geq 2$, by the dimension formula (5.3), it
suffices to show that the restriction map $Aa:\ker Ac|_{{\cal H}^{1}\otimes
S^{2}}\to{\cal H}^{0}\otimes S^{1}$ is surjective. Let
$u_{k\ell}=H_{k}e^{k}e^{\ell}-H_{\ell}e^{k}e^{k}$. Then

$Ac(u_{k\ell})=0\ \ \ \ \text{and}\ \ \ \ Aa(u_{k\ell})=e^{\ell}.$

This shows the restriction map $Aa:\ker Ac|_{{\cal H}^{1}\otimes
S^{2}}\to{\cal H}^{0}\otimes S^{1}$ is surjective. ∎

######  Lemma 5.5.

$\mathbb{K}^{2,2}=0$ for $n=1,2$, and if $n\geq 3$, we have

$\dim\mathbb{K}^{2,2}={n+1\choose 2}{n+1\choose 2}-n{n+2\choose 3}-{n+1\choose
2}.$ (5.4)

Proof. For notational simplicity, let $Ac=Ac|_{{\cal H}^{2}\otimes S^{2}}$.
Since the map

$Ca:{\cal H}^{1}\otimes S^{3}\to{\cal H}^{2}\otimes S^{2}$

is injective by (5.2), its adjoint $Ac$ is surjective; hence we have

$\dim\ker Ac={n+1\choose 2}{n+1\choose 2}-n{n+2\choose 3}.$

This shows $\mathbb{K}^{2,2}=0$ when $n=1$. Let $n\geq 2$. One may choose, by
inspection, linearly independent vectors in $\ker Ac$ as follows. For distinct
$i,j,k,\ell$, let

$\displaystyle u_{k\ell}$
$\displaystyle=2H_{k\ell}e^{k\ell}-H_{\ell\ell}e^{kk}-H_{kk}e^{\ell\ell}\
(n\geq 2),$ $\displaystyle v_{k\ell j}$ $\displaystyle=H_{kk}e^{\ell
j}-H_{k\ell}e^{kj}-H_{kj}e^{k\ell}+H_{\ell j}e^{kk}\ (n\geq 3),$
$\displaystyle w_{ijk\ell}$
$\displaystyle=H_{ij}e^{k\ell}-2H_{ik}e^{j\ell}+H_{i\ell}e^{jk}+H_{jk}e^{i\ell}-2H_{j\ell}e^{ik}+H_{k\ell}e^{ij}\
(n\geq 4).$

Here $H_{k\ell}=H_{(k,\ell)}$ and $e^{k\ell}=e^{(k,\ell)}$. These vectors form
a basis for $\ker Ac$ because they span a subspace of dimension

${n\choose 2}+3{n\choose 3}+2{n\choose 4}={n+1\choose 2}{n+1\choose
2}-n{n+2\choose 3}=\dim\ker Ac.$

Observing that

$Aa(u_{k\ell})=2H_{\ell}e^{\ell}+2H_{k}e^{k},\ \ Aa(v_{k\ell
j})=-H_{\ell}e^{j}-H_{j}e^{\ell},\ \ Aa(w_{ijk\ell})=0,$

we define a subspace of ${\cal H}^{1}\otimes S^{1}$ by

$S=\langle H_{k}e^{k}+H_{\ell}e^{\ell}\,|\,k\neq\ell\rangle\oplus\langle
H_{k}e^{\ell}+H_{\ell}e^{k}\,|\,k\neq\ell\rangle\subset{\cal H}^{1}\otimes
S^{1}.$

Obviously, $Aa(\ker Ac)\subset S$ and the restriction map $Aa:\ker Ac\to S$ is
injective for $n=2$, and surjective for $n\geq 3$. Therefore, the lemma
follows from counting dimension of $S$: for $n\geq 3$,

$\dim\langle H_{k}e^{k}+H_{\ell}e^{\ell}\,|\,k\neq\ell\rangle=n\ \ \
\text{and}\ \ \ \dim\langle
H_{k}e^{\ell}+H_{\ell}e^{k}\,|\,k\neq\ell\rangle={n\choose 2},$

where the first equality follows since $Aa(u_{k\ell}+u_{kj}-u_{\ell
j})=2H_{k}e^{k}$. ∎

## References

  * [1] C. Barbance, Sur les tenseurs symètriques, C.R. Acad. Sc. paris, 276, Sér. A, 387–389.
  * [2] N. Berline, E. Getzler, and M. Vergne, Heat kernels and Dirac operators, Springer-Verlag, Berlin 1992.
  * [3] M. Cahen, S. Gutt, L. Gravy, and J. Rawmsley, On Mpc-structures and symplectic Dirac operators, J. Geom. Phys. 86 (2014), 434–466.
  * [4] H. Cycon, R. Froese, W. Kirsch, and B. Simon, Schrödinger operators with application to quantum mechanics and global geometry, Texts Monogr. Phys. Springer Study Ed. Springer-Verlag, Berlin, 1987.
  * [5] R. Delong, Killing tensors and the Hamilton-Jacobi equation, Ph.D. thesis, University of Minnesota, 1982.
  * [6] G. Folland, Harmonic analysis in phase space, Ann. of Math. Stud., 122 Princeton University Press, Princeton, NJ, 1989.
  * [7] K. Heil, A. Moroianu and U. Semmelmann, Killing and conformal Killing tensors, J. Geom. Phys. 106 (2016), 383–400.
  * [8] A. Kimaczyńska and A. Pierzchalski, Elliptic operators in the bundle of symmetric tensors, Banach Center Publ., 113 Polish Academy of Sciences, Institute of Mathematics, Warsaw, 2017, 193–218.
  * [9] J. Lee and T. Parker, Spin Hurwitz numbers and the Gromov-Witten invariants of Kähler surfaces, Comm. Anal. Geom. 21 (2013), no. 5, 1015–1060.
  * [10] M. Maridakis, Spinor pairs and the concentrating principle for Dirac operators, Trans. Amer. Math. Soc. 369 (2017), no. 3, 2231–2254.
  * [11] R. McLenaghan, R. Milson, and R. Smirnov, Killing tensors as irreducible representations of the general linear group, C. R. Math. Acad. Sci. Paris 339 (2004), no. 9, 621–624.
  * [12] L. Nicolaescu, An invitation to Morse theory, Universitext Springer, New York, 2011.
  * [13] L. Nicolaescu, Lectures on Geometry of Manifolds, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2021.
  * [14] I. Prokhorenkov and K. Richardson, Perturbations of Dirac operators, J. Geom. Phys. 57 (2006), no. 1, 297–321.
  * [15] J.H. Sampson, On a theorem of Chern, Trans. Amer. Math. Soc. 177 (1973), 141–153.
  * [16] M. Shubin, Semiclassical asymptotics on covering manifolds and Morse inequalities, Geom. Funct. Anal. 6 (1996), no. 2, 370–409.
  * [17] M. Shubin, Novikov inequalities for vector fields, Gelfand Math. Sem. Birkhäuser Boston, Inc., Boston, MA, 1996, 243–274.
  * [18] B. Simon, Semiclassical analysis of low lying eigenvalues. I. Nondegenerate minima: asymptotic expansions, Ann. Inst. H. Poincaré Sect. A (N.S.) 38 (1983), no. 3, 295–308.
  * [19] M. Takeuchi, Killing tensor fields on spaces of constant curvature, Tsukuba J. Math. 7 (1983) 233–255.
  * [20] C. Taubes, Counting pseudo-holomorphic submanifolds in dimension 4, J. Diff. Geom., 44 (1996), no. 4, 818–893.
  * [21] G. Thompson, Killing tensors in spaces of constant curvature, J. Math. Phys. 27 (1986) 2693–2699.
  * [22] Edward Witten, Supersymmetry and Morse theory, J. Diff. Deom., 17 (1982), no. 4, 661–692.
  * [23] W. Zhang, Lectures on Chern-Weil theory and Witten deformations, Nankai Tracts Math., 4 World Scientific Publishing Co., Inc., River Edge, NJ, 2001.

Department of Mathematics, Dongguk University, Seoul 04620, Korea

<EMAIL_ADDRESS>
Department of Mathematics, University of Central Florida, Orlando, FL 32816

Korea Institute for Advanced Study, Seoul 02455, Korea

<EMAIL_ADDRESS>


# No Sub-Saturn-mass Planet Desert in the CORALIE/HARPS Radial-velocity Sample

David P. Bennett Laboratory for Exoplanets and Stellar Astrophysics, NASA
Goddard Space Flight Center, Greenbelt, MD 20771, USA  Department of
Astronomy, University of Maryland, College Park, MD 20742, USA  Clément Ranc
Sorbonne Université, CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98 bis
bd Arago, 75014 Paris, France  Rachel B. Fernandes  Lunar and Planetary
Laboratory, The University of Arizona, Tucson, AZ 85721, USA  Earths in Other
Solar Systems Team, NASA Nexus for Exoplanet System Science

###### Abstract

We analyze the CORALIE/HARPS sample of exoplanets (Mayor et al., 2011) found
by the Doppler radial velocity method for signs of the predicted gap or
“desert” at 10-$100{M_{\oplus}}$ caused by runaway gas accretion at semimajor
axes of $<3\,$AU. We find that these data are not consistent with this
prediction. This result is similar to the finding by the MOA gravitational
microlensing survey that found no desert in the exoplanet distribution for
exoplanets in slightly longer period orbits and somewhat lower host masses
(Suzuki et al., 2018). Together, these results imply that the runaway gas
accretion scenario of the core accretion theory does not have a large
influence on the final mass and semimajor axis distribution of exoplanets.

planetary systems — solar system: formation

##  1 Introduction

The runaway gas accretion process has long been considered to be a major step
in the formation of gas giant planets in the core accretion theory (Pollack et
al., 1996). Several authors (Ida & Lin, 2004a; Mordasini et al., 2009a, b)
have predicted that this runaway gas accretion process should produce a gap or
“desert” in the mass distribution of planets. The abstract of Ida & Lin
(2004a) states “Since planets’ masses grow rapidly from 10 to $100M_{\earth}$,
the gas giant planets rarely form with asymptotic masses in this intermediate
range.” They specifically indicate that this desert should exist at orbital
major axes of 0.2–$3\,$AU. The Bern population synthesis group (Mordasini et
al., 2009a, b) predict a lower amplitude desert, a factor of 2–3 at masses
between 30 to $100M_{\earth}$. The paper presenting the data that we analyze
in this paper (Mayor et al., 2011) claims that these data tend to confirm this
prediction. They state: “After correction of detection biases, we see even
more clearly the importance of the population of low-mass planets on tight
orbits, with a sharp decrease of the distribution between a few Earth masses
and $\sim 40{M_{\oplus}}$. We note that the planet population synthesis models
by Mordasini et al. (2009a, b) predicted such a minimum in the mass-
distribution at precisely this mass range.” As we shall see, our
interpretation differs somewhat from this, and it strongly contradicts a more
recent version of the Bern group’s population synthesis code (Emsenhuber et
al., 2020a, a), which predicts a deeper desert, with a planet occurrence rate
drop of a factor of $\sim 10$.

In this paper, we present the a re-analysis of the Mayor et al. (2011) Doppler
radial velocity detected planet sample, that was recently analyzed by
Fernandes et al. (2019). While Fernandes et al. (2019) focused primarily on
the orbital period dependence of the planet distribution, our analysis is
focused primarily on the dependence on planetary mass, so that we can
investigate the proposed sub-Saturn-mass planetary desert. Before considering
the Mayor et al. (2011) sample, however, we review recent results from
gravitational microlensing that seem to contradict this cold sub-Saturn-mass
desert prediction in Section 2. Then, in Section 3, we present the binned,
completeness corrected planet distribution with an estimate of the error bars
for each bin, and we show that this data shows no significant sign of a sub-
Saturn-mass planet desert. In Section 4, we present a more powerful forward
modeling method to search for evidence of a planet desert in the Mayor et al.
(2011) data using a planet occurrence model that employs a Gaussian desert
feature, and again, we find no evidence of a sub-Saturn-mass planet desert. We
also compare our results to the Fernandes et al. (2019) analysis and recover
their main result. We then consider the latest generation of planet population
synthesis models from the Bern group (Emsenhuber et al., 2020a, b) in Section
5, and we find that our forward Gaussian desert model does a reasonable job at
characterizing the planet desert feature in their simulated planet population
that is contradicted by the Mayor et al. (2011) CORALIE/HARPS data. Emsenhuber
et al. (2020c) have also pointed out that deep planet desert at
$20{M_{\oplus}}<M<200{M_{\oplus}}$ seen in their models is contradicted by the
Mayor et al. (2011) data. Next, we show that a possible sub-Jupiter “radius
desert” that may be seen in the Kepler data for wider orbit planets is likely
to be due to the planet mass-radius relation instead of a desert in the mass
distribution. In Section 6, we present our conclusions and discuss our
expectations for future developments in the study of wide orbit planet
demographics.

Note that the planet desert that we investigate in this paper differs from
other “deserts” that have been observed in or suggested by Kepler data. These
are considered in Appendix A.

##  2 Search for a Planet Desert Beyond the Snow Line

A recent comparison of theoretical population synthesis calculations, based on
the core accretion theory, with the distribution of planets detected by
gravitational microlensing (Suzuki et al., 2018), orbiting beyond the snow
line, found that the theoretical models predicted a strong desert at mass
ratios of $10^{-4}\leq q\leq 4\times 10^{-4}$ that was not reproduced by the
microlensing results. Instead, the exoplanet mass ratio distribution found by
microlensing indicated a smooth, power-law distribution down to a peak at a
mass ratio of $q\sim 10^{-4}$, with fewer planets at lower mass ratios
(Udalski et al., 2018; Jung et al., 2019). This discrepancy between theory and
observation was thought to be likely due to deficiencies in the assumptions
and approximations in the population synthesis calculations instead of
problems with the fundamental concepts of the core accretion theory. One
possibility that has been suggested for our own Solar System is that captured
planetesimals could heat the accreting gaseous envelope, which would act to
slow accretion, and other possibilities are that the scale height or viscosity
of the protoplanetary disk may be smaller than the values adopted for these
population synthesis calculations (Dobbs-Dixon et al., 2007; Fung et al.,
2014; Garaud & Lin, 2007). A reduction of the assumed disk viscosity in the
Bern group’s population synthesis model did indicate that this sub-Saturn-mass
desert could disappear with lower viscosity (Suzuki et al., 2018), as earlier
work had indicated (Szulágyi et al., 2014). Annular structures in
protoplanetary disks observed by ALMA are thought by many to be due to planets
in the process of formation. A comparison of these structures in a number of
protoplanetary disks also indicates that the growth of gas giant planets is
slower than predicted by standard runaway gas accretion models (Nayakshin et
al., 2019), so runaway gas accretion is disfavored unless these annular
structures are caused by another mechanism.

##  3 Planet Occurrence for the CORALIE/HARPS Sample

Figure 1:  Planet occurrence rate histograms for the planets in the Mayor et
al. (2011) sample in four different orbital period ranges, each spanning a
factor of eight in period. The histograms without any completeness corrections
are shown in grey and the light grey shaded regions show the location of the
“planet desert” predicted by Ida & Lin (2004a).

The CORALIE/HARPS data set (Mayor et al., 2011) that we analyze in this paper
was recently analyzed by Fernandes et al. (2019), who considered the
distribution of planets in the 30–$6000{M_{\oplus}}$ mass range. They fit this
with a broken power-law model that has a break and peak at a period of
$1720\pm 430\,$days. When they add a power-law for the mass distribution of
planets in this mass range, the range of possible locations of the peak
expands to $2100\pm 1200\,$days.

The Fernandes et al. (2019) analysis of the Mayor et al. (2011) data set
included a tabulation of the survey completeness from Figure 6 of Mayor et al.
(2011), and we use this tabulation of the survey completeness for our
analysis. As a first step in our analysis, we present the binned $M\sin i$
distribution of the Mayor et al. (2011) sample in Figure 1, where $M$ is the
planet mass and $i$ is the inclination of the planet’s orbital plane to the
line-of-sight. We break up the Mayor et al. (2011) sample into 4 different
period ranges, each spanning a factor of eight in period: 3.75–$30\,$days,
30–$240\,$days, 240–$1920\,$days and 1920–$15260\,$days. The boundary between
the third and fourth period ranges is set to be $1920\,$days because this
corresponds to the approximate outer semi-major axis boundary of $3\,$AU of
the Ida & Lin (2004a) prediction for a $1M_{\odot}$ host star and is within
$0.5\sigma$ of the break in the giant planet occurrence rate found by
Fernandes et al. (2019).

In order to relate the planets found in the CORALIE/HARPS data set to the
intrinsic planetary occurrence rate, we correct for the survey completeness,
which gives the colored histograms with error bars shown in Figure 1. The
occurrence rate is given by

$\eta={1\over N_{\star}}\sum_{j}^{n_{p}}{1\over C(M_{j}\sin i_{j},P_{j})}\ ,$
(1)

where $N_{\star}=822$ is the number of stars in the CORALIE/HARPS sample,
$n_{p}$ is the number of detected planets, and $M_{j}$, $i_{j}$, and $P_{j}$
are the mass, orbital inclination, and the period of the $j$th planet in the
sample. The survey completeness for mass $M$, inclination $i$ and period $P$
is given by $C(M\sin i,P)$.

The grey histograms without error bars in Figure 1. indicate the raw,
uncorrected histograms. The grey and colored histograms merge for large $M\sin
i$ and small periods, as expected, because these planets can be detected with
nearly 100% completeness.

In order to present the planet occurrence rate as a function of mass (or
$M\sin i$) and period, we bin the data in both $M\sin i$ and period. Because
the CORALIE/HARPS data set is sensitive to a very wide range of planetary
orbital periods and masses, it is necessary to select bins that span a
relatively large range in both $M\sin i$ and period in order to have a
reasonable number of planets in each bin. This complicates the issue of
determining the error bars for the completeness corrected occurrence rate for
each bin, because the completeness can vary significantly across some of the
bins. Using the average completeness over the range of parameters for each bin
could induce an error if the planet occurrence rate also varies significantly
over the bin size. Therefore, we use the completeness for the $M\sin i$ and
$P$ for each planet in equation 1. This means that planets in the same bin
will have different contributions to the planet occurrence rate, $\eta$, so
that we cannot use Poisson statistics to calculate the error bars for $\eta$.
Therefore, we employ a bootstrap method to determine the error bars. We take
the number of planets in each bin as the expected number for a Poisson
distribution, and then we randomly select the planets from the list of planets
in the bin. We repeat this procedure 100,000 times for each bin to determine
the root-mean-square dispersion of the completeness corrected $\eta$ value,
which we use as the error estimate for the $\eta$ value for each bin. For bins
without any planets, we calculate the mean completeness for the bin, and then
use the one-sided 1–$\sigma$ upper limit from Poisson statistics for zero
detections to determine the 1–$\sigma$ upper limit.

The grey shaded regions in the upper three panels of Figure 1 indicate the
mass range of the planet desert predicted by Ida & Lin (2004a). As noted by
Mayor et al. (2011), there is a significant decrease in the exoplanet
occurrence rate, $\eta$, between 10 and $\sim 40{M_{\oplus}}$ at the shorter
orbital periods, particularly in the $3.75<P<30\,$days range. (With our
binning, this minimum appears at $\sim 30{M_{\oplus}}$.) So, planets at the
bottom edge of the predicted “desert” are clearly more common, but to claim a
desert, the exoplanet occurrence must increase at higher masses, and such a
feature is not seen in Figure 1. Also, as mentioned in the introduction, Mayor
et al. (2011) point to a minimum in the planet distribution at $\sim
40{M_{\oplus}}$ as being consistent with the Mordasini et al. (2009a, b)
prediction of a “moderate” planet desert in the range 30–$100{M_{\oplus}}$,
with a reduction of a factor of 2–3 from the number of Jupiter-mass planets.
Our analysis does indicate a minimum at $\sim 30{M_{\oplus}}$ in the
$3.75<P<30\,$days bin, but our error bar estimates indicate that this is not
significant even at the 1–$\sigma$ level. It is possible that this minimum
might appear more significant with a different binning, but we will address
this issue in Section 4 with a forward modeling approach that involves no
binning. Also, note that there is no evidence for the predicted factor of 2–3
deficit of planets at 30–$100{M_{\oplus}}$ compared to Jupiter-mass ($\sim
318{M_{\oplus}}$) planets (Mordasini et al., 2009a, b) in any of the period
ranges shown in Figure 1. There does seem to be a marginal increase of
Jupiter-mass planets compared to 100–$200{M_{\oplus}}$ planets in the
240–$1920\,$days period range, but this is only marginally significant and is
not duplicated at shorter or longer periods.

In this discussion, we have avoided the distinction between planet masses,
$M$, and the $M\sin i$ values that are measured by radial velocity monitoring.
Because small inclinations, $i\ll 1$, are disfavored geometrically, this is
generally expected to have a small effect, but it could conceivably be
important if the occurrence rate, $\eta$, is strongly dependent on planet
mass. Fortunately, this complication disappears with the forward modeling
approach presented in the next section.

##  4 Forward Modeling Analysis

We now turn to a more rigorous forward modeling analysis to search for
evidence of this predicted desert, attributed to the runaway gas accretion
process of the core accretion theory. Forward modeling allows us to
investigate the possibility of a desert with somewhat different parameters
than proposed by Ida & Lin (2004a) and Mordasini et al. (2009a, b), and it
enables us to avoid the uncertainties relating to the binning of data and
error bar estimates mentioned in the previous section.

###  4.1 Exoplanet Desert Occurrence Rate Models

Following earlier Doppler radial velocity analyses (Cumming et al., 2008;
Fernandes et al., 2019; Wittenmyer et al., 2020), we base our model on power-
laws in both mass, $M$, and period, $P$. The previous analysis of the
CORALIE/HARPS sample (Fernandes et al., 2019) found a break in the planet
occurrence rate at a period of $P=2100\pm 1200\,$days, which is consistent
with the $P=1920\,$days boundary between the third and fourth period ranges
presented in Figure 1. (Fernandes et al. (2019) indicate that this broken
power-law applies to masses in the range 30-6000$\,{M_{\oplus}}$, but Figure 1
makes it clear that the Mayor et al. (2011) sample has little sensitivity for
$M<100{M_{\oplus}}$ at periods longer than the break.) Since the $M\sin i$
dependence does not appear to follow the same power law in these four period
ranges, we elect to fit these four period ranges separately without any power-
law break.

In order to describe a planetary desert, we add the following desert reduction
factor, $D(P,M_{\rm des})$ to the power law model:

$D(M,M_{\rm des})=1-B\exp\left(-{\left[\log M-\log M_{\rm des}\right]^{2}\over
2\sigma_{\rm des}^{2}}\right)\ ,$ (2)

where the planet desert is centered at period $M_{\rm des}$ with a Gaussian
dispersion of $\sigma_{\rm des}$ in $\log P$. (Throughout this paper, $\log$
refers to the base-10 logarithm.) This desert reduction factor multiplies the
power-law factors to yield our planet desert mass function model:

$\displaystyle f_{\rm des}(M,P;A,B,\log M_{\rm des},\sigma_{\rm des},p,m)$
$\displaystyle\equiv\frac{d^{2}N_{\rm pl}}{d\,\log M\ d\,\log P}$
$\displaystyle=A\left(\frac{P}{1995\,{\rm days}}\right)^{p}D(M,M_{\rm
des})\left({M\over 100{M_{\oplus}}}\right)^{m}.$ (3)

We select the Gaussian function to describe the desert because we think that a
smooth function is more plausible than a step function, and we do find that it
provides a reasonably good description of the planet distribution predicted by
the most recent version of the Bern group’s planetary population synthesis
model (Emsenhuber et al., 2020a, b) in some period ranges, as discussed in
Section 5. Note that we have selected the pivot point for the period power-law
to be $10^{3.3}=1995\,$days in order to match the break in the power-law when
we include all periods longer than $3.75\,$days, as discussed in Section 4.2.

One complication for our forward modeling effort is that the completeness,
$C(M\sin i,P)$, depends on the inclination angle, $i$, while intrinsic planet
properties will depend only on the planet mass, $M$. Fortunately, it is well
established that there is no physical process that is likely to produce a non-
random distribution of planetary orbital planes as seen from Earth, so we can
treat the inclination angles, $i$, as a random variable to integrate over. We
introduce the variable $\mathfrak{M}\equiv M\sin i$ to refer to the observable
parameter measured by the Doppler radial velocity method.

Figure 2:  The desert power-law planet occurrence rate models compared to the
histograms from the top three period range panels of Figure 1. The black
curves are the best fit models, and the dark and light grey shaded regions
indicate the central 68.3% and 95.4% distribution of planet occurrence rates
as a function of mass from the MCMC distributions. These indicate no
significant evidence of a sub-Saturn-mass planet desert in the Mayor et al.
(2011) CORALIE/HARPS data set.

We use a standard forward modeling Bayesian likelihood analysis (Alcock et
al., 1996; Suzuki et al., 2016) to determine the posterior distribution of
models that are consistent with the data. The likelihood function is given by

$\displaystyle\mathcal{L}_{\rm des}($ $\displaystyle A,B,\log M_{\rm
des},\sigma_{\rm des},p,m)=e^{-N_{\rm exp}}$
$\displaystyle\times\prod_{j}^{N_{\star}}\int di\,\sin i\,f_{\rm
des}(\mathfrak{M}_{j}/\sin i,P_{j};A,B,\log M_{\rm des},\sigma_{\rm
des},p,m)\,C(\mathfrak{M}_{j},P_{j})\ ,$ (4)

where

$N_{\rm exp}=\iiint dM\,dP\,di\,\sin i\,f_{\rm des}(M,P;A,B,\log M_{\rm
des},\sigma_{\rm des},p,m)\,C(M\sin i,P)\ ,$ (5)

is the expected number planets detected for the model with the given
parameters: $A$, $B$ ,$\log M_{\rm des}$, $\sigma_{\rm des}$, $p$, and $m$. We
evaluate this likelihood function with a Markov Chain Monte Carlo (MCMC) for
each of the period ranges: 3.75–$30\,$days, 30–$240\,$days, 240–$1920\,$days
and 1920–$15260\,$days, with chains of 140,000-330,000 steps after 2000 burn-
in steps. The results of these calculations are summarized in Figure 2. The
best fit models are shown as the black curves, while the dark and light grey
shading indicates the central 68.3% and 95.4% likelihood distributions for the
planet occurrence rate as a function of mass from the MCMC calculations. These
plots are shown for the median $\log P$ values for each of the period ranges.
These central periods are $10.6\,$days for $3.75<P<30\,$days, $84.9\,$days for
$30<P<240\,$days, and $679\,$days for $240<P<1920\,$days.

The histograms of the exoplanet occurrence from the CORALIE/HARPS data set
from the top three panels of Figure 1 are also shown in Figure 2. However, the
$x$-axis for Figure 2 is the planet mass, $M$, instead of the Doppler radial
velocity observable $\mathfrak{M}=M\sin i$. So, we shift the $x$-axis from
$\mathfrak{M}$ to $M\approx\mathfrak{M}/\sin i_{\rm med}=1.1547$, where
$i_{\rm med}=\cos^{-1}(0.5)$ is the median inclination value, assuming a
random orientation of the orbital planes. This has no effect on the Bayesian
likelihood analysis, since we integrate over all inclination angles, $i$,
which uses the measured period and $\mathfrak{M}$ values without any binning.
The parameters of the best fit and average desert power-law models are given
in Table 1. The scatter in the normalization parameter $A$ is large because
the pivot point of the period power law is set to a fairly large period of
$1995\,$days, and this results in a relatively large scatter in the
normalization parameter, $A$. Therefore, we also include a parameter, $A_{P\rm
m}=A(P_{\rm m}/1995\,{\rm days})^{p}$, where $P_{\rm m}$ logarithmic median of
the period range being considered. For the different period ranges, we have
$P_{\rm m}=10.607$, 84.853, and $678.820\,$days for the 3.75–30, 30–240, and
240–$1920\,$days ranges, respectively. These $A_{P\rm m}$ values reflect the
variation in the normalization in the middle of the period ranges under
consideration.

Table 1: Exoplanet Desert Power-Law Model Parameters parameter | Period ranges
---|---
| 3.75–30 days | 30–240 days | 240–1920 days
| best | average | best | average | best | average
$A$ | 0.745 | $4.03\pm 2.44$ | 0.104 | $0.38\pm 0.32$ | 0.105 | $0.178\pm 0.087$
$A_{P\rm m}$ | 0.0177 | $0.020\pm 0.009$ | 0.0561 | $0.068\pm 0.022$ | 0.0920 | $0.125\pm 0.052$
$B$ | 0.977 | $0.52\pm 0.26$ | 0.767 | $0.48\pm 0.24$ | 0.954 | $0.51\pm 0.27$
$\log M_{\rm des}$ | 1.529 | $2.08\pm 0.48$ | 1.698 | $2.08\pm 0.42$ | 2.183 | $1.61\pm 0.44$
$\sigma_{\rm des}$ | 0.155 | $0.45\pm 0.28$ | 0.137 | $0.55\pm 0.28$ | 0.110 | $0.53\pm 0.29$
$p$ | 0.714 | $0.97\pm 0.18$ | 0.194 | $0.473\pm 0.23$ | 0.123 | $0.30\pm 0.24$
$m$ | $-1.229$ | $-1.18\pm 0.15$ | $-1.008$ | $-0.94\pm 0.12$ | $-0.407$ | $-0.49\pm 0.14$

Figure 2 indicates that the best fit models do have narrow exoplanet deserts,
and the deserts in the 3.75–$30\,$days and 30–$240\,$days period ranges are
centered in the range of 10–$100\,$days predicted by Ida & Lin (2004a). The
depth of the desert features in the best fit models, represented by the $B$,
is fairly deep at their centers, reducing the occurrence rate by over 95% for
the 3.75–$30\,$day and 240–$1920\,$day ranges and by 76.7%in the 30–$240\,$day
range. However, the distribution of models consistent with the data, as
represented by the grey shaded regions in Figure 2 tell a different story. The
central 68.3% and 95.4% model occurrence rate distributions generally do not
follow the best fit models, although the lowest extent of the central 95.4%
occurrence rate distribution for the 3.75–$30\,$days period range does have a
dip at the location of the best fit desert model, but there is virtually no
such feature in the central 68.3% occurrence rate distribution. These features
are also narrower than the Ida & Lin (2004a) prediction of a desert spanning a
factor of 10 in mass. The Gaussian $\sigma_{\rm des}$ values (in $\log M$)
range from $\sigma_{\rm des}=0.110$ for the 240–$1920\,$days period range to
$\sigma_{\rm des}=0.155$ for the 3.75–$30\,$days period range. These imply
full-width half-max values in the range 1.8–2.3 in $M$. The most likely
explanation for this is that the very narrow desert features of the best fit
models are due to statistical noise in the observed exoplanet distribution.
This would certainly explain why the central 68.3% range of the occurrence
rate distribution does not show these features.

We do note that the location of the best fit desert in the 3.75–$30\,$day and
30–$240\,$day period ranges do roughly correspond to the minimum at $M\sin
i=40{M_{\oplus}}$ mentioned by Mayor et al. (2011), particularly in their
sample with periods $<100\,$days. However, our analysis indicates that this
minimum is not likely to be significant. In fact, this conclusion seems
consistent with the completeness corrected histogram shown in Figure 12 of
Mayor et al. (2011), where the bin centered near $M\sin i=40{M_{\oplus}}$ is
consistent at 1–$\sigma$ with the next two larger bins. The large rise at
lower masses is significant, as is the case with our histogram in the top
panel of Figure 2, although our error bars are more conservative. However,
this is does not support the runaway accretion prediction that planets should
grow rapidly across the 10–$100\,{M_{\oplus}}$ mass interval.

###  4.2 Power-Law Occurrence Rate Models

Table 2: Pure Power-Law Model Parameters parameter | Period ranges
---|---
| 3.75–30 days | 30–240 days | 240–1920 days
| best | average | best | average | best | average
$A$ | 0.542 | $1.77\pm 0.62$ | 0.081 | $0.22\pm 0.16$ | 0.099 | $0.113\pm 0.035$
$A_{P\rm m}$ | 0.0148 | $0.0148\pm 0.0044$ | 0.0481 | $0.048\pm 0.009$ | 0.0812 | $0.083\pm 0.016$
$p$ | 0.701 | $0.90\pm 0.12$ | 0.165 | $0.40\pm 0.23$ | 0.179 | $0.26\pm 0.24$
$m$ | $-1.251$ | $-1.26\pm 0.13$ | $-0.972$ | $-0.98\pm 0.11$ | $-0.367$ | $-0.38\pm 0.09$

Since we have concluded that there is no significant indication of a planet
desert in the CORALIE/HARPS data set, it is sensible to consider simpler
models that do not include the desert. The model parameters for such models
are presented in Table 2, which includes both the best fit models and the MCMC
averages. These power-law models are given by equation 3 with $B\equiv 0$ and
therefore $D(M,M_{\rm des})\equiv 1$. As was the case for the planet desert
models, the parameter $A_{P\rm m}$ is not a model parameter. Instead, it
reflects the normalization of the occurrence rate function at the logarithmic
center of the period range under consideration. So, the scatter in $A_{P\rm
m}$ is a better indication of the uncertainty in the occurrence rate at the
center of the period range under consideration than $A$.

We can compare these power-law models to the corresponding power-law plus
desert models to see if the addition of the desert feature significantly
improves the fit. We can do this by comparing the $\log{\cal L}$ values for
the models without and with the desert features. We find that the desert
features improve the $\log{\cal L}$ values by 1.526, 0.722, and 0.944 for the
3.75–30, 30–240, and 240–$1920\,$days period ranges, respectively. Using
frequentist statistics, these would imply that the pure power-law models were
disfavored by $p=0.030$, 0.189, and 0.115, for these three period ranges.
Following equation 27 of Trotta (2008), we can compute upper bounds on the
Bayes factor of 3.52, 1.17, and 1.48 for the 3.75–30, 30–240, and
240–$1920\,$days period ranges. These correspond to weak evidence favoring the
power-law plus desert model for the 3.75–$30\,$days period range, and no
significant evidence in favor of the desert additions to the power-law models
at the longer period ranges. While there may be marginal evidence that the
desert model is an improvement over the pure power-law model for the
3.75–$30\,$days period range, we do not consider this to be evidence in favor
of the desert feature itself based upon the MCMC results.

Table 3: Broken Power-Law Model Parameters for $0.1M_{\rm Jup}\leq M\sin i\leq 20M_{\rm Jup}$ parameter |  $P$ range: 1–11000 days |
---|---|---
| best | average | EPOS
$A$ | 0.167 | $0.165\pm 0.035$ | $0.157^{+0.036}_{-0.030}$
$\log(P_{\rm br})$ | 3.32 | $3.20\pm 0.17$ | $3.32^{+0.19}_{-0.24}$
$p_{1}$ | 0.553 | $0.592\pm 0.085$ | $0.70^{+0.32}_{-0.16}$
$p_{2}$ | $-0.75$ | $-0.99\pm 0.53$ | $-1.20^{+0.92}_{-1.26}$
$m$ | $-0.462$ | $-0.475\pm 0.075$ | $-0.46\pm 0.06$

Finally, to compare with Fernandes et al. (2019), we also consider a model
with a broken power law in the period with masses constrained to the range
$0.1M_{\rm Jup}\leq M\sin i\leq 20M_{\rm Jup}$. This mass range is considered
(Fernandes et al., 2019) to be the range of giant planets. Our broken power-
law mass function model is given by

$\displaystyle f_{\rm bpl}(M,P;A,P_{\rm br},$ $\displaystyle
p_{1},p_{2},m)\equiv\frac{d^{2}N_{\rm pl}}{d\,\log M\ d\,\log P}$
$\displaystyle=A\left[\left(\frac{P}{P_{\rm br}}\right)^{p_{1}}\Theta(P-P_{\rm
br})+\left(\frac{P}{P_{\rm br}}\right)^{p_{2}}\Theta(P_{\rm
br}-P)\right]M^{m},$ (6)

where $\Theta$ is the Heavyside theta function (or step function), $P_{\rm
br}$ is the period break and $p_{1}$, $p_{2}$, and $m$ are the power-law
exponents for small periods, large periods and the planet mass, $M$,
respectively. The likelihood function then becomes

$\displaystyle\mathcal{L}_{\rm bpl}(A,P_{\rm br},$ $\displaystyle
p_{1},p_{2},m)=e^{-N_{\rm exp}}$ $\displaystyle\times\prod_{j}^{N_{\star}}\int
di\,\sin i\,f_{\rm bpl}(\mathfrak{M}_{j}/\sin i,P_{j};A,P_{\rm
br},p_{1},p_{2},m)\,C(\mathfrak{M}_{j},P_{j})\ ,$ (7)

and the number of expected planet detections becomes

$N_{\rm exp}=\iiint dM\,dP\,di\,\sin i\,f_{\rm bpl}(M,P;A,P_{\rm
br},p_{1},p_{2},m)\,C(M\sin i,P)\ .$ (8)

Evaluating this with a MCMC, as we did with the planet desert models, we find
the best fit model, as wll the range of models from the MCMC distribution. The
model parameter results are shown in the first three columns of Table 3
indicating a power law break at $P_{\rm br}=10^{3.20\pm
0.17}=1600^{+800}_{-500}\,$days. Note that we have used $\log(P_{\rm br})$
instead of $P_{\rm br}$ as the parameter to indicate the location of the break
in the period distribution in our MCMC calculations. This is equivalent to
assuming that the prior distribution of the $P_{\rm br}$ variable in uniform
in $\log(P_{\rm br})$ rather than uniform in $P_{\rm br}$. These results are
consistent with the Fernandes et al. (2019) results for the same model as
determined by the EPOS program, which are reproduced in the fourth column of
Table 3. (Note that the normalization parameter, $A$, has been modified from
the value in the Fernandes et al. (2019) paper to account for the fact that
they define the broken power law mass function over $d\,\ln M\ d\,\ln P$,
whereas we use $d\,\log_{10}M\ d\,\log_{10}P$.) Thus, we confirm the main
conclusion of Fernandes et al. (2019) that the giant planet occurrence rate
has a maximum in the vicinity of the snow line.

##  5 Comparison of CORALIE/HARPS Data and Gaussian Desert Model to Generation
III Bern Models

Figure 3:  Comparison of Bern New Generation Planetary Population Synthesis
calculations (red histograms) to Mayor et al. (2011) data (blue, cyan and
green histograms) and our Gaussian desert model (black curves and grey shaded
regions) using the same three period ranges used in our previous analyses:
3.75–30 days, 30–240 days, and 240-1920 days. Results from simulations
starting with 20 $0.01{M_{\oplus}}$ planetary embryos are shown in the left
panel and results from 100 $0.01{M_{\oplus}}$ planetary embryo simulations are
shown on the right.

In this section, we compare the CORALIE/HARPS exoplanet sample (Mayor et al.,
2011) to the results Generation III Bern population synthesis model of
planetary formation and evolution (Emsenhuber et al., 2020a, b), and we also
fit these model results with our Gaussian desert, power-law mass function
model for the three period ranges used in this paper: 3.75–30 days, 30–240
days, and 240-1920 days. Emsenhuber et al. (2020b) present results for
different numbers of $0.01{M_{\oplus}}$ planetary embryos in the initial
conditions, and they consider their runs with 100 planetary embryos to be
their “nominal” or most realistic results. However, they also indicate that
their runs with 10, 20, and 50 planetary embryos yield results that are
similar to the 100 embryo results, and Fernandes et al. (2019) compared the
Mayor et al. (2011) CORALIE/HARPS sample to an early run of the Emsenhuber et
al. (2020b) simulations with 20 planetary embryos. So, we have presented
comparisons to both the Emsenhuber et al. (2020b) 20 and 100 embryo
simulations in Figure 3. The Emsenhuber et al. (2020a, b) simulations consider
only solar mass host stars, but this is similar to the Mayor et al. (2011)
sample which is dominated by solar type stars in the mass range
$0.7M_{\odot}\leq M_{\rm host}\leq 1.3M_{\odot}$. The Gaussian desert model
parameters for these models are given in Table 4. Each Emsenhuber et al.
(2020b) simulation includes 1000 simulated host stars, so it is comparable to
the 822 stars in the Mayor et al. (2011) CORALIE/HARPS sample. Our analysis of
the simulated planet sample is done in exactly the same way as our analysis of
the CORALIE/HARPS sample, with orbital inclinations assigned randomly.

Table 4: Exoplanet Desert Power-Law Model Parameters for New Bern Models | Period Ranges for 20 Embryos
---|---
parameter | 3.75–30 days | 30–240 days | 240–1920 days
| best | average | best | average | best | average
$A$ | 0.04216 | $0.16\pm 0.12$ | 6.644 | $11.8\pm 5.4$ | 0.1789 | $0.132\pm 0.058$
$A_{P\rm m}$ | 0.08073 | $0.097\pm 0.024$ | 0.3563 | $0.302\pm 0.075$ | 0.3754 | $0.255\pm 0.098$
$B$ | 0.9913 | $0.84\pm 0.10$ | 0.9934 | $0.970\pm 0.026$ | 0.9591 | $0.87\pm 0.011$
$\log M_{\rm des}$ | 1.726 | $1.83\pm 0.13$ | 1.707 | $1.726\pm 0.070$ | 1.509 | $1.55\pm 0.18$
$\sigma_{\rm des}$ | 0.245 | $0.41\pm 0.17$ | 0.850 | $0.80\pm 0.12$ | 0.960 | $0.74\pm 0.17$
$p$ | $0.124$ | $0.03\pm 0.17$ | 0.927 | $1.12\pm 0.19$ | $-0.687$ | $-0.63\pm 0.14$
$m$ | $-0.887$ | $-0.859\pm 0.067$ | $-0.621$ | $-0.583\pm 0.073$ | $-0.334$ | $-0.22\pm 0.11$
| Period Ranges for 100 Embryos
parameter | 3.75–30 days | 30–240 days | 240–1920 days
| best | average | best | average | best | average
$A$ | 0.00580 | $0.96\pm 0.65$ | 12.86 | $8.9\pm 2.3$ | 0.3044 | $0.28\pm 0.10$
$A_{P\rm m}$ | 0.0877 | $0.074\pm 0.028$ | 0.3941 | $0.274\pm 0.066$ | 0.4348 | $0.39\pm 0.12$
$B$ | 0.9993 | $0.84\pm 0.17$ | 0.9820 | $0.946\pm 0.028$ | 0.9996 | $0.97\pm 0.028$
$\log M_{\rm des}$ | 2.608 | $2.48\pm 0.24$ | 1.727 | $1.758\pm 0.082$ | 1.708 | $1.696\pm 0.082$
$\sigma_{\rm des}$ | 0.848 | $0.75\pm 0.20$ | 1.104 | $0.80\pm 0.12$ | 0.757 | $0.74\pm 0.14$
$p$ | $-0.081$ | $0.42\pm 0.20$ | 1.104 | $1.10\pm 0.11$ | $-0.331$ | $-0.34\pm 0.14$
$m$ | $-0.896$ | $-0.918\pm 0.087$ | $-0.733$ | $-0.661\pm 0.078$ | $-0.388$ | $-0.32\pm 0.10$

In contrast to the case with the CORALIE/HARPS sample, the Emsenhuber et al.
(2020b) simulations show very significant desert features for all period
ranges with maximum depths ($B$ values) ranging from $0.838\pm 0.168$ to
$0.974\pm 0.028$. The widths of these features (in $\log M$) are also larger
with $\sigma_{\rm des}$ values ranging from $0.414\pm 0.167$ to $0.804\pm
0.116$. The one peculiar result is for the 100 embryo simulations in the
3.75–30 day period range, where center of the desert occurs above
$100{M_{\oplus}}$ at $\log M_{\rm des}=2.475\pm 0.236$, but this is largely
because this data set has only 2 of the 98 simulated planets with $M\sin
i>100{M_{\oplus}}$.

Since the Mayor et al. (2011) CORALIE/HARPS sample shows no evidence for a
planet desert, it is not surprising that the population synthesis results do
not match the CORALIE/HARPS data. The simulations indicate the desert features
that are not seen in the data, and for the 30-240 day and 240–1920 day period
ranges the simulations predict a large excess of $\sim 1000{M_{\oplus}}$
planets that are not seen in the observed sample. These discrepancies were
also noted by Emsenhuber et al. (2020c). However, they did their comparison
using observed and simulated histograms of detected planets instead of the
detection efficiency corrected histograms that we present here. The problem
with their approach is that a steep rise in planet occurrence at $M\sin
i\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}30{M_{\oplus}}$ with a flat
distribution above $30{M_{\oplus}}$ can be converted into an apparent desert
due to the lower detection efficiency for 30–100$\,{M_{\oplus}}$ planets
compared to $>100{M_{\oplus}}$ planets.

##  6 Discussion and Conclusions

Our analysis has demonstrated that there is no evidence for the predicted sub-
Saturn-mass planetary desert (Ida & Lin, 2004a) in the Mayor et al. (2011)
radial velocity sample. Previously, the Suzuki et al. (2016) study
demonstrated that there was no such planetary desert in the exoplanet mass
ratio distribution of planets found by microlensing. These results are
complementary because microlensing is very sensitive to low-mass planets
(Bennett & Rhie, 1996) orbiting beyond the snow line (Gould & Loeb, 1992),
while the radial velocity method is currently not sensitive to planets near
the lower mass limit of the predicted desert, beyond the snow line. However,
the sensitivity of the radial velocity method improves for shorter period
orbits, and the sensitivity of the Mayor et al. (2011) sample extends below
the predicted planet desert lower mass limit of $\sim 10{M_{\oplus}}$ for
orbital periods in our first and second period bins with $P<240\,$days. In our
third period bin, at $240\,{\rm days}<P<1920\,$days, the sensitivity extends
down to only $\sim 12{M_{\oplus}}$, but this period range overlaps with the
sensitivity of the microlensing method. So, the combination of the radial
velocity analysis presented here and the earlier microlensing work (Suzuki et
al., 2016) appears to rule out this predicted planet desert over a wide range
of planetary orbits, out to $\sim 10\,$AU.

While the population synthesis models are often referred to as theoretical
models, it might be more reasonable to consider them to be methods that
attempt to interpret observational constraints in terms of the many
complicated physical processes that are thought to be involved in the
formation of planets. As mentioned in the introduction, Suzuki et al. (2018)
have presented a number of possible modifications to the Ida & Lin (2004a, b,
2005, 2008a, 2008b) and Bern group (Mordasini et al., 2009a, b, 2015). These
include a number of processes that could slow gas accretion or terminate it
well before the planet reaches $100\,{M_{\oplus}}$. These include heating of
the gaseous envelope by the accretion of planetesimals, a low disk viscosity,
low disk scale height, or early formation of a circumplanetary disk. Runaway
growth is often thought to terminate at lower masses at wider orbits, so
gravitational interactions between the planets could transport lower mass
planets to Jupiter-like orbits from the wider orbits where they formed.

Another possibility is that individual planetary systems might form in
protoplanetary disks that allow for runaway gas accretion growth that
disfavors an intermediate final planetary mass range that could be considered
a desert. However, if the properties of these disks have a significant
variation, then it could be that these disfavored planetary mass ranges do not
line up to create a desert in the combined mass or mass ratio distribution of
a large sample of planetary systems. So, it could be that the planet desert
feature is simply smoothed out by large variations in the properties of
protoplanetary disks. A wide variation in assumed planet accretion rates plays
a role in theoretical framework developed by Adams et al. (2021). This
approach is more phenomenological than the population synthesis approach, and
does not seek to describe most of the detailed physical processes involved in
the formation of planets. One of the basic assumptions of the Adams et al.
(2021) is that the exoplanet mass function has an approximate power-law form
that appears to describe the planet distribution above $\sim 10{M_{\oplus}}$
according to studies using a variety of different methods (Cumming et al.,
2008; Suzuki et al., 2016; Schlaufman, 2018; Nielsen et al., 2019; Wagner et
al., 2019). Adams et al. (2021) find that approximately power-law mass
functions can be achieved through a combination of exponential decay rate for
protoplanetary disk gas and a gas accretion rate roughly proportional to
planet mass if they also include a random gas accretion efficiency factor that
can vary significantly between different systems.

One issue that has yet to be explored in much detail is the host mass
dependence of the wide orbit planet distribution. It is possible that the
predicted planet desert does exist for some particular range of host masses,
although the results presented here tend to disfavor this idea since the Mayor
et al. (2011) sample consists of FGK stars, while the Suzuki et al. (2016)
sample extends to lower mass host stars. Nevertheless, Kepler observations
have established that planets in short period orbits have a higher occurrence
rate around low-mass stars than around solar type stars (Mulders et al.,
2015), so it would not be surprising to find a host mass dependence of the
wide orbit planet occurrence rate. However, most of the existing demographic
studies do not allow this. The exiting microlensing studies use only the
planet-star mass ratio (Gould et al., 2010; Cassan et al., 2012; Suzuki et
al., 2016), and the largest of the radial velocity studies (Mayor et al.,
2011) does not provide enough information for such a study. The recent radial
velocity analyses of Wittenmyer et al. (2020) and (Fulton et al., 2021) do
provide the individual detection efficiencies for both the stars with and
without detected planets, which is what is needed. However the Wittenmyer et
al. (2020) study is smaller, and both surveys are not as sensitive to low
planet masses as the Mayor et al. (2011) study. So, this data set is less
suitable for a study of the predicted sub-Saturn-mass desert for wide orbit
planets.

While most of the parameters of planetary microlensing events can only be
determined with data taken during the event, it is possible to determine the
masses of the host stars and planets and their separation with high angular
resolution follow-up observations from the Hubble Space Telescope (Bennett et
al., 2015; Bhattacharya et al., 2017) or adaptive optics (AO) observations on
8–10m class telescopes. The Keck telescopes have been the most successful at
this, due an AO system that is more effective with faint guide stars (Batista
et al., 2015; Beaulieu et al., 2016; Bhattacharya et al., 2018, 2020; Bennett
et al., 2020; Vandorou et al., 2020; Terry et al., 2021). A good AO correction
is important for the interpretation of these follow-up observations because it
is generally necessary to confirm the identification of the host star by
confirming the relative proper motion between the candidate host star and the
background source star (Koshimoto et al., 2020). The MOA Collaboration is in
the process of extending the Suzuki et al. (2016) statistical sample, while
the Keck and Hubble follow-up program is expanding with observations and host
and planet mass determinations for many more microlens planetary systems in
this extended statistical sample. New radial velocity results have recently
been reported from the California Legacy Survey (Fulton et al., 2021). So, we
expect significant developments in the understanding of wide orbit exoplanet
demographics in the next few years. Ultimately, the Nancy Grace Roman
Telescope’s Galactic Exoplanet survey (Bennett & Rhie, 2002; Penny et al.,
2018) will perform a much more comprehensive statistical study of the
occurrence rate of wide orbit planets at semi-major axes of $a\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}0.5\,$AU with masses down to the
mass of Mars ($0.1\,{M_{\oplus}}$).

We thank Fred Adams, B. Scott Gaudi, Gijs Mulders, and Ilaria Pascucci for
helpful suggestions and comments on a draft version of this paper. DPB and CR
were supported by NASA through grant NASA-80NSSC18K0274 and award number
80GSFC17M0002. CR is supported by the ANR project COLD-WORLDS of the French
_Agence Nationale de la Recherche_ with the reference ANR-18-CE31-0002. The
results reported herein have benefited from collaborations and/or information
exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research
coordination network sponsored by NASA’s Science Mission Directorate.

##  Appendix A Planet Deserts in Kepler Data?

There have been extensive discussions in the literature of planet deserts seen
in Kepler data. As mentioned in the introduction, there is the “hot Neptune
desert”, also known as the the “sub-Jovian desert” (Szabó & Kiss, 2011;
Lundkvist et al., 2016; Mazeh et al., 2016; Owen, 2019). This was an
observational discovery, based primarily on Kepler data, and it is thought to
be due to heating of the atmospheres of planets that come very close to their
host stars. There is also the Fulton gap (Fulton et al., 2017), which is
sometimes referred to as the ”sub-Neptune” desert and is thought to be caused
by photo-evaporation. These deserts occur for only very short period orbits
and are thought to be caused by heating or photoevaporation due to the close
proximity of the host star, so they are not relevant to the desert that the
topic of this paper, which is predicted by the runaway gas accretion scenario
at much wider orbits.

There has also been some, mostly informal, suggestions that the Kepler data
shows a sub-Jupiter radius desert for cooler planets at wider separations. For
example, this shows up, somewhat weakly, in Figure 7 of Thompson et al. (2018)
at orbital periods $\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}20\,$days, and it also can be seen
in Figure 3 of Berger et al. (2020) at an incident flux of $\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}10\times$ that received by the
Earth. In this section, we show that this sub-Jupiter radius gap is likely to
be real, but that it does not imply a sub-Jupiter mass gap.

Figure 4:  The predicted radius distribution for planets in the three period
bins shown in Figure 2 assuming the best fit power-law models with no desert
feature. The small planet desert feature at a radius of $R\simeq 9R_{\oplus}$
is due to the mass-radius relation that reflects the pile-up of planets at
radii of $\sim 14\,R_{\oplus}$ due to the dominance of electron degeneracy
pressure for objects in the mass range of $0.4M_{\rm Jup}$ to $0.08M_{\odot}$,
a factor of $\lower 2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}200$ in mass.

In order to relate the masses and radii of planets, we use the mass-radius
relation of Chen & Kipping (2017), and to describe the mass distribution of
exoplanets we employ the best fit pure power-law models from Table 2 for our
three different period ranges: 3.75–30 days, 30–240 days and 240-1920 days.
The results of this calculation are shown in Table 4. This shows a clear
minimum in the planetary radius distribution at $R\simeq 9R_{\oplus}$ and a
maximum at $R\simeq 14R_{\oplus}$. This feature is dominated by planets in the
longest period bin of $240\,{\rm days}<P\leq 1920\,{\rm days}$ because the
assumed long period mass function has the shallowest mass dependence: $M^{m}$
with $m=-0.367$. Note that this exercise is intended to show how a smooth mass
function can still imply a “planet radius desert”, and we do not intend to
claim that these exoplanet mass functions can provide a good description of
the Kepler data with its much larger exoplanet sample.

The reason that the smooth, monotonic mass function is converted into a radius
function with a dip and a peak is clear from Figure 3 of Chen & Kipping
(2017), which displays their mass-radius relation. This relation is nearly
flat between masses of $0.4M_{\rm Jup}$ to $0.08M_{\odot}$, a factor of 209 in
mass. This is a well known feature caused by the fact that gas giant planets,
brown dwarfs and the lowest mass stars are all largely supported by electron
degeneracy pressure (which is also responsible for the decreasing radius with
increasing mass for white dwarfs). Thus, the observed potential sub-Jupiter
desert in the Kepler planetary radius distribution is likely to be explained
by the mass-radius relation rather than a feature in the exoplanet mass
distribution, since the feature appears for a model that has no desert in the
mass distribution.

## References

  * Adams et al. (2021) Adams, F. C., Meyer, M. R., & Adams, A. D. 2021, ApJ, 909, 1
  * Alcock et al. (1996) Alcock, C., Allsman, R. A., Axelrod, T. S., et al. 1996, ApJ, 461, 84
  * Batista et al. (2015) Batista, V., Beaulieu, J.-P., Bennett, D. P., et al. 2015, ApJ, 808, 170
  * Beaulieu et al. (2016) Beaulieu, J.-P., Bennett, D. P., Batista, V., et al. 2016, ApJ, 824, 83
  * Bennett et al. (2020) Bennett, D. P., Bhattacharya, A., Beaulieu, J.-P., et al. 2020, AJ, 159, 68
  * Bennett & Rhie (1996) Bennett, D. P., & Rhie, S. H. 1996, ApJ, 472, 660
  * Bennett & Rhie (2002) Bennett, D. P., & Rhie, S. H. 2002, ApJ, 574, 985
  * Bennett et al. (2015) Bennett, D. P., Bhattacharya, A., Anderson, J., et al. 2015, ApJ, 808, 169
  * Berger et al. (2020) Berger, T. A., Huber, D., Gaidos, E., et al. 2020, AJ, 160, 108
  * Bhattacharya et al. (2018) Bhattacharya, A., Beaulieu, J.-P., Bennett, D. P., et al. 2018, AJ, 156, 289
  * Bhattacharya et al. (2017) Bhattacharya, A., Bennett, D. P., Anderson, J., et al. 2017, AJ, 154, 59
  * Bhattacharya et al. (2020) Bhattacharya, A., Bennett, D. P., Beaulieu, J. P., et al. 2020, arXiv:2009.02329
  * Cassan et al. (2012) Cassan, A., Kubas, D., Beaulieu, J.-P., et al. 2012, Nature, 481, 167
  * Chen & Kipping (2017) Chen, J. & Kipping, D. 2017, ApJ, 834, 17
  * Cumming et al. (2008) Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531
  * Dobbs-Dixon et al. (2007) Dobbs-Dixon, I., Li, S. L., & Lin, D. N. C. 2007, ApJ, 660, 791
  * Emsenhuber et al. (2020a) Emsenhuber, A., Mordasini, C., Burn, R., et al. 2020a, arXiv:2007.05561
  * Emsenhuber et al. (2020b) Emsenhuber, A., Mordasini, C., Burn, R., et al. 2020b, arXiv:2007.05562
  * Emsenhuber et al. (2020c) Emsenhuber, A., Mordasini, C., Mayor, M., et al. 2020c, European Planetary Science Congress
  * Fernandes et al. (2019) Fernandes, R. B., Mulders, G. D., Pascucci, I., et al. 2019, ApJ, 874, 81
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109
  * Fulton et al. (2021) Fulton, B. J., Rosenthal, L. J., Hirsch, L. A., et al. 2021, ApJS, 255, 14. doi:10.3847/1538-4365/abfcc1
  * Fung et al. (2014) Fung, J., Shi, J.-M., & Chiang, E. 2014, ApJ, 782, 88
  * Garaud & Lin (2007) Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606
  * Gaudi et al. (2020) Gaudi, B. S., Christiansen, J. L., & Meyer, M. R. 2020, arXiv:2011.04703
  * Gould & Loeb (1992) Gould, A., & Loeb, A. 1992, ApJ, 396, 104
  * Gould et al. (2010) Gould, A., Dong, S., Gaudi, B. S., et al. 2010, ApJ, 720, 1073
  * Ida & Lin (2004a) Ida, S., & Lin, D. N. C. 2004a, ApJ, 604, 388
  * Ida & Lin (2004b) Ida, S., & Lin, D. N. C. 2004b, ApJ, 616, 567
  * Ida & Lin (2005) Ida, S., & Lin, D. N. C. 2005, ApJ, 626, 1045
  * Ida & Lin (2008a) Ida, S., & Lin, D. N. C. 2008a, ApJ, 673, 487
  * Ida & Lin (2008b) Ida, S., & Lin, D. N. C. 2008b, ApJ, 685, 584
  * Jung et al. (2019) Jung, Y. K., Gould, A., Zang, W., et al. 2019, AJ, 157, 72
  * Koshimoto et al. (2020) Koshimoto, N., Bennett, D. P., & Suzuki, D. 2020, AJ, 159, 268 a
  * Koshimoto et al. (2017) Koshimoto, N., Shvartzvald, Y., Bennett, D. P., et al. 2017, AJ, 154, 3
  * Lundkvist et al. (2016) Lundkvist, M. S., Kjeldsen, H., Albrecht, S., et al. 2016, Nature Communications, 7, 11201
  * Mayor et al. (2011) Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv:1109.2497
  * Mazeh et al. (2016) Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75
  * Mordasini et al. (2009a) Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139
  * Mordasini et al. (2009b) Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1161
  * Mordasini et al. (2015) Mordasini, C., Mollière, P., Dittkrist, K.-M., Jin, S., & Alibert, Y. 2015, International Journal of Astrobiology, 14, 201
  * Mulders et al. (2015) Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 798, 112
  * Nayakshin et al. (2019) Nayakshin, S., Dipierro, G., & Szulágyi, J. 2019, MNRAS, 488, L12
  * Nielsen et al. (2019) Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13
  * Owen (2019) Owen, J. E. 2019, Annual Review of Earth and Planetary Sciences, 47, 67
  * Penny et al. (2018) Penny, M. T., Gaudi, B. S., Kerins, E., et al. 2018, arXiv:1808.02490
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
  * Raymond & Izidoro (2017) Raymond, S. N., & Izidoro, A. 2017, Icarus, 297, 134
  * Schlaufman (2018) Schlaufman, K. C. 2018, ApJ, 853, 37
  * Suzuki et al. (2018) Suzuki, D., Bennett, D. P., Ida, S., et al. 2018, ApJ, 869, L34
  * Suzuki et al. (2016) Suzuki, D., Bennett, D. P., Sumi, T., et al. 2016, ApJ, 833, 145
  * Szabó & Kiss (2011) Szabó, G. M. & Kiss, L. L. 2011, ApJ, 727, L44
  * Szulágyi et al. (2014) Szulágyi, J., Morbidelli, A., Crida, A., et al. 2014, ApJ, 782, 65
  * Terry et al. (2021) Terry, S. K., Bhattacharya, A., Bennett, D. P., et al. 2021, AJ, 161, 54.
  * Thompson et al. (2018) Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, ApJS, 235, 38
  * Trotta (2008) Trotta, R., 2008, Contemp. Phys., 49, 71
  * Udalski et al. (2018) Udalski, A., Ryu, Y.-H., Sajadian, S., et al. 2018, Acta Astron., 68, 1
  * Vandorou et al. (2020) Vandorou, A., Bennett, D. P., Beaulieu, J.-P., et al. 2020, AJ, 160, 121
  * Wagner et al. (2019) Wagner, K., Apai, D., & Kratter, K. M. 2019, ApJ, 877, 46
  * Wittenmyer et al. (2020) Wittenmyer, R. A., Wang, S., Horner, J., et al. 2020, MNRAS, 492, 377