capbtabboxtable[][3.5in] capbfigboxfigure[][2.3in]  FERMILAB-PUB-22-391-T

# Constraining the Local Pulsar Population with the Cosmic-Ray Positron
Fraction

Olivia Meredith Bitter  and Dan Hooper

###### Abstract

Observations of the TeV halos associated with nearby pulsars indicate that
these objects inject significant fluxes of very high-energy electron-positrons
pairs into the interstellar medium (ISM), thereby likely providing the
dominant contribution to the cosmic-ray positron flux. In this paper, we use
the cosmic-ray positron fraction as measured by the AMS-02 Collaboration to
constrain the characteristics of the local pulsar population. For reasonable
model parameters, we find that we can obtain good agreement with the measured
positron fraction up to energies of $E_{e}\sim 300\,{\rm GeV}$. At higher
energies, the positron fraction is dominated by a small number of pulsars,
making it difficult to reliably predict the shape of the expected positron
fraction. The low-energy positron spectrum supports the conclusion that
pulsars typically transfer approximately $\eta\sim 5-20\%$ of their total
spindown power into the production of very high-energy electron-positron
pairs, producing a spectrum of such particles with a hard spectral index,
$\alpha\sim 1.5-1.7$. Such pulsars typically spindown on a timescale on the
order of $\tau\sim 10^{4}\,{\rm years}$. Our best fits were obtained for
models in which the radio and gamma-ray beams from pulsars are detectable to
28% and 62% of surrounding observers, respectively.

##  1 Introduction

Recent observations by the HAWC Collaboration [8, 1, 2, 7] indicate that young
and middle-aged pulsars (and also possibly millisecond pulsars [15, 18])
produce significant fluxes of very high-energy gamma rays through the inverse
Compton scattering of electrons and positrons [19, 23, 29, 31]. This emission
leads to the appearance of “TeV halos” surrounding most, if not all, pulsars.
Furthermore, the observed characteristics of these TeV halos support the
conclusion that these sources are likely responsible for the cosmic-ray
positron excess [19, 16] reported by the AMS-02 Collaboration [32], and
previously by PAMELA [3], AMS-01 [4], and HEAT [9] (for earlier work, see
Refs. [10, 28]).

In previous work, discussions of TeV halos in relation to the cosmic-ray
positron fraction tended to focus on the contributions from a small number of
nearby and high-spindown power pulsars, such as Geminga and Monogem [19, 22,
14]. It has since become clear, however, that the positrons detected by AMS do
not originate from only a few pulsars. Instead, the cosmic-ray positron
fraction receives substantial contributions from a significant number of such
objects found throughout the surrounding kiloparsecs of the Milky Way [12,
27].

In this paper, we revisit TeV halos as a source of cosmic-ray positrons, using
the measured characteristics of the cosmic-ray positron fraction to constrain
the properties of the local pulsar population. Toward this goal, we have
calculated the cosmic-ray positron flux from the known pulsars contained
within the Australia Telescope National Facility (ATNF) pulsar catalog [25],
taking into account the effects of diffusion, inverse Compton scattering, and
synchrotron. We then accounted for those nearby pulsars which have thus far
gone undetected, employing a Monte Carlo to simulate this population of
sources. Using this Monte Carlo, we have calculated the contribution to the
cosmic-ray positron flux from these pulsars, and compared these results to the
positron fraction as reported by AMS-02, varying a variety of parameters
including the spectral shape of the injected positrons, the timescale for
pulsar spindown, the efficiency for electron-positron production, the
magnitude of the magnetic field in the interstellar medium (ISM), and the
duration of these objects’ pulsar wind nebula phase.

For reasonable choices of these parameters, we have found that it is possible
to obtain good agreement with the measured positron fraction, at least up to
energies of $E_{e}\sim 300\,{\rm GeV}$. This result supports the conclusion
that pulsars and their associated TeV halos are likely to be the dominant
source of the observed cosmic-ray positron spectrum. At energies above a few
hundred GeV, the positron fraction is dominated by a small number of pulsars,
making it difficult to reliably predict the shape of the expected positron
fraction or to use this high-energy information to draw reliable conclusions
pertaining to the Milky Way’s broader pulsar population. The low-energy
positron spectrum, on the other hand, supports the conclusion that pulsars
typically transfer about $\eta\sim 5-20\%$ of their total spindown power into
the production of very high-energy electron-positron pairs, and produces a
spectrum of such particles with a hard spectral index of $\alpha\sim 1.5-1.7$.

##  2 Pulsars as a Source of High-Energy Positrons

Once a high-energy electron or positron has escaped from a TeV halo, their
transport through the ISM is described by the following diffusion-loss
equation:

$\displaystyle\frac{\partial}{\partial
t}\frac{dn_{e}}{dE_{e}}(E_{e},r,t)=\vec{\mathbf{\nabla}}\cdot\bigg{[}D(E_{e})\vec{\mathbf{\nabla}}\frac{dn_{e}}{dE_{e}}(E_{e},r,t)\bigg{]}+\hskip
5.69054pt\frac{\partial}{\partial
E_{e}}\bigg{[}\frac{dE_{e}}{dt}(E_{e},r)\frac{dn_{e}}{dE_{e}}(E_{e},r,t)\bigg{]}+\delta(r)Q(E_{e},t),$
(2.1)

where $r$ is the distance from the pulsar, $Q$ characterizes the spectrum and
time profile of the injected particles, $dn_{e}/dE_{e}$ is the electron and/or
positron number density per unit energy, and $E_{e}$ is the energy of the
electrons and positrons. The normalization and energy dependence of the
diffusion coefficient, $D$, is constrained by the measurements of various
secondary-to-primary ratios in the cosmic-ray spectrum. Throughout this study
we will adopt a parameterization given by $D=(2\times 10^{28}\,{\rm
cm}^{2}/{\rm s})\,(E_{e}/{\rm GeV})^{0.4}$ [33, 17]. The energy losses of
high-energy electrons and positrons are dominated by a combination of inverse
Compton scattering and synchrotron processes, resulting in the following rate:

$\displaystyle-\frac{dE_{e}}{dt}(E_{e},r)$
$\displaystyle=\sum_{i}\frac{4}{3}\sigma_{T}\rho_{i}(r)S_{i}(E_{e})\bigg{(}\frac{E_{e}}{m_{e}}\bigg{)}^{2}+\frac{4}{3}\sigma_{T}\rho_{\rm
mag}(r)\bigg{(}\frac{E_{e}}{m_{e}}\bigg{)}^{2}$ (2.2) $\displaystyle\approx
1.02\times 10^{-16}\,{\rm GeV/s}\bigg{(}\frac{E_{e}}{\rm
GeV}\bigg{)}^{2}\times\bigg{[}\sum_{i}\frac{\rho_{i}(r)}{\rm
eV/cm^{3}}S_{i}(E_{e})+0.224\,\bigg{(}\frac{B(r)}{3\,\mu{\rm
G}}\bigg{)}^{2}\bigg{]},$

where the $\sigma_{T}$ is the Thomson cross section, $\rho_{\rm mag}\approx
0.224\,{\rm eV/cm}^{3}\times(B_{\rm ISM}/3\,\mu{\rm G})^{2}$ is the energy
density of the magnetic field, and the sum is performed over different
components of the radiation field, consisting of the cosmic microwave
background (CMB), infrared emission (IR), starlight (star), and ultraviolet
emission (UV). Lastly, the quantity $S_{i}(E_{e})$ is a factor that accounts
for the Klein-Nishina suppression of the inverse Compton scattering cross
section, and is given by

$\displaystyle
S_{i}(E_{e})\approx\frac{45m^{2}_{e}/64\pi^{2}T^{2}_{i}}{(45m^{2}_{e}/64\pi^{2}T^{2}_{i})+(E^{2}_{e}/m^{2}_{e})}.$
(2.3)

Throughout our analysis, we will follow Ref. [19] in adopting the following
parameters for the radiation model: $\rho_{\rm CMB}=0.260\,{\rm eV/cm}^{3}$,
$\rho_{\rm IR}=0.60\,{\rm eV/cm}^{3}$, $\rho_{\rm star}=0.60\,{\rm
eV/cm}^{3}$, $\rho_{\rm UV}=0.10\,{\rm eV/cm}^{3}$, $T_{\rm CMB}=2.7\,{\rm
K}$, $T_{\rm IR}=20\,{\rm K}$, $T_{\rm star}=5000\,{\rm K}$, and $T_{\rm
UV}=20,000\,{\rm K}$. At very high energies ($E_{e}\gg m^{2}_{e}/T$), inverse
Compton scattering is strongly Klein-Nishina suppressed. This is particularly
important for the scattering of cosmic-ray electrons and positrons with UV and
starlight photons.

For the time profile of the source term in Eq. 2.1, we assume that the rate at
which high-energy positrons and electrons are injected from a pulsar into the
ISM is proportional to the rate at which that pulsar is losing rotational
kinetic energy (i.e. its spindown power). For the case of magnetic-dipole
breaking, a pulsar’s spindown power is given by

$\displaystyle\dot{E}(t)$
$\displaystyle=-\frac{8\pi^{4}B^{2}R^{6}}{3c^{3}P(t)^{4}}$ (2.4)
$\displaystyle\approx 1.0\times 10^{35}\,{\rm
erg/s}\times\bigg{(}\frac{B}{1.6\,\times\hskip 2.84526pt10^{12}\,{\rm
G}}\bigg{)}^{2}\bigg{(}\frac{R}{15\,{\rm
km}}\bigg{)}^{6}\bigg{(}\frac{0.23\,{\rm sec}}{P(t)}\bigg{)}^{4},$

where $B$, $R$, and $P$ are the pulsar’s magnetic field, radius, and period.
Solving this differential equation, we obtain the following expression for the
time evolution of the pulsar’s period:

$\displaystyle P(t)=P_{0}\,\bigg{(}1+\frac{t}{\tau}\bigg{)}^{\frac{1}{2}},$
(2.5)

where $P_{0}$ is the pulsar’s initial period and $\tau$ is the spindown
timescale, given by

$\displaystyle\tau$
$\displaystyle=\frac{3c^{3}IP_{0}^{2}}{4\pi^{2}B^{2}R^{6}}$ (2.6)
$\displaystyle\approx 9.1\,\times 10^{3}\hskip 2.84526pt{\rm
yr}\times\bigg{(}\frac{1.6\times 10^{12}\,{\rm
G}}{B}\bigg{)}^{2}\bigg{(}\frac{M}{1.4\,M_{\odot}}\bigg{)}\bigg{(}\frac{15\,{\rm
km}}{R}\bigg{)}^{4}\bigg{(}\frac{P_{0}}{0.040\,{\rm sec}}\bigg{)}^{2},$

where $I$ and $M$ are the moment of inertia and mass of the neutron star.

During a young neutron star’s pulsar wind nebulae phase, positrons and
electrons can be confined to the region surrounding the object. With this in
mind, we set the positron and electron injection rate to zero prior to a time,
$t_{\rm PWN}$, which we expect to be on the order of $\mathcal{O}(10^{4}\,{\rm
yr})$.

To obtain the total source term, $Q$, we multiply the spindown power by a
constant efficiency factor, $\eta$, and parameterize the spectral shape using
a power-law with an exponential cutoff:

$\displaystyle
Q(E_{e},t)=\begin{cases}0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,t<t_{\rm
PWN}\\\\[10.0pt]
\eta\,\dot{E}(t)\,Q_{0}\,E^{-\alpha}\,e^{-E_{e}/E_{c}},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,t>t_{\rm
PWN}\end{cases}$ (2.7)

where the normalization coefficient, $Q_{0}$, is set such that $\eta$ is equal
to the fraction of the pulsar’s spindown power that goes into the production
of positrons and electrons with $E_{e}>10\,{\rm GeV}$ (for all $t>t_{\rm
PWN}$):

$\displaystyle\eta=\frac{\int_{10\,{\rm
GeV}}Q(E_{e},t)\,E_{e}\,dE_{e}}{\dot{E}(t)}.$ (2.8)

Putting this together, we solve the diffusion-loss equation as described in
Ref. [19] to find the spectrum of electrons and positrons from a given pulsar
at the location of the Solar System. The resulting positron spectrum depends
on many parameters, including the distance to the pulsar, the pulsar’s age,
spindown timescale, initial period, and efficiency, as well as the duration of
the pulsar wind nebula phase, the magnitude of the magnetic field in the ISM,
and the spectral index and cutoff energy of the injected spectrum.

Once we have calculated the spectrum of positrons and electrons from the
pulsars under consideration, we can combine this information with the spectrum
of additional primary electrons (such as those from supernova remnants), as
well as the spectra of secondary positrons and electrons that are produced
through the propagation of primary cosmic rays. Taken together, this can be
used to calculate the cosmic-ray positron fraction, defined as the positron
flux divided by the combined flux of electrons and positrons. In calculating
the positron fraction, we use the secondary spectra as obtained using the
publicly available code GALPROP111http://galprop.stanford.edu/., and fix the
primary spectrum using the measured cosmic-ray electron spectrum, assuming
that pulsars produce an identical spectrum of positrons and electrons [26,
33]. To account for the finite energy resolution of AMS-02, we convolve the
spectrum we calculate in each case with a Gaussian of width
$\sigma=0.02\,E_{e}$ [32].

##  3 Positrons From Known Pulsars

We begin our analysis by calculating the contribution to the cosmic-ray
positron flux from known pulsars; namely those contained within the Australia
Telescope National Facility (ATNF) pulsar catalog [25]. To this end, we
consider the 103 ANTF pulsars that are located within 3 kpc of the Solar
System and that are younger than $\tau_{c}<10^{6}\,{\rm years}$.222The
characteristic age of a pulsar is defined by $\tau_{c}\equiv P/2\dot{P}$. For
the case of magnetic dipole braking, this is related to a pulsar’s true age,
$t$, and spindown timescale according to $t=\tau_{c}-\tau$.

Figure 1: The cosmic-ray positron fraction individually calculated from each
of the 103 pulsars in the ATNF catalog that are located within 3 kpc of the
Solar System and younger than $\tau_{c}<10^{6}\,{\rm years}$ (thin colored
lines), and from the sum of all 103 of these sources (thick black line). For
the color key, see Table 1. In this figure, we have adopted the following
parameters: $\alpha=1.5$, $E_{c}=50\,{\rm TeV}$, $\eta=0.15$,
$\tau=10^{4}\,{\rm yr}$, $B_{\rm ISM}=3\,\mu{\rm G}$, and $t_{\rm PWN}=3\times
10^{4}\,{\rm yr}$. For each individual pulsar, we adopt values of $\tau_{c}$,
$r$, and $P$ as reported in the ATNF catalog.

In Fig. 1, we plot the contribution to the cosmic-ray positron fraction from
each of these 103 pulsars (thin colored lines), and from the sum of all 103 of
these sources (thick black line), and compare this to the spectrum measured by
the AMS-02 Collaboration [5, 6, 13]. Here, we have adopted the following
parameters: $\alpha=1.5$, $E_{c}=50\,{\rm TeV}$, $\eta=0.15$,
$\tau=10^{4}\,{\rm yr}$, $B_{\rm ISM}=3\,\mu{\rm G}$, and $t_{\rm PWN}=3\times
10^{4}\,{\rm yr}$. For each individual pulsar, we adopt values of $\tau_{c}$,
$r$, and $P$ as reported in the ATNF catalog (see Table 1). Keep in mind that
this figure only reflects the contributions to the cosmic-ray positron flux
from those pulsars contained in the ATNF catalog, which is largely limited to
relatively young sources with radio or gamma-ray beams pointed in the
direction of the Solar System.

Pulsar Name | $P\,({\rm s})$ | $d\,({\rm kpc})$ | $\tau_{c}\,({\rm yr})$ | Color (Fig. 1)
---|---|---|---|---
J1302-6350 | 0.048 | 2.63 | $3.22\times 10^{5}$ | Light Blue/Cyan
J1057-5226 | 0.197 | 0.093 | $5.25\times 10^{5}$ | Orange
J0117+5914 | 0.101 | 1.77 | $2.65\times 10^{5}$ | Red
J0633+1746 | 0.237 | 0.190 | $3.32\times 10^{5}$ | Light Green
J2032+4127 | 0.143 | 1.33 | $1.91\times 10^{5}$ | Violet
J0908-4913 | 0.107 | 1.00 | $1.02\times 10^{5}$ | Yellow
J2030+4415 | 0.227 | 0.720 | $5.45\times 10^{5}$ | Brown
J1846+0919 | 0.226 | 1.53 | $3.50\times 10^{5}$ | Pink
J1745-3040 | 0.367 | 0.200 | $5.36\times 10^{5}$ | Blue
J1530-5327 | 0.279 | 1.12 | $9.34\times 10^{5}$ | Green
Table 1:  The 10 pulsars contained in the ATNF catalog which contribute the
most to the local high-energy cosmic-ray positron flux. For each pulsar, we
provide the reported period, distance, and characteristic age. Also given are
the colors assigned to each of these pulsars in Fig. 1.

The locations of the peaks that appear in the spectrum shown in Fig. 1
correspond to the energies to which very high-energy positrons can cool in an
amount of time equal to the age of the pulsar in question [19]. For this
reason, the older pulsars in this sample tend to exhibit spectral features
which peak at relatively low energies. This opens the possibility that
individual pulsars could generate distinctive spectral features in the
positron spectrum. Conversely, the lack of such features could be used to
constrain the characteristics of such pulsars [11].

It is interesting to note that the Geminga and Monogem pulsars are responsible
for only a modest fraction of the total local positron flux. More
specifically, while the contribution from Geminga (PSR J0633+1746, shown in
light green in Fig. 1) is among the largest from these pulsars, the positron
flux from Monogem does not even make the list of sources presented in Table 1.

##  4 Positrons From Unknown Pulsars

A large fraction of the pulsars in the Milky Way, including many young pulsars
located near the Solar System, have not yet been detected, and thus are not
contained within the ATNF catalog. To estimate the contribution to the cosmic-
ray positron flux from these yet-to-be discovered pulsars, we have employed a
Monte Carlo with parameters tuned to reproduce the observed characteristics of
the ATNF catalog’s pulsar population.

For the locations of the pulsars modelled by our Monte Carlos, we draw from
the following spatial distribution [24]:

$\displaystyle n_{\rm pulsar}\propto R^{2.35}\,e^{-R/1530\,{\rm
pc}}\,e^{-|z|/300\,\rm pc},$ (4.1)

where $R$ and $z$ describe the location of a given pulsar in cylindrical
coordinates (we take the distribution to be symmetric with respect to the
angular coordinate). We take the solar system to be located at R$=8.25\,{\rm
kpc}$ and $z=0$. For the age of each pulsar, we simply draw from a flat
distribution.

In the calculations performed by our Monte Carlo, we focus on those pulsars
which are located within 3 kpc of the Solar System and that are younger than
$10^{6}$ years. While pulsars older than $10^{6}$ years do contribute to the
local positron flux, this portion of the spectrum has a negligible degree of
variation from realization-to-realization. Consequently, we use an average
output from our Monte Carlo to account for the pulsars in this age range. We
also find that pulsars located more than 3 kpc from the Solar System
contribute negligibly to the local positron flux, at a level that is smaller
than the line width used in our plots.

For each realization of our Monte Carlo, we adopt a value for the pulsar birth
rate, $R_{\rm PSR}$, as well as for the spindown timescale, $\tau$. For the
spatial distribution described by Eq. 4.1, 5.2% of the pulsar population is
located within 3 kpc, so in each realization, we simulate a total of $R_{\rm
PSR}\times 0.052\times 10^{6}\simeq 520\times(R_{\rm PSR}/0.01\,{\rm
yr^{-1}})$ pulsars [21]. Notice that for reasonable values of $R_{\rm PSR}$,
the number of these pulsars is significantly larger than the 103 discussed in
the previous section, demonstrating the significant incompleteness of the ATNF
catalog.

In order to compare the output of a given Monte Carlo realization to the
observed population of young and nearby pulsars, we must establish under what
conditions a given pulsar is likely to be detected. In the case of radio
detection, the main limiting factor is whether the radio beam is pointed in
the direction of the Solar System. To model this, we simply adopt a value for
the beaming fraction, $f_{\rm radio}$, which is defined as the probability
that a given pulsar produces radio emission in our direction. For this radio
emission to be detected, we require both that we reside within the geometry of
the radio beam, and that the pulsar’s spindown flux is larger than a minimum
detectable value, $F_{\rm SD}\equiv\dot{E}/d^{2}>F_{\rm radio}^{\rm min}$. In
considering the detection of the pulsed gamma-ray emission from pulsars, the
corresponding beaming fraction is significantly larger, $f_{\gamma}>f_{\rm
radio}$. In this case, however, Fermi is generally only able to detect those
pulsars with a relatively high spindown flux, corresponding to
$F_{\rm\gamma}^{\rm min}>F_{\rm radio}^{\rm min}$.

In Fig. 2, we show an example of one realization of our Monte Carlo, plotting
the spin-down flux and age of each detectable pulsar. In this figure, we have
adopted the following parameters: $\tau=10^{4}\,{\rm yr}$, $P_{0}=0.042\,{\rm
s}$, $R_{\rm PSR}=0.01\,{\rm yr}^{-1}$, $f_{\rm radio}=0.28$,
$f_{\gamma}=0.62$, $F_{\gamma}^{\rm min}$ = 4.22 $\times 10^{42}\,{\rm
erg/{kpc}^{2}/{\rm yr}}$, and $F_{\rm radio}^{\rm min}$ = 2.66 $\times
10^{41}\,{\rm erg/{kpc}^{2}/{\rm yr}}$. The symbols used indicate whether a
given pulsar is deemed to be detectable by their radio or gamma-ray emission.

To constrain the values of $f_{\gamma}$, $f_{\rm radio}$, $F_{\gamma}^{\rm
min}$, and $F_{\rm radio}^{\rm min}$, we compare the observed distribution of
detected pulsar ages to the predictions of our Monte Carlo for various
combinations of parameter values. These distributions are shown in Fig. 3 for
the case of $\tau=10^{4}\,{\rm yr}$ and $R_{\rm PSR}=0.01\,{\rm yr}^{-1}$,
using the best-fit parameter values of $F_{\gamma}^{\rm min}=4.22\times
10^{42}$ ergs/kpc2/yr, $F_{\rm radio}^{\rm min}=2.66\times 10^{42}$
ergs/kpc2/yr, $f_{\gamma}=0.62$, and $f_{\rm radio}=0.28$. The histograms in
these figures represent the number of pulsars detected by Fermi (left frame)
and by radio telescopes (right frame), which are then compared to the results
of our Monte Carlo, averaged over 1000 realizations. For each choice of
parameters, we calculate the Poisson likelihood of obtaining a distribution of
pulsars (in terms of age, and in the three distance increments shown) that is
consistent with the measured histograms. After repeating this exercise for
various parameter values, we identify the best-fit values of $F_{\gamma}^{\rm
min}$, $F_{\rm radio}^{\rm min}$, $f_{\gamma}$, and $f_{\rm radio}$ for each
value of $\tau$ and $R_{\rm PRB}$ as given in Table 2. Throughout the
remainder of our analysis, we will adopt the best-fit values of
$F_{\gamma}^{\rm min}$, $F_{\rm radio}^{\rm min}$, $f_{\gamma}$, and $f_{\rm
radio}$ for any given choice of $\tau$ and $R_{\rm PRB}$.

Figure 2: An example of one realization of our Monte Carlo, showing the spin-
down flux and age of each detectable pulsar. In this figure, we have adopted
the following parameters: $\tau=10^{4}\,{\rm yr}$, $P_{0}=0.042\,{\rm s}$,
$R_{\rm PSR}=0.01\,{\rm yr}^{-1}$, $f_{\rm radio}=0.28$, $f_{\gamma}=0.62$,
$F_{\gamma}^{\rm min}=4.22\times 10^{42}\,{\rm erg/kpc}^{2}/{\rm yr}$, and
$F_{\rm radio}^{\rm min}=2.66\times 10^{41}\,{\rm erg/kpc}^{2}/{\rm yr}$. The
symbols used reflect whether a given pulsar is detectable by their radio or
gamma-ray emission.

Figure 3: The distribution of ATNF pulsars detected in gamma-rays (left) and in radio emission (right), according to their age and distance from Earth. The curves represent the distributions predicted by our Monte Carlo, averaged over 1000 realizations, for the case of $\tau=10^{4}\,{\rm yr}$, $R_{\rm PSR}=0.01\,{\rm yr}^{-1}$, $F_{\gamma}^{\rm min}=4.22\times 10^{42}$ ergs/kpc2/yr, $F_{\rm radio}^{\rm min}=2.66\times 10^{42}$ ergs/kpc2/yr, $f_{\gamma}=0.62$, and $f_{\rm radio}=0.28$. $\tau$ (yr) |  $F_{\gamma}^{\rm min}$ (erg/kpc2/yr) | $f_{\gamma}$ |  $F_{\rm radio}^{\rm min}$ (erg/kpc2/yr) | $f_{\rm radio}$
---|---|---|---|---
$5\times 10^{3}$ | $2.37\times 10^{42}$ | 0.63 | $1.33\times 10^{41}$ | 0.27
$1\times 10^{4}$ | $4.22\times 10^{42}$ | 0.62 | $2.66\times 10^{41}$ | 0.28
$2\times 10^{4}$ | $9.44\times 10^{43}$ | 0.63 | $3.98\times 10^{42}$ | 0.26
Table 2: The best-fit values of $F_{\gamma}^{\rm min}$, $F_{\rm radio}^{\rm
min}$, $f_{\gamma}$, and $f_{\rm radio}$ for three values of $\tau$, and for
$R_{\rm PSR}$ = 0.01 $\,{\rm yr}^{-1}$.

In Fig. 4, we show the cosmic-ray positron fraction from the pulsars in the
ATNF catalog, from unknown pulsars as calculated by our Monte Carlo, and from
the sum of these contributions, adopting the best-fit values of
$F_{\gamma}^{\rm min}$, $F_{\rm radio}^{\rm min}$, $f_{\gamma}$, and $f_{\rm
radio}$ and otherwise using the same parameters as adopted in the previous
section. In the following section, we will consideration variations of these
parameters and assess their impact on our results.

Figure 4: The cosmic-ray positron fraction from the pulsars in the ATNF
catalog (black), from unknown pulsars as calculated by our Monte Carlo (blue),
and from the sum of these contributions (orange). In this figure, we have
adopted the following parameters: $\alpha=1.5$, $E_{c}=50\,{\rm TeV}$,
$\eta=0.15$, $\tau=10^{4}\,{\rm yr}$, $B_{\rm ISM}=3\,\mu{\rm G}$, and $t_{\rm
PWN}=3\times 10^{4}\,{\rm yr}$, and values of $\tau_{c}$, $r$, and $P$ as
reported in the ATNF catalog. The solid (shaded) bands around the blue and
orange curves reflect the variation observed across 68% (all) of the
realizations of our Monte Carlo.

##  5 Parameter Variations

Figure 5: The impact on the positron fraction of varying the injected spectral
index, $\alpha$, the injected spectral cutoff, $E_{c}$, the duration of the
pulsars wind nebula phase, $t_{\rm PWN}$, the spindown timescale, $\tau$, or
the magnetic field of the ISM, $B$.

The results shown in Fig. 4 were obtained using our default set of parameter
values: $\alpha=1.5$, $E_{c}=50\,{\rm TeV}$, $t_{\rm PWN}=3\times 10^{4}\,{\rm
yr}$, $\tau=10^{4}\,{\rm yr}$, $B=3\,\mu{\rm G}$, and $\eta=0.15$. Any changes
to these parameters will impact the predicted magnitude and spectral shape of
the cosmic-ray positron fraction. In this section, we explore these variations
and their connection to the observed features of the positron fraction.

In Fig. 5, we illustrate the impact of a variety of parameter variations,
modifying one of the above mentioned parameters at a time while keeping the
others fixed to their default values (with the exception of $\eta$, which is
selected in each case to normalize the positron fraction to its measured value
at $E_{e}=20\,{\rm GeV}$). Focusing first on the positron fraction as measured
at energies below a few hundred GeV, the positrons in this energy range arise
from many pulsars, and thus the positron fraction is largely insensitive to
the location, age, or other characteristics of any specific pulsars. From this
data, it is clear that an injected spectral index of approximately $\alpha\sim
1.5-1.6$ is preferred, although the precise value of this parameter also
depends on the choices of $t_{\rm PWN}$, $\tau$, and $B$ that are adopted. For
some choices of these parameters (such as $B=6\,\mu{\rm G}$, for example), the
shape of the predicted positron spectrum is more concave than is measured by
AMS-02, leading us to disfavor such scenarios.

In some cases, there are substantial degeneracies between the parameters we
have considered. In particular, we find that the value of the injected cutoff
energy, $E_{c}$, is approximately degenerate with the electron-positron
production efficiency, $\eta$. This is demonstrated in the upper right frame
of Fig. 5, where we show that the low-energy positron fraction is not
significantly impacted by variations in $E_{c}$, other than in terms of the
overall normalization. If we increase the value of $E_{c}$, we must compensate
by increasing the value of the efficiency, $\eta$.

Upon reviewing these results, one might reasonably be concerned that when we
select model parameters which provide a good fit to the positron fraction at
low energies, we consistently predict a positron flux at high energies that is
larger than that observed by AMS-02. This could be an indication that some
young pulsars lose their kinetic energy more rapidly than we have assumed (in
this study, we have limited our calculations to the case of magnetic dipole
breaking). This could also be impacted by the stochastic nature of Inverse
Compton scattering in the Klein-Nishina limit, which we have thus far treated
as continuous [20]. See our definitions in section 2 and in particular
equation 2 for details. Furthermore, in light of the fact that the positron
fraction is dominated at the highest measured energies by only a small number
of sources, this could simply be the consequence of a small number of young
and nearby pulsars which have characteristics that are substantially different
from the average of the local population. We will explore this possibility
further in the following section.

##  6 Pulsar-to-Pulsar Variations

Throughout this study, we have adopted parameter values which are intended to
reflect the average characteristics of the local pulsar population. At low and
moderate energies, we can expect this approach to be adequate, as the observed
positron fraction is the result of non-negligible contributions from several
dozen individual pulsars. At higher energies, however, the measured positron
fraction could potentially be dominated by the emission from only a few
individual pulsars. Judging from Fig. 1, for example, the single most
important pulsar (PSR J1057-5226) contributes less than 15% of the observed
positrons at 100 GeV, but is responsible for producing roughly half of these
particles at 500 GeV. The specific parameters for this individual pulsar are
thus quite important for our prediction of the positron fraction at high
energies. More generally speaking, stochastic variations in the parameters of
individual pulsars make it difficult to reliably predict the positron fraction
at energies above a few hundred GeV. It is, after all, very plausible that the
pulsar parameters being considered in this study are not universal, but rather
vary significantly from pulsar-to-pulsar.

Figure 6: The cosmic-ray positron fraction as calculated for $\alpha=1.55$
(top) or $\alpha=1.50$ (bottom), $E_{c}=10\,{\rm TeV}$, $t_{\rm PWN}=5\times
10^{4}\,{\rm yr}$, $\tau=10^{4}\,{\rm yr}$ and $B=4\,\mu{\rm G}$. In each
frame, we show the positron spectrum that results when we include all of the
pulsars (both those in the ATNF catalog, and those simulated by our Monte
Carlo), as well as that obtained when we neglect the contributions from the
first four pulsars listed in Table 1. This illustrates the importance of
pulsar-to-pulsar parameter variations, in particular at energies above a few
hundred GeV.

We explore some examples of the possible impact of pulsar-to-pulsar variations
in Fig. 6 where we show the positron fraction for the case of $\alpha=1.55$
(top) or 1.50 (bottom), $E_{c}=10\,{\rm TeV}$, $t_{\rm PWN}=5\times
10^{4}\,{\rm yr}$, $\tau=10^{4}\,{\rm yr}$ and $B=4\,\mu{\rm G}$. In each
frame, we show the positron spectrum that results when we include all of the
pulsars (both those in the ATNF catalog, and those simulated by our Monte
Carlo), as well as that obtained when we neglect the contributions from the
first four pulsars listed in Table 1. At low to moderate energies, these few
pulsars do not impact the positron fraction much, and we consistently obtain a
result which is in reasonably good agreement with the measured spectrum. At
energies above a few hundred GeV, in contrast, the result changes
significantly depending on the detailed properties of these few individual
pulsars. This makes it impossible for us to use our approach to reliably
predict the positron fraction at the highest energies measured by AMS-02.

We acknowledge that our assumed choice of sources, especially those used to
normalize to data at low energies, does lead to a unique set of best fit
parameters. In fact, it is possible that one could, in principle, find a
another set of ”best fit” free parameters which are chosen based upon other
additional sources which could reduce the aforementioned overshoot at high
energies [30]. Also, we note that the possible impacts of spatial variations
in the diffusion constant used could further play a role in how the overall
contribution is modeled, particularly at high energies. However, from earlier
work, it was found that lowering the diffusion coefficient too much, or indeed
varying it, would negate any meaningful contributions from local pulsars [17].

Finally, we did explore possible further model variations, especially in
comparisons of the known source contributions and the Monte Carlos. We find
that while our assumptions and approximations hold up at low energies, our
high energy discrepancies could be due to the fact that we do not fully model
the specific pulsar to pulsar spatial variations. In particular, our Monte
Carlo modeling does not take into account the variations in pulsar initial
period, spin-down time, or spectral index. Future studies might benefit from
further exploring and amending the pulsar modeling to account for this.

##  7 Summary and Conclusions

Observations by the HAWC Collaboration indicate that pulsars are accompanied
by TeV halos which produce a significant contribution to the cosmic-ray
positron flux. In this study, we have used the measured characteristics of the
cosmic-ray positron fraction to constrain the properties of the local
population of pulsars and their associated TeV halos. More specifically, we
have calculated the contribution to the cosmic-ray positron flux from the
pulsars contained within the Australia Telescope National Facility (ATNF)
pulsar catalog [25], and used a Monte Carlo to estimate the contribution from
those pulsars which have thus far gone undetected. After varying a number of
parameters, including the spectral shape of the injected positrons, the
timescale for pulsar spindown, the efficiency for electron-positron
production, the magnitude of the magnetic field in the ISM, and the duration
of the pulsar wind nebula phase, we have compared our results to the positron
fraction as measured by AMS-02. For reasonable parameter values, we find that
it is possible to obtain good agreement with the measured positron fraction,
at least up to energies of $E_{e}\sim 300\,{\rm GeV}$. This result further
supports the conclusion that TeV halos are the main source of the observed
cosmic-ray positron spectrum. At energies above a few hundred GeV, the
positron fraction is dominated by a small number of pulsars, limiting our
ability to reliably predict the shape of the expected positron fraction or to
use this high-energy information to draw reliable conclusions pertaining to
the Milky Way’s broader pulsar population.

Overall, the observed positron spectrum supports a picture in which pulsars
transfer $\eta\sim 5-20\%$ of their total spindown power into the production
of high-energy electron-positron pairs, generating a spectrum of such
particles with a hard spectral index, $\alpha\sim 1.5-1.7$. Our best fits were
obtained for models in which the radio beams of pulsars cover $\sim 28\%$ of
the surrounding solid angle, while the gamma-ray beams are significantly
wider, covering $\sim 62\%$ of the total solid angle.

## Acknowledgments

We would like to thank Tim Linden for helpful discussions. DH is supported by
the Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the
U.S. Department of Energy, Office of Science, Office of High Energy Physics.

## References

  * [1] A. U. Abeysekara  “Extended gamma-ray sources around pulsars constrain the origin of the positron flux at Earth”  In _Science_ 358.6365, 2017, pp. 911–914  DOI: 10.1126/science.aan4880
  * [2] A. U. Abeysekara  “Multiple Galactic Sources with Emission Above 56 TeV Detected by HAWC”  In _Phys. Rev. Lett._ 124.2, 2020, pp. 021102  DOI: 10.1103/PhysRevLett.124.021102
  * [3] O. Adriani et al.  “Cosmic-Ray Positron Energy Spectrum Measured by PAMELA”  In _Phys. Rev. Lett._ 111 American Physical Society, 2013, pp. 081102  DOI: 10.1103/PhysRevLett.111.081102
  * [4] M. Aguilar et al.  “Cosmic-ray positron fraction measurement from 1 to 30 GeV with AMS-01”  In _Physics Letters B_ 646.4 Elsevier BV, 2007, pp. 145–154  DOI: 10.1016/j.physletb.2007.01.024
  * [5] M. Aguilar et al.  “Towards Understanding the Origin of Cosmic-Ray Electrons”  In _Phys. Rev. Lett._ 122 American Physical Society, 2019, pp. 101101  DOI: 10.1103/PhysRevLett.122.101101
  * [6] M. Aguilar et al.  “Towards Understanding the Origin of Cosmic-Ray Positrons”  In _Phys. Rev. Lett._ 122 American Physical Society, 2019, pp. 041102  DOI: 10.1103/PhysRevLett.122.041102
  * [7] A. Albert  “3HWC: The Third HAWC Catalog of Very-High-Energy Gamma-ray Sources”  In _Astrophys. J._ 905.1, 2020, pp. 76  DOI: 10.3847/1538-4357/abc2d8
  * [8] A. Albert  “Evidence that Ultra-high-energy Gamma Rays Are a Universal Feature near Powerful Pulsars”  In _Astrophys. J. Lett._ 911.2, 2021, pp. L27  DOI: 10.3847/2041-8213/abf4dc
  * [9] S. W. Barwick et al.  “Measurements of the Cosmic-Ray Positron Fraction from 1 to 50 G[CLC]e[/CLC]V”  In _The Astrophysical Journal_ 482.2 American Astronomical Society, 1997, pp. L191–L194  DOI: 10.1086/310706
  * [10] Pasquale Blasi  “Origin of the Positron Excess in Cosmic Rays”  In _Physical Review Letters_ 103.5 American Physical Society (APS), 2009  DOI: 10.1103/physrevlett.103.051104
  * [11] Ilias Cholis, Tanvi Karwal and Marc Kamionkowski  “Features in the Spectrum of Cosmic-Ray Positrons from Pulsars”  In _Phys. Rev. D_ 97.12, 2018, pp. 123011  DOI: 10.1103/PhysRevD.97.123011
  * [12] Ilias Cholis, Tanvi Karwal and Marc Kamionkowski  “Studying the Milky Way pulsar population with cosmic-ray leptons”  In _Phys. Rev. D_ 98.6, 2018, pp. 063008  DOI: 10.1103/PhysRevD.98.063008
  * [13] V. Di Felice et al.  “Looking for cosmic ray data? The ASI Cosmic Ray Database”  In _35th International Cosmic Ray Conference (ICRC2017)_ 301, International Cosmic Ray Conference, 2017, pp. 1073
  * [14] Kun Fang, Xiao-Jun Bi and Peng-fei Yin  “Discriminating local sources of high-energy cosmic electrons and positrons by current and future anisotropy measurements”  In _Mon. Not. Roy. Astron. Soc._ 478.4, 2018, pp. 5660–5670  DOI: 10.1093/mnras/sty1463
  * [15] Dan Hooper and Tim Linden  “Evidence of TeV Halos Around Millisecond Pulsars”, 2021  arXiv:2104.00014 [astro-ph.HE]
  * [16] Dan Hooper and Tim Linden  “Measuring the local diffusion coefficient with H.E.S.S. observations of very high-energy electrons”  In _Physical Review D_ 98.8 American Physical Society (APS), 2018  DOI: 10.1103/physrevd.98.083009
  * [17] Dan Hooper and Tim Linden  “Measuring the Local Diffusion Coefficient with H.E.S.S. Observations of Very High-Energy Electrons”  In _Phys. Rev. D_ 98.8, 2018, pp. 083009  DOI: 10.1103/PhysRevD.98.083009
  * [18] Dan Hooper and Tim Linden  “Millisecond Pulsars, TeV Halos, and Implications For The Galactic Center Gamma-Ray Excess”  In _Phys. Rev. D_ 98.4, 2018, pp. 043005  DOI: 10.1103/PhysRevD.98.043005
  * [19] Dan Hooper, Ilias Cholis, Tim Linden and Ke Fang  “HAWC observations strongly favor pulsar interpretations of the cosmic-ray positron excess”  In _Physical Review D_ 96.10 American Physical Society (APS), 2017  DOI: 10.1103/physrevd.96.103013
  * [20] Isabelle John and Tim Linden  “Pulsars Do Not Produce Sharp Features in the Cosmic-Ray Electron and Positron Spectra”  arXiv, 2022  DOI: 10.48550/ARXIV.2206.04699
  * [21] Simon Johnston, D A Smith, A Karastergiou and M Kramer  “The Galactic population and properties of young, highly energetic pulsars”  In _Monthly Notices of the Royal Astronomical Society_ 497.2 Oxford University Press (OUP), 2020, pp. 1957–1965  DOI: 10.1093/mnras/staa2110
  * [22] Jagdish C. Joshi and Soebur Razzaque  “A self-consistent model of cosmic-ray fluxes and positron excess: Roles of nearby pulsars and a sub-dominant source population”  In _JCAP_ 09, 2017, pp. 029  DOI: 10.1088/1475-7516/2017/09/029
  * [23] Tim Linden et al.  “Using HAWC to discover invisible pulsars”  In _Phys. Rev. D_ 96.10, 2017, pp. 103016  DOI: 10.1103/PhysRevD.96.103016
  * [24] Dunc R. Lorimer  “The galactic population and birth rate of radio pulsars”  In _IAU Symp._ 218, 2004, pp. 105  arXiv:astro-ph/0308501
  * [25] R. N. Manchester, G. B. Hobbs, A. Teoh and M. Hobbs  “The Australia Telescope National Facility Pulsar Catalogue”  In _The Astronomical Journal_ 129.4 American Astronomical Society, 2005, pp. 1993–2006  DOI: 10.1086/428488
  * [26] I. V. Moskalenko and A. W. Strong  “Production and propagation of cosmic ray positrons and electrons”  In _Astrophys. J._ 493, 1998, pp. 694–707  DOI: 10.1086/305152
  * [27] Luca Orusa, Silvia Manconi, Fiorenza Donato and Mattia Di Mauro  “Constraining positron emission from pulsar populations with AMS-02 data”  In _JCAP_ 12.12, 2021, pp. 014  DOI: 10.1088/1475-7516/2021/12/014
  * [28] Stefano Profumo  “Dissecting cosmic-ray electron-positron data with Occam’s Razor: the role of known Pulsars”  In _Central Eur. J. Phys._ 10, 2011, pp. 1–31  DOI: 10.2478/s11534-011-0099-z
  * [29] Stefano Profumo, Javier Reynoso-Cordova, Nicholas Kaaz and Maya Silverman  “Lessons from HAWC pulsar wind nebulae observations: The diffusion constant is not a constant; pulsars remain the likeliest sources of the anomalous positron fraction; cosmic rays are trapped for long periods of time in pockets of inefficient diffusion”  In _Phys. Rev. D_ 97.12, 2018, pp. 123008  DOI: 10.1103/PhysRevD.97.123008
  * [30] Nir J. Shaviv, Ehud Nakar and Tsvi Piran  “Inhomogeneity in Cosmic Ray Sources as the Origin of the Electron Spectrum and the PAMELA Anomaly”  In _Physical Review Letters_ 103.11 American Physical Society (APS), 2009  DOI: 10.1103/physrevlett.103.111302
  * [31] Takahiro Sudoh, Tim Linden and Dan Hooper  “The Highest Energy HAWC Sources are Likely Leptonic and Powered by Pulsars”, 2021  DOI: 10.1088/1475-7516/2021/08/010
  * [32] V. Vagelli  “Results from AMS-02 on the ISS after 6 years in space”  In _Nuovo Cim. C_ 42.4, 2019, pp. 173  DOI: 10.1393/ncc/i2019-19173-y
  * [33] Andrey E. Vladimirov et al.  “GALPROP WebRun: an internet-based service for calculating galactic cosmic ray propagation and associated photon emissions”  In _Comput. Phys. Commun._ 182, 2011, pp. 1156–1161  DOI: 10.1016/j.cpc.2011.01.017


* [46] K. Ecker and G. Huisken, _Mean curvature evolution of entire graphs._ Ann. of Math. (2) 130 (1989), no.3, 453-471.
  * [47] J.H. Eschenburg and E. Heintze, _Comparison theory for Riccati equations._ Manuscripta Math. 68 (1990), 209-214.
  * [48] A. Farina, _Liouville-type theorems for elliptic problems._ In: M. Chipot (Ed.), Handbook of Differential Equations - vol. 4, Stationary Partial Differential Equations, Elsevier, 2007, pp. 60-116.
  * [49] A. Farina and J. Serrin, _Entire solutions of completely coercive quasilinear elliptic equations._ J. Differential Equations 250 (2011), no.12, 4367-4408.
  * [50] A. Farina and J. Serrin, _Entire solutions of completely coercive quasilinear elliptic equations, II._ J. Differential Equations 250 (2011), no.12, 4409-4436.
  * [51] P. Felmer, M. Montenegro, A. Quaas, _A note on the strong maximum principle and the compact support principle._ J. Differential Equations 246 (2009), 39-49.
  * [52] P. Felmer, A. Quaas and B. Sirakov, _Solvability of nonlinear elliptic equations with gradient terms._ J. Differential Equations 254 (2013), 4327-4346.
  * [53] R. Filippucci, _Nonexistence of positive weak solutions of elliptic inequalities._ Nonlinear Anal. 8 (2009), 2903-2916.
  * [54] R. Filippucci, P. Pucci and M. Rigoli, _Non–existence of entire solutions of degenerate elliptic inequalities with weights._ Archive Rat. Mech. Anal. 188 (2008), 155-179; Erratum, 188 (2008), 181.
  * [55] R. Filippucci, P. Pucci and M. Rigoli, _On weak solutions of nonlinear weighted $p$–Laplacian elliptic inequalities._ Nonlinear Anal. 70 (2009), 3008-3019.
  * [56] R. Filippucci, P. Pucci and M. Rigoli, _On entire solutions of degenerate elliptic differential inequalities with nonlinear gradient terms._ J. Math. Anal. Appl. 356 (2009), no.2, 689-697.
  * [57] R. Filippucci, P. Pucci and M. Rigoli, _Nonlinear weighted $p$–Laplacian elliptic inequalities with gradient terms._ Comm. Cont. Math. 12 (2010), no.3, 501-535.
  * [58] D. Gilbarg and N. Trudinger, _Elliptic Partial Differential Equations of Second Order._ Third ed., Springer-Verlag, 1998.
  * [59] E. Giusti, _Minimal surfaces and functions of bounded variation._ Monographs in Mathematics, 80. Birkhäuser Verlag, Basel, 1984. xii+240 pp.
  * [60] M. Ghergu, C. Niculescu and V. Radulescu, _Explosive solutions of elliptic equations with absorption and nonlinear gradient term._ Proc. Indian Acad. Sci. Math. Sci. 112 (2002), 441-451.
  * [61] A. Greco, _On the existence of large solutions for equations of prescribed mean curvature._ Nonlinear Anal. 34 (1998), no.4, 571-583.
  * [62] R. Greene and H.H. Wu, Function theory on manifolds which possess a pole, Lecture Notes in Math. 699, Springer-Verlag, 1979\.
  * [63] A. Grigor’yan, _Analytic and geometric background of recurrence and non-explosion of the Brownian motion on Riemannian manifolds._ Bull. Amer. Math. Soc. 36 (1999), 135-249.
  * [64] Y. Haitao, _A compact support principle for a class of elliptic differential inequalities._ Nonlinear Anal. 58 (2004), 103-119.
  * [65] E. Hebey, _Optimal Sobolev inequalities on complete Riemannian manifolds with Ricci curvature bounded below and positive injectivity radius._ Amer. J. Math. 118 (1996), no.2, 291-300.
  * [66] E. Hebey, _Nonlinear analysis on manifolds: Sobolev spaces and inequalities._ Courant Lecture Notes in Mathematics, 5. New York University, Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, 1999. x+309 pp.
  * [67] J. Heinonen and T. Kilpelainen, _$\mathscr{A}$ -superharmonic functions and supersolutions of degenerate elliptic equations._ Ark. Mat. 26 (1988), 87-105.
  * [68] J. Heinonen, T. Kilpeläinen and O. Martio, _Nonlinear potential theory of degenerate elliptic equations._ Unabridged republication of the 1993 original, Dover Publications Inc., Mineola, NY (2006), xii+404.
  * [69] D. Hoffman and J. Spruck, _Sobolev and isoperimetric inequalities for Riemannian submanifolds._ Comm. Pure Appl. Math. 27 (1974), 715-727.
  * [70] I. Holopainen, _Nonlinear potential theory and quasiregular mappings on Riemannian manifolds._ Ann. Acad. Sci. Fenn. Ser. A I Math. Dissertationes 74 (1990), 45pp.
  * [71] I. Holopainen, _Positive solutions of quasilinear elliptic equations on Riemannian manifolds._ Proc. London Math. Soc. 65 (1992), 651-672.
  * [72] I. Holopainen, _Volume growth, Green’s functions, and parabolicity of ends._ Duke Math. J. 97 (1999), no.2, 319-346.
  * [73] E. Hopf, _Elementäre Bemerkungen über die Lösungen partieller Differentialgleichungen zweiter Ordnung vom elliptischen Typus._ Sitzungsberichte Preussische Akademie der Wissenschaften (Berlin), 1927, 147-152.
  * [74] E.-P. Hsu, _Heat semigroup on a complete Riemannian manifold._ Ann. Probab. 17 (1989), no.3, 1248-1254.
  * [75] E.-P. Hsu, _Stochastic Analysis on Manifolds._ Grad. Stud. Math., vol. 38, Amer. Math. Soc., 2002.
  * [76] M. Kanai, _Rough isometries, and combinatorial approximations of geometries of non-compact Riemannian manifolds._ J. Math. Soc. Japan 37 (1985), 391-413.
  * [77] L. Karp, _Differential inequalities on complete Riemannian manifolds and applications._ Math. Ann. 272 (1985), no.4, 449-459.
  * [78] R.Z. Khas’minskii, _Ergodic properties of recurrent diffusion processes and stabilization of the solution of the Cauchy problem for parabolic equations._ Teor. Verojatnost. i Primenen., Akademija Nauk SSSR. Teorija Verojatnosteĭ i ee Primenenija 5 (1960), 196-214.
  * [79] T. Kato _Schrödinger operators with singular potentials._ Israel J. Math. 13 (1972), 135-148
  * [80] J.B. Keller, _On solutions of $\Delta u=f(u)$._ Comm. Pure Appl. Math. 10 (1957), 503-510.
  * [81] S. Kichenassamy and L. Véron, _Singular solutions of the $p$-Laplace equation_, Math. Ann 275 (1986), no.4, 599-615.
  * [82] T. Kura, _The weak supersolution-subsolution method for second order quasilinear elliptic equations._ Hiroshima Math. J. 19 (1989), no.1, 1-36.
  * [83] A. Lair and A. Wood, _Large solutions of semilinear elliptic equations with nonlinear gradient terms._ Int. J. Math. Sci. 22 (1999), no.4, 869?883.
  * [84] J.M. Lasry and P.L. Lions, _Nonlinear elliptic equations with singular boundary conditions and stochastic control with state constraints. I. The model problem._ Math. Ann. 283 (1989), no.4, 583-630.
  * [85] V.K. Le, _On some equivalent properties of sub- and supersolutions in second order quasilinear elliptic equations._ Hiroshima Math. J. 28 (1998), 373-380.
  * [86] X. Li and F. Li, _Nonexistence of solutions for singular quasilinear differential inequalities with a gradient nonlinearity._ Nonlinear Anal. 75 (2012), no.5, 2812-2822.
  * [87] P. Li and L. Karp, _The heat equation on complete Riemannian manifolds._ unpublished.
  * [88] P. Li and S.T. Yau, _On the parabolic kernel of the Schrödinger operator._ Acta Math. 156 (1986), no.3-4, 153-201.
  * [89] M. Magliaro, L. Mari, P. Mastrolia and M. Rigoli, _Keller–Osserman type conditions for differential inequalities with gradient terms on the Heisenberg group._ J. Differential Equations 250 (2011), 2643-2670.
  * [90] C. Mantegazza, G. Mascellani and G. Uraltsev, _On the distributional Hessian of the distance function._ Pacific J. Math. 270 (2014), no.1, 151-166.
  * [91] C. Mantegazza, _Lecture notes on mean curvature flow._ Progress in Mathematics, 290. Birkhäuser/Springer Basel AG, Basel, 2011. xii+166 pp.
  * [92] C. Mantegazza and A.C. Mennucci, _Hamilton-Jacobi equations and distance functions on Riemannian manifolds._ Applied Math. Opt. 47 (2002), no.1, 1-25.
  * [93] L. Mari and L.F. Pessoa, _Duality between Ahlfors-Liouville and Khas’minskii properties for nonlinear equations._ To appear on Comm. Anal. Geom., Available at arXiv:1603.09113.
  * [94] L. Mari, M. Rigoli and A.G. Setti, _Fake distances and gradient estimates for the Green kernel of the $p$-Laplacian under Ricci lower bounds._ In preparation.
  * [95] L. Mari, M. Rigoli and A.G. Setti, _Keller-Osserman conditions for diffusion-type operators on Riemannian Manifolds._ J. Funct. Anal. 258 (2010), no.2, 665-712.
  * [96] L. Mari and D. Valtorta, _On the equivalence of stochastic completeness, Liouville and Khas’minskii condition in linear and nonlinear setting._ Trans. Amer. Math. Soc. 365 (2013), no.9, 4699-4727.
  * [97] O. Martio and G. Porru, _Large solutions of quasilinear elliptic equations in the degenerate case._ Complex analysis and differential equations (Uppsala, 1997) Acta Univ. Upsaliensis Skr. Uppsala Univ. C Organ. Hist. 64 (1999), 225-241.
  * [98] E. Mitidieri and S.I. Pokhozhaev, _A priori estimates and the absence of solutions of nonlinear partial differential equations and inequalities._ Tr. Mat. Inst. Steklova, 234 (2001), 1-384.
  * [99] K. Motomiya, _On functions which satisfy some differential inequalities on Riemannian manifolds._ Nagoya Math. J. 81 (1981), 57-72.
  * [100] Y. Naito and H. Usami, _Entire solutions of the inequality $\mathrm{div}(A(|\nabla u|)\nabla u)\geq f(u)$._ Math. Z. 225 (1997), 167-175.
  * [101] L. Ni, _Mean value theorems on manifolds._ Asian J. Math. 11 (2007), no.2, 277-304.
  * [102] M. Obata, _Conformal transformations of compact Riemannian manifolds._ Illinois J. Math. 6 (1962) 292-195.
  * [103] H. Omori, _Isometric immersions of Riemannian manifolds._ J. Math. Soc. Japan 19 (1967), 205-214.
  * [104] R. Osserman, _On the inequality $\Delta u\geq f(u)$._ Pacific J. Math. 7 (1957), 1641-1647.
  * [105] R. Osserman, _A survey of minimal surfaces._ $2^{\mathrm{nd}}$ ed., Dover Publications, Inc., New York, 1986, vi+207.
  * [106] P. Petersen, _Riemannian Geometry._ Springer-Verlag, 1997.
  * [107] S. Pigola, M. Rigoli and A.G. Setti, _Maximum principles and singular elliptic inequalities._ J. Funct. Anal. 193 (2002), no.2, 224-260.
  * [108] S. Pigola, M. Rigoli and A.G. Setti, _A remark on the maximum principle and stochastic completeness._ Proc. Amer. Math. Soc. 131 (2003), no.4, 1283-1288.
  * [109] S. Pigola, M. Rigoli and A.G. Setti, _Volume growth, “a priori” estimates, and geometric applications._ Geom. Funct. Anal. 13 (2003), no.6, 1303-1328.
  * [110] S. Pigola, M. Rigoli and A.G. Setti, _Maximum principles on Riemannian manifolds and applications._ Mem. Amer. Math. Soc. 174 (2005), no.822.
  * [111] S. Pigola, M. Rigoli and A.G. Setti, _Maximum principles at infinity on Riemannian manifolds: an overview._ Matemática Contemporânea 31 (2006), 81-128.
  * [112] S. Pigola, M. Rigoli and A.G. Setti, _Some non-linear function theoretic properties of Riemannian manifolds._ Rev. Mat. Iberoam. 22 (2006), no.3, 801-831.
  * [113] S. Pigola, M. Rigoli and A.G. Setti, _Vanishing and finiteness results in Geometric Analysis. A generalization of the Böchner technique._ Progress in Mathematics 266, Birkäuser, 2008, xiv+282 pp.
  * [114] S. Pigola and A.G. Setti, _The Feller property on Riemannian manifolds._ J. Funct. Anal. 262 (2012), no.5, 2481-2515.
  * [115] S. Pigola, A.G. Setti and M. Troyanov, _The topology at infinity of a manifold supporting an ${L^{q,p}}$-Sobolev inequality._ Expo. Math. 32 (2014), no.4, 365-383.
  * [116] M.H. Protter and H.F. Weinberger, _Maximum principles in differential equations._ Corrected reprint of the 1967 original, Springer-Verlag, 1984.
  * [117] P. Pucci, M. García-Huidobro, R. Manásevich, J. Serrin, _Qualitative properties of ground states for singular elliptic equations with weights._ Ann. Mat. Pura Appl. 185 (2006), no.4, suppl., 205-243.
  * [118] P. Pucci and V. Radulescu, _The maximum principle with lack of monotonicity_ , in preparation.
  * [119] P. Pucci and M. Rigoli, _Entire solutions of singular elliptic inequalities on complete manifolds._ Disc. Cont. Dyn. Sys. 20 (2008), no.1, 115-137.
  * [120] P. Pucci, M. Rigoli and J. Serrin, _Qualitative properties for solutions of singular elliptic inequalities on complete manifolds._ J. Differential Equations 234 (2007), no.2, 507-543.
  * [121] P. Pucci and J. Serrin, _Dead cores and bursts for quasilinear singular elliptic equations._ SIAM J. Math. Anal. 38 (2006), 259-278.
  * [122] P. Pucci and J. Serrin, _The Strong Maximum Principle Revisited._ J. Differ. Equations, 196 (2004), 1-66; Erratum, J. Differ. Equations, 207 (2004), 226-227.
  * [123] P. Pucci and J. Serrin, _The maximum principle._ Progress in Nonlinear Differential Equations and their Applications, 73, Birkhäuser Verlag, Basel, 2007, x+235 pp.
  * [124] P. Pucci and J. Serrin, _A remark on entire solutions of quasilinear elliptic equations._ J. Differential Equations 250 (2011), no.2, 675-689.
  * [125] P. Pucci, J. Serrin and H. Zou, _A strong maximum principle and a compact support principle for singular elliptic inequalities._ J. Math. Pures Appl. 78 (1999), 769-789.
  * [126] V. Radulescu, _Singular phenomena in nonlinear elliptic problems: from boundary blow-up solutions to equations with singular nonlinearities._ In: Michel Chipot (Ed.), Handbook of Differential Equations: Stationary Partial Differential Equations (2007), vol. 4, 483-591.
  * [127] A. Ratto and M. Rigoli, _Elliptic differential inequalities with applications to harmonic maps._ J. Math. Soc. Japan 45 (1993), 321-337.
  * [128] A. Ratto and M. Rigoli and A.G. Setti, _On the Omori-Yau maximum principle and its application to differential equations and geometry._ J. Funct. Anal. 134 (1995), 486-510.
  * [129] R. Redheffer, _On the inequality $\Delta u\geq f(u,\,|\,{\rm grad}\,u|)$._ J. Math. Anal. Appl. 1 (1960), 277-299.
  * [130] R. Redheffer, _Maximum principles and duality._ Monatsh. Math. 62 (1958), 56-75.
  * [131] R. Redheffer, _Nonlinear differential inequalities and functions of compact support._ Trans. Amer. Math. Soc. 220 (1976) 133-157.
  * [132] M. Rigoli, M. Salvatori and M. Vignati, _On the compact support principle on complete manifolds._ J. Diff.Eq. 246 (2009), 870-894.
  * [133] M. Rigoli, M. Salvatori and M. Vignati, _Volume growth and $p$-subharmonic functions on complete manifolds._ Math. Z. 227 (1998), no.3, 367-375.
  * [134] M. Rigoli, M. Salvatori and M. Vignati, _Some remarks on the weak maximum principle._ Rev. Mat. Iberoamericana 21 (2005), no.2, 459-481.
  * [135] M. Rigoli, M. Salvatori and M. Vignati, _A Liouville type theorem for a general class of operators on complete manifolds._ Pacific J. Math. 194 (2000), no.2, 439-453.
  * [136] M. Rigoli, M. Salvatori and M. Vignati, _The compact support principle for differential inequalities with gradient terms._ Nonlinear Anal. 72 (2010), 4360-4376.
  * [137] H. Rosenberg, F. Schulze and J. Spruck, _The half-space property and entire positive minimal graphs in $M\times\mathbb{R}$._ J. Diff. Geom. 95 (2013), 321-336.
  * [138] R. Schoen and S.-T. Yau. _Lectures on differential geometry, volume 1._ International press Cambridge, 1994.
  * [139] J. Serrin, _Local behaviour of solutions of quasilinear equations._ Acta Math. 111 (1964), 247-302.
  * [140] J. Serrin, _Isolated singularities of solutions of quasilinear equations._ Acta Math. 113 (1965), 219-240.
  * [141] J. Serrin, _Entire solutions of quasilinear elliptic equations._ J. Math. Anal. Appl. 352 (2009), no.1, 3-14.
  * [142] J. Serrin, _The Liouville theorem for homogeneous elliptic differential inequalities._ Problems in mathematical analysis. no.61. J. Math. Sci. (N.Y.) 179 (2011), no.1, 174-183.
  * [143] J. Serrin, _Entire solutions of nonlinear Poisson equations._ Proc. London. Math. Soc. 24 (1972), no.3, 348-366.
  * [144] V.G. Tkachev, _Some estimates for the mean curvature of nonparametric surfaces defined over domains in $\mathbb{R}^{n}$._ Ukr. Geom. Sb. 35 (1992), 135-150; translated in: J. Math. Sci. (N. Y.) 72 (1994), 3250-3260.
  * [145] P. Tolksdorf, _Regularity of a more general class of quasilinear elliptic equations._ J. Differential Equations 51 (1984), 126-150.
  * [146] M. Troyanov, _Parabolicity of manifolds._ Siberian Adv. Math. 9 (1999), no.4, 125-150.
  * [147] M. Troyanov, _Solving the $p$-Laplacian on manifolds._ Proc. Amer. Math. Soc. 128 (2000), no.2, 541-545.
  * [148] H. Usami, _Nonexistence of positive entire solutions for elliptic inequalities of the mean curvature type_. J. Differential Equations 111 (1994) 472-480.
  * [149] N.T. Varopoulos, _The Poisson kernel on positively curved manifolds._ J. Funct. Anal. 44 (1981), no.3, 359-380.
  * [150] L. Véron, _Singularities of solutions of second order quasilinear equations_ , Pitman Research Notes in Mathematics Series 353 (1996), viii+377.
  * [151] X.-J. Wang, _Convex solutions to the mean curvature flow._ Ann. of Math. (2) 173 (2011), no.3, 1185-1239.
  * [152] L. Wang, _A Bernstein type theorem for self-similar shrinkers._ Geom. Dedicata 151 (2011), 297-303.
  * [153] X. Wang and L. Zhang, _Local gradient estimate for $p$-harmonic functions on Riemannian manifolds._ Comm. Anal. Geom. 19 (2011), no.4, 759-771.
  * [154] S.T. Yau, _A general Schwarz lemma for Kähler manifolds._ Amer. J. Math. 100 (1978), no. 1, 197-203.
  * [155] S.T. Yau, _Remarks on conformal transformations._ J. Diff. Geom. 8 (1973), 369-381.
  * [156] S.T. Yau, _Some function-theoretic properties of complete Riemannian manifold and their applications to geometry._ Indiana Univ. Math. J. 25 (1976), no.7, 659-670.


# Automatic Code Generation using Pre-Trained Language Models

Luis Perez
Department of Computer Science
Stanford University
<EMAIL_ADDRESS>Lizi Ottens
Department of Computer Science
Stanford University
<EMAIL_ADDRESS>
Sudharshan Viswanathan
Department of Computer Science
Stanford University
<EMAIL_ADDRESS>
###### Abstract

Recent advancements in natural language processing [1] [2] have led to near-
human performance in multiple natural language tasks. In this paper, we seek
to understand whether similar techniques can be applied to a highly structured
environment with strict syntax rules. Specifically, we propose an end-to-end
machine learning model for code generation in the Python language built on-top
of pre-trained language models. We demonstrate that a fine-tuned model can
perform well in code generation tasks, achieving a BLEU score of 0.22, an
improvement of 46% over a reasonable sequence-to-sequence baseline. All
results and related code used for training and data processing are available
on GitHub. 111See shared repository located at
https://github.com/kandluis/code-gen.

##  1 Introduction

Automating even small parts of software development is an active research area
[3], with multiple approaches proposed methods (See Section 1). Succeeding in
the automation of even small tasks can save time for countless software
engineers, which translates to saved resources across multiple industries.
Furthermore, as software continues to eat the world 222See article by Marc
Andreessen. and demand for experienced software developers continues to
outpace supply, automatic code generation will become increasingly important.

In this paper, we propose a machine learning model to automate the task of
writing code by assisting developers in writing individual units of
functionality (or “functions”). Automating code generation can take on many
forms, from auto-completing lines of source code to generating lines of source
code from comments, generating source code from UI images, or generating unit
tests from source code. In this project, we aim to take the initial lines of
code (a function signature) along with a doc-string (function documentation)
and generate the corresponding function body. In order to do this, we use a
pre-trained language model and fine-tune it on a canonical corpus of Python
code scraped from GitHub [4].

##  2 Background

A primary challenge in code generation is that it is still an active area of
research, with many possible solutions and ongoing investigation [5]. State of
the art solutions have not yet come close to automating basic tasks software
engineers perform on a daily basis.

###  2.1 Traditional Code Completion

The most traditional and well-known approach used by multiple IDEs across a
range of languages simply consists of token completion based on structured
information obtained from static analysis of code. For example, when a
developer types a sequence of characters, the system will attempt to find
near-matching strings corresponding to function definitions and propose
completing these function calls. Similarly, for object methods, on the typing
of the accessor token (such as “-¿” or “.”), the IDE will propose auto-
completing different methods belonging to the object.

The biggest drawback of these approaches is that they lack true understanding
of the programmers intent, and also lack context relating to the surrounding
code other than that from heuristics by the tool’s developers.

###  2.2 Using Machine Learning for Code Search

Another approach taken in multiple papers in the literature [4] involves
framing the problem as a code search problem. Rather than trying to generate
code or complete the code that the developer is making, we can re-frame the
problem as one of searching for relevant pre-existing snippets. This is the
primary approach we take in three of our baseline models.

###  2.3 Using Machine Learning for Code Generation

Other more novel approaches from literature [5] are typically applied to
restricted language domains, and have massive complexity in evaluation
results, etc. Specifically, while pre-trained models are trained on free-form
language data, programming languages often utilize non-natural variable names,
function names, and syntax with more structure [5]. Work in this area has
focused on creating more structured models that take advantage of specific
architectures [6]. In [7], the authors work to first decompose the input
sequence of text tokens for the context into a tree-like structure. Other
approaches involve restricting the output of the model to a context-free
grammar (CFG) or domain-specific language (DSL) [8]. A code generation model’s
output must adhere to a very specific form in order to be syntactically
correct.

In this paper, we instead focus on taking a different approach. As has been
demonstrated by ever-increasing sizes of language models, we focus on
improving the performance on the code prediction task by making use of pre-
trained language models that are then fine-tuned on code.

###  2.4 Dataset and Feature

In this project, we are leveraging the CodeSearchNet dataset [4]. The dataset
consists of 2 million (comment, code) pairs from open source libraries,
ranging in languages from Python to Javascript, PHP, Java, Go and Ruby. Median
code-length consists of 60-100 text tokens, with 95% code-length of up to 350
tokens. Median documentation length consists of 10 text tokens. The
distributions of methods and (comment, code) pairs across programming language
are visualized in Figure 3.

We restrict our dataset to samples in the Python programming language rather
than the others available. Focusing on Python, there are over 1M methods and
approximately 500k (comment, code) pairs that make up our dataset. We make
this decision both for practical and modeling reasons. From a practical
perspective, restricting to a reasonably-sized dataset focused on a single-
language domains permits for more thorough ablation studies. From a modeling
perspective, we belief that transfer learning from natural language to a
programming language such as Python is an easier task to accomplish.

##  3 Methodology

In this section, we explain our methodology for multiple experiments and
baselines proposed as well as details on the training data and distribution.

###  3.1 CodeSearchNet Models

Figure 4 explains the general architecture of the baseline models from the
CodeSearchNet task. We successfully trained and evaluated two baselines:
Neural-Bag-Of-Words and an RNN-based baseline. See Section 4.

Generally speaking, the baselines models take as input examples of (comments,
code) pairs and learn to retrieve a specific code snippet. Each programming
language has its own encoder network (see three columns to the right in Figure
4), which are tasked with encoding a set of candidate code snippets. They are
then combined through a dot product operation with the embedding generated by
the query (docstring) encoder to produce a matrix comparison.

The matrix diagonal serves as the scores of each query doc string/code
snippet. Through this methodology, these baseline models are able to extract
meaningful information and learn a joint distribution over the query and
comment pairs. We train these models as a baseline since we believe they will
be useful in the downstream task of code generation. The models are trained on
the following loss function:

$\displaystyle-\frac{1}{N}\sum_{i}\log\left(\frac{\exp(E_{c}(\textbf{c}_{i}^{T})E_{q}(\textbf{d}_{i}))}{\sum_{j}\exp(E_{c}(\textbf{c}_{j}^{T})E_{q}(\textbf{d}_{j}))}\right)$
(1)

###  3.2 From Scratch RNN Models

The above baseline is useful only in the sense that it would allow our system
to find pre-existing code snippets which might be relevant to the developer.
Since our goal is rather to make novel code, we propose a different baseline
based on a more traditional sequence-to-sequence model.

In this case, we use a traditional RNN architecture which takes as input
individual characters. The reason we take this approach is to circumvent the
need to learn word-level embeddings. Furthermore, we hypothesize that making
use of entire words, from NLP models, will actually harm the performance of
the model for code generation. The primary reason for this being that most of
the syntax involved in writing code does not generally map directly to the
English language. Concretely, we encode each character present in the training
data as a 1-of-k encoding (one-hot encoding) and feed them into an RNN one at
a time. Our output will be a k-dimensional output vector corresponding to a
probability distribution over the entire set of characters.

For the model architecture, we sweep over multiple types of RNN cells,
including LSTM, RNN, and GRU. We find the best performing model to consists of
an LSTM-based model using a hidden state size of $128$ with two hidden layers
in the internal RNN cell. Our training takes place using sequences of $50$
characters, sampled at random from our input code. Given a sequence from $i$
to $i+50$, the model is trained to predict the sequence from $i+1$ to $i+51$.
This means we have a many-to-many sequence model (See Figure 6.2.1). We use
batch size of $50$ and train for a total of $50$ epochs.

To avoid issues with gradient explosion and stabilize training, we make
liberal use of gradient-clipping. In particular, we clip all gradients to an
absolute size of $5$.

We sweep over learning rates and find that a started learning rate of $0.002$
with an exponentially decaying schedule appears to perform best as measured by
a held-out validation set. We use a decay rate of $0.97$ per epoch. We also
experiment with the use of dropout, but find little impact on final
performance.

###  3.3 Fine-Tuned Pre-Trained Large Language Models

Our final approach relies on the use of pre-trained language models. We fine
tune our code generator using the small GPT-2 model with 117 million
parameters. Using such a large backbone and continuing to fine tune allows us
to generate synthetic code samples with even higher quality, treating
programming languages as another specific domain alongside encyclopedia
articles, news or books.

Figure 1: Decoder-Only Architecture used by GPT-2.

The general architecture of the GPT-2 model consists of a sequence-to-sequence
predictive task based on the transformer architecture [9] [1]. However, it
consists solely of the 12-layer decoder-only, as visualized in Figure 1. Each
layer has 12 independent attention heads, leading to 144 distinct attention
patterns. By making use of an attention-based framework, the model is more
adept at dealing with long-range dependencies. This is because the attention
mechanism allows the model to focus on the encoding of any of the input
sequence tokens.

##  4 Results

CodeSearchNet provides a good starting point as we are able to train different
models on the input code streams. We trained a simple LSTM model as well as a
neural bag of words model on a combination of all the available (code,
documentation) pairs. For details on these simple baselines, please see
Appendix Section 6.1.

###  4.1 Code Generation with Char-RNN

As both of the above baselines focus on understanding and extracting useful
embeddings for our overall task, our primary baseline consists of a straight-
forward sequence-to-sequence model. Given that code typically does not consist
of English words and can instead have quite a varied syntax, our baseline
model is a model which uses character level embedding, so it is character
aware [10].

Due to computational constraints, we train only on the Python subset of the
data and only on 10% of the total data available. For the char-rnn model [10],
this corresponds to around 50MB of raw text, or 78,357,395 characters with
1,618 distinct symbols. Figure 9 shows the training and validation losses on
the model. The loss is simply a softmax loss on the 1,618 characters for a
sequence of length 128 (the model is trained on sequences of length 128 by
default). Figure 10 shows the perplexity, or the amount of meaningful
information encoded.

We include a sample generated from the best performing model for reference
(See Section 2 in Appendix). A hyper-parameter tuning of learning rate and
batch side for a total of 20 epochs has final measured performance as shown in
Table 4.1.

| Batch
---
Size
| Starter
---
Learning
Rate
| Regularization
---
Weight
| BLEU Score
---
on Train
| BLEU Score
---
on Eval
64 | 0.02 | 0.1 | 0.022 | 0.012
64 | 0.02 | 0.01 | 0.023 | 0.019
64 | 0.002 | 0.1 | 0.034 | 0.028
64 | 0.002 | 0.01 | 0.037 | 0.012
64 | 0.0002 | 0.1 | 0.09 | 0.073
64 | 0.0002 | 0.01 | 0.094 | 0.014
128 | 0.02 | 0.1 | 0.024 | 0.021
128 | 0.02 | 0.01 | 0.021 | 0.013
128 | 0.002 | 0.1 | 0.033 | 0.029
128 | 0.002 | 0.01 | 0.038 | 0.011
128 | 0.0002 | 0.1 | 0.117 | 0.093
128 | 0.0002 | 0.01 | 0.113 | 0.034

###  4.2 Code Generation with GPT-2

We have been working with the publicly available small GPT-2 model with 117
million parameters. We trained using the small GPT-2 model for 100,000 mini-
batch iterations with a batch size of 2. We have included some sample code
that our model generated directly in the report. Qualitatively, our model
generates code which is far more reasonable than our baseline. The generated
code is novel, as verified by doing n-gram overlap analysis between the
generated code and the training dataset. We also note that the model learns
appropriate understanding of Python syntax, with uses of if-statements,
function and method calls, as well as regularly commented code. For full
output, see Appendix Section 6.2.

We observed that the idea of using Byte Pair encoding as used in GPT-2 is a
much better strategy to generate code than just using characters, while of
course the size of the models itself has a very observable effect in
generating Python-like code.

Overall, the GPT-2 model quickly achieves performance that’s much better than
the baseline. Continued training of the model shows that our BLEU score
performance will continue to increase, as seen in Figure 2

Figure 2: BLEU Score During Training of GPT-2 Based Model for Python Code
Generation

##  5 Conclusions

In this paper, we explore the problem of automatically completing a function
from the given function signature and human-readable documentation. We find
the best performing model to be a fine-tuned version GPT-2, a transformer-
based NLP model which is trained to generate natural text on an extremely
large dataset. Despite the fact that our model focuses specifically on code
rather than natural language, we hypothesize that it is able to treat
programming language as another specific domain alongside the encyclopedia
articles, news or books that its backbone has been trained on. We are able to
achieve a BLEU score of 0.22, improving our baseline by ¿40%.

##  6 Contributions

All team member contributed equally to this project. Baselines from the
CodeSearchNet models for code search were trained and tuned by Luis Perez and
Sudharshan Viswanathan. Data analysis and understanding of the features
(including histograms, distribution of tokens, and other data insights) was
primarily performed by Lizi Ottens.

Training of the baseline char-rnn model, as well as analysis of results and
discussion was contributed primarily by Luis Perez. Fine-tuning and training
with the small and medium GPT-2 models was primarily explored and analyzed by
Lizi Ottens and Sudharshan Viswanathan.

All written submissions were co-written by all three authors.

## References

  * [1] Alec Radford, Jeff Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever.  Language models are unsupervised multitask learners.  None, 2019.
  * [2] Kenton Lee Kristina Toutanova Jacob Devlin, Ming-Wei Chang.  BERT: pre-training of deep bidirectional transformers for language understanding.  CoRR, abs/1810.04805, 2018.
  * [3] Miltiadis Allamanis, Earl T. Barr, Premkumar T. Devanbu, and Charles A. Sutton.  A survey of machine learning for big code and naturalness.  CoRR, abs/1709.06182, 2017.
  * [4] Hamel Husain, Ho-Hsiang Wu, Tiferet Gazit, Miltiadis Allamanis, and Marc Brockschmidt.  Codesearchnet challenge: Evaluating the state of semantic code search, 2019.
  * [5] Yasir Hussain, Zhiqiu Huang, Senzhang Wang, and Yu Zhou.  Codegru: Context-aware deep learning with gated recurrent unit for source code modeling.  CoRR, abs/1903.00884, 2019.
  * [6] Zeyu Sun, Qihao Zhu, Yingfei Xiong, Yican Sun, Lili Mou, and Lu Zhang.  Treegen: A tree-based transformer architecture for code generation, 2019\.
  * [7] Xinyue Liu, Xiangnan Kong, Lei Liu, and Kuorong Chiang.  Treegan: Syntax-aware sequence generation with generative adversarial networks.  CoRR, abs/1808.07582, 2018.
  * [8] Zeyu Sun, Qihao Zhu, Lili Mou, Yingfei Xiong, Ge Li, and Lu Zhang.  A grammar-based structural CNN decoder for code generation.  CoRR, abs/1811.06837, 2018.
  * [9] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  CoRR, abs/1706.03762, 2017.
  * [10] Yoon Kim, Yacine Jernite, David A. Sontag, and Alexander M. Rush.  Character-aware neural language models.  CoRR, abs/1508.06615, 2015.
  * [11] KyungHyun Cho, Bart van Merrienboer, Dzmitry Bahdanau, and Yoshua Bengio.  On the properties of neural machine translation: Encoder-decoder approaches.  CoRR, abs/1409.1259, 2014.

## Appendix and Figures

###  6.1 CodeSearchNet Results

The Neural Bag of Words and LSTM CodeSearchNet baselines both report metrics
in the same fashion. Below, we show the training curves, which correspond to
the loss in Equation (1).

Additionally, given that the baselines models for CodeSearchNet focus on code
snippet retrieval, we also report the achieved mean reciprocal rank. The MRR
is a statistic measure for evaluating any process that produces a list of
possible responses to a sample of queries, ordered by probability of
correctness. The reciprocal rank of a query response is the multiplicative
inverse of the rank of the first correct answer: 1 for first place,
$\frac{1}{2}$ for second place, $\frac{1}{3}$ for third place and so on. The
mean reciprocal rank is the average of the reciprocal ranks of results for a
sample of queries, as in Equation (2).

$\displaystyle\text{MRR}=\frac{1}{|Q|}\sum_{i=1}^{|Q|}\frac{1}{\text{rank}_{i}}$
(2)

####  6.1.1 Neural Bag of Words Baselines

This baseline consists of a simple encoder architecture which takes as input
bag-of-words representation of the code and using a single neural network
encodes these token representation into an embedding [4]. This baseline
actually performs the best, achieving the lowest overall training and
validation losses (see Figure 5) as well as the highest MRR on the validation
set (See Figure 6).

####  6.1.2 Bi-directional RNN Model

In this model, we employ the GRU cell [11] to summarize the input sequence.
This baseline performs significantly worse, suffering from what appears to be
obvious over-fitting. In Figure 7, we can see that while the training loss
appears to plateau, the validation loss begins quickly climbing. While this
behavior does not appear to affect the overall MRR achieved on the validation
set, it is still clear that the model performs worse than the bag of words
baseline as per Figure 8.

###  6.2 Example Code

Listing 1: Sample generated using Char-RNN model.

⬇

Downloads Dailymotion videos by URL. (x or pillares (if macks), style as a
bool to you extner method instruxError,

to oio instance of that Seir to be two categorical String to mandation
:attr:‘ColumnserverName(

zr, j)

dimensed_source

also = “axis“. Element. This repr_dures a substimcle of

”””

␣␣␣␣␣␣␣␣code␣=␣item

␣␣␣␣return␣self.filename

␣␣␣␣␣␣␣␣␣␣␣␣restroxig,␣self.get_channels():

␣␣␣␣”””Get the function as firmap_list {1} and a :attracted_coordint: value of
Time

end_aes:

target_ext = int(cmd_dict):

# In data applicate horinad(string): the channel for decoded collfile

Runnee network number of element. What
you’s␣associates␣py␣ch␣of␣List␣does␣request.␣Channels␣represented␣should␣be␣a␣TargetExcept␣the␣more␣notes

␣␣␣␣␣␣␣␣“index“.

def␣vert_event_notify_channel(self,␣infna=1)

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣elif␣trying.format(

␣␣␣␣␣␣␣␣comps␣+␣’,’␣%␣(random)␣+␣1,

␣␣␣␣␣␣␣␣␣␣␣␣H.data:

␣␣␣␣␣␣␣␣␣␣␣␣if␣list␣+␣lookbing.get(’get_dict.get’]}’

Listing 2: Sample generated using small (117M) GPT-2 model after 40000
iterations.

⬇

input_snippet = def
semilocal_linear_trend_transition_matrix(autoregressive_coef):

”””Build␣the␣transition␣matrix␣for␣a␣semi-local␣linear␣trend␣model.”””,

output_code =

# Create the matrix

# (autoregressive_linear_trend_matrix) where :param
semilocal_linear_trend_matrix:

# (n_samples) is length of the semi-sparse linear regression matrix;

# i, j = semilocal_linear_trend

# Create the columns

columns = np.asarray(columns)

if auto_correlated_trend_matrix:

column = np.asarray(columns)

# Create the indices for the columns

i = column.min()

j = column.max()

# Ensure the columns have uniformity and the axis position

assert hasattr(

column, ”sigma”)

assert hasattr(axes, ”axis”)

# Create the columns

i = columns.copy()

j = columns.copy()

# Ensure the indices for the indices for the columns

ax = ColumnSetAxes(i, j)

axes.set_axes(ax)

# Create the columns

ax.set_xlabel(self.default_label)

ax.set_ylabel(self.default_ylabel)

axes.set_title(self._title)

axes.text(u”Subsim”]

# Create the column index

self.get_index(axis = columns, cols = cols)

# Create the columns if not already provided by this function

if not columns:

if auto_correlated_trend_matrix:

ax.columns = columns

else:

if i < 0 or j < 0 or i+1 <= i

####  6.2.1 Figures

Figure 3: Histogram of the the number of (comment, code) pairs available in
our dataset, as well as the number of unique function methods for each
language.

Figure 4: General CodeSearchNet architecture for all of our baselines. Each
language is processed through different encoder mechanisms. The query encoder
is shared (an NLP encoder), and the purpose of the CodeSearchNet tasks is to
retrieve the most relevant code snippets subject to the natural language
query.

Figure 5: Training and Validation losses for the Neural Bag of Words model in
CodeSearchNet. Figure 6: MRR on validation set for the baseline neural bag of
words model in the CodeSearchNet Challenge.

Figure 7: Training and Validation losses for the RNN model in CodeSearchNet.
Figure 8: MRR on validation set for the baseline RNN in the CodeSearchNet
Challenge.

Figure 9: Training and Validation Losses on the Baseline Char-RNN Model. This
is the cross-entropy loss over 128 predicted character sequence.

Figure 10: Training and Validation Perplexity on the Baseline Char-RNN Model.
This is the cross-entropy loss over 128 predicted character sequence.