# Charge kinetics across a negatively biased semiconducting plasma-solid
interface

K. Rasek, F. X. Bronold and H. Fehske  Institut für Physik, Universität
Greifswald, 17489 Greifswald, Germany

(August 27, 2024)

###### Abstract

An investigation of the selfconsistent ambipolar charge kinetics across a
negatively biased semiconducting plasma-solid interface is presented. For the
specific case of a thin germanium layer with nonpolar electron-phonon
scattering, sandwiched between an Ohmic contact and a collisionless argon
plasma, we calculate the current-voltage characteristic and show that it is
affected by the electron microphysics of the semiconductor. We also obtain the
spatially and energetically resolved fluxes and charge distributions inside
the layer, visualizing thereby the behavior of the charge carriers responsible
for the charge transport. Albeit not quantitative, because of the crude model
for the germanium band structure and the neglect of particle-nonconserving
scattering processes, such as impact ionization and electron-hole
recombination, which at the energies involved cannot be neglected, our results
clearly indicate (i) the current through the interface is carried by rather
hot carriers and (ii) the perfect absorber model, often used for the
description of charge transport across plasma-solid interfaces, cannot be
maintained for semiconducting interfaces.

##  I Introduction

Low-temperature gas discharges are bound by solid objects, acting either as
confining walls or as electrodes. To maintain the discharge, an electric
current has to flow across the plasma-electrode interface. Since the plasma of
the discharge contains electrons, ions, and radicals, the charge transfer
across the interface is ambipolar, consisting of electrons deposited into the
electrode and extracted electrons (holes) arising from the neutralization of
ions and the de-excitation of radicals. Hence, inside the electrode, a flux of
electrons and holes builds up whose fate depends on the electron microphysics
of the electrode material.

The transport scenario just described is in some sense obvious, hardly
addressed in textbooks on plasma discharges Franklin (1976); Lieberman and
Lichtenberg (2005), and of course qualitatively known since the beginning of
gaseous electronics Langmuir and Mott-Smith (1924). On a fundamental level,
however, it implies a subtle interplay of gaseous and solid state transport
processes in any man-made gas discharge. Its investigation may thus perhaps
bear novel possibilities for controlling discharges by manipulating processes
inside the electrodes. The electric breakdown in dielectric barrier
discharges, for instance, depends on the charge distribution inside the
dielectric Nemschokmichal  _et al._ (2018); Peeters  _et al._ (2016); Wagner
_et al._ (2005); Massines  _et al._ (2003), and hence on the transport
processes to which they give rise to. Revealing how gaseous and solid-based
charge transport merge at the plasma-solid interface may thus allow to control
the breakdown by a judicious choice of the dielectric. It may also suggest
optimization strategies for large scale industrial barrier discharges
Kogelschatz (2003). However, a quantitative description of charge transport
across the interface will be most beneficial for the further development of
microdischarges embedded in semiconducting substrates Chiang  _et al._ (2019);
Eden  _et al._ (2013); P. A. Tchertchian, C. J. Wagner, T. J. Houlahan Jr., B.
Li, D. J. Sievers, and J. G. Eden (2011); Dussart  _et al._ (2010), where the
time and length scales of electron transport and energy relaxation are no
longer well separated (see the Introduction of Ref. Bronold and Fehske (2017)
for a discussion of this point). It is thus the purpose of this work to
provide first steps towards a selfconsistent kinetic description of the
ambipolar charge transport across a biased semiconducting plasma-solid
interface.

In an attempt to model the whole electric double layer forming at a plasma-
solid interface, and not only the positive space charge on the plasma side
(that is, the plasma sheath Schwager and Birdsall (1990); Riemann (1991);
Franklin (2003); Brinkmann (2009); Robertson (2013)), we recently developed a
model for floating dielectric interfaces that treats the electrons and holes
in the solid on the same kinetic footing as the electrons and ions in the
plasma Rasek  _et al._ (2020); Bronold and Fehske (2017). The model links the
electron-ion plasma to the electron-hole plasma inside the dielectric by
allowing electrons to cross the interface in both directions. Electrons from
the plasma may thus not only enter the conduction band of the solid by
traversing the surface potential but also leave it due to internal
backscattering and subsequent traversal of the surface potential in the
reverse direction. In addition, electrons can be extracted from the valence
band by the neutralization of ions. Tracking the charge distributions by a set
of Boltzmann equations, the resulting charge imbalance can be determined and
used as a source in the Poisson equation, determining the selfconsistent
electric field, which in turn influences the kinetics of the charge carriers
on both sides of the interface. The selfconsistent solution of the Boltzmann-
Poisson system for the distribution functions and the electric potential is
thus at the core of our approach. For the floating interface, where no net
flux is flowing through the interface, we have solved this set of equations
under simplifying assumptions Rasek  _et al._ (2020). We now extend the model
to an interface which carries a net current and remedy also some of the
limiting simplifications used before.

As a first step towards a realistic treatment of the current-carrying
interface, we consider a planar semiconductor of finite thickness, in which
electrons and holes loose or gain energy by scattering on optical phonons,
sandwiched between an Ohmic contact and a collisionless plasma. At the
interface between the plasma and the semiconductor, electrons may be reflected
when impinging on the interface from either side, while ions are neutralized
with unit probability. At the other end of the semiconductor, the Ohmic metal
serves as a sink for electrons and holes. The resulting setup resembles thus a
Langmuir probe Mott-Smith and Langmuir (1926); Cherrington (1982); Lam (1965)
coated with a semiconducting layer. A negatively biased probe attracts fluxes
of electrons and ions, the magnitude of which depends on the bias voltage. At
large enough negative bias, electrons cannot reach the probe anymore,
resulting in a negative net flux due solely to ions. Otherwise the flux is
dominated by the electrons, due to their lower mass and higher temperature. By
calculating the functional dependence of the net flux on the bias voltage,
that is, the current-voltage characteristic, and comparing it with the
characteristic of a perfectly absorbing interface, which assumes that any
charge carriers hitting the interface get instantaneously absorbed and never
enter the plasma again, we can determine the influence the electron
microphysics of the semiconducting layer has on the electric current flowing
through the device. Our simulations show that the assumptions of the perfect
absorber model cannot be maintained. The charge kinetics inside the
semiconductor is an essential part of the kinetics of the gas discharge and
should thus be included in its modeling.

The outline of the paper is as follows. In the two parts of Sec. II we present
the equations of the kinetic model and discuss its numerical implementation
and solution, focusing on aspects which differ from our treatment of the
floating interface. Using for illustration germanium in contact with an argon
plasma, numerical results are given in Sec. III. It is divided into three
parts, corresponding to the three perspectives from which one may consider the
interface. First, in subsection III.1, it is regarded as an electric device,
discussing thus the current-voltage characteristic that results from the
kinetic theory. Spatially resolved macroscopic properties, such as density
distributions and potential profiles, are discussed in subsection III.2, while
a microscopic view, based on spatially and energetically resolved distribution
functions is presented in subsection III.3. The paper concludes in Sec. IV
with an outlook to what could be the next steps.

##  II Theoretical framework

Figure 1: Illustration of the electric potential across a negatively biased
semiconducting plasma-solid interface (not to scale). The bias voltage $U_{B}$
is applied between the bulk plasma and the Ohmic contact to the left of the
semiconducting layer, as sketched in the bottom right. Measuring the potential
with respect to its value in the plasma bulk, $U_{B}=U_{I}-U_{SC}$, with
$U_{I}$ the sheath (interface) potential at $z=0$ and $U_{SC}$ the potential
drop across the semiconducting layer. In the Schwager-Birdsall approach
Schwager and Birdsall (1990) to the collisionless plasma, the presheath
potential $U_{W}$ stretches over the entire region $z>z_{p}$, where $z_{p}$ is
the location of the bulk plasma, which is effectively infinitely far away from
the interface.

###  II.1 Model equations

The basis for the kinetic modeling of the current-carrying interface is a set
of Boltzmann equations, describing the dynamics of the distribution functions
of the charge carriers on both sides of the interface, augmented by the
Poisson equation for the electric potential and matching as well as boundary
conditions for the distribution functions and the potential Bronold and Fehske
(2017). Due to the particle flux through the interface, the boundary
conditions differ from the ones used for a floating interface. To make the
modeling more realistic, we improve in this work also the matching conditions
at the plasma-semiconductor interface, considering now realistic injection
energies for holes and the possibility for electrons to be quantum-
mechanically reflected. If not noted otherwise, all equations are written in
atomic units, measuring energy in Rydbergs, length in Bohr radii, and masses
in electron masses.

Figure 1 shows the electric potential energy $U(z)$ across the interface in
the manner it is implemented in the model of this paper. Within the plasma
($z>0$), we employ the Schwager-Birdsall approach Schwager and Birdsall (1990)
for a collisionless plasma to model the merging of the plasma sheath with the
bulk plasma. The sheath potential $U_{I}=U(0)$ is thus the potential
difference between the interface at $z=0$ and a point $z_{p}$, where the bulk
plasma is established and which also serves as the reference point from which
electric potentials are measured. The presheath potential $U_{W}=U(z_{w})$
accelerates the ions to make the Bohm criterion satisfiable at $z_{p}$. Since
$z_{p}$ tends to infinity, the presheath does not belong to the physically
relevant part of the plasma-solid interface. In the following, we adopt the
sign convention that $U_{B}=U(z_{0})$ and $U_{I}$ are negative, whereas
$U_{W}$ and $U_{SC}$ are positive. Hence,

$\displaystyle U_{B}=U_{I}-U_{SC}\leavevmode\nobreak\ .$ (1)

When no collisions are considered within the plasma, the particle densities
can readily be expressed as functions of $U(z)$. It is thus practical to use
the once integrated Poisson equation to calculate the electric field in the
form

$\frac{\mathrm{d}U}{\mathrm{d}z}=\mathcal{E}(z)=\left(16\pi\int_{U(z_{p})}^{U(z)}\mathrm{d}U^{\prime}n(U^{\prime})\right)^{1/2}\leavevmode\nobreak\
,$ (2)

with $U(z_{p})=0$ and also $\mathcal{E}(z_{p})=0$ since the bulk plasma is
field-free.

The semiconducting material, having a dielectric constant $\varepsilon$,
stretches from $z=z_{0}<0$ to $z=0$. It has thus a fixed width and it is more
reasonable to keep $z$ as the spatial variable. Hence,

$\mathcal{E}(0^{-})-\mathcal{E}(z)=\frac{8\pi}{\varepsilon}\int_{z}^{0}\mathrm{d}z^{\prime}n(z^{\prime})\leavevmode\nobreak\
.$ (3)

Using the matching condition for the electric potential at the interface,

$\varepsilon\mathcal{E}(0^{-})=\mathcal{E}(0^{+})\leavevmode\nobreak\ ,$ (4)

the magnitude of the electric field, and thus the electric potential, across
the entire interface can be determined for a given charge density

$n(z)=\left[n_{e}(z)-n_{i}(z)\right]\Theta(z)\\\
+\left[n_{*}(z)-n_{h}(z)+n_{A}-n_{D}\right]\Theta(-z)\leavevmode\nobreak\ ,$
(5)

where $n_{s}(z)$ denotes the density of electrons $(s=e)$, ions $(s=i)$,
conduction band electrons $(s=*)$, and valence band holes $(s=h)$. For an
intrinsic semiconductor, the acceptor ($n_{A}$) and donor densities ($n_{D}$)
are absent. Inside the metal, the electric potential is constant and the
electric field vanishes.

The Boltzmann equations governing the distribution functions
$F_{s}^{\gtrless}(z,E,T)$ for the charge carriers on either side of the
interface are given by Rasek  _et al._ (2020)

$\pm v_{s}\frac{\partial}{\partial
z}F^{\gtrless}_{s}=\Phi^{\gtrless}_{s}-\gamma^{\gtrless}_{s}F^{\gtrless}_{s}\leavevmode\nobreak\
,$ (6)

with $\Phi^{\gtrless}_{s}$ the in-scattering part of the collision integral,
$\gamma^{\gtrless}_{s}$ the scattering rate, both will be specified below, and

$v_{s}(z,E,T)=2\sqrt{m_{s}^{-1}(E-U_{s}(z)-T)}\leavevmode\nobreak\ $ (7)

the modulus of the velocity in $z-$direction.

For brevity, the independent variables, which are the total energy $E$, the
lateral kinetic energy $T$, and the spatial variable $z$ are suppressed in Eq.
(6) and the distribution functions for left- and right-moving particles are
distinguished by the superscripts $<$ and $>$. In the lateral directions the
interface is isotropic. Equation (7) holds only for parabolic dispersions,
which are of course valid for the free charge carriers on the plasma side of
the interface, but for the free carriers inside the solid it is an
approximation specified by effective masses. Anticipating to use for
illustration germanium, we take for the effective electron mass the density-
of-state effective mass, $m_{*}=\left(9m_{l}m_{t}^{2}\right)^{\frac{1}{3}}$,
and for the hole mass the average of the masses of light and heavy holes. With
the numerical values from Ref. Jacoboni and Reggiani (1983), we then obtain
the masses given in Table 1. However, since the injection energies for
electrons and holes are rather high, the parabolic dispersion is only a crude
approximation to the band structure. More than one valley as well as
nonparabolicities should be actually considered. But it is beyond the scope of
this exploratory work. The species potential $U_{s}(z)$, finally, takes each
species’ charge and energy offset into account, relating thus to the electric
potential $U(z)$ via $U_{i}=U$, $U_{e}=-U$, $U_{*}=-U-\chi$, and
$U_{h}=U+E_{g}+\chi$, with the electron affinity $\chi$ and band gap $E_{g}$.

Table 1: Material parameters for the germanium layer Jacoboni and Reggiani (1983) and the argon plasma facing it.  $k_{B}T_{SC}(\text{eV})$ | 0.025 | $m_{*}(m_{e})$ | 0.2
---|---|---|---
$k_{B}T_{i}(\text{eV})$ | 0.025 | $m_{h}(m_{e})$ | 0.34
$k_{B}T_{e}(\text{eV})$ | 2 | $m_{i}(m_{e})$ | 73551
$D_{t}K(10^{8}\textrm{eV/cm})$ | 9.5 | $n_{i}(10^{13}\textrm{cm}^{-3})$ | 2
$n_{p}(\textrm{cm}^{-3})$ | $10^{13}$ | $\rho(\textrm{g/cm}^{3})$ | 5.32
$\varepsilon$ | 16.2 | $\hbar\omega_{0}(\text{eV})$ | 0.037
$E_{g}(\text{eV})$ | 0.67 | $\chi(\textrm{eV})$ | 4
$I_{\textrm{inj}}(\text{eV})$ | 15.76 | $\Gamma_{\mathrm{inj}}(\textrm{eV})$ | 3

While no collisions are considered on the plasma side, implying
$\Phi_{e,i}=\gamma_{e,i}=0$, within the semiconductor, we include collisions
with optical phonons. For nonpolar materials, such as silicon or germanium,
the phonon collision integral is isotropic. Hence, no distinction between
left- and right-moving particles must be taken into account, implying
$\Phi_{s}^{>}=\Phi_{s}^{<}=\Phi_{s}$ and
$\gamma_{s}^{>}=\gamma_{s}^{<}=\gamma_{s}$. The collision rates entering the
Boltzmann equations for conduction band electrons and valence band holes are
then Ridley (1999); Roth (1992)

$\gamma_{s}(z,E)=\frac{(D_{t}K)^{2}m_{s}}{4\pi\rho\omega_{0}}\sqrt{m_{s}}\left[n_{b}\sqrt{E+\hbar\omega_{0}-U_{s}}\right.\\\
\left.+(n_{b}+1)\sqrt{E-\hbar\omega_{0}-U_{s}}\right]\leavevmode\nobreak\ ,$
(8)

where the second term in the square brackets only applies if the argument of
the root is positive, while the in-scattering parts of the collision integrals
read

$\Phi_{s}(z,E)=\frac{(D_{t}K)^{2}}{8\pi\rho\omega_{0}}\left[n_{b}N_{s}(z,E-\hbar\omega_{0})\right.\\\
\left.+(n_{b}+1)N_{s}(z,E+\hbar\omega_{0})\right]\leavevmode\nobreak\ ,$ (9)

with the optical deformation potential $D_{t}K$, the mass density $\rho$, the
optical phonon frequency $\omega_{0}$, the phonon occupation number
$n_{b}=1/\left(\exp(\hbar\omega_{0}/k_{B}T_{SC})-1\right)$, and

$N_{s}(z,E)=m_{s}\sum_{\gtrless}\int_{0}^{E-U_{s}(z)}\mathrm{d}T\frac{F_{s}^{\gtrless}(z,E,T)}{v_{s}(z,E,T)}\leavevmode\nobreak\
,$ (10)

the spatially and energetically resolved density of the species $s$, from
which the densities entering the Poisson equation (3) follow by one more
integration,

$n_{s}(z)=\int_{U_{s}(z)}^{\infty}\frac{\mathrm{d}E}{8\pi^{2}}N_{s}(z,E)\leavevmode\nobreak\
.$ (11)

For germanium, the material parameters required for $\gamma_{*,h}$ and
$\Phi_{*,h}$ are given in Table 1.

As in our previous work Rasek  _et al._ (2020), we solve the equations on the
plasma side analytically and use the iterative approach by Grinberg and Luryi
Grinberg and Luryi (1992) inside the semiconductor. The boundary conditions at
the outer limits of the kinetically modeled interval are essential to
determine the solution. On the right boundary, at $z=z_{w}$, the Schwager-
Birdsall model Schwager and Birdsall (1990) prescribes half Maxwellian
distributions

$F_{s}^{\mathrm{LM}}(z_{w})=n_{s}^{\mathrm{LM}}\left(\frac{4\pi}{k_{B}T_{s}m_{s}}\right)^{3/2}\exp\left(-\frac{E-U_{s}(z_{w})}{k_{B}T_{s}}\right)\leavevmode\nobreak\
,$ (12)

with densities $n_{i,e}^{\textrm{LM}}$ such that at $z=z_{p}$ the densities
for both electrons and ions are equal to the plasma density $n_{p}$. On the
left boundary, at $z=z_{0}$, that is, at the metal-semiconductor interface,
the boundary condition is not strictly known. Following standard semiconductor
device modeling Moglestue (1993), we assume an Ohmic contact, implying at
$z=z_{0}$ Maxwellian distributions with densities and temperatures applying to
the bulk of the semiconductor. For an undoped, intrinsic semiconductor, for
instance, the boundary condition at $z=z_{0}$ is thus given by (12) with
$z_{w}$ replaced by $z_{0}$, $T_{s}=T_{SC}$, and $n_{s}^{\mathrm{LM}}=n_{\rm
int}$, where $n_{\rm int}$ is the intrinsic density of the semiconductor. Put
together, we thus have to enforce at the system boundaries,

$\displaystyle F_{s}^{>}(z_{0})=F_{s}^{\mathrm{LM}}(z_{0})\hskip
30.0pt\text{for }$ $\displaystyle s=h,*\leavevmode\nobreak\ ,$ (13a)
$\displaystyle F_{s}^{<}(z_{w})=F_{s}^{\mathrm{LM}}(z_{w})\hskip
30.0pt\text{for }$ $\displaystyle s=i,e\leavevmode\nobreak\ .$ (13b)

We also need to match distribution functions at the interface between the
plasma and the semiconductor. For ions, we assume perfect neutralization, that
is, an impinging ion extracts (injects) with unit probability an electron
(hole) from the valence band. For electrons, the interface is quantum-
mechanically reflecting with a reflection coefficient $R(E,T)$. Thus, the
matching conditions at $z=0$ read

$\displaystyle F_{i}^{>}(0,E,T)$ $\displaystyle=0$ (14a) $\displaystyle
F_{h}^{<}(0,E,T)$ $\displaystyle=F_{h}^{>}(0,E,T)+S_{h}^{<}(E,T)$ (14b)
$\displaystyle F_{*}^{<}(0,E,T^{\prime})$
$\displaystyle=R(E,T)F_{*}^{>}(0,E,T^{\prime})+S_{*}^{<}(E,T^{\prime})$ (14c)
$\displaystyle F_{e}^{>}(0,E,T)$ $\displaystyle=R(E,T)F_{e}^{<}(0,E,T)$
$\displaystyle+(1-R(E,T))F_{*}^{>}(0,E,T^{\prime})$ (14d)

with source terms $S^{<}_{s}$, describing injection of holes and electrons
into the semiconductor, given by

$S_{h}^{<}(E,T)=n_{h}^{\mathrm{inj}}\left(\frac{4\pi}{k_{B}T_{h}m_{h}}\right)^{\frac{3}{2}}\\\
\times\exp\left(-\frac{(E-U_{h}(0)-I_{\mathrm{inj}})^{2}}{\Gamma_{\mathrm{inj}}^{2}}\right)\leavevmode\nobreak\
,$ (15a) which is in fact independent of $T$, and
$S_{*}^{<}(E,T^{\prime})=(1-R(E,T))F_{e}^{<}(0,E,T)\leavevmode\nobreak\ .$
(15b)

In the matching conditions for the electron distribution function,
$T^{\prime}=T/m_{*}$. The change in lateral energy from $T$ to $T^{\prime}$
arises from the conservation of lateral momentum. Since, the effective mass
$m_{*}<m_{e}$, the electron gains (looses) lateral energy while passing
through the interface from the plasma (solid) side. The mass mismatch leads
also to total reflection for electrons coming from the plasma when
$E-U_{*}(0)<T^{\prime}$.

Figure 2: (color online) Energy-resolved source functions $S_{*}^{<}$ (left)
and $S_{h}^{<}$ (right) in arbitrary units for the floating interface, that
is, the situation where electron and ion fluxes are equal. The material
parameters are taken from Table 1, the thickness of the germanium layer is 1
$\mu\mathrm{m}$, $U_{W}=0.97\,\mathrm{eV}$, and $U_{B}=-5.125\,\mathrm{eV}$.
Note, $T^{\prime}>E-U_{*}(0)$ and $T>E-U_{h}(0)$ are energetically not
allowed. Since the injection of holes spreads over a larger energy range, the
absolute values of $S_{h}^{<}$ are smaller than the values of $S_{*}^{<}$.
Integrated over energy, however, the source functions ensure flux equality, as
required for the floating interface.

Both source terms are illustrated in Fig. 2 for the parameters given in Table
1, which are used for the numerical calculations described in the next
section. The normalization density $n_{h}^{\textrm{inj}}$ in the source term
for the holes is chosen such that the flux is conserved across the interface.
With this source term, holes are injected at the ionization energy
$I_{\textrm{inj}}$ of an argon atom, homogeneously distributed in lateral
direction. The width $\Gamma_{\textrm{inj}}$ accounts for an energy spread in
the neutralization process. The source for electrons depends on the reflection
coefficient $R(E,T)$. For a surface potential with the depth of the electron
affinity $\chi$ and a $1/z$ tail on the plasma side due to the image charge
(Schottky effect), the coefficient becomes, adapting results from Ref. MacColl
(1939),

$R(E,T)=\left|\frac{v_{e}(0,E,T)-y\,v_{*}(0,E,T^{\prime})}{v_{e}(0,E,T)+y^{*}\,v_{*}(0,E,T^{\prime})}\right|^{2}\leavevmode\nobreak\
$ (16)

with

$y=-2\frac{W^{\prime}_{\lambda,\frac{1}{2}}(\xi_{0})}{W_{\lambda,\frac{1}{2}}(\xi_{0})}\leavevmode\nobreak\
,$ (17)

where $W_{\lambda,\frac{1}{2}}(x)$ denotes the Whittaker function,
$W^{\prime}_{\lambda,\frac{1}{2}}(x)$ its first derivative with respect to
$x$,

$\displaystyle\lambda$
$\displaystyle=-\mathrm{i}\frac{\varepsilon-1}{\varepsilon+1}\frac{1}{\sqrt{8}\sqrt{E-T-
U_{I}}}\leavevmode\nobreak\ ,$ (18) $\displaystyle\xi_{0}$
$\displaystyle=\frac{\mathrm{i}\sqrt{2}}{\chi}\frac{\varepsilon-1}{\varepsilon+1}\sqrt{E-T-
U_{I}}\leavevmode\nobreak\ ,$ (19)

and $y^{*}$ is the complex conjugate of $y$. As can be seen in Fig. 2, due to
the high temperature, electrons are injected into the conduction band of the
semiconductor over a wide range of energies, with most weight at the low
energy cutoff given by $E-U_{*}(0)=\chi$.

Since the plasma is treated collisionless, the solutions of the Boltzmann
equations on the plasma side are completely determined by the profile of the
electric potential and the boundary conditions at $z=0$ and $z=z_{w}$. The
latter are given by (13b), where the ions are restricted to energies above the
presheath potential $U_{W}$, which needs to be determined selfconsistently. It
is responsible for the acceleration of ions before reaching the sheath and can
be determined from the generalized Bohm criterion Riemann (1991),

$\displaystyle\partial_{U}(n_{e}-n_{i})|_{z_{p}}\geq 0\leavevmode\nobreak\ .$
(20)

In our formalism, it follows from Eq. (2), where the radicand needs to be
positive for $z\lesssim z_{p}$. As usual, we enforce marginal fulfillment.
Thus, $U_{W}$ is determined by Eq. (20) with the equal sign. In this work, we
will prescribe the plasma density $n_{p}$. Enforcing Eq. (20) instead of
$\mathcal{E}(z_{w})=0$, as in Rasek  _et al._ (2020), is then numerically
advantageous.

For the current-voltage characteristic we need the particle fluxes. As for the
densities given in Eq. (10), we initially define energy-resolved fluxes,

$J_{s}^{\gtrless}(z,E)=m_{s}\int_{0}^{E-U_{s}(z)}\mathrm{d}TF_{s}^{\gtrless}(z,E,T)\leavevmode\nobreak\
,$ (21)

in terms of which the macroscopic fluxes required for the characteristic
become

$j_{s}(z)=\int_{U_{s}(z)}^{\infty}\frac{\mathrm{d}E}{8\pi^{2}}\left[J_{s}^{>}(z,E)-J_{s}^{<}(z,E)\right]\leavevmode\nobreak\
.$ (22)

###  II.2 Numerical strategy

In the following, we give a sketch of the numerical strategy used for solving
the kinetic problem stated in the previous subsection. Due to the changes in
the boundary and matching conditions, the strategy differs somewhat from the
one used previously Rasek  _et al._ (2020). In particular, the procedure for
establishing selfconsistency between the plasma and the solid side of the
interface is different due to the reflectivity of the interface. The isotropy
of the collision integrals enables us moreover to discretize a much larger
energy domain.

Besides the distribution functions $F_{s}^{\gtrless}(z,E,T)$ and the electric
potential profile $U(z)$, three energy parameters $U_{W}$, $U_{I}$, and
$U_{SC}$ and two density parameters $n^{\textrm{LM}}_{e}$ and
$n^{\textrm{LM}}_{i}$ have to be selfconsistently determined for prescribed
plasma density $n_{p}$ and external voltage $U_{B}$. The five equations
required for it are, the charge neutrality condition at $z=z_{p}$, providing
two equations, $n_{p}=n_{e}(z_{p})=n_{i}(z_{p})$, the electric matching
condition (4), the generalized Bohm criterion (20), and the condition (1)
following from the definitions of the potential drops at the interface. The
net flux through the interface,

$\displaystyle j\equiv j_{*}(0)-j_{h}(0)=j_{e}(0)-j_{i}(0)\leavevmode\nobreak\
,$ (23)

due to flux conservation identical to the flux anywhere in the device, is then
obtained as a function of $U_{B}$. It will be however numerically advantageous
to specify $U_{SC}$ instead of $U_{B}$ and to initially determine $j(U_{SC})$
from which $j(U_{B})$ follows straight by applying Eq. (1).

For the numerical implementation of the selfconsistent calculation of the
current-voltage characteristic we rewrite the Boltzmann equations (6) for the
charge carriers inside the semiconductor in integral form (suppressing the
parametrical dependencies on $E$ and $T$) Rasek  _et al._ (2020),

$F_{s}^{>}(z)=\xi_{s}(z,z-\Delta)F_{s}^{>}(z-\Delta)\\\
+\int_{z-\Delta}^{z}\mathrm{d}z^{\prime}\frac{\Phi_{s}(z^{\prime})}{v_{s}(z^{\prime})}\xi_{s}(z,z^{\prime})$
(24a) and $F_{s}^{<}(z)=\xi_{s}(z+\Delta,z)F_{s}^{<}(z+\Delta)\\\
+\int_{z}^{z+\Delta}\mathrm{d}z^{\prime}\frac{\Phi_{s}(z^{\prime})}{v_{s}(z^{\prime})}\xi_{s}(z^{\prime},z)$
(24b)

with the integrating factor

$\xi_{s}(z,z^{\prime})=\exp\left(-\int_{z^{\prime}}^{z}\mathrm{d}\bar{z}\frac{\gamma_{s}(\bar{z})}{v_{s}(\bar{z})}\right)\leavevmode\nobreak\
$ (25)

and $s=*,h$.

To avoid the integrable divergences of $1/v_{s}(z,E,T)$ at $T=E-U_{s}(z)$,
that is, at the turning points for the perpendicular motion, where the charge
carriers move parallel to the interface, it is convenient to perform a
coordinate transformation from $z$ to $X_{s}=\sqrt{E-U_{s}(z)-T}$. Using
$\partial U_{h}/\partial z=\mathcal{E}$, $\partial U_{*}/\partial
z=-\mathcal{E}$, and $v_{s}=2X_{s}/\sqrt{m_{s}}$, we find

$\displaystyle\frac{\mathrm{d}z}{v_{h}}$
$\displaystyle=-\sqrt{m_{h}}\frac{\mathrm{d}X_{h}}{\mathcal{E}}\leavevmode\nobreak\
,$ (26a) $\displaystyle\frac{\mathrm{d}z}{v_{*}}$
$\displaystyle=\sqrt{m_{*}}\frac{\mathrm{d}X_{*}}{\mathcal{E}}\leavevmode\nobreak\
,$ (26b)

which enables us to rewrite the $z-$integrals as $X_{s}-$integrals. If the
semiconducting layer is not too thick, the electric field $\mathcal{E}$ is
finite in the whole integration domain. Thus, using this substitution we can
avoid diverging integrands in the numerical solution of the Boltzmann
equations.

After the transformation, the coordinates $E$, $T$, and $X_{s}$ are
discretized, with the discretization kept fixed during the iteration. Through
the domain of $X_{*,h}$, we thus prescribe the value $U_{SC}$, although it is
actually a parameter to be determined selfconsistently for given $U_{B}$.
Instead of $U_{SC}$, we use $U_{B}$ as a derived parameter, which allows a
better handling of the turning points, where $v_{*,h}=0$ and
$F_{*,h}^{>}=F_{*,h}^{<}$. Due to Eq. (1) this is permissible.

The remaining four parameters, $U_{I}$, $U_{W}$, $n_{e}^{\mathrm{LM}}$, and
$n_{i}^{\mathrm{LM}}$ are determined from the charge neutrality at $z=z_{p}$,
bringing in two equations and yielding $n_{e,i}^{\mathrm{LM}}$ in terms of
$U_{W}$ and $U_{I}$, the matching condition (4), and the generalized Bohm
criterion (20). The latter two provide at the end two coupled equations for
$U_{I}$ and $U_{W}$ which have to be solved selfconsistently with the
Boltzmann equations and the Poisson equation on both sides of the interface.

To get the two equations, we relate the potential drop $U_{SC}$ to
$\mathcal{E}(z)$ by the integral

$U_{SC}=\int_{z_{0}}^{0}\mathrm{d}z\mathcal{E}(z)\leavevmode\nobreak\ .$ (27)

Inserting $\mathcal{E}(z)$ from Eq. (3), setting $n_{A}=n_{D}=0$, since we
consider in the next section an undoped germanium layer, and solving for
$\mathcal{E}(0^{-})$ yields

$\mathcal{E}(0^{-})=\left[U_{SC}+\frac{8\pi}{\varepsilon}\int_{z_{0}}^{0}\mathrm{d}z\int_{z}^{0}\mathrm{d}z^{\prime}(n_{*}(z^{\prime})-n_{h}(z^{\prime}))\right]/z_{0}\leavevmode\nobreak\
,$ (28)

that is, the electric field at $z=0$ necessary to produce, with the net charge
distribution $n_{*}(z)-n_{h}(z)$, the potential drop $U_{SC}$ over the width
of the semiconductor.

The electric field $\mathcal{E}(0^{+})$, on the other side of the interface,
in turn can be obtained from Eq. (2). Inserting the electron and ion densities
arising from the electron and ion distribution functions $F_{e,i}^{\gtrless}$,
which can be largely worked out analytically Rasek  _et al._ (2020), we get

$\displaystyle\mathcal{E}(0^{+})=\sqrt{16\pi\left(G_{e}+G_{e}^{>}+G_{i}\right)}$
(29)

with

$\displaystyle G_{e}$
$\displaystyle=\frac{n_{p}-n_{e}^{>}}{1+\mathrm{erf}\left(\sqrt{-\tilde{\beta}_{e}\tilde{U}_{I}}\right)}\left[\exp\left(\tilde{\beta}_{e}\tilde{U}_{I}\right)\left(1+2\sqrt{-\frac{\tilde{\beta}_{e}\tilde{U}_{I}}{\pi}}\right)-1-\mathrm{erf}\left(\sqrt{-\tilde{\beta}_{e}\tilde{U}_{I}}\right)\right]\leavevmode\nobreak\
,$ (30) $\displaystyle G_{e}^{>}$
$\displaystyle=\frac{m_{e}^{3/2}}{32\pi^{2}}\int_{0}^{\infty}\mathrm{d}E\int_{0}^{E}\mathrm{d}TF_{e}^{>}(z_{p},E,T)\left(\sqrt{E+U_{I}-T}-\sqrt{E-T}\right)\leavevmode\nobreak\
,$ (31) $\displaystyle G_{i}$
$\displaystyle=\frac{n_{p}}{\mathrm{erfc}\left(\sqrt{\tilde{U}_{W}}\right)}\bigg{\\{}\mathrm{erfc}\left(\sqrt{\tilde{U}_{W}}\right)-\exp\left(-\tilde{U}_{I}\right)\mathrm{erfc}\left(\sqrt{\tilde{U}_{I}-\tilde{U}_{W}}\right)+2\exp\left(-\tilde{U}_{W}\right)\left(\sqrt{\frac{\tilde{U}_{W}}{\pi}}-\sqrt{\frac{\tilde{U}_{W}-\tilde{U}_{I}}{\pi}}\right)\bigg{\\}}\leavevmode\nobreak\
,$ (32)

where we used $U_{e}(z_{p})=0$, because of the choice of the reference point
for the potentials, set $\tilde{U}_{I,W}=U_{I,W}/k_{B}T_{i}$, and
$\tilde{\beta}_{e}=k_{B}T_{i}/k_{B}T_{e}$. The symbols $\mathrm{erf}(x)$ and
$\mathrm{erfc}(x)=1-\mathrm{erf}(x)$ denote the error and complementary error
function.

Multiplying Eq. (28) by the dielectric constant $\varepsilon$ and equating it
with Eq. (29), an equation is obtained relating $U_{SC}$ to $U_{W}$, $n_{p}$,
and $U_{I}$. Since $n_{p}$ and $U_{SC}$ are external parameters, we thus have
a relation between $U_{W}$ and $U_{I}$, as required. For a perfectly absorbing
interface $G_{e}^{>}=0$.

A second relation between $U_{W}$ and $U_{I}$ follows from the generalized
Bohm criterion. Expressing again $n_{e}$ and $n_{i}$ in terms of the
distribution functions $F_{e,i}^{\gtrless}$, and using the definitions
introduced above, we obtain

$\displaystyle\frac{f_{+}\left(-\tilde{\beta}_{e}\tilde{U}_{I}\right)}{k_{B}T_{e}}\left(1-\frac{n_{e}^{>}}{n_{p}}\right)+\frac{\Delta_{e}^{>}}{n_{p}}-\frac{f_{-}\left(\tilde{U}_{W}\right)}{k_{B}T_{i}}=0$
(33)

with

$\displaystyle f_{\pm}(x)$ $\displaystyle=\frac{e^{-x}}{\sqrt{\pi
x}(1\pm\mathrm{erf}(\sqrt{x}))}\pm 1\leavevmode\nobreak\ $ (34)

and

$\displaystyle\Delta_{e}^{>}$
$\displaystyle=-\frac{m_{e}^{3/2}}{32\pi^{2}}\int_{0}^{\infty}\mathrm{d}E\int_{0}^{E}\mathrm{d}T\frac{F_{e}^{>}(z_{p},E,T)}{(E-T)^{3/2}}\leavevmode\nobreak\
,$ (35) $\displaystyle n_{e}^{>}$
$\displaystyle=\frac{m_{e}}{8\pi^{2}}\int_{0}^{\infty}\mathrm{d}E\int_{0}^{E}\mathrm{d}T\frac{F_{e}^{>}(z_{p},E,T)}{\sqrt{E-T}}\leavevmode\nobreak\
,$ (36)

where $U(z_{p})=0$ has been used once again. We thus have a second equation
connecting $U_{W}$ and $U_{I}$. Note, through Eq. (14d), the electron
distribution function $F_{e}^{>}$ on the plasma side depends on the electron
distribution function $F_{*}^{>}$ inside the solid, providing an additional
feedback of the solid to the plasma, in addition to the electric matching (4).
For a perfectly absorbing interface this kind of feedback is absent.

With the parameters $U_{W}$ and $U_{I}$, the particle fluxes and hence the
source functions $S_{*,h}^{>}$ are fixed. Inserted in the boundary conditions
(14b) and (14c), the distribution functions for valence band holes and
conduction band electrons are obtained from Eqs. (24), with the substitutions
(26), from which follow also the densities $n_{*,h}(z)$ according to Eq. (11),
which in turn can be fed into the Poisson equation (3) to yield a new
potential $U(z)$ to be used in the next iteration step. For semiconductors
with high intrinsic densities, such as germanium, it turns out that the change
in densities caused by the injected charge carriers is negligible compared to
the Maxwellian background of intrinsic carriers. It is thus possible to
neglect it in the source term of the Poisson equation, reducing thereby the
calculational costs substantially. The densities plotted in Fig. 5 verify a
posteriori the validity of this simplification.

The iteration starts with constant Maxwellian distributions in accordance to
the boundary conditions (14), but for conduction band electrons continued to
energies which are not accessible at $z=z_{0}$ due to the higher value of
$U_{*}$. Successively, the distribution functions are updated by Eq. (24), as
are the parameters $U_{W}$ and $U_{I}$ by Eq. (4), expressed in terms of Eqs.
(28) and (29), and Eq. (33). In contrast to the perfectly absorbing interface
Rasek  _et al._ (2020), the plasma parameters now change slightly in each
iteration step, which in turn changes also the matching conditions at the
interface. The updating is repeated until convergence is reached. It should be
noted that the conduction band electron and valence band hole distribution
functions also enter the in-scattering collision integrals $\Phi_{s}$ in Eqs.
(24). Hence, even without the coupling to the plasma and changing boundary
conditions, an iteration is required to solve the Boltzmann equations inside
the germanium layer. This is also the case when the Grinberg-Luryi approach
Grinberg and Luryi (1992) is applied to semiconductor device modeling A. R.
St. Denis and D. L. Pulfrey (1998); Konistis and Hu (2002).

Figure 3: (color online) Current-voltage characteristic across a device
consisting of a germanium layer sandwiched between an argon plasma and an
Ohmic contact. The orange lines belong to the reflecting plasma-germanium
interface whereas the blue lines show data for an interface which perfectly
absorbs electrons from the plasma and lets them never return to it. Different
layer thicknesses are considered. For comparison, we also plot the
characteristic of a perfect absorber without a germanium layer, that is, for
$U_{SC}=0$ (solid black line, hardly seen in the main panel). The inset shows
for the perfectly absorbing interfaces the zero-crossings of the net flux. The
germanium layer shifts the crossings to lower voltages because of the
potential drops $U_{SC}$.

##  III Results

This section discusses the numerical results obtained for an argon plasma in
contact with a germanium layer. The material parameters are given in Table 1.
We split the discussion into three parts, depending on the perspective from
which the device shown in Fig. 1 is analyzed. First, regarding it as part of
an electric circuit, we present the current-voltage characteristic in
subsection III.1. Then, we proceed to discuss in subsection III.2 the spatial
profiles of the electric potential and the species’ densities as well as
fluxes. Finally, in subsection III.3, we turn to the energetically and
spatially resolved distribution functions of the charge carriers inside the
semiconducting layer.

###  III.1 Electric picture

Figure 4: (color online) As a function of $U_{B}$ and for the reflecting
plasma-solid interface, the electric field at the interface (top panel), the
potential drop $U_{SC}$ across the germanium layer (center panel), and the
electron and ion/hole fluxes flowing from the plasma towards the Ohmic contact
(bottom panel). Reflected and emitted electrons refer, respectively, to
electrons quantum-mechanically reflected at the potential step connected with
the plasma-solid interface and electrons emitted back to the plasma due to
collisions inside the semiconductor. While the electric field is identical for
the layer thicknesses considered, $U_{SC}$ changes with the thickness.
Electric field and $U_{SC}$ are effectively identical to what one obtains for
absorbing matching conditions. Fluxes are given only for the case of a
$1\,\mu$m thick germanium layer.

A biased plasma-solid interface can be considered as part of an electric
circuit and thus as an electric device characterized by a current-voltage
characteristic, that is, the net flux (current density) flowing through the
system as a function of the bias voltage $U_{B}$. For germanium layers of
different thicknesses in contact with an argon plasma, the characteristics are
shown in Fig. 3.

The germanium layers are terminated by an Ohmic contact, which is not further
characterized in the model. It merely acts as a sink for any flux carrying
particles reaching the semiconductor-metal interface. For $U_{SC}=0$, that is,
without the semiconducting layer, the system is thus a perfectly electron
absorbing Langmuir probe Mott-Smith and Langmuir (1926); Cherrington (1982);
Lam (1965). Current-voltage characteristics are given for different
thicknesses of the germanium layer, for both the electron reflecting and the
perfectly absorbing semiconductor-plasma interface. In the latter, the
matching conditions (14c) and (14d) are replaced by the assumption that every
electron impinging on the interface from the plasma is absorbed by the
germanium layer and every electron reaching the interface from the inside of
the layer is specularly reflected. The boundary conditions then take a form
similar to Eqs. (14a) and (14b), with the injection term (15b) normalized such
that the electron flux is conserved across the interface.

Within the perfect absorber model, there are hardly any differences in the
current-voltage characteristics of the Ohmic and the semiconducting plasma-
solid interface (black and blue lines in Fig. 3). The characteristics of the
latter are only shifted to lower voltages because of the voltage drops
$U_{SC}$ inside the semiconductor, as can be seen in the inset of the figure
for different film thicknesses. The shifts are a few tens of mV, with a
thicker layer giving rise to a larger shift. If the layers were thick enough
to host the whole negative space charge, the shifts would saturate, because of
the vanishing electric field deep inside the semiconductor. For the
thicknesses shown in Figs. 3 this is however not yet the case. Since, as
discussed in the next paragraph, the voltage drops $U_{SC}$ across reflecting
and absorbing germanium layers of the same thickness turn out to be
essentially identical, because $n_{e}^{>}$ and $\Delta_{e}^{>}$ are rather
small, and hence effectively negligible in Eq. (33), the large shift between
the blue and orange lines visible in the main panel is a consequence of the
reflected and emitted electron fluxes, which can be rather significant. From a
broader perspective, the results shown in Fig. 3 imply that modifications of
the surface of a Langmuir probe, for instance, by an oxide film, affects the
current-voltage characteristic only by a small amount if the perfect absorber
assumption holds. However, if this is not the case, the characteristic depends
on the film. In particular, the zero-crossing, that is, the floating
potential, would depend strongly on the emissive properties of the film.

To gain more insights, we plot in Fig. 4, as a function of the bias voltage
$U_{B}$, the electric field at the interface, the potential drop $U_{SC}$
across the semiconducting layer, and the electron and hole fluxes. The
electric field at the interface is independent of the thickness of the layers,
only the potential drop depends on it. The fluxes are shown representatively
only for a $1\,\mu\mathrm{m}$ thick semiconducting layer, with electron fluxes
split into the contributions arising from electrons reflected at the potential
step and electrons emitted from the semiconductor, that is, electrons which
made it into the layer but are backscattered by electron-phonon collisions
inside it. The two contributions arise, respectively, from the first and
second term of Eq. (14d).

While the reflected flux is independent of the thickness of the layer, since
it only depends on the reflection coefficient $R(E,T)$, the emitted flux
increases with the thickness of the layer, leading to smaller total net fluxes
towards the Ohmic contact, as can be also seen in Fig. 3. This may be an
artifact of the simple model we employ for the germanium layer. Flux-carrying
electrons, having typically large energy, as we will see in subsection III.3,
undergo in this model only a few collisions while passing through the layer.
The thicker the layer, the higher is thus the chance of an electron belonging
to this group for being backscattered. Hence, the emitted flux increases with
layer thickness.

The electron fluxes decay exponentially the more negative the bias voltage is,
since the fraction of the Maxwellian distribution of plasma electrons
contributing to the flux by overcoming the sheath potential $U_{I}$ decreases.
Only few electrons in its high energy tail contribute to the flux across the
interface and have thus a chance to get reflected or emitted. The emitted and
reflected electrons, encoded in $\Delta_{e}^{>}$ and $n_{e}^{>}$, modify thus
(33) only weakly and are hence negligible, as it is also the case in the
source term of the Poisson equation, where emitted and reflected electrons are
dominated by the plasma electrons belonging to the low energy part of the
Maxwellian.

###  III.2 Macroscopic picture

Besides regarding the interface as an electric device, characterized by a
current-voltage characteristic, it is also instructive to consider it as a
system of charged particles, characterized by density, flux, and potential
profiles. We limit the discussion in this and the next subsection to the
reflecting interface at the floating point, where the electron and ion/hole
fluxes are equal, noting, however, that for different bias voltages the
fundamental observations are similar, differing only by the numerical values.

Figure 5: (color online) Macroscopic properties for the reflecting
$1\,\mu\mbox{m}$ thick interface at the floating point, for which
$U_{I}=5.1\,$eV, $U_{SC}=0.025\,$eV, and $U_{W}=0.97\,$eV. Top to bottom:
Electric field, electric potential (shifted so that $U(z_{0})=0$), fluxes, and
densities. The last two are only given for the injected surplus electrons and
holes, the Maxwellian contributions due to the intrinsic carriers, giving rise
to the double layer shown in the inset of the top panel, are subtracted. In
addition to the net fluxes and densities per species, we also plot the
direction-resolved quantities.

Figure 5 shows from top to bottom the spatial profiles for the electric field,
the electric potential, the fluxes, and the densities of the injected carriers
for an interface with a $1\,\mu\mathrm{m}$ thick germanium layer. As it was
the case on the plasma side, the injected carriers, for which we show the
fluxes and densities, are also negligible in the source term of the Poisson
equation on the solid side of the interface. The solid-bound parts of the
electric field and potential shown in the top two panels are thus determined
by the intrinsic charge carriers. Hence, the solid-bound part of the electric
double layer plotted in the inset of the top panel, is the result of the
selfconsistent distribution of the intrinsic charge carriers in the electric
field arising due to the interface. For germanium, the intrinsic density at
room temperature is $2\cdot 10^{13}\,\textrm{cm}^{-3}$, and thus much larger
than the net density of the injected carriers, which, according to the bottom
panel of Fig. 5, is even at $z=0$ only only around $3\cdot
10^{11}\,\textrm{cm}^{-3}$, and hence about two orders of magnitude smaller.
The simplification we made for the solution of the Poisson equation is thus
valid. It not only stabilizes the iteration process, but allows also for an
efficient implementation, which reduces the computation time by about two
orders of magnitude.

The electric potential, shown in the second panel from the top of Fig. 5, is
shifted such that $U(z_{0})=0$. The potential drop across the germanium layer
$U_{SC}=25\,$meV can thus be read off directly at $z=0$. It can be also
immediately seen that the electric field is the derivative of the electric
potential. The fluxes, finally, are shown in the third panel from the top. In
addition to the net fluxes, $j_{*}$ and $j_{h}$, which are conserved and hence
independent of $z$, left- and right-moving fluxes are also plotted. They decay
exponentially with decreasing $z$. The net flux is the difference of the
right- and left-moving fluxes. It is negative, showing that more electrons and
holes are moving towards the Ohmic contact than in the other direction. The
conservation of the net flux is a consequence of the modeling of the germanium
layer, which contains only particle-conserving collisions with phonons. Had we
incorporated also particle non-conserving scattering processes, such as
radiative and non-radiative electron-hole recombination or impact ionization,
the net fluxes of each species would not be conserved. However, to have flux
conservation as an important indicator of the numerical accuracy of the
implementation of the selfconsistency scheme (it is not explicitly kept
constant in the iteration scheme), we neglected in this exploratory work
particle non-conserving scattering processes.

###  III.3 Microscopic picture

We now turn to the microscopic picture as it arises from the distribution
functions satisfying the two coupled sets of Boltzmann-Poisson equations. The
electric and macroscopic pictures contain this information only in an integral
manner. Now, we take an energy-resolved look at the kinetics of the charge
carriers across the interface. The dependence of the distribution functions on
the total energy $E$ and the lateral kinetic energy $T$ enables us to
visualize energy relaxation due to collisions with phonons as well as the
motion of the charge carriers parallel to the interface.

Let us first look at the energy-resolved net fluxes per species defined in Eq.
(21). For an interface with a $1\,\mu\mathrm{m}$ thick germanium layer at the
floating point, the fluxes flowing from the plasma towards the Ohmic contact
are plotted in Fig. 6. On the left (right), fluxes inside the germanium layer
(plasma) are shown. While the species’ energy-integrated net fluxes are
constant, the energy-resolved net fluxes display a rich behavior along the
germanium layer.

The net electron flux approaching the interface from the plasma side splits
inside the solid into a high- and a low-energy component. The former comprises
electrons remaining close to the energy where they have been initially
injected. Because the germanium layer is rather thin, a large number of
electrons reaching the Ohmic contact suffer only a few electron-phonon
collisions. Hence, they loose only a small amount of energy and remain
energetically high in the conduction band. Yet, some electrons scatter often
enough to end up at the band minimum. There is thus also a net flux at low
energy. In absolute numbers, this is even larger than the one at high energies
(note the different scales for high and low energies). However, its
contribution to the macroscopic net flux is rather small because the energy
range over which it is integrated is rather narrow. The main contribution to
the net total electron flux arises thus from high-energy electrons. In
numbers, we find $98\%$ of the total electron flux to be carried by the hot
electrons, with energies above $0.1\,\mathrm{eV}$, and only $2\%$ by the low
energy electrons.

Figure 6: (color online) In arbitrary units, the energy-resolved fluxes across
the interface with a $1\,\mu\mathrm{m}$ thick germanium layer. The interface
is at the floating point with potential drops given in the caption of Fig. 5.
On the left (right), the net electron and hole (electron and ion) fluxes
inside the germanium layer (plasma) are shown, with high and low energy
contributions plotted in separate panels. The ion flux, located at
$U_{I}+U_{W}+U_{SC}-\chi-E_{g}=1.39\,$eV, is shown by the red line (not
belonging to the color scale). Its energetic spread is negligible on the scale
of the other fluxes. The energy scale is chosen such that within the germanium
layer $E=0$ occurs at the minimum of $U_{*}$ and $U_{h}$, respectively.

The net flux of valence band holes, arising from the neutralization of ions,
hitting the interface with thermal kinetic energy, shows qualitatively the
same behavior. Holes, injected high up in the valence band at $E=I_{\rm inj}$,
relax due to collisions with phonons to the bottom of the band.
Quantitatively, however, the situation is different. Due to the larger
effective mass, holes loose energy due to collisions with phonons more
efficiently than electrons. The scattering rate $\gamma_{s}$, for instance,
defined in Eq. (8), is proportional to $m_{s}^{3/2}$. Hence, it increases with
the effective mass. The same holds for the in-scattering parts of the
collision integral $\Phi_{s}$. Thus, the larger the effective mass of the
charge carriers, the smaller is the inelastic mean free path, and hence the
spatial scale required to loose a substantial amount of energy. The
$1\,\mu\mathrm{m}$ thick germanium layer provides apparently enough space for
a significant number of holes to relax to the bottom of the band and to give
rise to a rather pronounced hole flux at low-energies. Integrated over energy,
the low-energy flux provides about $20\%$ to the total net flux. The high-
energy flux is hence still dominant but not as dominant as in the case of
electrons.

Energy relaxation of electrons and holes stops when their energies, measured
from the bottom of the conduction and valence band, respectively, are less
then the phonon energy. There is thus a threshold for energy relaxation,
leading to a kind of energy backlog in the fluxes up to a phonon energy above
the band minima. By absorbing phonons, carriers caught in the backlog can
reach higher energies. For holes, this is clearly visible in Fig. 5. It leads
to the feature in the bottom right panel around $E=0.08\,\mathrm{eV}$, which
is roughly twice the phonon energy. After the injected holes emitted several
hundred phonons, each carrying away $\hbar\omega_{0}=0.037\,$eV, their
energies fall below the threshold. Collisions are then less frequent, as the
emission (second) term in the scattering rate (8) disappears. The relatively
high population of the holes just below the phonon energy makes however the
absorption process operative in the collision integrals leading to the feature
in the hole flux around twice the phonon energy.

Figure 7: (color online) Distribution functions in arbitrary units for
electrons (top two rows) and holes (bottom two rows) at different spatial
locations of a $1\,\mu\mathrm{m}$ thick nearly floating germanium layer.
Regions of high and low energy are shown in separate panels. On the horizontal
axis, positive values of $|T_{z}|=E-U_{s}-T$ show $F^{>}_{s}$ and negative
values depict $F_{s}^{<}$. The triangular shape is due to the energy
restriction $0<T<E-U_{s}$. Note the logarithmic (linear) scale for low (high)
energies. At $z=-0.9\,\mu$m and $z=-0.5\,\mu$m, the high energy distributions
are multiplied by the factors given in the panels to utilize the same scale.
The potential drops are as in Figs. 5 and 6.

Finally, we discuss for the $1\,\mu\mathrm{m}$ thick germanium layer the
distribution functions $F_{*,h}^{\gtrless}(z,E,T)$. In addition to the effects
already present in the energy-resolved net fluxes, there are now also features
due to the lateral motion of the charge carriers. The distribution functions
are plotted in Fig. 7, with the two top (bottom) rows showing the data for
conduction band electrons (valence band holes). Three spatial locations are
considered and data for high and low energies are plotted in separate panels.
The lateral coordinate of the plots encodes by its sign also the direction of
motion, and is chosen such, that for $|T_{z}|=E-U_{s}-T=0$ it describes the
turning points. Note the linear scale for high and the logarithmic scale for
low energies. As for the macroscopic densities plotted in Fig. 5, the
Maxwellian background of intrinsic carriers is not included in the data. Only
the distribution functions for the surplus electrons and holes arising from
the plasma are shown.

In the high energy parts of the distribution functions right at the interface,
at $z=0$ (right column), the source functions $S_{*}^{<}$ and $S_{h}^{<}$,
describing the injection of electrons and holes, can be clearly identified.
Due to the collisions with phonons there is also a noticeable portion of
right-moving charge carriers, as can be seen from the rather large numerical
values for $F_{*,h}^{>}$. For electrons, this group of charge carriers
eventually leads to the flux of emitted electrons, if they also make it
through the potential barrier, and hence to secondary electron emission. As a
consequence of energy relaxation, Maxwellian distributions are established for
both species at low energies. Indeed, the low-energy parts of the distribution
functions are rather homogeneous in the lateral direction and exponentially
decaying with total energy, as it should be for Maxwellians. As a consequence
of the Maxwellian distribution at low energies, there is no low-energy flux
present at $z=0$, in accordance with the data shown in Fig. 6.

In the center of the germanium layer, at $z=-0.5\,\mu$m (center column), the
sharp features of the source functions are softened due to collisions with
phonons. A surplus motion to the left is also building up, most distinctly for
conduction band electrons, and weaker for valence band holes. At low energies,
the surplus of left-moving charge carriers remains, but it is now somewhat
more dominant for holes. Due to the imbalances, the distribution functions at
low energy start to deviate from Maxwellian distributions. As a result, low-
energy fluxes are building up, as can be also seen in Fig. 6. The behavior of
the energy-resolved net fluxes can thus be explained by the changes in the
distribution functions caused by the scattering processes.

The distribution functions at $z=-0.9\,\mu$m (left column), finally, are
already dominated by the boundary condition $F_{*,h}^{>}=0$, set for the
injected carriers at the interface to the Ohmic contact. For both species, the
motion is heavily biased to the left, irrespective of the energy. While right-
moving electrons and holes are strongly suppressed, there is a faint feature
in $F_{*}^{>}$ and $F_{h}^{>}$ slightly above the phonon energy which arises
from collisions with phonons. The features of the source functions are still
visible in the high-energy parts of $F_{*}^{<}$ and $F_{h}^{<}$, albeit
severely washed out, and at low energies $F_{*}^{<}$ and $F_{h}^{<}$ remain
Maxwellian. It should be noted that at high energies the absolute values of
the distribution functions decay with $z$ approaching $z_{0}$. The ratio,
however, of left to right moving distributions grows. Hence, the net flux at
high energies stays high, even though the individual distribution functions
$F_{*,h}^{\gtrless}$ for left- and right-moving carriers decay with decreasing
$z$, just like the direction-resolved fluxes $j_{*,h}^{\gtrless}$ shown in
Fig. 5.

##  IV Conclusion

We presented a kinetic description of ambipolar charge transport across a
biased plasma-solid interface consisting of a semiconducting germanium layer
sandwiched between an argon plasma and an Ohmic contact. The electron-hole
plasma within the semiconductor is coupled to the electron-ion plasma in front
of it through matching conditions for the distribution functions and the
electric field at the interface. Argon ions impinging on the germanium layer
create holes in the valence band, whereas electrons may be quantum-
mechanically reflected or transmitted. Electrons entering the semiconductor
from the plasma may be emitted back to it due to collisions with phonons
inside the solid, which also cause energy relaxation in its valence and
conduction band. To drive a current through the setup, a bias voltage is
applied between the bulk of the argon plasma, which we provide with a
prescribed plasma density, and the Ohmic contact used to collect the current.

From the distribution functions for the charge carriers on both sides of the
interface, we calculated the current-voltage characteristic. Due to the
quantum-mechanical reflection at the plasma-solid interface and the collisions
inside the solid, the characteristic differs significantly from the one
obtained for an interface which absorbs electrons from the plasma perfectly
and keeps them in the solid forever. Hence, the electron microphysics of the
semiconductor affects the electric properties of the interface and should thus
be considered by its theoretical description as well as its experimental
analysis.

We focused in this work on the implementation of a numerical scheme for the
selfconsistent calculation of the distribution functions and potential
profiles building up at the flux-carrying plasma-solid interface. For that
purpose, we kept the argon plasma collisionless and allowed electrons and
holes inside the germanium layer to scatter only on phonons. Particle
nonconserving collisions are not included. Within the simplified model, we
find flux-carrying conduction band electrons to remain at the high energies
set by the injection process. Due to the boundary and matching conditions, the
surplus electrons inside the semiconductor are thus rather hot and not in
thermal equilibrium with the intrinsic electrons. The same holds for injected
holes. With increasing layer thickness, we find secondary electron emission,
encoded in the emitted electron flux, to increase because a surplus electron
is then more likely to suffer collisions with phonons, which may bring it back
to the plasma, if it also successfully traverses the potential step at the
plasma-solid interface. The selfconsistent distribution functions enable us,
moreover, to visualize how the ambipolar gaseous charge transport in the argon
plasma merges with the solid-bound ambipolar charge transport inside the
germanium layer.

Albeit the description of the argon plasma is also somewhat crude, containing
an unspecified source and a collisionless sheath made consistent with the
germanium layer due to an equally unspecified presheath, the idealized
electronic structure of the germanium layer is more limiting. A quantitative
modeling of the biased plasma-solid interface has to be based on a realistic
band structure of the plasma-facing semiconductor. It should contain, over the
energy range set by the injection processes, the electronic structure of the
surface, entering the calculation of the electron reflectivity and the
modeling of the hole source function, as well as the electronic structure of
the bulk, which enters the collision integrals of the Boltzmann equations and
determines the velocity of the charge carriers. Of particular importance are
Bragg gaps preventing the transmission of electrons across the interface and
surface states trapping electrons and/or holes close to it. Since the
injection of electrons and holes occurs at rather high energies, impact
ionization, creating electron-hole pairs across the band gap, and its inverse,
the recombination of electron-hole pairs have to be moreover also included in
a modeling which attempts to provide more information than the insight that
current-carrying conduction band electrons and valence band holes are rather
hot.

## Acknowledgments

Support by the Deutsche Forschungsgemeinschaft through project BR-1994/3-1 is
greatly acknowledged.

## References

  * Franklin (1976) R. N. Franklin,  _Plasma Phenomena in Gas Discharges_ (Clarendon Press, Oxford, 1976).
  * Lieberman and Lichtenberg (2005) M. A. Lieberman and A. J. Lichtenberg,  _Principles of Plasma Discharges and Materials Processing_ (Wiley-Interscience, New York, 2005).
  * Langmuir and Mott-Smith (1924) I. Langmuir and H. Mott-Smith, Gen. Electr. Rev. 27, 449 (1924).
  * Nemschokmichal  _et al._ (2018) S. Nemschokmichal, R. Tschiersch, H. Höft, R. Wild, M. Bogaczyk, M. M. Becker, D. Loffhagen, L. Stollenwerk, M. Kettlitz, R. Brandenburg, and J. Meichsner, Eur. Phys. J. D 72, 89 (2018).
  * Peeters  _et al._ (2016) F. J. J. Peeters, R. F. Rumphorst, and M. C. M. van de Sanden, Plasma Sources Sci. Technol. 25, 03LT03 (2016).
  * Wagner  _et al._ (2005) H.-E. Wagner, Y. V. Yurgelenas, and R. Brandenburg, Plasma Phys. Control. Fusion 47, B641 (2005).
  * Massines  _et al._ (2003) F. Massines, P. Segur, N. Gherardi, C. Khamphan, and A. Ricard, Surface and Coatings Technology 174-175, 8 (2003).
  * Kogelschatz (2003) U. Kogelschatz, Plasma Chem. and Plasma Proc. 23, 1 (2003).
  * Chiang  _et al._ (2019) W.-H. Chiang, D. Mariotti, R. M. Sankaran, J. G. Eden, and K. Ostrikov, Adv. Materials 2019, 1905508 (2019).
  * Eden  _et al._ (2013) J. G. Eden, S.-J. Park, J. H. Cho, M. H. Kim, T. J. Houlahan, B. Li, E. S. Kim, T. L. Kim, S. K. Lee, K. S. Kim, J. K. Yoon, S. H. Sung, P. Sun, C. M. Herring, and C. J. Wagner, IEEE Trans. Plasma Science 41, 661 (2013).
  * P. A. Tchertchian, C. J. Wagner, T. J. Houlahan Jr., B. Li, D. J. Sievers, and J. G. Eden (2011) P. A. Tchertchian, C. J. Wagner, T. J. Houlahan Jr., B. Li, D. J. Sievers, and J. G. Eden, Contr. Plasma Phys. 51, 889 (2011).
  * Dussart  _et al._ (2010) R. Dussart, L. J. Overzet, P. Lefaucheux, T. Dufour, M. Kulsreshath, M. A. Mandra, T. Tillocher, O. Aubry, S. Dozias, P. Ranson, J. B. Lee, and M. Goeckner, Eur. Phys. J. D 60, 601 (2010).
  * Bronold and Fehske (2017) F. X. Bronold and H. Fehske, J. Phys. D: Appl. Phys. 50, 294003 (2017).
  * Schwager and Birdsall (1990) L. A. Schwager and C. K. Birdsall, Phys. Fluids B 2, 1057 (1990).
  * Riemann (1991) K. U. Riemann, J. Phys. D: Appl. Phys. 24, 493 (1991).
  * Franklin (2003) R. N. Franklin, J. Phys. D: Appl. Phys. 36, R309 (2003).
  * Brinkmann (2009) R. P. Brinkmann, J. Phys. D: Appl. Phys. 42, 194009 (2009).
  * Robertson (2013) S. Robertson, Plasma Phys. Control. Fusion 55, 093001 (2013).
  * Rasek  _et al._ (2020) K. Rasek, F. X. Bronold, and H. Fehske, Phys. Rev. E 102, 023206 (2020).
  * Mott-Smith and Langmuir (1926) H. M. Mott-Smith and I. Langmuir, Phys. Rev. 28, 727 (1926).
  * Cherrington (1982) B. E. Cherrington, Plasma Chem. and Plasma Proc. 2, 113 (1982).
  * Lam (1965) S. H. Lam, Phys. Fluids 8, 73 (1965).
  * Jacoboni and Reggiani (1983) C. Jacoboni and L. Reggiani, Rev. Mod. Phys. 55, 645 (1983).
  * Ridley (1999) B. K. Ridley,  _Quantum Processes in Semiconductors_ (Clarendon Press, Oxford, 1999).
  * Roth (1992) L. Roth, in  _Handbook on Semiconductors_ , edited by T. S. Moss (Elsevier Science Publishers B. V., Amsterdam, NL, 1992) p. 489.
  * Grinberg and Luryi (1992) A. A. Grinberg and S. Luryi, Solid-St. Electron. 35, 1299 (1992).
  * Moglestue (1993) C. Moglestue,  _Monte Carlo Simulation of Semiconductor Devices_ (Chapman & Hall, London, 1993).
  * MacColl (1939) L. A. MacColl, Phys. Rev. 56, 699 (1939).
  * A. R. St. Denis and D. L. Pulfrey (1998) A. R. St. Denis and D. L. Pulfrey, J. Appl. Phys. 84, 4959 (1998).
  * Konistis and Hu (2002) K. Konistis and Q. Hu, J. Appl. Phys. 91, 5400 (2002).


# Optimal angular window for observing Dark Matter annihilation from the
Galactic Center region: the case of gamma-ray lines.

Pasquale D. Serpico  Center for Particle Astrophysics, Fermi National
Accelerator Laboratory, Batavia, IL 60510-0500 USA  Gabrijela Zaharijas  HEP
Division, Argonne National Laboratory, 9700 Cass Ave., Argonne, IL 60439

###### Abstract

Although the emission of radiation from dark matter annihilation is expected
to be maximized at the Galactic Center, geometric factors and the presence of
point-like and diffuse backgrounds make the choice of the angular window size
to optimize the chance of a signal detection a non-trivial problem. We find
that the best strategy is to focus on an annulus around the Galactic Center of
$\sim 1^{\circ}$ to $\raise
1.29167pt\hbox{$\;>$\kern-7.5pt\raise-4.73611pt\hbox{$\sim\;$}}30^{\circ}$,
where the optimal size depends on the angular distribution of the signal and
the backgrounds. Although our conclusions are general, we illustrate this
point in the particular case of annihilation into two monochromatic photons in
the phenomenologically most interesting range of energy $45\,{\rm GeV}\raise
1.29167pt\hbox{$\;<$\kern-7.5pt\raise-4.73611pt\hbox{$\sim\;$}}E\raise
1.29167pt\hbox{$\;<$\kern-7.5pt\raise-4.73611pt\hbox{$\sim\;$}}80\,$GeV, which
is of great interest for the GLAST satellite. We find for example that Dark
Matter models with sufficiently strong line annihilation signals, like the
Inert Doublet Model, may be detectable without or with reasonable boost
factors.

###### pacs:

95.35.+d, 14.80.Cp, 98.70.Rz ANL-HEP-PR-08-16, FERMILAB-PUB-08-036-A

##  I Introduction

Although the existence of Dark Matter seems to be required by a wealth of
astrophysical and cosmological data (for reviews see Bertone:2004pz ;
Bergstrom:2000pn ), the identification of its physical nature remains still
elusive. The most popular candidates are the so-called weakly interacting
massive particles (or WIMPs), believed to be neutral and colorless thermal
relics from the early universe with electroweak-scale masses ($m_{X}\sim
10-1000\,$GeV) and couplings. These generic assumptions seem indeed to roughly
reproduce the amount of DM needed in the Universe and in our Galaxy in
particular, while arising at the same time naturally in many extensions of the
Standard Model of Particle Physics.

Among the different methods proposed to identify WIMP dark matter (production
in accelerators, direct detection via recoil of nuclei underground, indirect
detection via their annihilation or decay products Hooper:2008sn ), a
particularly promising technique is the detection of gamma rays from
annihilation in the galactic halo. However, it is nowadays clear that the DM
signal should be likely singled out of a larger astrophysical background,
either point-like (the case of the Galactic Center) or diffuse. In assessing
the capability of detecting a given WIMP candidate it is important to optimize
the search strategy based on our current expectation for the energy range, but
also the angular distribution of the signal. This information proves actually
to be at least as useful in guiding DM searches as the energy spectrum and the
choice of astrophysical targets. Here we illustrate this point with reference
to the inner Galaxy region. Although the emission of radiation per unit solid
angle from dark matter annihilation is expected to be maximized at the
Galactic Center, geometric factors and the presence of point-like and diffuse
backgrounds make the choice of the window size to optimize the chance of a
signal detection a non-trivial problem. To the best of our knowledge, this
issue has been only raised within over-simplified assumptions for the signal
and background Stoehr:2003hf  or, by converse, treated in some technical
studies by experimental groups devoted to dark matter searches (see for
example poster ), with the background profile determined e.g. via a numerical
model. In no case a comparative study for different profiles and background
assumptions, which illustrates the physical origin of the optimal cuts to
adopt, has been presented. This is the main task we address in this article.
We shall prove that—in all the commonly considered models for the DM
distribution in our Galaxy and given the present knowledge of point-like and
diffuse backgrounds—the optimal strategy for DM searches is never to focus on
the inner 0.1∘-0.5∘ around the Galactic Center. Unless astrophysical processes
dramatically enhance the DM distribution in the innermost region of our
Galaxy, producing effectively a bright point-source, an annulus of radius from
$\sim 1^{\circ}$ to $\raise
1.29167pt\hbox{$\;>$\kern-7.5pt\raise-4.73611pt\hbox{$\sim\;$}}30^{\circ}$ is
always to be preferred. The optimal size depends on the angular distribution
of the signal and the backgrounds, and the improvement in the sensitivity
reaches a factor of a few, thus affecting the prospects for the detectability
of specific models. To give an explicit example, we shall refer mostly to the
detection of gamma-ray lines from annihilation. Gamma-ray lines of energy
$E=m_{X}$ would be an unambiguous signature of WIMPs as DM particles
populating our Galactic Halo, but the two-photon channel is loop(or
radiatively)-suppressed with respect to the main WIMP annihilation channels:
fermion pairs, weak gauge bosons or higgses. Even models where the line
emission is enhanced, as the Inert Doublet Model (IDM) Barbieri:2006dq ;
Gustafsson:2007pc  are very challenging to detect via line emission. Yet, with
an optimized choice of the angular size, the IDM model may be detectable
without or with reasonable boost factors (i.e., enhancement of the flux with
respect to the smooth halo signal due to the “clumpiness” of DM). Similar
considerations apply for the more conventional continuum gamma-ray flux from
DM annihilation, usually considered a more promising channel.

This paper is structured as follows. In Sec. II we discuss the ingredients
contributing to the DM signal. In Sec. III we describe our model for the
point-like and diffuse backgrounds, while Sec. IV is devoted to present our
results. In Sec. V we conclude.

##  II The Signal

It is customary to write the number flux of gamma rays (per energy and solid
angle) produced in dark matter (DM) annihilations in our Galactic halo as

$\Phi^{\rm DM}(E,\psi)=\left[\frac{{\rm d}N_{\gamma}}{{\rm
d}E}\frac{\langle\sigma v\rangle}{8\pi m^{2}_{X}}\right]\,d_{\rm f}\,\rho_{\rm
f}^{2}\,J(\psi)\,,$ (1)

where the factor in square brackets depends only on the particle physics
model: $\langle\sigma v\rangle$ is the WIMP annihilation cross section
multiplied by the relative velocity of the two WIMPs (averaged over the WIMP
velocity distribution), $m_{X}$ is the mass of the WIMP, ${\rm
d}N_{\gamma}/{\rm d}E$ is the photon spectrum per annihilation. If the dark
matter is not its own antiparticle (as assumed here), Eq. (1) should be
multiplied further by a factor 1/2 (if $X$ and $\bar{X}$ are equally
abundant).

The factor outside the brackets only depends on the DM astrophysics, in
particular on the distribution of DM in the halo and on the position of the
Solar System within the halo. The quantities $d_{\rm f}$ and $\rho_{\rm f}$
are constants with dimensions of distance and energy density respectively,
whose only purpose is of making $J$—which encodes the angular dependence of
the DM signal—a dimensionless function. Following Bergstrom:1997fj , we fix
these parameters at $d_{\rm f}=8.5\,$kpc and $\rho_{\rm f}=0.3\,$GeV/cm3
(indicative the galactocentric distance of the sun and of the local density of
dark matter). The $J$-function is thus defined in terms of the DM density
$\rho$ as

$J(\psi)\equiv\frac{1}{d_{\rm f}}\frac{1}{\rho_{\rm
f}^{2}}\int_{\rm{los}}\rho^{2}[r(\ell,\psi)]{\rm d}\ell,$ (2)

where

$r(\ell,\psi)=\sqrt{r_{\odot}^{2}+\ell^{2}-2\,r_{\odot}\,\ell\cos\psi}\,.$ (3)

In the previous equations, $\psi$ is the angle between the direction in the
sky and the Galactic Center so that, in terms of galactic latitude $b$ and
longitude $l$, $\cos\psi=\cos b\cos l$; $r_{\odot}\approx 8.0\,$kpc Yao:2006px
is the solar distance from the Galactic Center; $\ell$ is the distance from
the Sun along the line-of-sight (l.o.s.). Note that we implicitly assumed that
the halo is spherically symmetric and that the angular dependence only comes
from the off-center position of the Sun in the halo, which implies for example
the azimuthal symmetry of the signal (see e.g. Hooper:2007be  and refs.
therein for a discussion of sub-leading effects determining the angular
distribution of the signal). A general class of spherically symmetric, smooth
halo distributions can be fitted as

$\rho(r)=\left(\frac{r_{s}}{r}\right)^{\gamma}\frac{\rho_{0}}{[1+(r/r_{s})^{\alpha}]^{(\beta-\gamma)/\alpha}},$
(4)

where $\rho_{0}$ is a normalization constant and $r_{s}$ is a characteristic
radius below which the profile scales as $r^{-\gamma}$. Two of the most well
known profiles have been proposed by Navarro, Frenk and White (NFW) NFW  and
Moore et al. Moore . Another fit, proposed by Kravtsov et al. Kravtsov:1997dp
mimics a flat core potential, which seems to be in better agreement with some
observations. More recent simulations (see e.g. Power:2002sw ; Navarro:2003ew
; Reed:2003hp ; Merritt:2005xc ) suggest that halo density profiles are better
represented by a function with a continuously-varying slope. At least above
the resolved radius ($\sim 1\,$kpc) and for Galactic-size halos, the newly
proposed fitting formulae provide only marginal improvement with respect to
the more traditional ones, although it seems that the inner slope of a DM halo
is very likely shallower that the Moore one, and probably shallower than the
NFW one, too. For illustrative purposes, we shall consider the Moore, NFW and
Kravtsov profiles with the parameters reported in Table I, since the former
and the latter bracket the different cases, with the NFW (as well as the
paramterization used in Merritt:2005xc ) falling in between. Similarly to Ref.
Yuksel:2007ac , the normalizations are chosen so that the mass contained
within the solar circle provides the appropriate DM contribution to the local
rotational curves, such that $\rho_{\odot}\equiv\rho(r_{\odot})$—the DM
density at the solar distance from the GC—yields to the canonical value
$\rho_{\odot}=0.3\,$GeV/cm3 for the NFW profile.

Table 1: Halo parameters, see text and Ref. Yuksel:2007ac . Model | $\alpha$ | $\beta$ | $\gamma$ | $\rho_{\odot}$ | $r_{s}$
---|---|---|---|---|---
|  |  |  | [GeV cm-3] | [kpc]
Moore | 1.5 | 3 | 1.5 | 0.27 | 28
NFW | 1.0 | 3 | 1.0 | 0.30 | 20
Kravtsov | 2.0 | 3 | 0.4 | 0.47 | 10

Note that in none of the above-mentioned models astrophysical processes
modifying the pure DM profile are included. For example, “adiabatic
compression” may lead to the formation of a DM spike around the central
supermassive black hole; mergers as well as scattering on the dense stellar
cusp around the central black hole may destroy density enhancements, etc.
While potentially relevant in determining the dark matter annihilation
brightness of the Galactic Center, they only affect regions too close to the
Galactic Center to be resolved angularly by present detectors. In general,
they might lead to a measurable gamma ray flux from the innermost angular bin
even in presence of relatively large astrophysical backgrounds, but in this
case only the energy information is relevant for the separation, for which the
present discussion is irrelevant. We not address these effects further;
details can be found in the refs. cited in Bertone:2004pz  or ref.
Fornasa:2007nr .

The spectrum per annihilation depends on the nature of the WIMP, although the
gamma-ray yield is often largely dominated by the result of the fragmentation
and hadronization of the WIMP annihilation products. Neutralinos in the MSSM,
for example, typically annihilate to final states consisting of heavy fermions
and gauge or Higgs bosons jungman . With the exception of the
$\tau^{+}\tau^{-}$ channel, each of these annihilation modes result in a very
similar spectrum of gamma rays. In most of the following considerations,
however, it is the angular distribution of the DM radiation that determine the
optimal angular window. To avoid unnecessary complications, we shall
illustrate our point considering mostly the case of annihilation in two
gammas, for which

$\frac{{\rm d}N_{\gamma}}{{\rm d}E}=2\,b_{\gamma\gamma}\delta(E-m_{X})\,,$ (5)

where $b_{\gamma\gamma}$ is the branching ratio for the 2-photon line.
Although it would constitute a striking signature of a DM annihilation
process, in general the photon line is considered hard to detect, since it is
loop-suppressed with respect to the main annihilation modes Bergstrom:1997fj .
This is not necessarily the case in some DM models, as for example the Inert
Doublet Model (IDM) Barbieri:2006dq ; Gustafsson:2007pc . In general, the line
feature may be enhanced whenever: (i) the DM almost decouples from fermions,
while being mostly coupled to gauge bosons; (ii) the DM candidate is in the
mass range $45\,{\rm GeV}\simeq m_{Z}/2<M_{X}<m_{W}\simeq{\rm 80}\,$GeV, the
lower limit coming from the bound on the $Z^{0}$-width measurement and the
upper limit to avoid the annihilation mode $X\,X\to W^{+}W^{-}$ by kinematical
constraints (other cases where a prominent line emission is possible at lower
or higher energies have been treated respectively in Pullen:2006sy  and
Bringmann:2007nk ). Interestingly, this is probably the best region where the
GLAST satellite GLASTurl  has a reasonable chance to detect a line emission
and at the same time, potentially interesting for the next generation of low-
threshold atmospheric cherenkov telescopes (ACTs), like CTA CTA  or AGIS AGIS
. In the following, for GLAST we shall use an integrated area exposure over
time of $10^{4}\,{\rm cm}^{2}\,{\rm yr}$, corresponding to about 5 years of
observations, and an energy resolution of $\pm 7\%$ at $1\sigma$. For ACTs, we
shall assume an effective array area of $A_{\rm eff}=1\,$km2 , energy
resolution of $15\%$ and 200 hour observations, with a $2^{\circ}\times
2^{\circ}$ field of view.

Whenever needed, we shall compare our results with two benchmark models for
IDM candidates among the ones considered in Gustafsson:2007pc , see Table 2.
Also, to get a feeling of how things change if the continuum gamma-ray
spectrum is considered instead, we shall also report the results obtained when
a fiducial case of a 100 GeV MSSM neutralino mainly annihilating into
$W^{+}W^{-}$ is considered.

Table 2: Benchmark models in the IDM scenario, see Gustafsson:2007pc . Model | $m_{X}$ | $\langle\sigma v\rangle$ | $b_{\gamma\gamma}$
---|---|---|---
| [GeV] | [cm3/s] |
$I$ | $70$ | $1.6\times 10^{-28}$ | $0.36$
$II$ | $50$ | $8.2\times 10^{-29}$ | $0.29$

##  III The Backgrounds

The Galactic Center is a complex region of the sky at all wavelengths, the
gamma-ray window being no exception. In this section, we discuss how, in our
analysis, we treat the backgrounds for dark matter searches due to known and
unknown astrophysical sources of gamma rays. Since we focus on an extended
signal, diffuse (or unresolved) backgrounds are our main concern. However, to
illustrate how effective is a large field-of-view search, we shall compare it
with the sensitivity achieved when focusing on the inner region around the
Galactic Center, when at least one astrophysical point-like background is
known to exist. This is the relatively bright, very high-energy gamma ray
source observed by HESS, MAGIC, WHIPPLE, and CANGAROO-II GCtev . This source
is consistent with point-like emission and is located at
$l=359^{\circ}56^{\prime}41.1^{\prime\prime}\pm 6.4^{\prime\prime}$ (stat.),
$b=-0^{\circ}2^{\prime}39.2^{\prime\prime}\pm 5.9^{\prime\prime}$ (stat.) with
a systematic pointing error of 28′′ van Eldik:2007yi . It appears to be
coincident with the position of Sgr A⋆, the black hole constituting the
dynamical center of the Milky Way. Following the measurements of HESS, we
describe the spectrum of this source as a power-law given by:

$\Phi^{\rm GC}=1.0\times 10^{-8}\left({{E}\over{{\rm GeV}}}\right)^{-2.25}{\rm
GeV}^{-1}\,{\rm cm}^{-2}\,{\rm s}^{-1}.$ (6)

At energies below $\sim 200$ GeV, the spectrum of this source has not yet been
measured. GLAST, however, will be capable of measuring the spectrum of this
source at energies below the thresholds of HESS and other ACTs.

In addition to the HESS source, a yet unidentified source has been detected by
EGRET in the GeV range, approximately 0.2∘ away from the dynamical center of
our galaxy dingus ; pohl . It has not been detected so-far at the energies
explored by HESS and other ACTs, which means that its power-law spectrum cuts-
off between $\sim$10 and $\sim$100 GeV. Although it should be included in more
realistic analyses of the DM signal from the inner galactic regions
Dodelson:2007gd , in the following we neglect it, thus overestimating the
diagnostic power of a narrow field of view search and providing a conservative
estimate of the improvement achieved by a large angle search.

The overall diffuse gamma-ray radiation can be qualitatively divided into a
galactic and an extragalactic contribution. Since the latter is not simply the
isotropic part of the flux, the separation of these two components can be done
at present only assuming a specific model for the production of secondaries by
cosmic rays in the galactic disk and halo. (However, a measurement of the
cosmological Compton-Getting effect that should be achievable for GLAST would
provide a model-independent way to separate the two contributions
Kachelriess:2006aq ). A significant fraction of the quasi-isotropic component,
especially in the GeV range, may be due to high-latitude galactic emission
coming from processes in the magnetized halo of the Milky Way. We employ here
a fit of the galactic diffuse flux proposed in Bergstrom:1997fj  and
calibrated on EGRET data around the GeV Hunter:1997 ,

$I_{\rm gal}(E)=N_{0}(l,b)\times 10^{-6}\left(\frac{E}{{\rm
GeV}}\right)^{-\alpha}\,{\rm cm}^{-2}{\rm s}^{-1}{\rm sr}^{-1}{\rm
GeV}^{-1}\,,$ (7)

where the arguments $l,b$ are in degrees, $-180^{\circ}\leq l\leq 180^{\circ}$
and $-90^{\circ}\leq b\leq 90^{\circ}$, $\alpha\simeq 2.7$, and

$N_{0}(l,b)=\left\\{\begin{array}[]{cc}\frac{85.5}{\sqrt{1+(l/35)^{2}}\sqrt{1+[b/(1.1+0.022\,|l|)]^{2}}}+0.5\,,&\>\>|l|\geq
30^{\circ}\\\
\frac{85.5}{\sqrt{1+(l/35)^{2}}\sqrt{1+(b/1.8)^{2}}}+0.5\,,&\>\>|l|\leq
30^{\circ}\end{array}\right..$ (8)

The EGRET collaboration derived the intensity of the extragalactic gamma-ray
flux as Sreekumar:1997un

$I_{\rm ex}(E)=k_{0}\times 10^{-6}\left(\frac{E}{0.451{\rm
GeV}}\right)^{-\beta}{\rm cm}^{-2}{\rm s}^{-1}{\rm sr}^{-1}{\rm GeV}^{-1}\,,$
(9)

with $\beta=2.10\pm 0.03$, $k_{0}=(7.32\pm 0.34)$ and the fits is valid from
$E\sim\,$10 MeV to $E\sim\,$100 GeV.

In the following, we shall consider two models of diffuse backgrounds: a
“conservative” model assuming $\alpha=2.7$ in Eq. (7) and $k_{0}=7.32$ in Eq.
(9), and an “optimistic” model, where a harder spectrum $\alpha=3.0$ is
assumed for the Galactic background and $k_{0}=3.66$. Based on presently
available models, one expects indeed that GLAST should resolve roughly half of
the extragalactic background Stecker:2001dk ; analogously, one expects that
unresolved sources contribute significantly to the diffuse Galactic spectrum
as measured by EGRET, and once detected and removed by GLAST the Galactic
background should be steeper (see e.g. discussion in Pullen:2006sy ).

Finally, for ACTs the dominant background is the isotropic cosmic ray one. At
energies well below $\sim 100\,$GeV one expects the cosmic-ray electron
background to play a dominant role. In the following, we account for these
backgrounds following the formulae reported in Bergstrom:1997fj .

##  IV Results

First, let us repeat a simple argument already sketched in Stoehr:2003hf . The
DM signal coming from integrating over a solid angle $\Delta\Omega$ is
proportional to $\langle J(\Omega)\rangle\Delta\Omega$, where

$\langle
J(\Omega)\rangle\equiv\frac{1}{\Delta\Omega}\int_{\Delta\Omega}J\,{\rm
d}\Omega\,.$ (10)

For an annulus around the GC with inner opening semi-angle $\psi_{\rm m}$ and
outer one $\psi_{\rm M}$, the previous formula reduces to

$\langle J\rangle=\frac{1}{\cos\psi_{\rm m}-\cos\psi_{\rm M}}\int_{\psi_{\rm
m}}^{\psi_{\rm M}}J(\psi)\sin\psi{\rm d}\psi\,.$ (11)

The significance $S/N$ of a Signal $S$ in presence of a background $B$ is
defined as

$\frac{S}{N}\equiv\frac{S}{\sqrt{S+B}}$ (12)

The dark matter signal goes like $S\propto\langle J\rangle\,\Delta\Omega$,
while for an isotropic background $B\propto\Delta\Omega$. In the limit $S\ll
B$ (which is likely to hold in all realistic DM searches), we have

$\frac{S}{N}\propto\langle J\rangle\sqrt{\Delta\Omega}\sim\langle
J\rangle\psi_{\rm M}\,.$ (13)

In Fig. 1, we plot this function vs. $\psi_{\rm M}$ for the three models
considered in Table II, assuming $\psi_{\rm m}=0.1^{\circ}$ (solid lines) or
$\psi_{\rm m}=0.3^{\circ}$ (dashed lines). Note the exact choice of the inner
bound is irrelevant for cored profiles: the peak of this function is above
20∘. Even for profiles as cusped as the NFW (or more) the peak of this
function is at a few degrees. Of course, the Moore profile gives a formally
divergent signal when $\psi_{\rm m}\to 0$, so the inner angle cutoff does
matter.

Figure 1:  Plot of the function $\langle J\rangle\psi_{\rm M}$ defined in Eq.
(11) (proportional to the significance of a DM signal for an isotropic
background in the limit $S\ll B$) for the three models considered in Table I:
the solid lines assume $\psi_{\rm m}=0.1^{\circ}$, the dashed lines $\psi_{\rm
M}=0.3^{\circ}$; the leftmost (red) lines are for the Moore profile, the
rightmost (green) lines for the Kravtsov, in the middle (blue) lines for the
NFW profile.

Figure 1 already shows that the $S/N$ is typically maximized for values of
$\psi_{\rm M}\gg\psi_{\rm m}$. However, the diffuse backgrounds are not
$\\{l,b\\}$-independent, but really decline away from the Galactic Plane. When
using the models for the galactic and extragalactic background introduced in
the previous section, we see how this effect is appreciable. In Fig. 2 , in
order to make the comparison among the different models easier, we plot the
rescaled $S/N$ vs. $b_{\rm max}$ (in such a way that their maximum values are
equal) for a region $0.4^{\circ}<|b|<b_{\rm max}$, $0^{\circ}<|l|<l_{\rm
max}=b_{\rm max}$ for the fiducial IDM model I of Tab. II, and for the
continuum spectrum of a 100 GeV MSSM neutralino annihilating into
$W^{+}W^{-}$. In the former case, we consider the 1-$\sigma$ smoothing of the
gamma-line due to the GLAST energy resolution, while for the continuum we
actually show the S/N in the most significant energy bin (which happens to
fall in the $\sim 1\,$GeV range). The inner cut has been chosen to
$0.4^{\circ}$ to avoid contamination from known point sources close to the GC.
Note how the peak of $S/N$ moves to larger angles in all cases compared with
Fig. 1, developing a maximum at finite angle also for the Moore profile. For
the continuum spectrum the peak of the significance moves to larger angles,
due to the fact that the isotropic component of the background is harder, and
thus at the lower energies relevant for this signal the overall background is
dominated by the Galactic background which falls rapidly with $b$. It is also
worth nothing that even for heavy dark matter particles (TeV mass scale) the
continuum photon spectrum peaks at or below about 30 GeV. Since the $S/N$ for
continuum spectra is dominated by its value at the peak, the angular window
results reported in Figs. 1 and 2 are really more general than the case they
were discussed for, reflecting the situation for the continuum spectrum from
DM in all the interesting mass range.

Figure 2:  relative Signal/Noise vs. annulus size $b_{\rm max}$ (we equate
their maximum value to ease the comparison) for a region
$0.4^{\circ}<|b|<b_{\rm max}$, $0^{\circ}<|l|<l_{\rm max}=b_{\rm max}$ for the
continuum spectrum of a 100 GeV MSSM neutralino annihilating into $W^{+}W^{-}$
(dot-dashed line, NFW profile) and for the fiducial IDM model I of Tab. II for
the Moore, NFW and Kravtsov profile (solid lines respectively from the left to
the right).

In Fig. 3 we show the curves with $S/N=3$ for the branching ratio into two
gammas, normalized in terms of the fiducial annihilation cross section value
$\langle\sigma v\rangle=3\times 10^{-26}\,$cm3/s in the phenomenologically
most interesting mass range below $m_{W}$. The triangles represent the
fiducial points of the IDM model reported in Table II, and for the window size
that maximizes $S/N$ (respectively $b_{\rm max}\sim 1.5^{\circ}$, 15∘ and 30∘,
for Moore, NFW and Kravtsov profiles of Tab. I). The solid lines are for the
conservative background estimate, the dashed lines for the optimistic one. It
is interesting to note that for the Moore profile, the fiducial models should
be detectable at the $3\,\sigma$ level without need for boost factors. For the
NFW case, boost factors of order 10 should suffice. Even for the unfortunate
case of a Kravtsov profile, a boost factor $\sim 20$ should be sufficient.
These values appear to be viable in cold DM cosmologies, see e.g.
Strigari:2006rd . Thus, our conclusions are more promising that what
envisioned in Gustafsson:2007pc , where boost factors $\raise
1.29167pt\hbox{$\;>$\kern-7.5pt\raise-4.73611pt\hbox{$\sim\;$}}100$ were
considered. Also, the top dot-dashed line shows the sensitivity achieved if
one only includes the inner $0.1^{\circ}$ region around the GC, where the
background is dominated by the central source. These limits, comparable with
the results reported in gabi , show an evident improvement.

It is also worth commenting on the fact that, according to our estimate,
although ACTs are not optimized to search for line signals in the range below
100 GeV, the performance of the next generation, km2 instruments presently
being considered may be comparable or better than GLAST. For example, in
accordance with existing estimates ACTwp , for the fiducial exposure
parameters we have used and the NFW profile we find a sensitivity a factor
$\sim$1.5 better in the parameter space of Fig. 3, despite their worse energy
resolution. It is also fortunate that the ground-based instruments, that have
necessarily a much smaller field of view with respect to satellite
observatories like GLAST, do not suffer too much from limited field of view,
at least for relatively cusped profiles. The dominant backgrounds of ACTs are
indeed isotropically distributed, and for this case the optimal annulus size
is $\lesssim 2^{\circ}\div 3^{\circ}$ (for cusped profiles), within the field
of view range considered by present designs for future ACT.

Figure 3:  Plot of the curves with $S/N=3$ in the $10-100$ GeV range, with
vertical lines in correspondence of $m_{Z}/2$ and $m_{W}$. The dot-dashed line
shows the sensitivity achieved if one only includes the inner $0.1^{\circ}$
region around the GC for the NFW profile. The couples of solid/dashed lines
(solid for conservative background, dashed for optimistic one) represent
instead the sensitivity for the optimal annulus size for each model: from top
to bottom, for Kravtsov, NFW and Moore profiles, respectively.

##  V Discussion and Conclusions

A large field of view instrument like GLAST promises to revolutionize the
field of high energy gamma ray astrophysics. But it could also turn into a
discovery instrument for DM, if it is made by WIMPs which annihilate into
standard model particles in the halo of our galaxy. Since the DM signal grows
like the square of its density, it is usually assumed that the Galactic
Center, where the DM density is expected to be the highest, is the most
promising region to focus on. However, two arguments conspire to modify this
naive expectation: i) the presence of astrophysical point-like backgrounds
close to the GC; ii) geometric considerations, in particular related to the
different angular shape of the DM signal and astrophysical backgrounds.

It is thus important to establish the level of the diffuse galactic and
extragalactic diffuse backgrounds, a task which will be accomplished with
exquisite precision by GLAST itself. Yet, by employing simple models based on
EGRET data, we showed that it is always to be preferred a window size
$>1^{\circ}$ around the GC, even for profiles as steep as the Moore one. The
flatter the profile, the larger the optimal window size: for example we find
an optimal galactic latitude cut $b_{\rm max}\raise
1.29167pt\hbox{$\;>$\kern-7.5pt\raise-4.73611pt\hbox{$\sim\;$}}10^{\circ}$ for
NFW, $b_{\rm max}\raise
1.29167pt\hbox{$\;>$\kern-7.5pt\raise-4.73611pt\hbox{$\sim\;$}}30^{\circ}$ for
Kravtsov profile. The discovery of additional point-like sources is not
expected to change things dramatically, since especially at high energies the
good angular resolution of GLAST will allow one to remove the quasi-totality
of the resolved point-source contaminations with only a minor loss in the
solid angle of the diffuse signal.

When these geometric effects are taken into account and the search strategy is
adjusted consistently, the prospects for DM detection improve considerably.
For example, we found that models with sufficiently strong line annihilation
signals, like the Inert Doublet Model, may be detectable without or with still
reasonable boost factors, of the order of 10 even for Kravtsov-like profiles.
Also, the perspectives to detect a DM signal become less sensitive to the
unknown profile of the DM halo. In particular the $S/N$ of quasi-cored
profiles is only a factor $\sim 2$ worse than for the NFW case. We also showed
with one example that our considerations remain true also for the widely
considered case of continuum gamma-ray spectrum from DM annihilation. Given
the generality of these considerations, we believe they should be taken into
account when estimating the reach of DM searches via gamma-ray telescopes.

Acknowledgments. Authors are grateful to James Buckley and Dan Hooper for
careful reading of the manuscript and useful comments. This work is supported
in part by the US Department of Energy, Division of High Energy Physics, under
Contract DE-AC02-06CH11357. P.S. is supported by the US Department of Energy
and by NASA grant NAG5-10842. Fermilab is operated by Fermi Research Alliance,
LLC under Contract No. DE-AC02-07CH11359 with the United States Department of
Energy.

## References

  * (1) G. Bertone, D. Hooper and J. Silk, Phys. Rept. 405, 279 (2005) [hep-ph/0404175].
  * (2) L. Bergstrom, Rept. Prog. Phys. 63, 793 (2000) [hep-ph/0002126].
  * (3) D. Hooper and E. A. Baltz, arXiv:0802.0702 [hep-ph].
  * (4) F. Stoehr, S. D. M. White, V. Springel, G. Tormen and N. Yoshida, Mon. Not. Roy. Astron. Soc. 345, 1313 (2003) [astro-ph/0307026].
  * (5) http://confluence.slac.stanford.edu/download/attachments/19303/HEAD_Poster.pdf?version=1
  * (6) R. Barbieri, L. J. Hall and V. S. Rychkov, Phys. Rev. D 74, 015007 (2006) [hep-ph/0603188].
  * (7) M. Gustafsson, E. Lundstrom, L. Bergstrom and J. Edsjo, Phys. Rev. Lett. 99, 041301 (2007) [astro-ph/0703512].
  * (8) L. Bergstrom, P. Ullio and J. H. Buckley, Astropart. Phys. 9, 137 (1998) [astro-ph/9712318].
  * (9) W. M. Yao et al. [Particle Data Group], J. Phys. G 33, 1 (2006).
  * (10) D. Hooper and P. D. Serpico, JCAP 0706, 013 (2007) [astro-ph/0702328].
  * (11) J. F. Navarro et al., Astrophys. J. 462 563 (1996).
  * (12) B. Moore et al., Phys. Rev. D 64 063508 (2001).
  * (13) A. V. Kravtsov, A. A. Klypin, J. S. Bullock and J. R. Primack, Astrophys. J. 502, 48 (1998) [astro-ph/9708176].
  * (14) C. Power et al., Mon. Not. Roy. Astron. Soc. 338, 14 (2003) [astro-ph/0201544].
  * (15) D. Reed et al., Mon. Not. Roy. Astron. Soc. 357, 82 (2005) [astro-ph/0312544].
  * (16) J. F. Navarro et al., Mon. Not. Roy. Astron. Soc. 349, 1039 (2004) [astro-ph/0311231].
  * (17) D. Merritt, J. F. Navarro, A. Ludlow and A. Jenkins, Astrophys. J. 624, L85 (2005) [astro-ph/0502515].
  * (18) H. Yuksel, S. Horiuchi, J. F. Beacom and S. Ando, Phys. Rev. D 76, 123506 (2007) [arXiv:0707.0196 [astro-ph]].
  * (19) M. Fornasa and G. Bertone, arXiv:0711.3148 [astro-ph].
  * (20) G. Jungman, M. Kamionkowski and K. Griest, Phys. Rept. 267, 195 (1996).
  * (21) A. R. Pullen, R. R. Chary and M. Kamionkowski, Phys. Rev. D 76, 063006 (2007) [astro-ph/0610295].
  * (22) T. Bringmann, L. Bergstrom and J. Edsjo, arXiv:0710.3169 [hep-ph].
  * (23) http://www-glast.slac.stanford.edu/
  * (24) http://www.mpi-hd.mpg.de/hfm/CTA/
  * (25) http://cherenkov.physics.iastate.edu/wp/
  * (26) F. Aharonian et al. [The HESS Collaboration], astro-ph/0408145; J. Albert et al. [MAGIC Collaboration], Astrophys. J. 638, L101 (2006) [astro-ph/0512469]; K. Kosack et al. [The VERITAS Collaboration], Astrophys. J. 608, L97 (2004) [astro-ph/0403422]; K. Tsuchiya et al. [CANGAROO-II Collaboration], Astrophys. J. 606, L115 (2004) [astro-ph/0403592].
  * (27) C. van Eldik, O. Bolz, I. Braun, G. Hermann, J. Hinton and W. Hofmann [for the HESS Collaboration], arXiv:0709.3729 [astro-ph].
  * (28) D. Hooper and B. L. Dingus, Phys. Rev. D 70, 113007 (2004) [astro-ph/0210617]; D. Hooper and B. Dingus, Proc. of the 34th COSPAR Scientific Assembly, Houston, Texas (2002), astro-ph/0212509.
  * (29) M. Pohl, Astrophys. J. 626, 174 (2005).
  * (30) S. Dodelson, D. Hooper and P. D. Serpico, Phys. Rev. D 77, 063512 (2008) [arXiv:0711.4621 [astro-ph]].
  * (31) M. Kachelrieß and P. D. Serpico, Phys. Lett. B 640, 225 (2006) [astro-ph/0605462].
  * (32) S. D. Hunter et al. [EGRET Collaboration], Astrophys. J. 481, 205 (1997).
  * (33) P. Sreekumar et al. [EGRET Collaboration], Astrophys. J. 494, 523 (1998) [astro-ph/9709257].
  * (34) F. W. Stecker and M. H. Salamon, astro-ph/0104368.
  * (35) L. E. Strigari, S. M. Koushiappas, J. S. Bullock and M. Kaplinghat, Phys. Rev. D 75, 083526 (2007) [astro-ph/0611925].
  * (36) G. Zaharijas and D. Hooper, Phys. Rev. D 73, 103501 (2006). [astro-ph/0603540].
  * (37) J. H. Buckley et al., whitepaper for “Dark Matter Searches with a Future VHE -Ray Observatory”, in preparation.