11institutetext:  INAF/Trieste Astronomical Observatory, via G.B. Tiepolo 11,
I-34100 Trieste, Italy
INFN-Trieste, via A. Valerio 2, I-34127 Trieste, Italy
11email<EMAIL_ADDRESS>22institutetext: School of Physics &
Astronomy, Tel Aviv University, Tel Aviv 69978, Israel
Center for Astrophysics and Space Sciences, University of California at San
Diego, La Jolla, CA 92093, USA
22email<EMAIL_ADDRESS>
# Diffuse non-thermal emission in the disks of the Magellanic Clouds

M. Persic 11 Y. Rephaeli 22

###### Abstract

Context. The Magellanic Clouds, two dwarf galaxy companions to the Milky Way,
are among the Fermi Large Area Telescope (LAT) brightest $\gamma$-ray sources.

Aims. Comprehensive modeling of the non-thermal electromagnetic and neutrino
emission in both Clouds.

Methods. We self-consistently model the radio and $\gamma$-ray spectral energy
distribution from their disks based on recently published Murchison Widefield
Array and Fermi/LAT data. All relevant radiative processes involving
relativistic and thermal electrons (synchrotron, Compton scattering, and
bremsstrahlung) and relativistic protons (neutral pion decay following
interaction with thermal protons) are considered, using exact emission
formulae.

Results. Joint spectral analyses indicate that radio emission in the Clouds
has both primary and secondary electron synchrotron and thermal bremsstrahlung
origin, whereas $\gamma$ rays originate mostly from $\pi^{0}$ decay with some
contributions from relativistic bremsstrahlung and Compton scattering off
starlight. The proton spectra in both galaxies are modeled as power laws in
energy with similar spectral indices, $\sim$2.4, and energy densities, $\sim$1
eV cm-3. The predicted 0.1–10 GeV neutrino flux is too low for detection by
current and upcoming experiments.

Conclusions. We confirm earlier suggestions of a largely hadronic origin of
the $\gamma$-ray emission in both Magellanic Clouds.

###### Key Words.:

galaxies: cosmic rays – galaxies: individual: Large Magellanic Cloud –
galaxies: individual: Small Magellanic Cloud – gamma rays: galaxies –
radiation mechanisms: non-thermal

##  1 Introduction

The spectral energy distributions (SEDs) of Cosmic Rays (CRs) outside their
galactic accelerators are important for determining basic properties of CR
populations and for assessing the impact of the particle interactions in the
magnetized plasma in galactic disks and halos. Knowledge of these
distributions is generally limited as it is usually based only on spectral
radio observations. When, in addition, non-thermal (NT) X/$\gamma$-ray
observations are available, reasonably detailed spectral modeling of the CR
electron (CRe) and proton (CRp) distributions in star-forming environments can
be very useful: Sampling the SED over these spectral regions yields important
insight on the emitting CRe and possibly also CRp, whose interactions with the
ambient plasma may dominate (via $\pi^{0}$ decay) high energy ($\sim$$>$ 100
MeV) emission.

The nearby Magellanic Clouds (MC), two irregular dwarf satellite galaxies of
the Milky Way, may exemplify the level of spectral modeling currently feasible
in a joint analysis of radio and $\gamma$-ray measurements. The Large MC (LMC)
is located at a distance $d=50$ kpc (Foreman et al. 2015); for the purpose of
this study its structure is modeled as a cylinder with radius $r=3.5$ kpc and
height $h=0.4$ kpc (Ackermann et al. 2016). Likewise, for the Small MC (SMC)
these quantities are $d=60$ kpc, $r=1.6$ kpc, and $h=4$ kpc (Abdo et al.
2010b).

Both MCs have been detected with the Fermi Large Area Telescope (LAT) as
extended sources (Abdo et al. 2010a,b). Based on recent data (Lopez et al.
2018), there is statistically significant emission along the SMC ”bar” and
”wing” where there is active star formation. A spectro-spatial analysis of the
LAT data suggests its integrated $\gamma$-ray emission to be mostly
($\geq$90%) hadronic. Similar analyses of the LMC with multi-source models
(Foreman et al. 2015; Ackermann et al. 2016; Tang et al. 2017) attempted to
extract the spectral content of these sources in a statistically viable way
(including details on the quality of the fits). The sources were found to be
point-like (pulsars, supernova remnants, plerions, unidentified background
AGNs) and extended. The latter included the large-scale disk denoted as 111
Emission components are labeled with ”E” in Ackermann et al. (2016) and ”G” in
Tang et al. (2017).  E0 $\sim$ G1, and smaller-scale components denoted as
E1+E3 $\sim$ G2 (the region west of 30 Doradus), E2 $\sim$ G3 (the LMC
northern region), and E4 $\sim$ G4 (the LMC western region) – all the latter
possibly encompassing different sources with different spectra. These
spectral/spatial analyses revealed a galaxy-scale component dominating the
emission, with a spectral shape suggesting a lepto-hadronic (Foreman et al.
2015) or hadronic (Ackermann et al. 2016; Tang et al. 2017) origin. 222 As to
the prominent 30 Dor massive-star–forming region, Ackermann et al. (2016)
emphasize that it shines in GeV $\gamma$ rays mainly because of the presence
of PSR J0540-6919 and PSR J0537-6910.

In this paper we focus on the emission in the MC disks in an attempt to
determine mean values of their magnetic fields and CRe and CRp energy
densities. In spite of their detailed LAT data analyses, the above-mentioned
studies were not based on self-consistent modeling of the broadband
radio/$\gamma$ SED of the two galaxies. SED modeling is important because a
firm proof of the (mostly) pionic nature of the $\gamma$-ray emission should
be based on a quantitative estimate of the CRe spectrum. Clearly, the latter
is primarily determined from measurements of synchrotron radio emission –
which was not accounted for in those earlier works. In an attempt to clarify,
and possibly remove, some uncertainties inherent in previous modeling work,
here we re-assess key aspects of (average) conditions in the MC large-scale
disks. Employing a one-zone model for the extended disk emission, we self-
consistently carry out detailed calculations of the emission by CRe and CRp.
The radio and $\gamma$-ray data used in our analyses are taken from the afore-
mentioned papers Lopez et al. (2018; SMC), Tang et al. (2017; LMC), Ackermann
et al. (2016; LMC), and For et al. (2018; SMC, LMC).

In Section 2 we briefly review the observations of extended NT emission from
the MC disks. In section 3 we review the IR/optical radiation fields
permeating the MCs. In Section 4 we describe calculations of the MC disk SEDs
and perform fits to the data. Prospects of neutrino detection are discussed in
Section 5, followed with a conclusion in Section 6.

##  2 Observations of extended emission

The MCs have been extensively observed over a wide range of radio/microwave,
(soft) X-ray, and $\gamma$-ray bands. As mentioned above, point-source, and
extended small and large scale emission has been detected from both galaxies.
In this paper we focus on the extended large-scale disk emission of each
galaxy because such emission traces the mean galactic properties of NT
particles and magnetic fields. The spectral data sets used in our analysis are
public (either tabulated or plotted) and are fully specified in Table 1 (SMC)
and Table 2 (LMC). In this section we briefly review the observations most
relevant to NT emission leaving out details that are discussed in the cited
papers.

$\bullet$ Radio. In a multifrequency radio continuum study of the MCs, For et
al. (2018) presented closely-sampled 76–227 MHz Murchison Widefield Array
data, supplemented by previous radio measurements at lower (for the LMC), and
higher (for both MCs) frequencies. Measured fluxes include emission from MC
(and background) point sources, which were estimated to contribute 11% and 23%
of the measured emission of, respectively, the SMC and the LMC. Based on
statistical analyses of a set of four fitting spectral models (power-law
[1PL], curved PL, and double PL [2PL]) with either free and fixed high-$\nu$
index), they found that:

(a) for the SMC the best-fitting model is a 1PL (in the 76–8550 MHz band) with
$\alpha($85.5 MHz – 8.55 GHz$)=0.82\pm 0.03$ (consistent with Haynes et al.
1991); whereas

(b) for the LMC the preferred model is a 2PL (19.7–8550 MHz) with a low-$\nu$
(19.7–408 MHz) index $\alpha_{0}=0.66\pm 0.08$ (consistent with Klein et al.
1989) and a (fixed) high-$\nu$ index $\alpha_{1}=0.1$ suggesting that
synchrotron radiation and thermal free-free (ff) emission dominate at,
respectively, low and high frequencies.

$\bullet$ X-rays. The observed diffuse X-ray emission in the MCs has been
determined to be of thermal origin (Wang et al. 1991: Einstein Observatory;
Points et al. 2001: Röntgen Satellit (ROSAT) Position Sensitive Proportional
Counters (PSPC); Nishiuchi et al. 1999, 2001: Advanced Satellite for Cosmology
and Astrophysics (ASCA)); as such, it is not directly relevant to our SED
analysis except for the estimated thermal plasma density, which is required to
calculate the pionic emission (see below).

$\bullet$ $\gamma$-rays. Both MCs were detected at $>$100 MeV $\gamma$-rays as
galaxy-scale extended sources whose emission is dominated by a diffuse
component. The latter has been interpreted as mainly originating from CRp
interacting with the interstellar gas.

(a) The SMC was first detected with 17 months of Fermi/LAT data (Abdo et al.
2010b) as a $\sim$3o-size source, in which emission was not strongly
correlated with prominent sites of star formation. More recently Lopez et al.
(2018), based on 105 months of LAT Pass8 data, produced maps of the extended
$>$2 GeV emission (no signal at $>$13 GeV) and found statistically significant
emission along the galaxy-scale quietly–star-forming ”bar” and ”wing”. Within
a set of single-component spectral models – i.e. PL, broken PL (with frozen
spectral shape and free normalization, representing pulsars), and
exponentially-cutoff PL (representing pionic emission) – the latter provides
the best fit to the total $\gamma$-ray spectrum, with only a marginal
improvement of the fit if the broken PL is added. Lopez et al.’s analysis
suggests that although pulsars may contribute $\leq$10% at $>$100 MeV, the
extended emission is mainly ($\geq$90%) of pionic origin, with a $\gamma$-ray
emissivity $\sim$$>$ 5 times smaller than the local (Galactic) one.

(b) The LMC was marginally (¿4.5 $\sigma$) detected with the Compton Gamma Ray
Observatory’s (CGRO) Energetic Gamma Ray Telescope (EGRET) (Sreekumar et al.
1992) and then confirmed (33 $\sigma$) with 11 months of Fermi/LAT data (Abdo
et al. 2010a). More recent studies have focused on the spatial and spectral
modeling of the $\gamma$-ray surface brightness distribution; these are
briefly reviewed below.
(i) The first LAT-based spectro-spatial analysis of the $\gamma$-ray surface
brightness distribution was performed by Foreman et al. (2015) using 5.5 years
of data, in the photon energy range 0.2$-$20 GeV. They modeled the CR
distribution and $\gamma$-ray production based on observed maps of the LMC
interstellar medium, star formation, radiation fields, and radio emission. The
$\gamma$-ray spectrum was described by means of analytical fitting functions
for the $\pi^{0}$-decay, bremsstrahlung, and inverse-Compton yields: these
processes were estimated to account for, respectively, 50%, 44%, and 6% of the
total 0.2$-$20 GeV emission. In particular they inferred the CRp spectral
index to be $q_{p}=2.4\pm 0.2$ and the equipartition (with CRp) magnetic field
to be $B_{\rm eq}=2.8\mu$G.
(ii) A subsequent spectro-spatial analysis (Ackermann et al. 2016; A+16) of 6
years of LAT (P7REP) data in the 0.2–100 GeV band reported extended and point-
like emissions. The (dominant) extended emission is in the form of a disk-
scale component, denoted as E0 in their paper, and additional degree-scale
emissions from several regions of enhanced star formation (including 30
Doradus): if pionic, the E0 emission component implies a population of CRp
with $\sim$1/3 the local Galactic density, whereas the superposed small-scale
emissions imply local enhancements of the CRp density by factors of least 2-6.
The spectrum of the E0 component (which does not include emission from the
degree-scale emissions from several regions) shows some curvature (see their
Fig.7-top): with no reference to the underlying emission mechanism, A+16
fitted the data with a tabulated function derived from the local (Galactic)
gas emissivity spectrum (as a result of its interactions with CRp), an
exponentially-cutoff PL, a broken PL, and a log-parabola, and concluded that
the best-fitting analytical model was the log-parabola – noting, however, that
the similarity between the log-parabola and tabulated-function models (see
their Fig.9) suggests a pionic nature of the large-scale disk emission.
(iii) More recently Tang et al. (2017; T+17) analyzed 8 years of LAT Pass8
data and deduced a deeper, more spectrally extended 0.08–80 GeV spectrum of
the large-scale disk component (identified as G1; this too did not include
emission from other regions) largely overlaps with E0 of A+16. The G1 shape,
now extending down to 60 MeV, was determined to be best described as a
$\pi^{0}$-decay hump from a CRp spectrum harder than that in our Galaxy –
although they could not rule out CRe relativistic bremsstrahlung completely.
In this analysis we use the T+17 G1 data as our reference set for the diffuse
large-scale LMC disk emission, and the earlier A+16 E0 data as an auxiliary
set for a consistency check.

Table 1: SMC radio and $\gamma$-ray data. Frequency | Energy Flux | Reference |  | Frequency | Energy Flux | Reference
---|---|---|---|---|---|---
log($\nu$/Hz) |  $10^{-12}$erg/(cm2s) |  |  | log($\nu$/Hz) |  $10^{-12}$erg/(cm2s) |
7.932 | $0.393\pm 0.171$ | Mills (1955, 1959) |  | 8.340 | $0.481\pm 0.125$ | For et al. (2018)
7.881 | $0.307\pm 0.080$ | For et al. (2018) |  | 8.356 | $0.489\pm 0.127$ | For et al. (2018)
7.924 | $0.245\pm 0.064$ | For et al. (2018) |  | 8.611 | $0.543\pm 0.041$ | Loiseau et al. (1987)
7.964 | $0.230\pm 0.060$ | For et al. (2018) |  | 9.146 | $0.588\pm 0.084$ | Loiseau et al. (1987)
7.996 | $0.260\pm 0.068$ | For et al. (2018) |  | 9.146 | $0.486\pm 0.028$ | For et al. (2018)
8.029 | $0.379\pm 0.099$ | For et al. (2018) |  | 9.362 | $0.713\pm 0.138$ | Mountfort et al. (1987)
8.061 | $0.311\pm 0.081$ | For et al. (2018) |  | 9.362 | $0.623\pm 0.161$ | For et al. (2018)
8.090 | $0.316\pm 0.082$ | For et al. (2018) |  | 9.389 | $0.637\pm 0.073$ | Haynes et al. (1991)
8.114 | $0.317\pm 0.083$ | For et al. (2018) |  | 9.677 | $0.902\pm 0.190$ | Haynes et al. (1991)
8.155 | $0.464\pm 0.121$ | For et al. (2018) |  | 9.932 | $1.282\pm 0.342$ | Haynes et al. (1991)
8.176 | $0.387\pm 0.101$ | For et al. (2018) |  | $22.860^{\scriptscriptstyle+0.125}_{\scriptscriptstyle-0.167}$ | $7.047\pm 0.860$ | Lopez et al. 2018
8.199 | $0.417\pm 0.109$ | For et al. (2018) |  | $23.161^{\scriptscriptstyle+0.125}_{\scriptscriptstyle-0.167}$ | $7.998\pm 0.853$ | Lopez et al. 2018
8.220 | $0.360\pm 0.094$ | For et al. (2018) |  | $23.462^{\scriptscriptstyle+0.125}_{\scriptscriptstyle-0.167}$ | $7.194\pm 0.658$ | Lopez et al. 2018
8.241 | $0.620\pm 0.161$ | For et al. (2018) |  | $23.781^{\scriptscriptstyle+0.125}_{\scriptscriptstyle-0.167}$ | $5.284\pm 0.563$ | Lopez et al. 2018
8.258 | $0.512\pm 0.133$ | For et al. (2018) |  | $24.055^{\scriptscriptstyle+0.125}_{\scriptscriptstyle-0.167}$ | $5.598\pm 0.683$ | Lopez et al. 2018
8.276 | $0.473\pm 0.123$ | For et al. (2018) |  | $24.383^{\scriptscriptstyle+0.125}_{\scriptscriptstyle-0.167}$ | $3.357\pm 0.754$ | Lopez et al. 2018
8.294 | $0.486\pm 0.127$ | For et al. (2018) |  | $24.684^{\scriptscriptstyle+0.125}_{\scriptscriptstyle-0.167}$ | $<2.972$ | Lopez et al. 2018
8.310 | $0.594\pm 0.155$ | For et al. (2018) |  | $24.985^{\scriptscriptstyle+0.125}_{\scriptscriptstyle-0.167}$ | $<0.879$ | Lopez et al. 2018
8.326 | $0.489\pm 0.127$ | For et al. (2018) |  |  |  |
Table 2: LMC radio and $\gamma$-ray data. Frequency | Energy Flux | Reference |  | Frequency | Energy Flux | Reference
---|---|---|---|---|---|---
log($\nu$/Hz) |  $10^{-12}$erg/(cm2s) |  |  | log($\nu$/Hz) |  $10^{-12}$erg/(cm2s) |
7.294 | $1.038\pm 0.208$ | Shain (1959) |  | 9.389 | $9.555\pm 0.490$ | Haynes et al. (1991)
7.653 | $1.349\pm 0.203$ | Alvarez et al. (1987) |  | 9.677 | $17.243\pm 1.425$ | Haynes et al. (1991)
7.881 | $1.410\pm 0.240$ | For et al. (2018) |  | 9.677 | $14.250\pm 1.900$ | For et al. (2018)
7.924 | $1.492\pm 0.254$ | For et al. (2018) |  | 9.932 | $23.085\pm 2.993$ | Haynes et al. (1991)
7.932 | $3.154\pm 0.342$ | Mills (1959) |  |  | $=========$ |
7.964 | $1.449\pm 0.246$ | For et al. (2018) |  |  |  $10^{-11}$erg/(cm2s) |
7.986 | $2.748\pm 0.581$ | Mills (1959) |  |  | ——————— |
7.996 | $1.437\pm 0.245$ | For et al. (2018) |  | $22.2946\pm 0.1336$ | $0.93\pm 0.31$ | Tang et al. (2017)
8.029 | $1.956\pm 0.333$ | For et al. (2018) |  | $22.5635\pm 0.1353$ | $1.64\pm 0.33$ | Tang et al. (2017)
8.061 | $1.872\pm 0.318$ | For et al. (2018) |  | $22.8325\pm 0.1336$ | $1.91\pm 0.19$ | Tang et al. (2017)
8.090 | $2.021\pm 0.344$ | For et al. (2018) |  | $23.1004\pm 0.1343$ | $2.08\pm 0.13$ | Tang et al. (2017)
8.114 | $2.043\pm 0.347$ | For et al. (2018) |  | $23.3691\pm 0.1344$ | $2.17\pm 0.15$ | Tang et al. (2017)
8.155 | $2.379\pm 0.405$ | For et al. (2018) |  | $23.6376\pm 0.1342$ | $1.57\pm 0.15$ | Tang et al. (2017)
7.653 | $1.349\pm 0.203$ | Alvarez et al. (1987) |  | $23.9061\pm 0.1343$ | $1.44\pm 0.17$ | Tang et al. (2017)
8.176 | $2.175\pm 0.370$ | For et al. (2018) |  | $24.1747\pm 0.1342$ | $1.16\pm 0.21$ | Tang et al. (2017)
8.199 | $2.134\pm 0.363$ | For et al. (2018) |  | $24.4432\pm 0.1343$ | $0.52\pm 0.24$ | Tang et al. (2017)
8.199 | $2.743\pm 0.774$ | Mills (1959) |  | $24.7117\pm 0.1342$ | $1.09\pm 0.31$ | Tang et al. (2017)
8.220 | $1.999\pm 0.340$ | For et al. (2018) |  | $24.9801\pm 0.1342$ | $<0.73$ | Tang et al. (2017)
8.241 | $2.324\pm 0.397$ | For et al. (2018) |  | $25.2486\pm 0.1342$ | $0.63\pm 0.36$ | Tang et al. (2017)
8.258 | $2.257\pm 0.384$ | For et al. (2018) |  | $22.781\pm 0.100$ | $3.524\pm 0.912$ | Ackermann et al. (2016)
8.276 | $2.313\pm 0.393$ | For et al. (2018) |  | $22.985\pm 0.100$ | $3.357\pm 0.349$ | Ackermann et al. (2016)
8.294 | $2.186\pm 0.372$ | For et al. (2018) |  | $23.182\pm 0.100$ | $3.041\pm 0.309$ | Ackermann et al. (2016)
8.310 | $2.520\pm 0.429$ | For et al. (2018) |  | $23.383\pm 0.100$ | $2.958\pm 0.256$ | Ackermann et al. (2016)
8.326 | $2.378\pm 0.404$ | For et al. (2018) |  | $23.587\pm 0.100$ | $2.685\pm 0.232$ | Ackermann et al. (2016)
8.340 | $2.261\pm 0.385$ | For et al. (2018) |  | $23.781\pm 0.100$ | $1.811\pm 0.254$ | Ackermann et al. (2016)
8.356 | $2.315\pm 0.394$ | For et al. (2018) |  | $23.985\pm 0.100$ | $0.859\pm 0.240$ | Ackermann et al. (2016)
8.611 | $3.774\pm 0.122$ | Klein et al. (1989) |  | $24.161\pm 0.100$ | $1.167\pm 0.326$ | Ackermann et al. (2016)
9.146 | $5.367\pm 0.420$ | For et al. (2018) |  | $24.383\pm 0.100$ | $0.980\pm 0.315$ | Ackermann et al. (2016)
9.146 | $7.406\pm 0.420$ | Klein et al. (1989) |  | $24.587\pm 0.100$ | $<0.706$ | Ackermann et al. (2016)
9.362 | $9.476\pm 1.150$ | Mountfort et al. (1987) |  | $25.183\pm 0.504$ | $<1.119$ | Ackermann et al. (2016)

##  3 Radiation fields

A reasonably precise determination of the ambient radiation field is needed
for predicting the level of $\gamma$-ray emission from Compton scattering of
the radio-emitting electrons (and positrons). The total radiation field
includes cosmic (background) and local (foreground) components.

Relevant cosmic radiation fields include the Cosmic Microwave Background (CMB)
and the Extragalactic Background Light (EBL). The CMB is a pure Planckian
described by a temperature $T_{\rm CMB}=2.735\,(1+z)$ K and energy density
$u_{\rm CMB}=0.25\,(1+z)^{4}$ eV cm-3. The EBL originates from direct and
dust-reprocessed starlight integrated over the star formation history the
Universe. It shows two peaks, corresponding respectively to the Cosmic
Infrared Background (CIB, at $\sim$100 $\mu$m) that originates from dust-
reprocessed starlight integrated over the star formation history of galaxies,
and the Cosmic Optical Background (COB, at $\sim$1 $\mu$m) that originates
from direct starlight integrated over all stars that formed (e.g. Cooray
2016). The two peaks are described as diluted Planckians, characterized by a
temperature $T$ and a dilution factor $C_{\rm dil}$. The latter is the ratio
of the actual energy density, $u$, to the energy density of an undiluted
blackbody at the same temperature $T$, i.e. $u=C_{\rm dil}\,aT^{4}$, where $a$
is the Stefan-Boltzmann constant. The dilution factors of the cosmic fields
are, $C_{\rm CMB}=1$, $C_{\rm CIB}=10^{-5.629}$, and $C_{\rm
COB}=10^{-13.726}$. A recent updated EBL model, based on accurate galaxy
counts in several spectral bands, is due to Franceschini & Rodighiero (2017);
locally ($z=0$) it can be numerically approximated as a combination of diluted
Planckians,

$\displaystyle n_{\rm EBL}(\epsilon)\leavevmode\nobreak\ =\leavevmode\nobreak\
\sum_{j=1}^{8}A_{j}\,\frac{8\pi}{h^{3}c^{3}}\,\frac{\epsilon^{2}}{e^{\epsilon/k_{B}T_{j}}-1}\hskip
14.22636pt{\rm cm^{-3}\leavevmode\nobreak\ erg^{-1}}$ (1)

with: $A_{1}=10^{-5.629}$, $T_{1}=29$ K; $A_{2}=10^{-8.522}$, $T_{2}=96.7$ K;
$A_{3}=10^{-10.249}$, $T_{3}=223$ K; $A_{4}=10^{-12.027}$, $T_{4}=580$ K;
$A_{5}=10^{-13.726}$, $T_{5}=2900$ K; $A_{6}=10^{-15.027}$, $T_{6}=4350$ K;
$A_{7}=10^{-16.404}$, $T_{7}=5800$ K; $A_{8}=10^{-17.027}$, $T_{8}=11600$ K.

Local radiation fields in the MCs arise from their intrinsic stellar
populations and (given its close proximity) also from the Milky Way; we refer
to this as the Galactic Foregound Light (GFL). Similar to the EBL by shape and
origin, the GFL is dominated by two thermal humps, IR and optical. The
corresponding energy densities (in eV cm-3) are: (i) LMC: $u_{\rm IR}=0.12$
and $u_{\rm opt}=0.20$, estimated from Foreman et al. 2015 ($u_{\rm
IR+opt}=0.32$) and the IR/opt SED (from NED 333 The NASA/IPAC Extragalactic
Database (NED) is operated by the Jet Propulsion Laboratory (JPL), Caltech,
under contract with the National Aeronautic and Space Administration (NASA).
); (ii) SMC: $u_{\rm IR}=0.026$, $u_{\rm opt}=0.062$, scaling down the LMC
values taking (from the SED, cf. NED) the optical and IR peaks to be factors
of 0.215 and 0.144, respectively, of the corresponding LMC values.

Figure 1:  The MC $\gamma$-ray spectra. Symbols are as in Fig.4.

Figure 2:  CRe spectra in the Magellanic Clouds. Curves denote secondary
spectra (dashed), primary spectra (dotted), and the latter’s local PL
representations (solid). In each panel, the normalization of the primary-CRe
spectrum to its PL counterpart is chosen such that they nearly overlap. The
primary-CRe injection indices, $q_{i}$, are 2.28 (SMC) and 2.23 (LMC). The two
LMC CRe spectra panels refer, respectively, to the T+17 and A+16 Fermi/LAT
datasets, as secondary spectra are derived from the pionic fits to the
$\gamma$-ray data.

Figure 3:  SMC radio spectrum: data points (dots) are from Table 3 of For et
al. (2018). Left: the ”preferred” 1PL model 1 of For et al. (2018) is shown as
a solid line, and the best fitting 3rd-order polynomial described in Table 3
is shown as a dashed curve. Right: the two-component model (solid curve) that
mimicks the polynomial includes synchrotron radiation (dotted/short-dashed
curve; dashed, dotted curves denote primary, secondary contributions) and
thermal-ff emission (dotted/long-dashed curve).  Figure 4:  The broad-band (76
MHz – 53.7 GeV) NT SED of the SMC disk. Data points (Table 1) are shown by
dots, model predictions by curves. Emission components are plotted by the
following line types: synchrotron, dotted/long-dashed; thermal-ff,
dotted/short-dashed; total radio, solid; Compton/CMB, long-dashed;
Compton/starlight, short-dashed; NT bremsstrahlung, dotted; pionic, solid. For
each type of NT leptonic yield secondary, primary, and total emissions are
denoted as curves of progressively higher flux.

##  4 SED models

Our main objective in this study is to determine CR electron and proton
spectra in the MC disks by spectral modeling of NT emission in all relevant
energy range accessible to observations. The particle, gas density, magnetic
and radiation field distributions clearly vary significantly across the disk,
obviating the need for a spectro-spatial treatment. Indeed, a modeling
approach based on a solution to the diffusion-advection equation has been
applied in the study of the Galaxy and several nearby galaxies. However, even
for (just) a diffusion-based treatment to be feasible and meaningful, the
spatial profiles of key quantities, such as the particle acceleration sources,
gas density, magnetic field, and generally also the diffusion coefficient have
to be specified. Given the generally very limited observational basis for
reliably determining these profiles, a parameter-intensive spectro-spatial
modeling approach is rarely warranted. An example is the very approximate
diffusion-based approach adopted in modeling NT emission in the disks and
halos of the star-forming galaxies NGC4631 and NGC4666 (Rephaeli & Sadeh 2019)
for which reasonably detailed radio spectra and spatial profiles are
available; even so, the results of these analyses are not conclusive given the
substantial uncertainty in the values of key parameters.

Radio emission in the MC disks is electron synchrotron in a disordered
magnetic field whose mean value $B$ is taken to be spatially uniform, and
thermal-ff from a warm ionized plasma. A significant CRp component could yield
additional radio emission by secondary $e^{\pm}$ produced by $\pi^{\pm}$
decays, and $\gamma$-ray emission from $\pi^{0}$ decay. In addition,
$\gamma$-ray emission is produced by CRe scattering off photons of local and
background radiation fields. The calculations of the emissivities from all
these processes are well known and standard.

We assume the CR spectral distributions to be time-independent and locally
isotropic with spectral PL form. The primary CRe spectral density (per cm3 and
per unit of the electron Lorentz factor $\gamma$) is,
$N_{e}(\gamma)=N_{e,0}\,\gamma^{-q_{e}}$ for
$\gamma_{min}<\gamma<\gamma_{max}$ with $\gamma_{min}>>1$. As discussed later,
this spectrum proves to be a good approximation to the actual primary steady-
state spectrum in the relevant electron energy range. The secondary CRe
spectrum can be analytically approximated as a smoothed-2PL, exponentially
cutoff at high energies (see below). The assumed CRp spectrum (in units of
cm-3 GeV-1) as a function of $E_{p}$/GeV is,
$N_{p}(E_{p})=N_{p,0}\,E_{p}^{-q_{p}}$ for $m_{p}<E<E_{p}^{max}$.

###  4.1 SMC

The emission spectrum from $\pi^{0}$-decay has more constraining power than a
generic PL; thus, our modeling procedure begins with fitting a pionic emission
profile to the $\gamma$-ray data with free normalization, slope, and high-
energy cutoff. We then use the deduced (essentially, fully-determined)
secondary-CRe spectrum together with an assumed primary-CRe spectrum to
calculate the combined synchrotron emission using an observationally estimated
value of $B$. The fit to the radio data includes, at high frequencies, also a
thermal bremsstrahlung component computed with previously determined values of
the gas density and temperature. Finally, the full Compton X/$\gamma$-ray
yield is determined.

The $\pi^{0}$-decay $\gamma$-ray flux is computed using the emissivity in Eq.
(15) of Persic & Rephaeli (2019a), where the total (thermal) gas proton
density is $n_{g}=n_{\rm HI}+2\,n_{\rm H_{2}}$ with the average neutral-
hydrogen density $n_{\rm HI}=0.54$ cm-3, inferred from $M_{\rm
HI}=4.23\,10^{8}M_{\sun}$ (this value of $M_{\rm HI}$ was deduced from direct
determination of HI column density from the 21 cm emission line; Stanimirović
et al. 1999), and the molecular-hydrogen density $n_{\rm H_{2}}=0.02$ cm-3,
inferred from $M_{\rm H_{2}}=3.2\,10^{7}M_{\sun}$ (which, in turn, is derived
from modeling Spitzer, Cosmic Background Explorer (COBE), InfraRed
Astronomical Satellite (IRAS), and Infrared Space Observatory (ISO) far-
infrared data; Leroy et al. 2007). The CRp spectrum derived from our fit to
the Fermi/LAT data is $N_{p0}=2.4\,10^{-10}$ cm-3, $q_{p}=2.40$, with
$E_{p}^{max}=30$ GeV. The model spectrum fully reproduces the 0.2–50 GeV
Fermi/LAT spectrum (Fig.1). With these values, the estimated CRp energy
density $u_{p}=\int_{m_{p}}^{E_{p}^{max}}E_{p}\,N_{p}(E_{p})\,dE_{p}$ is $\sim
0.5$ eV cm-3.

The closely related $\pi^{\pm}$-decay secondary CRe spectrum,
$N_{se}(\gamma)$, has no free parameters once the CRp spectrum is determined
444 Denoting the total cross-section for inelastic $pp$ collisions
$\sigma_{pp}(E_{p})$ (see Eq.[79] of Kelner et al. 2006), the secondary CRe
injection spectrum is
$Q_{se}(\gamma)=(8/3)\,m_{e}\,(c/k_{\pi^{\pm}})\,n_{g}\,N_{p0}\,[m_{p}+(4m_{e}/\kappa_{\pi^{\pm}})\,\gamma]^{-q_{p}}\,\sigma_{pp}(m_{p}+E_{\pi^{\pm}}/k_{\pi^{\pm}})$
for $\gamma_{thr}<\gamma<\gamma_{se}^{max}$, with
$\gamma_{se}^{thr}=m_{\pi^{\pm}}/(4m_{e})=68.5$ and
$\gamma_{se}^{max}$$=$$m_{\pi^{\pm}}^{max}/(4m_{e})\simeq(E_{p}^{max}-3/2m_{p})/(4m_{e})$.
The corresponding steady-state distribution is
$N_{se}(\gamma)=1/b(\gamma)\cdot\int_{\gamma}^{\gamma_{se}^{max}}Q_{se}(\gamma)\,d\gamma$,
where $b(\gamma)$ is the radiative energy loss term. In a magnetized medium
consisting of ionized, neutral and molecular gas, it is
$b(\gamma)=\sum_{j=0}^{2}b_{j}(\gamma)$ where $b_{j}(\gamma)$ are loss terms
appropriate to each gas phase (see Eqs.[4]-[6] of Rephaeli & Persic 2015).
$N_{se}(\gamma)$ is analytically approximated by $N_{se}^{\rm
fit}(\gamma)=N_{se,0}\gamma^{-q_{1}}(1+\gamma/\gamma_{b1})^{q_{1}-q_{2}}{\rm
exp}[{-(\gamma/\gamma_{b2})^{\eta}}]$, where $q_{1}$, $q_{2}$ and
$\gamma_{b1}$, $\gamma_{b2}$ are the low-/high-energy spectral indices and
breaks; $\eta$ gauges the steepness of the high-end cutoff. . We use
$N_{se}^{\rm fit}(\gamma)$, with parameter values reported in Table 4, to
compute the corresponding leptonic yields. However, the uncertainty in the
total (HI and H2) gas density clearly affects the energy loss rate $b(\gamma)$
and the resulting spectral normalization and shape of $N_{se}(\gamma)$, hence
also of $N_{e}(\gamma)$, and their yields. To compute the synchrotron emission
we assume $B=3.5\mu$G, based on estimates of the ordered (1.7 $\mu$G) and
random (3 $\mu$G) fields in the SMC (Mao et al. 2008).

Once the secondary-CRe synchrotron yield has been computed, its low-frequency
residuals from the data are modeled using a CRe spectrum, $N_{e}(\gamma)$,
with $N_{e0}=3.7\,10^{-8}$ cm-3, $q_{e}=2.23$, $\gamma_{max}=10^{4}$. The
latter is the 1PL approximation, for $100\raise-2.0pt\hbox{\hbox
to0.0pt{\hbox{$\sim$}\hss}\raise 5.0pt\hbox{$<$}\
}\gamma\raise-2.0pt\hbox{\hbox to0.0pt{\hbox{$\sim$}\hss}\raise
5.0pt\hbox{$<$}\ }\gamma_{max}$, of the actual steady-state primary-CRe
spectrum. Neglecting diffusion and advection losses, this primary spectrum is
$N_{pe}(\gamma)=k_{i}\gamma^{-(q_{i}-1)}/[b(\gamma)\,(q_{i}-1)]$ cm-3 (units
of $\gamma$)-1 where $\dot{N}_{i}(\gamma)=k_{i}\,\gamma^{-q_{i}}$ is the CRe
spectral injection rate. We find that, over the mentioned electron energy
range, $q_{i}=2.28$ provides a decent match between the shapes of the curved
spectrum and the 1PL spectrum – this value is typical for $\gamma$-ray
emission from Galactic CR accelerators (e.g. SNRs, the Crab nebula). The
relative normalization between the two spectra can be based on the measured
flux density at some radio frequency or on imposing the same CRe energy
density on the two spectra over the electron energy range of interest (see
discussion in Rephaeli & Persic 2015). In our case here the normalization of
$N_{pe}(\gamma)$, and hence $k_{i}$, can be found by using both
$N_{pe}(\gamma)$ and $N_{se}(\gamma)$ to compute the total (primary plus
secondary) synchrotron yield, which in turn is fitted to the radio synchrotron
data: in doing this, the acceptable match (in the relevant energy range)
between the curved and PL spectra allows us to use $N_{e}(\gamma)$ without
substantial loss of accuracy. Given this procedure, the uncertainties in the
gas density ultimately affect $N_{e}(\gamma)$ as well.

The spectra of the synchrotron-emitting CRe are shown in Fig. 2. The resulting
primary electron component is quite dominant, with $\sim$70% of the total
synchrotron and Compton yields (Fig. 3), but the normalization, $N_{e0}$, is
not strongly constrained due to its dependence on the assumed magnetic field.
If measurements of NT–X-ray and $\sim$1-30 MeV $\gamma$-ray fluxes –
corresponding to, respectively, the CMB and CIB peaks – become available, most
of the uncertainty could be removed (e.g. Persic & Rephaeli 2019a). In
addition to this modeling uncertainty, the considerable coupling between the
CRe spectral parameters implies an appreciable range of the deduced values of
$N_{e0}$ and $q_{e}$: This range can be roughly bracketed by a flatter slope
($q_{e}=2.2$) and lower normalization ($N_{e0}=2.9\,10^{-8}$ cm-3) or,
alternatively, a steeper slope ($q_{e}=2.3$) and higher normalization
($N_{e0}=6.7\,10^{-8}$ cm-3) – both, with the same $\gamma_{max}$ as the
reference model (Table 4).

At high frequencies a relatively flat, $\propto\nu^{-0.1}$, component
represents diffuse thermal-ff emission (Spitzer 1978). The latter flux may be
gauged to the H$\alpha$ flux if both emissions come from the same emitting
volume (HII regions), because in this case the relevant warm-plasma parameters
(temperature, density, filling factor) are the same. So the measured (optical)
H$\alpha$ flux may be used to predict the (radio) ff emission. We model the
thermal-ff emission combining Eqs.(3),(4a) of Condon (1992) and using
$T_{e}=1.1\,10^{4}$ ${}^{\circ}K$ (Toribio San Cipriano et al. 2017) and
$F(H_{\alpha})=1.6\,10^{-8}$ erg cm-2 s-1 (Kennicutt et al. 1995). The
(primary) synchrotron and thermal-ff normalizations are spectrally quite apart
so they do not significantly affect each other. The radio model is shown
alongside data in Fig. 3-right.

Although subdominant, secondary CRe contribute appreciably to the total CRe
population of the SMC. Whereas the primary spectrum is approximately PL, the
secondary spectrum is clearly curved (see Fig. 2). The resulting curvature of
the total CRe spectrum is reflected in the total synchrotron spectrum
(computed using Eq. (9) of Persic & Rephaeli (2019a) for $N_{e}(\gamma)$ and
its straightforward generalization for $N_{se}(\gamma)$). Therefore the total
radio spectrum, which consists of synchrotron and thermal-ff emission at low
and high frequencies, is not a smooth 2PL but shows some extra structure. A
3rd-order polynomial (in log units) outperforms For et al.’s (2018)
”preferred” 1PL model ($\Delta$BIC$>0$; Table 3). This polynomial (4 free
parameters; Fig. 3-left) is matched by the physical model described above,
i.e. a low-frequency component representing synchrotron emission (3 free
parameters) plus a high-frequency $\propto\nu^{-0.1}$ PL representing thermal-
ff emission (1 free parameter); this model is shown in Fig. 3-right.

Table 3: Summary of fits to the SMC radio spectrum.

| ${\chi^{2}}^{\bullet}$ | $\chi^{2}_{\rm red}$ | BIC‡ |  $\Delta$BIC | dof | log$(S_{0})$ | $\alpha$ | $c_{0}$ | $c_{1}$ | $c_{2}$ | $c_{3}$
---|---|---|---|---|---|---|---|---|---|---|---
PL⋆ | 32.3 | 1.2 | 298.0 |  | 27 | 2.284 | $-0.82$ |  |  |  |
P${}_{3}^{+}$ | 22.3 | 0.9 | 294.7 | 3.3 | 25 |  |  | 1.763 |  $-$0.968 | 0.09 | 0.315

⋆ Power law: ${\rm log}(S_{\nu})={\rm log}(S_{0})+\alpha\,{\rm
log}(\nu/0.2\,{\rm GHz})$ from Eq.(4) and Table 4 of For et al. (2018).

+ 3rd-order polynomial: ${\rm log}(S_{\nu})=c_{0}+c_{1}\,x+c_{2}\,x^{2}+c_{3}\,x^{3}$, with $x={\rm log}(\nu/{\rm GHz})$.

∙ $\chi^{2}$ calculations use actual fluxes (Table 3 of For et al. 2018), not
their logarithms.

‡ See definitions and discussion in Sect. 4.2.1 of For et al. 2018.

With the CRe spectra essentially determined, we can calculate the Compton and
NT-bremsstrahlung yields from CRe scattering off, respectively, CMB/EBL/GFL
photons and thermal-plasma nuclei using, respectively, Eqs. (2.42) and (3.1)
of Blumenthal & Gould (1970). The photon fields are described in section 3; as
to the densities of plasma nuclei we assume $n_{i}=0.033$ cm-3 (Mao et al.
2008) and metallicity $Z^{2}=1.2$ (Sasaki et al. 2002). Although the shape of
the $\gamma$-ray spectrum is distinctly pionic, NT-bremsstrahlung peaks near
the pionic peak ($\nu\sim 10^{23}$ Hz) where it contributes $\sim$3% of the
measured flux; Compton/starlight emission, that peaks at $\nu\sim 10^{22}$ Hz,
contributes $\sim$$<$ 1% of the flux at the lowest LAT data point. The
predicted Compton/CMB spectral flux, too, is shown in Fig. 4. We found no need
to introduce a further, broken-PL component in the $\gamma$-ray band
representing pulsars, in agreement with Lopez et al. (2018).

The resulting SED model, overlaid on data, is plotted in Fig.4.

An upper limit on the mean magnetic field, $B_{max}$, can be estimated by
assuming that the measured radio emission is fully produced by secondary CRe.
By matching the measured and predicted emission levels we infer a limiting
value $B_{lim}\simeq 7\,\mu$G. This upper limit, which is twice the
observationally deduced value, is unrealistically high because, most notably
at low frequencies, the shapes of the measured and predicted synchrotron
spectra do not match well, and because primary CRe obviously contribute to the
measured radio flux. A lower limit, $B_{min}$, can estimated from the fact
that if $B$ decreases the secondary synchrotron yield decreases too, so to
keep the synchrotron yield unchanged, the primary yield must increase – by
increasing the spectral normalization ($N_{e0}$) and possibly also varying the
spectral shape ($q_{e}$, $\gamma_{max}$) of the primary CRe spectrum. At low
energies the primary spectrum is steeper than the secondary, so at some point
the resulting synchrotron spectrum becomes progressively and systematically
steeper than the $<$230 GHz radio data; this yields $B_{min}\simeq 2\,\mu$G.
Based on our SED analysis, the value $B=3.5\mu$G suggested by Mao et al.
(2008) lies comfortably between the estimated bounds.

Figure 5:  LMC radio spectrum: data points (dots) are from Table 2 of For et
al. (2018). Left: the 2PL model 3 of For et al. (2018) is shown as a solid
curve. Right: the two-component model (solid curve) includes synchrotron
radiation (dotted/short-dashed curve; dashed and dotted curves denote,
respectively, primary and secondary contributions - the latter corresponding
to T+17 $\gamma$-ray data) and thermal-ff emission (dotted/long-dashed curve).

Figure 6:  The broad-band (radio/$\gamma$, 19.7 MHz – 100 GeV) SED of the LMC
disk. The $\gamma$-ray data are from T+17 (top; reference set) and A+16
(bottom; auxiliary set). Symbols are as in Fig. 4

###  4.2 LMC

Our fitting approach differs from that adopted for the SMC. The fact the radio
spectrum is best fitted by a 2PL model interpreted as a combination of
(individually 1PL) synchrotron and thermal-ff components (For et al. 2018),
indicates that synchrotron emission is generated by a 1PL CRe spectrum. This
is unlike the case of the SMC for which the interplay of primary and secondary
CRe with different spectra results in a more structured, non-PL form. With the
primary CRe spectrum determined from fitting the radio data, the Compton
$\gamma$-ray yield is calculated. To this a modeled pionic component is added
and fit to the $\gamma$-ray data. Doing so enables determination of the the
secondary CRe spectrum and the resulting yields, which turn out to be very
minor in comparison with those of primary CRe; thus, no iteration of the
fitting process is required.

The radio flux is dominated by synchrotron emission with $\alpha=0.66$ at the
lowest measured frequencies, and by thermal-ff emission with $\alpha=-0.1$ at
the highest frequencies (For et al. 2018). Adopting these spectral indices for
the two components, we fit the full radio database with $T_{e}=1.3\,10^{4}$
${}^{\circ}K$ and $F(H\alpha)=2\,10^{-7}$ erg cm-2 s-1 (Toribio San Cipriano
et al. 2017; Kennicutt et al. 1995) for the thermal-ff component, and
$B=4.3\mu$G (Gaensler et al. 2005) 555 This field value, from a Faraday
rotation study of the LMC, is consistent with Mao et al.’s (2012) estimates,
$B_{\rm eq}\sim 2\mu$G from equipartition with CRp assuming $q_{p}=2$ and
$B<7\mu$G from using a (deduced) lower limit on $N_{e0}$ in fitting the 1.4
GHz radio flux. , $N_{e0}\simeq 3.5\,10^{-7}$ cm-3 and $\gamma_{\rm max}\simeq
2.4\,10^{4}$ for the synchrotron. (The cutoff is needed for consistency with
the measured flux at $\nu\raise-2.0pt\hbox{\hbox
to0.0pt{\hbox{$\sim$}\hss}\raise 5.0pt\hbox{$>$}\ }3\,10^{23}$ Hz.) We should
mention at this point that subsequent results (see below) suggest that the
secondary CRe are only $\sim$10% of the total. So the (quasi-)1PL primaries
dominate the CRe population (Fig. 2) and the emerging synchrotron emission.
Therefore, the combined synchrotron and thermal-ff radio spectrum matches the
smooth 2PL profile of For et al.’s (2018) fit. The radio model is shown in
Fig. 5.

Next, we calculate the X/$\gamma$-ray Compton and NT-bremsstrahlung yields
from (primary) CRe scattering off CMB/EBL/GFL photons and thermal-plasma
nuclei: to this aim we assume the photon fields described in section 3 and,
respectively, a plasma characterized by an average density $n_{i}=0.012$ cm-3
(Sasaki et al. 2002; Cox et al. 2006) and metallicity $Z^{2}=1.2$ (Sasaki et
al. 2002). Both emissions fall short of reproducing the $\gamma$-ray data.

The pionic emission is computed with $n_{\rm HI}=1$ cm-3 (from $M_{\rm
HI}=3.8\,10^{8}M_{\sun}$: Staveley-Smith et al. 2003) and $n_{\rm
H_{2}}=6.6\,10^{-2}$ cm-3 (from $M_{\rm H_{2}}=5\,10^{7}M_{\sun}$, measured
from CO emission: Fukui et al. 2008). From our fitting, the matching CRp
spectrum has $q_{p}=2.38$ and $E_{p}^{max}=80$ GeV for the T+17 data (Fig.
6-top), and $q_{p}=2.6$ and $E_{p}^{max}=25$ GeV for the A+16 data (Fig.
6-bottom). The corresponding CRp energy densities are, respectively, 1 and 1.5
eV cm-3 – the difference being largely due to the somewhat higher flux level
in the A+16 data. In spite of the observational uncertainties in the two
datasets, we see that for either the pionic yield largely dominates the
$\gamma$-ray emission, so the pionic nature of the measured $\gamma$-ray
spectrum appears well established. In addition to the dominant pionic
component a subdominant, 7% (5%), leptonic contributions peaks at 500 (400)
MeV, similarly to the pionic one: Compton/COB and NT bremsstrahlung emissions
contribute, respectively, 6.5% (4%) and 0.5% (1%) to the peak of the observed
emission in the T+16 (A+16) database.

The secondary CRe spectrum is derived and fitted analytically as described in
section 4.1: the fitting parameters are reported in Table 4 for both sets of
LAT data. As mentioned, in both cases the radiative contribution of the
secondary CRe is minor: this means that the corresponding primary PL spectra
are essentially the same 666 For the adopted nominal densities of the gas and
of the magnetic and radiation fields, the inferred primary CRe injection index
in the LMC is $q_{i}=2.22$ for both the T+17 and A+16 fits.  and (for each
data set) no fitting iteration is required. The primary and secondary CRe
spectra are shown in Fig. 2.

Primary CRe largely dominate the CRe population in the LMC, even more so than
in the SMC. However, their density is effectively unconstrained because it
depends on assuming a magnetic field value (even though deduced from, as in
our case, Faraday-rotation measurements) rather than using a value derived by
modeling the NT–X-ray emission as Compton/CMB (e.g., Persic & Rephaeli
2019a,b) or the hard–X-ray/MeV emission (when it becomes available) as
Compton/IR. Also, appreciable modeling uncertainty stems from the inter-
dependence among primary-CRe parameters which results in a range of acceptable
values of the spectral index and normalization, $q_{e}=2.26$ with
$N_{e0}=2.22\,10^{-7}$ cm-3 and $\gamma_{max}=2.0\,10^{4}$, as well as
$q_{e}=2.40$ with $N_{e0}=6.21\,10^{-7}$ cm-3 and $\gamma_{max}=2.7\,10^{4}$.
777 Values referring to the T+17 LAT data. A similar covariance, with only
slightly different values, is found when the A+16 LAT data are used.

Table 4: SED model parameters.

| $N_{e0}$ | $q_{e}$ | $\gamma_{max}$ | $u_{p}$ | $q_{p}$ | $E_{p}^{max}$ | $N_{se0}$ | $q_{1}$ | $q_{2}$ | $\gamma_{b1}$ | $\gamma_{b2}$ | $\eta$ | F(H$\alpha$) | $T_{e}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| cm-3 |  | $10^{4}$ | eV/cm3 |  | GeV |  $10^{-11}$cm-3 |  |  | $10^{2}$ | $10^{4}$ |  | erg/(cm2s) | $10^{4\,\circ}K$
SMC | $3.70\leavevmode\nobreak\ 10^{-8}$ | 2.23 | $1.0$ | 0.45 | 2.40 | 30 | $0.14$ | 0.10 | 2.60 | 1.5 | $0.95$ | 2.65 | $1.6\leavevmode\nobreak\ 10^{-8}$ | $1.1$
LMC+ [1] | $3.50\leavevmode\nobreak\ 10^{-7}$ | 2.32 | $2.4$ | 1.00 | 2.38 | 80 | $0.50$ | 0.05 | 2.67 | 1.2 | $2.75$ | 3.20 | $2.0\leavevmode\nobreak\ 10^{-7}$ | $1.3$
[2] | $3.30\leavevmode\nobreak\ 10^{-7}$ | 2.32 | $2.6$ | 1.55 | 2.60 | 25 | $1.55$ | 0.13 | 2.75 | 1.0 | $0.85$ | 3.20 | $2.0\leavevmode\nobreak\ 10^{-7}$ | $1.3$

+ Fermi/LAT $\gamma$-ray data from T+17 [1] and A+16 [2].

##  5 Neutrino emission

With an apparent dominant $\pi^{0}$-decay origin of the $\gamma$ rays produced
in the interstellar medium of both MCs it is clear that $\pi^{\pm}$-decay
neutrinos are also produced. Our calculations (following Kelner et al. 2006)
of the predicted muon- and electron-neutrino spectra of both galaxies indicate
(Fig. 7) that the broadly-peaked ($\sim$0.1–10 GeV) neutrino flux is too low
for detection by current and upcoming $\nu$ projects. This conclusion is based
on the estimated observation time needed to detect the LMC with an experiment
with detection sensitivity comparable to the Antarctica-based IceCube+DeepCore
Observatory, the most sensitive current/planned $\nu$-detector at neutrino
energies $10<E_{\nu}/{\rm GeV}<100$ (e.g. Bartos et al. 2013). The latter’s
effective area (Abbasi et al. 2012) is $A_{\rm
eff}(E_{\nu})=40\left(\frac{E_{\nu}}{100\,{\rm GeV}}\right)^{2}$ cm2 (for
$\nu_{\mu}$, twice smaller for $\nu_{e}$) in this energy range (Bartos et al.
2013). Only in the narrow energy range 10-50 GeV do the IceCube+DeepCore
sensitivity and the LMC predicted diffuse spectral $\nu$-flux effectively
overlap (cf. the T+17 LAT dataset). In this band the latter can be
approximated as $\frac{dN_{\nu}}{dE_{\nu}}\sim
10^{-10}\left(\frac{E_{\nu}}{100\,{\rm GeV}}\right)^{-2.5}{\rm cm}^{-2}{\rm
s}^{-1}{\rm GeV}^{-1}$. The corresponding number of detected neutrinos,
$N_{\nu}=t_{\rm obs}\int_{10\,{\rm GeV}}^{50\,{\rm
GeV}}\frac{dN_{\nu}}{dE_{\nu}}A_{\rm eff}(E_{\nu})dE_{\nu}$, is $N_{\nu}\sim
10(t_{\rm obs}/{\rm yr})$. The detector background (for up-going events) is
dominated by atmospheric neutrinos produced by cosmic rays in the northern
hemisphere; its energy spectrum is approximately flat in the relevant energy
range, at a level $\sim$$10^{2}(d\Omega_{\rm LMC}/1.6\,{\rm sr})$ Gev-1 yr-1
(Bartos et al. 2013 and references therein). So the net 10-50 GeV background
rate is $\sim$$4\times 10^{3}(\Delta\Omega_{\rm LMC}/1.6\,{\rm sr})$
yr${}^{-1}\sim 60$ yr-1 assuming the LMC angular radius to be 5 deg). Based on
this crude estimate, observation of diffuse GeV neutrinos from the LMC disk
would imply $S/N<0.2$.

Figure 7:  Neutrino spectral emissions in the Magellanic Clouds. The LMC panel
refers to the T+17 Fermi/LAT dataset.

##  6 Conclusion

The SMC and the LMC, star-forming Galactic satellite galaxies, are among the
brightest sources in the Fermi/LAT $\gamma$-ray sky. We self-consistently
modeled the radio/$\gamma$ SED of both galaxies using the latest available
radio and LAT data, using exact emissivity formulae for the relevant emission
processes. Both SEDs were modeled with the radio data interpreted as a
combination of NT electron synchrotron emission and thermal electron
bremsstrahlung, and the $\gamma$-ray data as a combination of $\pi^{0}$-decay
emission from CRp interacting with the ambient gas plus leptonic emission. The
CRp spectra appear similar in the two galaxies with $q_{p}\sim 2.4$, and
$u_{p}\sim 1$ eV cm-3. Our spectral findings are qualitatively and
quantitatively new for the SMC, and confirm and strengthen previous results
for the LMC. In detail, we quantify Lopez et al.’s (2018) suggestion of a
mostly pionic origin of the LAT-measured $\gamma$-ray emission of the SMC,
only $\sim$$<$ 3% of the peak flux being by NT bremsstrahlung; as to the LMC,
our analysis suggests that 93%, 6.5%, and 0.5% of the LMC peak $\gamma$-ray
flux (at 0.5 GeV) are accounted for by pionic, Compton/(EBL+GFL), and NT-
bremsstrahlung emission – in broad agreement with Foreman et al. (2015). As we
have noted, for the SMC disk we used Lopez et al.’s (2018) LAT dataset which
is only mildly contaminated by an estimated $\sim$10% contribution from
unresolved pulsars, whereas for the LMC disk we used LAT datasets that are
based on maps that do not include emission from local gas and CR
inhomogeneities associated with actively star-forming regions (A+16; T+17).
Thus, our pionic emission modeling results are unlikely to be appreciably
affected by emission from individual sources.

As stated in Section 4, our treatment is based on determining the particle
steady-state spectral distributions in the disk by fits to the radio and
$\gamma$-ray data, using previously deduced mean disk values of the gas
density, magnetic and radiation fields. An assessment of the combined
uncertainties in the predicted particle radiative (and neutrino) yields is
based firstly on the precision level of the observationally determined values
of $n_{g}$ and $B$. Whereas the deduced CRe normalization is (essentially)
independent of $n_{g}$, its dependence on $B$ is significant, $N_{e0}\propto
B^{-(\alpha+1)}$, which is roughly $B^{-1.7}$. The deduced CRp normalization
is $N_{p0}\propto 1/n_{g}$. Given the substantial uncertainties in both
$n_{g}$ and $B$, it is clear that these impact also the overall normalization
of the particle spectra, which can be uncertain by a factor of 2-3. However,
the constraining power of the fits to the measurements is mostly in the well
known spectral shapes of the predicted synchrotron and $\pi^{0}$-decay
processes. As explained in the previous section, the spectral fits are quite
good, with a typical level of uncertainty of a few $10\%$. Even though
substantial, the overall level of uncertainty is not large enough to question
our qualitative conclusion on the pionic nature of the $\gamma$-ray emission.

The SMC radio data have large error bars below 230 MHz where $q_{e}$ would be
best determined, are unsampled (but for one point) between 230 and 1400 MHz
where the effect of the spectral truncation would be observed, and are only
sparsely sampled above 1400 MHz. This leaves both $q_{e}$ and $\gamma_{max}$
poorly determined. The LMC radio spectral data permit determination only of
$q_{e}$, as no truncation is apparent; we estimated $\gamma_{max}$ by weighing
its lower limit implied by the radio spectrum and its upper limit implied by
the $\gamma$ spectrum. Indeed, a larger value would shift the
Compton/(EBL+FGL) blue peak to higher frequencies and to a higher flux level,
deteriorating the fit to the $\gamma$ spectrum at frequencies log$(\,\nu)>24$.
Our joint analysis of the broad-band SED of the MCs using exact emissivity
formulae, and the leptonic radio and $\gamma$ emissions are coupled and rest
on independent measures (independent of particle/field energy equipartition
considerations) of the magnetic field, so the pionic emission is formally
determined by modeling the residuals of the LAT data after the leptonic yield
is accounted for. In general, uncertainties on the CRe spectrum affect the
determination of the CRp spectrum because the latter is obtained by fitting
the $\gamma$ data once the leptonic yields have been properly accounted for.
However, we feel our CRp results here are reasonably safe because of the clear
dominance of pionic emission in the $\gamma$ spectrum.

For the SMC, the LAT data are sufficient to define the CRp spectral
parameters, $q_{p}=2.4$ and $E_{p}^{max}=30$ GeV. For the LMC, however, some
inconsistency between the T+17 and A+16 LAT datasets implies the derived CRp
spectrum to be, respectively, flatter ($q_{p}=2.4$) with no clearly
discernible high-energy truncation ($E_{p}^{max}\raise-2.0pt\hbox{\hbox
to0.0pt{\hbox{$\sim$}\hss}\raise 5.0pt\hbox{$<$}\ }80$ GeV), or steeper
($q_{p}=2.6$) with a clear truncation ($E_{p}^{max}=25$ GeV) and a $\sim$50%
higher normalization. (Both $q_{p}$ values are consistent within errors with
Foreman et al.’s (2015) earlier estimate.) The reason for the discrepancy
between the two data sets is not obvious, as the T+17 and A+16 data extraction
regions largely overlap with each other and with the LMC disk, and both have
been cleaned for identified (point-like or extended) sources. It should be
pointed out that in the T+17 dataset the quality of the four highest-frequency
points does not allow a definite estimate of $E_{p}^{max}$, whereas the four
lowest-frequency points enable a clear view of the rising portion of the
pionic hump. Ongoing Fermi/LAT measurements would likely improve the quality
of the spectrum and clarify this issue. For both galaxies we obtained
truncation energies systematically higher for CRp than for CRe – probably,
because energy losses are much less efficient for the former than for the
latter in the relatively low-density MC interstellar gas.

The secondary CRe spectra are at best as well determined as those of the
parent CRp. Our model secondary spectra (with no free parameters) are
relatively well defined in the MCs. The primary CRe spectra, determined by
fitting radio data after secondary CRe have been accounted for, are however
less well constrained. First, the measured radio fluxes include emission from
individual background and MC-disk point sources that contribute $\sim$20% of
the total emission and are spectrally similar to the overall emission (For et
al. 2018), so there is some intrinsic uncertainty in the predicted extended
emission’s spectral shape and normalization ($N_{e0}$ is biased high).
Secondly, because of the $q_{e}$-$N_{e0}$ degeneracy the uncertainty on
$N_{e0}$ affects the determination of the spectral slope. Thus, primary CRe
are rather poorly constrained in both MCs. In spite of this, our assumption of
a (truncated) 1PL representing primary CRe is sound. The PL is a fair
approximation, over the relevant CRe energy range, to the actual spectrum
computed accounting for realistic energy losses in the MC disk environments.
The approximation is particularly good for the LMC.

A possible underestimate of the gas content in the MC due to the presence of
”dark” (unseen) neutral gas, i.e. HI gas optically thick in the 21 cm line
and/or H2 has with no associated CO emission, could increase the uncertainty
on the total, atomic and molecular, gas density. The CO (J=1-0) transition,
often used to detect the total H2 mass, may be hard to detect in low-$Z$
galaxies such as the MC (Rolleston et al. 2002; Requena-Torres et al. 2016),
so the relatively bright [C II]$\lambda$158 $\mu$m line can be used instead:
some $>$70% may not have been traced by CO (1-0) in dwarf galaxies but is well
traced by [C II]$\lambda$158 $\mu$m (Madden et al. 2020). Another important
aspect to investigate is the mixing of dust into the interstellar medium (ISM)
and the spatial variations of their properties (gas clumpiness affects
emission levels, see above). Comparing the distributions of the IR dust
emission and several gas tracers (H$\alpha$, 21 cm, CO emission), their
different emission processes highlight the distribution of gas under different
conditions. In the LMC, this analysis reveals that: (i) dust emission, sampled
in the Multiband Imaging Photometer for Spitzer (MIPS) 70 and 160 $\mu$m
bands, is well mixed with the large-scale HI 21 cm emission; (ii) H$\alpha$
from star-forming H II regions is confined to the massive star formation
region where also warmer ($\sim$120 K) 24 $\mu$m emission is found; and (iii)
the Spitzer InfraRed Array Camera (IRAC) 8 $\mu$m band, that traces polycyclic
aromatic hydrocarbons, correlates well with the HI gas but is absent from the
massive-star–forming H II regions. All in all, the dust emission revealed by
the combined IRAC 8 $\mu$m and MIPS 24, 70 and 160 $\mu$m bands traces all
three phases of the ISM gas (Meixner et a. 2006). Underestimating the gas
content would affect the particle spectra in the following ways: (i)
overestimate the CRp spectral normalization; (ii) underestimate the Coulomb
and bremsstrahlung energy losses, $b_{0}(\gamma)$ and $b_{1}(\gamma)$
(Rephaeli & Persic 2015), which would cause the computed steady-state
secondary-CRe spectrum to be, respectively, higher and steeper; (iii) bias the
derived primary-CRe spectrum lower and flatter.

Uncertainties in the 3D structure of the MCs (mainly for the SMC, Abdo et al.
2010b) directly affect only the CRe and CRp spectral normalizations, to match
the spectral flux when the emitting volume changes – this holds for the
relevant leptonic emissivities, which depend linearly on the CRe
normalization. Thus, model SEDs are nearly unaffected by variations of the MC
structure parameters. Another source of uncertainty in $N_{e0}$ stems from
assuming a magnetic field (for both galaxies), in our case field values
derived by Faraday rotation studies of extragalactic polarized sources seen
through the MCs. These field strengths are measured along the line of sight to
the background sources; thus, they are line-averaged fields, whereas the field
in the synchrotron emissivity expression is a volume average. The two values
may or may not be the same: the potential mismatch introduces another
systematic uncertainty in the determination of $N_{e0}$ from fitting the radio
spectrum. For example, in principle if $B$ is biased high, $N_{e0}$ and all
leptonic yields are biased low, so the resulting $N_{p0}$ is biased high.

The main breakthrough needed in MC SED studies is the measurement of diffuse
NT X-ray emission from their disks. As exemplified in our studies of radio-
lobe SEDs (Persic & Rephaeli 2019a,b), the NT 1 keV flux, interpreted as
Compton/CMB emission, sets the normalization of the CRe spectrum: in the MC
case, given that the excellent pionic fit to the $\gamma$-ray emission firmly
defines the secondary spectrum, the 1 keV flux would effectively measure the
primary normalization. Modeling the radio spectrum as synchrotron radiation
would then provide the primary spectral shape – and, importantly, the magnetic
field. A spectral determination of $B$ would bypass the need to assume
particles/field equipartition, and would return a volume-averaged value better
representing the mean field than the line averages yielded by Faraday-
rotation–based measurements. Deeper X-ray observations of both MCs, with
better spatial/spectral resolution than currently available, will be highly
beneficial to searching the diffuse disk emission for a NT component.

Acknowledgement. We acknowledge insightful comments by an anonymous referee
that improved clarity and presentation of our work.

## References

  * (1) Abbasi R., Abdou Y., Abu Zayyad T., et al. (IceCube Collaboration), 2012, Astropart Phys, 35, 615
  * (2) Abdo A.A., Ackermann M., Ajello M., et al., 2010a, A&A, 512, A7
  * (3) Abdo A.A., Ackermann M., Ajello M., et al., 2010b, A&A, 523, A46
  * (4) Ackermann M., Albert A., Atwood W.B., et al., 2016, A&A, 586, A71 [A+16]
  * (5) Alvarez H., Aparici J., May J., 1987, A&A, 176, 25
  * (6) Bartos I., Beloborodov A.M., Hurley K., Márka S. 2013, Phys Rev Lett, 110, 241101
  * (7) Blumenthal G.R., Gould R.J., 1970, RvMP, 42, 237
  * (8) Condon J.J., 1992, ARAA, 30, 575
  * (9) Cooray A., 2016, R. Soc. Open Sci., 3, 150555
  * (10) Cox N.L.J., Cordiner M.A., Cami J., et al., 2006, A&A, 447, 991
  * (11) For B.-Q., Staveley-Smith L., Hurley-Walker N., et al., 2018, MNRAS, 480, 2743
  * (12) Foreman G., Chu Y.H., Gruendl L., et al. 2015, ApJ, 808:44
  * (13) Franceschini A., Rodighiero G., 2017, A&A, 603, A 34 (and Erratum: 2018, A&A, 614, C 1)
  * (14) Fukui Y., Kawamura K., Minamidami T., et al. 2008, ApJS, 178, 56
  * (15) Gaensler B.M., Haverkorn M., Staveley-Smith L., et al., 2005, Science, 307, 1610
  * (16) Haynes R.F., Klein U., Wayte S.R., et al., 1991, A&A, 252, 475
  * (17) Kelner S.R. et al., 2006, Phys. Rev. D, 74, 034018 [and Erratum in 2009, Phys. Rev. D, 79, 039901(E)]
  * (18) Kennicutt R.J. Jr., Bresolin F., Bomans D.J., Bothun G.D., Thompson I.B., 1995, AJ, 109, 594
  * (19) Klein U., Wielebinski R., Haynes R.F., Malin D.F., 1989, A&A, 211, 280
  * (20) Loiseau N., Klein U.,Greybe A., Wielebinski R., Haynes R.F., 1987, A&A, 178, 62
  * (21) Lopez L.A., Auchettl K., Linden T., et al., 2018, ApJ, 867:44
  * (22) Madden S.C., Cormier D., Hony S., et al. 2020, A&A, 643, A 141
  * (23) Mao S.A., Gaensler B.M., Stanimirovic S., et al., 2008, ApJ, 688, 1029
  * (24) Mao S.A., McClure-Griffiths N.M., Gaensler B.M., et al. 2012, ApJ, 759:25
  * (25) Meixner M., Gordon K.D., Indebetouw R., et al. 2006, ApJS, 132, 2268
  * (26) Mills B.Y., 1955, Aust. J. Phys., 8, 368
  * (27) Mills B.Y., 1959, Handbuch der Physik, 53, 239
  * (28) Mountfort P.I., Jonas J.L., De Jager G., Baart E.E., 1987, MNRAS, 226, 917
  * (29) Nishiuchi M., Imanishi K., Tsujimoto M., Yokogawa J., Koyama K., 1999, in ”Star Formation 1999” (ed. T. Nakamoto, Nobeyama Radio Observatory), p.54
  * (30) Nishiuchi M., Yokogawa J., Koyama K., 2001, in ”New Century of X-ray Astronomy” (ASP Conference Series, Vol. 251; eds. H. Inoue, H. Kunieda), p.264
  * (31) Persic M., Rephaeli Y., 2019a, MNRAS, 485, 2001 (and Erratum: 2019, MNRAS, 486, 950)
  * (32) Persic M., Rephaeli Y., 2019b, MNRAS, 490, 1489
  * (33) Points S.D., Chu Y.H., Snowden S.K., Smith R.C., 2001, ApJS, 136, 99
  * (34) Rephaeli Y., Persic M., 2015, ApJ, 805:111
  * (35) Rephaeli Y., Sadeh S., 2019, MNRAS, 486, 2496
  * (36) Requena-Torres M.A., Israel F.P., Okada Y., et al. 2016, A&A, 589, A 28
  * (37) Rolleston W.R.J., Trundle C., Dufton P.L. 2002, A&A, 396, 53
  * (38) Sasaki M., Haberl F., Pietsch W., 2002, A&A, 392, 103
  * (39) Shain C.A., 1959, in IAU Symp. 9, Observations of Extragalactic Radio Emission (Bracewell R.N., ed.; Stanford Univ. Press, Stanford, CA), p.328
  * (40) Spitzer L., Jr, 1978, Physical Processes in the Interstellar Medium (Wiley)
  * (41) Sreekumar P., Bertsch D.L., Dingus B.L., et al., 1992, ApJ, 400, L67
  * (42) Stanimirović S., Staveley-Smith L., Dickey J.M., Sault R.J., Snowden S.L. 1999, MNRAS, 302, 417
  * (43) Staveley-Smith L., Kim S., Calabretta M.R., Haynes R.F., Kesteven M.J., 2003, MNRAS, 339, 87
  * (44) Tang Q.W., Peng F.K., Liu R.Y., Tam P.H.T., Wang X.Y., 2017, ApJ, 843:42 [T+17]
  * (45) Toribio San Cipriano L., Domínguez-Guzmán G., Esteban C., et al., 2017, MNRAS, 467, 3759
  * (46) Wang Q., Hamilton T., Helfand D.J., Wu X., 1999, ApJ, 374, 475


# Deploying a sharded MongoDB cluster as a queued job on a shared HPC
architecture

Aaron Saxton<EMAIL_ADDRESS>National Center For Super Computing
Applications
Blue Waters Project OfficeUrbanaIllinois and  Stephen Squaire
<EMAIL_ADDRESS>National Center For Super Computing Applications
Integrated Cyber InfrastructureUrbanaIllinois

###### Abstract.

Data stores are the foundation on which data science, in all its variations,
is built upon. They provide a queryable interface to structured and
unstructured data. Data science often starts by leveraging these query
features to perform initial data preparation. However, most data stores are
designed to run continuously to service disparate user requests with little or
no downtime. Many HPC architectures process user requests by job queue
scheduler and maintain a shard filesystem to store a jobs persistent data. We
deploy a MongoDB sharded cluster with a run script that is designed to run a
data science workload concurrently. As our test piece, we run data ingest and
data queries to measure the performance with different configurations on the
Blue Waters supper computer.

distributed datastore, Mongodb, high performentce computing, shard filesystem

††ccs: Information systems Data management systems

##  1\. Introduction

Queryable data stores have been a crucial part of modern data science. They
provide a relatively small set of primitive operations by which a large
variety of data analysis can be performed. Machine learning and other
sophisticated statistical analysis are becoming more popular across all
disciplines (Wilson et al., 2016). These algorithms require more flexibility
than the average SQL or NoSQL data store and can take advantage of a wider
variety of compute resources than the average datastore host architecture.
However, it is common for a data scientist to computationally balance their
process between a datastore and compute resource. Data and algorithms are
becoming large enough and computationally intensive enough that they can take
advantage of a high performance architecture. The 2018 Gordon Bell prize
”Employing Deep Learning Methods to Understand Weather Patterns” is an
example. In this paper we present a performance profile of the NoSQL
distributed data store MongoDB running as transient queued batch job on Blue
Waters super computer. This will provide guidance about how to balance tasks
in a data science pipelines on a HPC architectures.

##  2\. Background

The major primitive operations of a traditional SQL data store are CRUD
(Create, Read, Update, Delete), groupby, and join. NoSQL data stores, such as
MongoDB, still maintain a core set of primitives that are analogous to
traditional SQL: insert, find, update, remove, aggregate, lookup. This work
will not examine all of these primitives, but we still wish to speak of them.

Finds and inserts naturally lend themselves to distributed computing. A
program can dispatch them to many workers witch report back either nothing is
found or the results of their task. The more sophisticated aggregate and
lookup can not in general be executed without some synchronization or
interprocess communication. Take for example a SQL inner join on column $C$
between tables $T_{A}$ and $T_{B}$. The naive algorithm goes as: Take first
item in $C$ of table $T_{A}$, $C_{A0}$, perform a find on column $C$ in table
$T_{B}$. In a distributed datastore, this will require each worker to look at
$C_{AN}$ of $T_{A}$ hosted locally, search through the portion of $T_{B}$
hosted locally and pass $C_{AN}$ to every other worker. If the columns are
indexed, there is some optimization you can do, but every worker will have to
touch every other work at some point. Groupby haves similar problems. A map-
reduce execution model can provide a single stage inter-worker communication
to perform a groupby. The algorithms for performing joins and groupbys on a
distributed data store are rebranded to give the user more flexibility. As
such, this paper will only be concerned with insert and find type operations.

##  3\. System Setup

Blue Waters is a mixture of Cray XE and XK blades with 4 Interlogos AMD
processors XE nodes and two Interlogs AMD, K20 Nvidia GPU pairs on XK nodes.
These nodes are all connected by Cray’s Gemini interconnect. Shared filesystem
storage is hosted on Cray’s Synexion racks with a luster file system. Blue
Waters has 4 login nodes that users can use to access the shared file system
and submit run scripts for Moab and Torque schedule on the requested compute
resources. (Bode et al., 2013) (Kramer et al., 2015)

Data stores are often intended to be persistent and service tasks that demand
zero downtime. For instance, users data on the backend of a web app, hospital
patient data, or product order data from a retailer. Compute resources on an
HPC system that are available to a typical user are often ephemeral and are
contrary to the execution model of common data stores. Fortunately, the final
destination of a data stores data is on an underlying filesystem. MongoDB by
default uses the WiredTiger storage engine to manage data underlying files on
the host filesystem. For Blue Waters, that is a massively parallel and
redundant file system(Bode et al., 2013) (Kramer et al., 2015).

###  3.1. MongoDB Basics

A sharded MongoDB cluster runs with 3 types of worker processes that we host
one on each processing element of a Blue Waters job: configuration server,
shard server, and router. The configuration server manages the metadata of the
collections (Mon, [n. d.]d). Config servers store the metadata for a sharded
cluster. The metadata reflects state and organization for all data and
components within the sharded cluster. The metadata includes the list of
chunks on every shard and the ranges that define the chunks(Mon, [n. d.]a). A
shard contains a subset of sharded data for a sharded cluster. Together, the
cluster?s shards hold the entire data set for the cluster(Mon, [n. d.]c). Each
shard will be assigned a unique filesystem path to deposit data too. MongoDB
routers route queries and write operations to shards in a sharded cluster. A
router provide the only interface to a sharded cluster from the perspective of
applications. Applications never connect or communicate directly with the
shards(Mon, [n. d.]b).

Figure 1. Sharded MongoDB cluster diagram (Mon, [n. d.]d)

###  3.2. User Execution Model

For a typical Blue Waters user to deploy a MongoDB cluster, they must
construct a run-script that assigns to each processing element which roll it
will be taking (config, shard, router). MongoDB is nativly deployed on a
TCP/IP network therefore worker reference each other by hostname and port. In
our run script we chose the number of processing elements to be equal number
of node hosts. Once the assignments have been made, the run scrip will
configure each worker to reference the config server. The config server will
then share the settings with all members of the cluster. In addition to taking
advantage of distributed compute resources, a sharded MongoDB cluster will
natively leverage the distributed nature of Luster on a Cray Synexian. When
each shard worker is assigned a directory to place files, luster will
distribute those files to an object storage server that should optimize
further I/O.

This describes a run scripts execution for the datastore portion. In practice,
once the datastore is active, a run script will continue making queries and
processing. The runsript makes available through environment variables or a
shared file a list of host names of the MongoBD clusters router servers. With
this list, a run script may use either the mongo shell command or the Python
package pymongo as the API to perform queries.

##  4\. Performance Data

We perform 2 tasks, ingest data with insertMany with ordered=False and a
conditional find on two indexed fields. These tasks are written in python
using pymongo. The data we choose to ingest was time series metric data of
Blue Waters compute nodes collected by OVIS. The data spans 5 years, sample
each node independently once every minute, and includes about 75 distinct
metrics (e.g. memory use, cpu activity, network activity). The totality of
this data is about 70 billion rows. Storage for this data in flat csv file on
Blue Waters Luster filesystem is about 200 terabytes. We index on timestamp
and node id.

We perform scaling analysis by increasing the number of shard-router pairs.
For example, a job of 32 nodes is scheduled. 2 nodes will be for the
configuration server, 7 shards, and 7 routers. This leaves 16 nodes to run the
ingest script. Ingest is run with 4 processing elements per node, thus 64
insertMany will be processed concurrently across 7 MongoDB routers. A job of
64 nodes would have 2 for configuration, 15 shards, 15 router servers and so
on. The larger the cluster, the more data we upload.

Table 1.  Nodes | Days of Data
---|---
32 | 3
64 | 7
128 | 14
256 | 14

A instertMany is performed by collecting a list of python dictionaries from
the metric data csv file. The keys of the dictionaries become analogous to SQL
table column names. A list of dictionaries is allowed to have different keys.
However, ingesting this dataset allows the same keys for each document.
Indexes were created for the timestamp key and node id key. In Figure 2
MongoDB scales close to linear between 32, 64, and 128 nodes. We are still
investigating the limitations at 256 nodes and beyond.

Figure 2.

The query test was done by doing a conditional find. The query is constructed
by reading user jobs metadata for time run, duration, and which nodes were
assigned. The candidate user jobs were selected from a time period starting
January 1, 2018 until the number of days described in Table 1. The total
number of documents returned by a query is number of user job nodes times
duration of user job in minutes. Indeed, Figure 3 indicates cluster size
maintains a similar query performance for various MongoDB cluster sizes. It is
important to point out that each cluster size is servicing more concurrent
quarries. Cluster with size 32 was servicing between 16 and 64 concurrent find
queries, cluster size 64 was between 32 and 128 queries, and so on.

Figure 3. Sharded MongoDB cluster diagram (Mon, [n. d.]d)

###### Acknowledgements.

This research is part of the Blue Waters sustained-petascale computing
project, which is supported by the National Science Foundation (awards
OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint
effort of the University of Illinois at Urbana-Champaign and its National
Center for Supercomputing Applications.

## References

  * (1)
  * Mon ([n. d.]a) [n. d.]a.  MongoDB Manual, Config Server.  Retrieved April 9, 2019 from https://docs.mongodb.com/manual/core/sharded-cluster-config-servers/
  * Mon ([n. d.]b) [n. d.]b.  MongoDB Manual, Mongos.  Retrieved April 9, 2019 from https://docs.mongodb.com/manual/core/sharded-cluster-query-router/
  * Mon ([n. d.]c) [n. d.]c.  MongoDB Manual, Shard.  Retrieved April 9, 2019 from https://docs.mongodb.com/manual/core/sharded-cluster-shards/
  * Mon ([n. d.]d) [n. d.]d.  MongoDB Manual, Sharding.  Retrieved April 9, 2019 from https://docs.mongodb.com/manual/sharding/
  * Bode et al. (2013) Brett Bode, Michelle Butler, Thom Dunning, Torsten Hoefler, William Kramer, William Gropp, and Wen-mei Hwu. 2013.  The Blue Waters Super-System for Super-Science.  In _Contemporary High Performance Computing_. Chapman and Hall/CRC, 339–366.  https://www.taylorfrancis.com/books/e/9781466568358
  * Kramer et al. (2015) William Kramer, Michelle Butler, Gregory Bauer, Kalyana Chadalavada, and Celso Mendes. 2015\.  Blue Waters Parallel I/O Storage Sub-system.  In _High Performance Parallel I/O_ , Prabhat and Quincey Koziol (Eds.). CRC Publications, Taylor and Francis Group, 17–32.
  * Wilson et al. (2016) H. James Wilson, Allan Alter, and Prashant Shukla. 2016\.  Companies Are Reimagining Business Processes with Algorithms.  Retrieved April 9, 2019 from https://hbr.org/2016/02/companies-are-reimagining-business-processes-with-algorithms


with $B$ exactly $a-1$ times. The second term corrects for this over-counting,
akin to an inclusion-exclusion argument. The following lemmas bound each of
the two terms:

######  Lemma A.12.

For all $i\in B$, we have

$\displaystyle\frac{1}{M_{i}}\left|\sum_{S\in\mathcal{S}_{i},i^{*}\notin
S}\left(g(X_{S})-g(X^{\prime}_{S})\right)\right|\leq\left(\xi+\frac{20kCL}{n}\right)\epsilon+\frac{12kC}{n}+\frac{14k^{2}CL}{n^{2}}.$

###### Proof.

We have

$\displaystyle\sum_{\begin{subarray}{c}S\in\mathcal{S}_{i}\\\ i^{*}\notin
S\end{subarray}}\left(g(X_{S})-g(X^{\prime}_{S})\right)=M_{i}\left(\hat{g}_{X}({i})-\hat{g}_{X^{\prime}}({i})\right)-\sum_{\begin{subarray}{c}S\in\mathcal{S}_{i,i^{*}}\end{subarray}}\left(g(X_{S})-g(X^{\prime}_{S})\right).$
(A.59)

By Lemmas A.10 and A.11, we have

$\displaystyle\left|\hat{g}_{X}({i})-\hat{g}_{X^{\prime}}({i})\right|\leq\left(\xi+\frac{kC(14+6L)}{n}\right)\epsilon+\frac{12kC}{n}+\frac{4k^{2}C}{n^{2}}(1+2L).$

Moreover, the second term in equation (A.59) is clearly upper-bounded by

$2CM_{i,i^{*}}\leq\frac{4kC}{n}M_{i}\leq\frac{8k^{2}C}{n^{2}}M_{i}.$

Summing these two bounds yields the result. ∎

For the second term of equation (A.58), we use the following lemma:

######  Lemma A.13.

We have

$\displaystyle\sum_{a=2}^{\min(k,|B|)}\sum_{\begin{subarray}{c}S\in\mathcal{S}\\\
|S\cap B|=a\\\ i^{*}\notin
S\end{subarray}}(a-1)\left|g(X_{S})-g(X^{\prime}_{S})\right|\leq\frac{36k^{2}L^{2}}{n^{2}}\left(2C\epsilon+\frac{6kC}{n}\right)\min(k,2L)M.$

Moreover, if $\mathcal{S}=\mathcal{I}_{n,k}$, with $m=\binom{n}{k}$, we have
the stronger inequality

$\displaystyle\sum_{a=2}^{\min(k,|B|)}\sum_{\begin{subarray}{c}S\in\mathcal{S}\\\
|S\cap B|=a\\\ i^{*}\notin
S\end{subarray}}(a-1)\left|g(X_{S})-g(X^{\prime}_{S})\right|\leq\frac{9k^{2}L^{2}}{n^{2}}\left(2C\epsilon+\frac{6kC}{n}\right)m.$

###### Proof.

For any $S$ not containing $i^{*}$, we have

$\displaystyle\left|g(X_{S})-g(X^{\prime}_{S})\right|$
$\displaystyle=\left|(h(X_{S})-A_{n})\text{wt}(S)+A_{n}-(h(X^{\prime}_{S})-A_{n}^{\prime})\text{wt}^{\prime}(S)-A_{n}^{\prime}\right|$
$\displaystyle\leq\left|(h(X_{S})-A_{n})(\text{wt}(S)-\text{wt}^{\prime}(S))\right|+\left|(h(X_{S})-h(X^{\prime}_{S}))\text{wt}^{\prime}(S)\right|+\left|A_{n}-A_{n}^{\prime}\right|$
$\displaystyle\leq 2C\epsilon+\frac{6kC}{n},$

using Lemma A.10 and the fact that the second term is zero. Moreover, we have

$\displaystyle\sum_{a=2}^{\min(k,|B|)}\sum_{\begin{subarray}{c}S\in\mathcal{S}\\\
|S\cap B|=a\\\ i^{*}\notin S\end{subarray}}(a-1)$
$\displaystyle\leq\sum_{a=2}^{k}\sum_{\begin{subarray}{c}S\in\mathcal{S}\\\
|S\cap
B|=a\end{subarray}}\min(k,|B|)=\sum_{\begin{subarray}{c}S\in\mathcal{S}\\\
|S\cap B|\geq 2\end{subarray}}\min(k,|B|)$
$\displaystyle\leq\sum_{\begin{subarray}{c}i,j\in B\\\ i\neq
j\end{subarray}}\sum_{S\in\mathcal{S}_{ij}}\min(k,|B|)\leq\frac{9k^{2}}{n^{2}}\min(k,|B|)|B|^{2}m.$

The last inequality follows from Lemma A.9, which implies that
$\frac{M_{ij}}{M}=O\left(\frac{k^{2}}{n^{2}}\right)$.

In the case when $\mathcal{S}=\mathcal{I}_{n,k}$, we have

$\displaystyle\sum_{a=2}^{\min(k,|B|)}\sum_{\begin{subarray}{c}S\in\mathcal{S}\\\
|S\cap B|=a\\\ i^{*}\notin S\end{subarray}}(a-1)$
$\displaystyle\leq\sum_{a=2}^{\min(k,|B|)}(a-1)\binom{|B|}{a}\binom{n-|B|}{k-a}$
$\displaystyle=\binom{n}{k}\left(\frac{k|B|}{n}-1\right)+\binom{n-|B|}{k}$
$\displaystyle\leq\binom{n}{k}\left(\frac{k|B|}{n}-1\right)+\binom{n}{k}\left(1-\frac{|B|}{n}\right)^{k}$
$\displaystyle\leq\binom{n}{k}\left(\frac{k|B|}{n}-1\right)+\binom{n}{k}\frac{1}{\frac{k|B|}{n}+1}\leq\frac{k^{2}|B|^{2}}{n^{2}}\binom{n}{k},$

where the first equality used the identities
$\sum_{a=0}^{k}\binom{n}{a}\binom{m}{k-a}=\binom{n+m}{k}$ and
$\sum_{a=0}^{k}a\binom{n}{a}\binom{m}{k-a}=\frac{nk}{m+n}\binom{n+m}{k}$, and
the third inequality used the fact that $(1-x)^{k}\leq\frac{1}{1+kx}$ for all
$x\in[0,1]$. The statement in the lemma follows because $|B|\leq 2L+1\leq 3L$.
∎

Combining the results of Lemma A.12 and A.13 yields Lemma A.14.

######  Lemma A.14.

We have

$\left|T_{1}\right|\leq\frac{kM}{n}\left(\left(\xi+\frac{20kCL}{n}\right)2\epsilon+\frac{24kC}{n}+\frac{28k^{2}CL}{n^{2}}+\frac{9kL^{2}}{n}\left(2C\epsilon+\frac{6kC}{n}\right)\gamma\right),$

where $\gamma=4\min(k,2L)$. If $\mathcal{S}=\mathcal{I}_{n,k}$, then the bound
also holds for $\gamma=1$.

It remains to bound the third term, $T_{3}$, of equation (A.4.1), which we do
in the following lemma:

######  Lemma A.15.

We have

$\left|T_{3}\right|\leq\frac{2k}{n}\left(2\xi+\frac{kC(11+18L)}{n}+\frac{12kC}{n\epsilon}\right)M.$

###### Proof.

Using Lemmas A.10 and A.11, we have

$\displaystyle\frac{1}{M_{i^{*}}}\left|T_{3}\right|$
$\displaystyle=\left|\hat{g}_{X}({i^{*}})-\hat{g}_{X^{\prime}}({i^{*}})\right|\leq\left|\hat{g}_{X}({i^{*}})-A_{n}\right|+\left|\hat{g}_{X^{\prime}}({i^{*}})-A_{n}^{\prime}\right|+\left|A_{n}-A_{n}^{\prime}\right|$
$\displaystyle\leq\left(\xi+\frac{9kCL}{n}+\frac{6kC}{n\epsilon}\right)+\left(\xi+\frac{9kCL^{\prime}}{n}+\frac{6kC}{n\epsilon}\right)+\frac{2kC}{n}$
$\displaystyle\leq 2\xi+\frac{kC(11+18L)}{n}+\frac{12kC}{n\epsilon}.$

The lemma follows after using the fact that $M_{i^{*}}\leq\frac{2k}{n}M$. ∎

Combining the bounds on $T_{1}$ and $T_{3}$ from Lemmas A.14 and A.15 in
equation (A.4.1), the local sensitivity of $\tilde{A}_{n}$ at $\mathbf{X}$ is
then bounded as

$\displaystyle
LS_{\tilde{A}_{n}}\left(\mathbf{X}\right)=O\left(\frac{k}{n}\left(\xi+\frac{kCL}{n}\right)\left(1+\epsilon
L\right)+\frac{k^{2}CL^{2}\min(k,L)}{n^{2}}\left(\epsilon+\frac{kC}{n}\right)+\frac{k^{2}C}{n^{2}\epsilon}\right).$

Let $g(\xi,L,n)$ denote the upper bound, where to simplify the following
argument, we assume the constant prefactor is $1$, i.e.,

$\displaystyle
g(\xi,L,n):=\frac{k}{n}\left(\xi+\frac{kCL}{n}\right)\left(1+\epsilon
L\right)+\frac{k^{2}CL^{2}\min(k,L)}{n^{2}}\left(\epsilon+\frac{k}{n}\right)+\frac{k^{2}C}{n^{2}\epsilon}.$

Note that $g$ is strictly increasing in $L$. Also define

$\displaystyle S(\mathbf{X})=\max_{\ell\in\mathbb{Z}_{\geq
0}}e^{-\epsilon\ell}g(\xi,L_{\mathbf{X}}+\ell,n).$ (A.60)

######  Lemma A.16.

The function $S(\mathbf{X})$ is an $\epsilon$-smooth upper bound on
$LS_{\tilde{A}_{n}}(\mathbf{X})$.

$S(\mathbf{X})=O\left(\frac{k}{n}\left(\xi+\frac{kC(1/\epsilon+L)}{n}\right)\left(1+\epsilon
L\right)+\frac{k^{2}C(1/\epsilon+L)^{2}\min(k,1/\epsilon+L)}{n^{2}}\left(\epsilon+\frac{k}{n}\right)+\frac{k^{2}C}{n^{2}\epsilon}\right).$

###### Proof of Lemma A.16.

Clearly, we have $S(\mathbf{X})\geq g(\xi,L_{\mathbf{X}},n)\geq
LS_{\tilde{A}_{n}}(\mathbf{X})$, and for any two adjacent $\mathbf{X}$ and
$\mathbf{X}^{\prime}$, we have

$\displaystyle S(\mathbf{X}^{\prime})$
$\displaystyle=\max_{\ell\in\mathbb{Z}_{\geq
0}}e^{-\epsilon\ell}g(\xi,L_{\mathbf{X}^{\prime}}+\ell,n)\leq\max_{\ell\in\mathbb{Z}_{\geq
0}}e^{-\epsilon\ell}g(\xi,L_{\mathbf{X}}+\ell+1,n)$
$\displaystyle=\max_{\ell\in\mathbb{Z}_{>0}}e^{-\epsilon\left(\ell-1\right)}g(\xi,L_{\mathbf{X}}+\ell,n)\leq
e^{\epsilon}\max_{\ell\in\mathbb{Z}_{\geq
0}}e^{-\epsilon\ell}g(\xi,L_{\mathbf{X}}+\ell,n)=e^{\epsilon}S(\mathbf{X}).$

This shows that $S$ is indeed a $\epsilon$-smooth upper bound on the local
sensitivity. As for the upper bound on $S$, for any $\ell\geq 0$, we have

$\displaystyle
e^{-\epsilon\ell}g(\xi,L_{\mathbf{X}}+\ell,n)=\frac{k}{n}\left(\xi
e^{-\epsilon\ell/2}+\frac{kC(\ell
e^{-\epsilon\ell/2}+Le^{-\epsilon\ell/2})}{n}\right)\left(e^{-\epsilon\ell/2}+\epsilon(\ell
e^{-\epsilon\ell/2}+Le^{-\epsilon\ell/2})\right)$
$\displaystyle+\frac{k^{2}}{n^{2}}C(\ell
e^{-\epsilon\ell/3}+Le^{-\epsilon\ell/3})^{2}\min(ke^{-\epsilon\ell/3},\ell
e^{-\epsilon\ell/3}+Le^{-\epsilon\ell/3})\left(\epsilon+\frac{k}{n}\right)+\frac{k^{2}Ce^{-\epsilon\ell}}{n^{2}\epsilon}$
$\displaystyle\leq\frac{k}{n}\left(\xi+\frac{kC(1/\epsilon+L)}{n}\right)\left(1+\epsilon(1+L)\right)+\frac{k^{2}}{n^{2}}C\left(\frac{2}{\epsilon}+L\right)^{2}\min\left(k,\frac{2}{\epsilon}+L\right)\left(\epsilon+\frac{k}{n}\right)+\frac{k^{2}C}{n^{2}\epsilon},$

where we used in multiple places the inequalities $e^{-\epsilon\ell}\leq 1$
and $\ell e^{-\ell/c}\leq c/e$, for any $c>0$. ∎

By Lemma A.16, it is clear that the term $S(\mathbf{X})$ added to
$\tilde{A}_{n}$ in Algorithm 3 is the smoothed sensitivity defined in equation
(A.60).

Therefore, $\tilde{A}_{n}+\frac{S(\mathbf{X})}{\epsilon}\cdot Z$, where $Z$ is
sampled from the distribution with density $h(z)\propto 1/(1+|z|^{4})$, is
$O(\epsilon)$-differentially private, by Lemma 4. Moreover, if
$\mathcal{S}=\mathcal{I}_{n,k}$, the above bound on the smooth sensitivity
holds without the $\min(k,1/\epsilon+L)$ term, due to Lemma A.13.

#### Utility.

By Chebyshev’s inequality, we have

$\displaystyle\left|A_{n}-\theta\right|\leq\frac{1}{\sqrt{\gamma}}\sqrt{\text{var}(A_{n})}=\frac{1}{\sqrt{\gamma}}\left(\sqrt{\text{var}(U_{n})}+\sqrt{\frac{\zeta_{k}}{m}}\right),$

with probability at least $1-\gamma$. Moreover, with probability at least
$1-\gamma$, each of the Hájek projections is within $\xi$ of $A_{n}$. This
implies that every index $i$ has weight $1$, which further implies that
$g(X_{S})=h(X_{S})$ for all $S\in\mathcal{S}$, and consequently,
$\tilde{A}_{n}=A_{n}$. Also, $L=1$ for such an $\mathbf{X}$. Finally, with
probability at least $1-\gamma$, we have $Z\leq\frac{3}{\sqrt{\gamma}}$.
Combining these inequalities and using Lemma A.16, we have

$\displaystyle\left|\mathcal{A}(\mathbf{X})-\theta\right|\leq\left|\tilde{A}_{n}-A_{n}\right|+\left|A_{n}-\theta\right|+\left|S(\mathbf{X})/\epsilon\cdot
Z\right|$
$\displaystyle=\frac{1}{\sqrt{\gamma}}O\left(\sqrt{\text{var}(U_{n})}+\sqrt{\frac{\zeta_{k}}{m}}+\frac{k\xi}{n\epsilon}+\left(\frac{k^{2}C}{n^{2}\epsilon^{2}}+\frac{k^{3}C}{n^{3}\epsilon^{3}}\right)\min\left(k,\frac{1}{\epsilon}\right)\right),$
(A.61)

with probability at least $1-4\gamma$, recalling that $\epsilon=O(1)$ when
simplifying the expression. Algorithm 3 uses a constant failure probability of
$\gamma=0.01$, which ensures a success probability of at least $0.75$. This is
further boosted by Wrapper 1. Now, an application of Lemma A.4 gives the
stated result.

###  A.4.2 Proof of Theorem 2

The proof of this theorem proceeds nearly identically to that of Theorem 4
with some exceptions. If $\mathcal{S}=\mathcal{I}_{n,k}$, the smoothed
sensitivity bound has no $\min(k,1/\epsilon)$ term, owing to Lemmas A.13 and
A.14, which gives the bound

$\displaystyle\left|\mathcal{A}(\mathbf{X})-\theta\right|\leq\frac{1}{\sqrt{\gamma}}O\left(\sqrt{\text{var}(A_{n})}+\frac{k\xi}{n\epsilon}+\frac{k^{2}C}{n^{2}\epsilon^{2}}+\frac{k^{3}C}{n^{3}\epsilon^{3}}\right).$
(A.62)

Furthermore, since $\mathcal{S}=\mathcal{I}_{n,k}$, we have $A_{n}=U_{n}$ with
probability at least $1-3\gamma$. Algorithm 3 uses a constant failure
probability of $\gamma=0.01$. This ensures a success probability of at least
$0.75$, which is further boosted by Wrapper 1. An application of Lemma A.4
gives the stated result.

###  A.4.3 Concentration for local Hájek projections

###### Proof of Lemma 5.

We first show that if $Y_{1},Y_{2},\dots,Y_{t}$ are random variables such that
each $Y_{j}$ is $\tau_{j}$-sub-Gaussian, the sum ${Y_{1}+\dots+Y_{t}}$ is
$\left(\sqrt{\tau_{1}}+\sqrt{\tau_{2}}+\dots+\sqrt{\tau_{t}}\right)^{2}$-sub-
Gaussian.

Define $p_{j}=\frac{\sum_{i=1}^{t}\sqrt{\tau_{i}}}{\sqrt{\tau_{j}}}$. Clearly,
we have $\sum_{j=1}^{t}1/p_{j}=1$. By Hölder’s inequality, for any
$\lambda>0$, we have

$\displaystyle\mathbb{E}\left[\exp\left(\lambda\sum_{i=1}^{t}Y_{i}\right)\right]$
$\displaystyle=\mathbb{E}\left[\prod_{i=1}^{t}\exp(\lambda
Y_{i})\right]\leq\prod_{i=1}^{t}\mathbb{E}\left[\exp(\lambda
Y_{i})^{p_{i}}\right]^{1/p_{i}}\leq\prod_{i=1}^{t}\left(\exp\left(\frac{\lambda^{2}p_{i}^{2}\tau_{i}}{2}\right)\right)^{1/p_{i}}$
$\displaystyle=\prod_{i=1}^{t}\exp\left(\frac{\lambda^{2}p_{i}\tau_{i}}{2}\right)=\exp\left(\frac{\lambda^{2}(\sqrt{\tau_{1}}+\sqrt{\tau_{2}}+\dots+\sqrt{\tau_{t}})^{2}}{2}\right).$

Now, $h(X_{S})$ is sub-Gaussian$(\tau)$ for all $S\in\mathcal{I}_{n,k}$. Since
$\hat{h}({i})$ is the average of $t=\binom{n-1}{k-1}$ such quantities, it is
clear from the previous claim that $\hat{h}({i})$ sub-Gaussian with parameter

$\frac{1}{t^{2}}(\underbrace{\sqrt{\tau}+\sqrt{\tau}+\dots+\sqrt{\tau}}_{t\text{
terms}})^{2}=\tau.$

∎

###### Proof of Lemma 6.

Define $\sigma_{i}^{2}:=\text{var}\left(h(X_{S_{i}})|X_{i}=x_{i}\right)$.
First, conditioned on $X_{i}$, the projection $\hat{h}_{X}({i})$ can be viewed
as a U-statistic on the other $n-1$ data points. First, for this lemma, we use
$\mathcal{S}=\mathcal{I}_{n,k}$, and

$\displaystyle\hat{h}_{X}({i})=\widehat{\mathbb{E}}\left[h(X_{S})|X_{i}\right]=\frac{\sum_{S\in\mathcal{S}_{i}}h(X_{S})}{{n-1\choose
k-1}}$

Since $h$ is bounded, the random quantity
$\hat{h}_{X}({i})-\mathbb{E}\left[h(X_{S})|X_{i}\right]\in[-2C,2C]$ satisfies
the Bernstein moment condition and also the Bernstein tail inequality (cf.
Proposition 2.3 in [55]). By Bernstein’s inequality for U-statistics (see
inequality (A.32)), for all $t>0$, we have

$\displaystyle\mathbb{P}\left(\left.\left|\hat{h}_{X}({i})-\mathbb{E}\left[h(X_{S})|X_{i}\right]\right|\geq
t\right|X_{i}\right)\leq
2\exp\left(\frac{-{\left\lfloor\frac{n-1}{k-1}\right\rfloor}t^{2}}{2\sigma_{i}^{2}+4Ct/3}\right).$

We can check that the choice
$t=\sigma_{i}\sqrt{\frac{4k}{n}}\sqrt{\log\left(\frac{2n}{\beta}\right)}+\frac{8Ck}{n}\log\left(\frac{2n}{\beta}\right)$
makes this latter probability at most $\frac{\beta}{n}$. ∎

##  Appendix A.5 U-statistic applications

###  A.5.1 Uniformity testing

To motivate the test, consider the expectation
$\theta:=\mathbb{E}[h(X_{i},X_{j})]$ and the variance $\text{var}(U_{n})$:

######  Lemma A.17.

We have $\mathbb{E}[h(X_{1},X_{2})]=\frac{1}{m}+\frac{\left\lVert
a\right\rVert^{2}}{m^{2}}$. In particular, the means under the two hypothesis
classes differ by at least $\frac{\delta^{2}}{2m}$.

###### Proof.

We have

$\displaystyle\mathbb{E}[h(X_{1},X_{2})]=\sum_{i=1}^{m}p_{i}^{2}=\sum_{i=1}^{m}\frac{1+2a_{i}+a_{i}^{2}}{m^{2}}=\frac{1}{m}+\frac{\left\lVert
a\right\rVert^{2}}{m^{2}}.$

Under approximate uniformity, this is at most $\frac{\delta^{2}}{2m}$; and
under the alternative hypothesis, this quantity is at least
$\frac{\delta^{2}}{m}$. ∎

######  Lemma A.18.

The variance of $U_{n}$ is

$\displaystyle\text{var}(U_{n})=\frac{2}{n(n-1)}\left(2(n-2)\sum_{i<j}p_{i}p_{j}(p_{i}-p_{j})^{2}+\sum_{i=1}^{m}p_{i}^{2}-\left(\sum_{i=1}^{m}p_{i}^{2}\right)^{2}\right).$

###### Proof.

The conditional variances $\zeta_{1}$ and $\zeta_{2}$ can be written as

$\displaystyle\zeta_{1}$
$\displaystyle=\text{cov}(h(X_{1},X_{2}),h(X_{1},X_{3}))$
$\displaystyle=\mathbb{E}[{\mathbbm{1}[X_{1}=X_{2}]\mathbbm{1}[X_{1}=X_{3}]}]-\mathbb{E}[{\mathbbm{1}[X_{1}=X_{2}]}]\mathbb{E}[{\mathbbm{1}[X_{1}=X_{3}]}]$
$\displaystyle=\sum_{i}p_{i}^{3}-\left(\sum_{i}p_{i}^{2}\right)^{2}=\sum_{i<j}(p_{i}^{3}p_{j}+p_{i}p_{j}^{3}-2p_{i}^{2}p_{j}^{2})=\sum_{i<j}p_{i}p_{j}(p_{i}-p_{j})^{2}\geq
0,\text{ and}$ $\displaystyle\zeta_{2}$
$\displaystyle=\text{cov}(h(X_{1},X_{2}),h(X_{1},X_{2}))=\sum_{i=1}^{m}p_{i}^{2}-\left(\sum_{i=1}^{m}p_{i}^{2}\right)^{2}.$

Also from equation (3), we have

$\displaystyle\text{var}(U_{n})=\binom{n}{2}^{-1}\left(2(n-2)\zeta_{1}+\zeta_{2}\right).$

Combining the above bounds with equation (3) shows the result. ∎

#### Proof of Theorem 5:

Recall that $\tilde{\theta}$ denotes the private test statistic, which is
thresholded at the value $\frac{1+3\delta^{2}/4}{m}$ to determine the output
of the hypothesis test. We claim that the validity of the test is established
if we can show that

$\mathbb{P}\left(|\tilde{\theta}-\mathbb{E}[U_{n}]|\leq\frac{\delta^{2}}{4m}\right)\geq
1-O(\gamma)$ (A.63)

under both hypotheses. Indeed, it would then hold that:

  * (i)

Under approximate uniformity,

$\tilde{\theta}<\frac{1}{m}+\frac{\delta^{2}}{2m}+\frac{\delta^{2}}{4m}=\frac{1+3\delta^{2}/4}{m}.$

  * (ii)

Under the alternative hypothesis,

$\tilde{\theta}\geq\frac{1}{m}+\frac{\delta^{2}}{m}-\frac{\delta^{2}}{4m}=\frac{1+3\delta^{2}/4}{m}.$

To establish inequality (A.63), we further write

$\mathbb{P}\left(|\tilde{\theta}-\mathbb{E}[U_{n}]|>\frac{\delta^{2}}{4m}\right)\leq\mathbb{P}\left(|\tilde{\theta}-U_{n}|>\frac{\delta^{2}}{8m}\right)+\mathbb{P}\left(|U_{n}-\mathbb{E}[U_{n}]|>\frac{\delta^{2}}{8m}\right).$
(A.64)

The second term can be controlled using an argument in Diakonikolas et al.
[17], which further develops the variance bound in Lemma A.18 for the two
hypothesis classes and then uses Chebyshev’s inequality. It is shown that the
second probability term in inequality (A.64) can be bounded by $\alpha$ if
$n=\Omega\left(\frac{\sqrt{m}}{\gamma\delta^{2}}\right)$.

To bound the first probability term in inequality (A.64), we study the
concentration parameter $\xi$ for the local Hájek projection
$\hat{h}({i})=\frac{1}{n-1}\sum_{j\neq i}h(X_{i},X_{j})$. We have the
following result:

######  Lemma A.19.

If $\xi=\frac{6}{m}+\frac{8\log(4n/\gamma)}{n}$ and $n\geq\frac{16}{\gamma}$,
then $|\hat{h}({i})-U_{n}|\leq\xi$ for all $i$, with probability at least
$1-\gamma$.

###### Proof.

By the triangle inequality, we have

$\displaystyle|\hat{h}({i})-U_{n}|\leq|\hat{h}({i})-\mathbb{E}[h(X_{1},X_{2})|X_{1}]|+|\mathbb{E}[h(X_{1},X_{2})|X_{1}]-\theta|+|U_{n}-\theta|.$
(A.65)

We will provide a bound on each of these three terms. Note that

$h(X_{1},X_{2})|X_{1}\sim\operatorname{Ber}(p_{X_{1}}),$

which has variance

$\displaystyle\sigma_{i}^{2}$
$\displaystyle=\text{var}(h(X_{i},X_{j})|X_{i})=p_{X_{i}}(1-p_{X_{i}})\leq\frac{2}{m}.$

Hence, with probability at least $1-\frac{\gamma}{2}$, the first term in
inequality (A.65) can be bounded as

$\displaystyle|\hat{h}({i})-\mathbb{E}[h(X_{1},X_{2})|X_{1}]|$
$\displaystyle\leq
2\sqrt{\frac{2}{mn}\log\left(\frac{4n}{\gamma}\right)}+\frac{16}{3n}\log\left(\frac{4n}{\gamma}\right)\leq\frac{2}{m}+\frac{7}{n}\log\left(\frac{4n}{\gamma}\right),$

where we have used the AM-GM inequality. The second term in inequality (A.65)
can be bounded as

$\displaystyle\left|\frac{1+a_{X_{i}}}{m}-\left(\frac{1}{m}+\frac{\left\lVert
a\right\rVert^{2}}{m^{2}}\right)\right|\leq\frac{a_{X_{i}}}{m}+\frac{\left\lVert
a\right\rVert^{2}}{m^{2}}\leq\frac{2}{m}.$

Finally, by Chebyshev’s inequality, the third term can be bounded as
$|U_{n}-\theta|\leq\sqrt{\frac{2\text{var}(U_{n})}{\gamma}}$ with probability
at least $1-\frac{\gamma}{2}$. It remains to find the variance of $U_{n}$.

By Lemma A.18, if $\left|a_{i}\right|\leq 1$ for all $i$, we have

$\displaystyle\text{var}(U_{n})$
$\displaystyle=\frac{2}{n(n-1)}\left(2(n-2)\sum_{i<j}\frac{(1+a_{i})(1+a_{j})(a_{i}-a_{j})^{2}}{m^{4}}+\sum_{i}\frac{(1+a_{i})^{2}}{m^{2}}-\left(\sum_{i}\frac{(1+a_{i})^{2}}{m^{2}}\right)^{2}\right)$
$\displaystyle\leq\frac{2}{n(n-1)}\left(2(n-2)\sum_{i<j}\frac{(1+a_{i})(1+a_{j})(a_{i}-a_{j})^{2}}{m^{4}}+\sum_{i}\frac{(1+a_{i})^{2}}{m^{2}}\right)$
$\displaystyle\leq\frac{2}{n(n-1)}\left(2(n-2)\frac{4}{m^{4}}\binom{m}{2}+\frac{4m}{m^{2}}\right)\leq\frac{8}{m^{2}n}+\frac{8}{mn^{2}}.$
(A.66)

Combining the three bounds into inequality (A.65), with probability at least
$1-\gamma$, we have

$\displaystyle|\hat{h}({i})-U_{n}|$
$\displaystyle\leq\left(\frac{2}{m}+\frac{7}{n}\log\left(\frac{4n}{\gamma}\right)\right)+\frac{2}{m}+\frac{\sqrt{2}}{\sqrt{\gamma}}\left(\frac{2\sqrt{2}}{m\sqrt{n}}+\frac{2\sqrt{2}}{\sqrt{m}n}\right)$
$\displaystyle=\frac{4}{m}+\frac{7\log(4n/\gamma)}{n}+\frac{4/\sqrt{\gamma}}{m\sqrt{n}}+\frac{4/\sqrt{\gamma}}{\sqrt{m}n}$
$\displaystyle\leq\frac{6}{m}+\frac{8\log(4n/\gamma)}{n},$

where in the second inequality, we used the AM-GM inequality and the
assumption $n\geq\frac{16}{\gamma}$. The statement of the lemma follows after
discarding lower-order terms. ∎

By Lemma A.19, with probability at least $1-\gamma$, the weights of all
projections in Algorithm 3 are equal to $1$ and $U_{n}=\tilde{A}_{n}$. Then
$|\tilde{\theta}-U_{n}|$ is simply the magnitude of the noise added in the
final step of Algorithm 3 (which uses a constant $\gamma=0.01$), which (cf.
the proof of Theorem 2) takes the form

$\displaystyle
O\left(\frac{\xi}{n\epsilon}+\frac{1}{n^{2}\epsilon^{2}}+\frac{1}{n^{3}\epsilon^{3}}\right)=O\left(\frac{\log
n}{n^{2}\epsilon}+\frac{1}{mn\epsilon}+\frac{1}{n^{2}\epsilon^{2}}+\frac{1}{n^{3}\epsilon^{3}}\right),$

with probability at least $0.75$. This is bounded by $\frac{\delta^{2}}{8m}$
as long as

$\displaystyle
n=\Omega\left(\frac{m^{1/2}}{\delta\epsilon^{1/2}}\log\left(\frac{m^{1/2}}{\delta\epsilon^{1/2}}\right)+\frac{m^{1/2}}{\delta\epsilon}+\frac{m^{1/3}}{\delta^{2/3}\epsilon}+\frac{1}{\delta^{2}}\right).$

This constant probability of success is further boosted by Wrapper 1. Now, an
application of Lemma A.4 gives the stated result.

###  A.5.2 Proof of Theorem 6

The privacy of the algorithm follows by composing (see Lemma 1) the
$\epsilon$-privacy of $\nu$ and the $O(\epsilon)$-privacy of $\tilde{\theta}$
conditioned on $\nu$. It remains to show the utility of the algorithm.

The kernel $h(x,y)=1(\|x-y\|\leq r_{n})$ is degenerate, since
$P(\|X_{i}-X_{j}\|\leq r_{n}|X_{i})$ does not depend on $X_{i}$. So
$\text{var}(\mathbb{E}[h(X_{i},X_{j})|X_{i}=x])=0$. We have
$\text{var}[h(X_{i},X_{j})]=O(r_{n}^{2})$, so the non-private error is
$O(r_{n}/n)$ (Eq 3). Using Proposition 2.3 from Arcones and Gine [5], there
exist universal constants $c_{1},c_{2},$ and $c_{3}$ such that

$\displaystyle
P\left(\left|\frac{n-1}{2}(U^{\prime}_{n}-r_{n}^{2}/4)\right|\geq t\right)\leq
c_{1}\exp\left(-\frac{c_{2}t}{c_{3}r_{n}+(t/n)^{1/3}}\right).$

Setting $t=nr_{n}^{2}/16$, we have, for large enough $n$, since
$r_{n}=\Omega(n^{-1/2})$,

$\displaystyle P\left(\left|U^{\prime}_{n}-r_{n}^{2}/4\right|\geq
r_{n}^{2}/8\right)\leq
c_{1}\exp\left(-\frac{c_{2}nr_{n}^{2}}{c_{3}r_{n}+r_{n}^{2/3}}\right)\leq
c_{1}\exp(-c^{\prime}nr_{n}^{4/3})=\tilde{O}\left(\exp(-n^{1/3})\right).$

Therefore, with probability $1-o(1)$, we have

$\displaystyle U^{\prime}_{n}$ $\displaystyle\in[r_{n}^{2}/8,3r_{n}^{2}/8]$
$\displaystyle\nu$
$\displaystyle:=U^{\prime}_{n}+\frac{Z}{n\epsilon}\in[r_{n}^{2}/9,r_{n}^{2}/2].$
(A.67)

In particular, the probability that $U^{\prime}_{n}+\frac{Z}{n\epsilon}$
computed in step 3 of Algorithm 5 is then positive.

From Lemma 6, we have

$\max_{i}|\hat{g}_{X}({i})-\mathbb{E}[g(X_{i},X_{j},X_{k})|X_{i}]|=\tilde{O}\left(\frac{r_{n}^{2}}{\sqrt{n}}+\frac{1}{n}\right),$

with probability at least $1-\gamma$. Moreover, since $g$ is degenerate, we
have $\mathbb{E}[g(X_{i},X_{j},X_{k})|X_{i}]=\theta_{n}$. Using the fact that
$|U_{n}-\theta_{n}|\leq\sqrt{\frac{2\text{var}(U_{n})}{\gamma}}\leq\frac{r_{n}^{2}}{n}\sqrt{\frac{2}{\gamma}}$
with probability at least $1-\frac{\gamma}{2}$, we have

$\displaystyle\max_{i}|\hat{g}_{X}({i})-U_{n}|$
$\displaystyle\leq\max_{i}|\hat{g}_{X}({i})-\theta_{n}|+|\theta_{n}-U_{n}|$
$\displaystyle=2\sigma_{i}\sqrt{\frac{2}{n}\log\left(\frac{2n}{\gamma}\right)}+\frac{16}{3n}\log\left(\frac{2n}{\gamma}\right)+\frac{r_{n}^{2}}{n}\sqrt{\frac{2}{\gamma}}$
$\displaystyle\leq
18\nu\sqrt{\frac{2}{n}\log\left(\frac{2n}{\gamma}\right)}+\frac{16}{3n}\log\left(\frac{2n}{\gamma}\right)+\frac{9\nu}{n}\sqrt{\frac{2}{\gamma}}=:\xi,$

with probability $1-O(\gamma)$. Hence, using Theorem 2, and noting that
$\xi=\tilde{O}\left(\frac{r_{n}^{2}}{\sqrt{n}}\right)$, we can ensure that the
estimate $\tilde{\theta}$ output by Algorithm 3 satisfies

$\displaystyle\left|\tilde{\theta}-\theta\right|$
$\displaystyle=O\left(\frac{r_{n}^{2}}{n}+\frac{\xi}{n\epsilon}+\frac{1}{n^{2}\epsilon^{2}}+\frac{1}{n^{3}\epsilon^{3}}\right)=\tilde{O}\left(\frac{r_{n}^{2}}{n}+\frac{r_{n}^{2}}{n^{3/2}\epsilon}+\frac{1}{n^{2}\epsilon}+\frac{1}{n^{2}\epsilon^{2}}+\frac{1}{n^{3}\epsilon^{3}}\right)$
$\displaystyle=\tilde{O}\left(\frac{r_{n}^{2}}{n}+\frac{1}{n^{2}\epsilon^{2}}\right),$

with probability at least $0.75$, where we used the AM-GM inequality to
eliminate the term $\frac{r_{n}^{2}}{n^{3/2}\epsilon}$. This constant
probability of success is further boosted by Wrapper 1. Now, an application of
Lemma A.4 gives the stated result.