# The Black–Hole Masses of High–Redshift QSOs

Andrew King1,2,3
1 School of Physics & Astronomy, University of Leicester, Leicester LE1 7RH UK
2 Astronomical Institute Anton Pannekoek, University of Amsterdam, Science
Park 904, NL-1098 XH Amsterdam, The Netherlands
3 Leiden Observatory, Leiden University, Niels Bohrweg 2, NL-2333 CA Leiden,
The Netherlands

###### Abstract

Observations of high–redshift quasars frequently promote suggestions of large
black hole masses, whose presence so early in cosmic time is not easily
explicable. I consider the parallel with ultraluminous X–ray sources (ULXs) –
now known to be stellar–mass black hole (and neutron star) binaries apparently
radiating far above their Eddington luminosities $L_{\rm Edd}$. The true
luminosity in ULXs is actually only of order $L_{\rm Edd}$, for stellar–mass
accretors, but has a very anisotropic (‘beamed’) component, plus a
near–isotropic component of similar luminosity but much lower specific
intensity. Observers viewing ULXs from within the beam but assuming spherical
symmetry deduce a luminosity $\gg L_{\rm Edd}$. These features appear because
the accretors are fed mass at highly super–Eddington rates, most of it
expelled in high–speed ($v>0.2c$) outflows from the accretion disc.

I show that in similarly–beamed AGN, emission–line properties would be
essentially the same as in unbeamed sources, but standard virial mass
indicators unusable because velocity widths are dominated by the outflows, not
bound motions about the black holes. In an ensemble of this kind the
apparently most luminous systems are always the most distant, but have the
lowest black hole masses. Interpreting observations of this ensemble without
knowing that they are beamed leads instead to very high black hole mass
estimates. The analogy with ULXs therefore suggests that high–redshift quasars
might actually have central black hole masses which could have grown from
stellar values within the lookback time. I consider how one might test these
ideas observationally.

###### keywords:

galaxies: active: supermassive black holes: black hole physics: X–rays:
galaxies

11footnotetext: E-mail<EMAIL_ADDRESS>
##  1 Introduction

Observations of quasars at redshift $\sim 6$ are often interpreted as
requiring black hole masses $M\gtrsim 3\times 10^{9}\rm M_{\odot}$ (see e.g.
Willott, McLure & Jarvis, 2003 for a very clear description). Depending on the
detailed accretion history, these lead to suggestions that these holes descend
from very large black hole seed masses $\gtrsim 10^{8}\rm M_{\odot}$ which
must have formed very rapidly at very high redshift, and there have been many
efforts to model what these might be (e.g. Begelman, Volonteri & Rees, 2006;
Begelman, 2009; Regan & Haehnelt, 2013; Johnson et al., 2013)111One can avoid
the need for massive seeds if the accretion process keeps the black hole spin
low. This in turn keeps the accretion efficiency low, so that the same
luminosity corresponds to a higher mass. This can happen if accretion on to
the hole is chaotic, i.e. from a succession of accretion events with
uncorrelated angular momenta. See King, Pringle & Hofmann, 2008 and Section
5.4 of King (2023) for a detailed discussion of this point..

An important part of the argument leading to this is the assumption that the
apparent luminosity222In this paper I will use the term ‘apparent luminosity’
and the symbol $L_{\rm sph}$ to denote cases where emission is assumed
isotropic, even though it may in reality be strongly anisotropic. $L\sim\nu
L_{\nu}$ (where $L_{\nu}$ is the monochromatic luminosity, assumed isotropic)
of the observed quasars cannot significantly exceed the Eddington limit

${L_{\rm Edd}}\simeq 1.3\times 10^{46}M_{8}\,{\rm erg\,s^{-1}}$ (1)

for black hole mass $M=10^{8}\rm M_{\odot}M_{8}$.

Accretion on to a black hole is an essentially scale–free process, since the
mass $M$ of the accretor defines typical lengthscales through the
gravitational radius $R_{g}=GM/c^{2}$, timescales through $R_{g}/c$, and
luminosities through ${L_{\rm Edd}}(M)$. Accordingly it is often helpful to
compare analogous accreting stellar–mass and supermassive black hole systems.
Here I ask what we can learn from the stellar–mass version of the argument
leading to very large black hole masses for high–redshift quasars.

##  2 Ultraluminous X–ray Sources

The stellar–mass analogues of quasars with very high deduced black hole masses
are the ultraluminous X–ray sources (ULXs). About 1800 systems are known by
now (see King, Lasota & Middleton, 2023 for a recent review).

ULXs are extragalactic systems which are not located in the centres of their
host galaxies, but whose apparent luminosities ($\gg 10^{39}\,{\rm
erg\,s}^{-1}$), if assumed isotropic, require (through (1)) accretor masses
significantly larger than likely for black holes descending from stellar
evolution. Colbert & Mushotzky (1999) identified ULXs as a class, and
suggested that they were systems with accretion on to hypothetical
‘intermediate–mass’ black holes (or IMBH), where the hole mass was
significantly larger ($\sim 10^{4}-10^{5}\rm M_{\odot}$) than could result
from stellar evolution.

But King et al. (2001) noted that mass transfer at super–Eddington rates in
otherwise standard X–ray binaries (with normal stellar–mass accretors) was
very likely to lead to strong ejection of dense winds from the accretion disc,
as pointed out by Shakura & Sunyaev (1973). That paper had shown that a disc
accreting at a rate which exceeds the local radiation pressure limit within a
certain disc radius would eject the excess at each radius further in,
ultimately producing total luminosity

$L=L_{\rm Edd}(1+\ln\dot{m}),$ (2)

where $\dot{m}=\dot{M}_{\rm tr}/\dot{M}_{\rm Edd}\gg 1$ is the ratio of the
rate $\dot{M}_{\rm tr}$ of mass transfer towards the accretor to the value
$\dot{M}_{\rm Edd}=L_{\rm Edd}/\eta c^{2}$ (with $\eta\simeq 0.1$ the
accretion efficiency, given by the specific binding energy of the black hole
ISCO radius), which is the rate at which the accretor mass actually grows.

King et al. (2001) suggested that this would make a significant part of the
luminosity strongly anisotropic, i.e. beamed333In this paper ‘beaming’ simply
means collimated by scattering, and not more exotic processes such as
relativistic beaming., as much of the accretion luminosity emitted from close
to the black hole finds the two narrow vacuum funnels along the rotational
axis of the accretion disc and escapes along them. The disc emission which
does not find these funnels produces a near–isotropic luminosity of a similar
order $\sim L_{\rm Edd}$. This is emitted from the outer photosphere of the
wind (cf King & Muldrew 2016) and spread over $\sim 4\pi$ steradians.

An observer viewing one of these systems from a location within the beam solid
angle $4\pi b$ (where $b\ll 1$) would measure a very high specific intensity,
and if they assumed isotropic emission, would assign an apparent luminosity
$L_{\rm sph}=(1/b)L\gg L$, where $L$ is the true total power output of the
source. In this case the deduction of unusually large black hole masses from
the assumption that $L_{\rm Edd}>L_{\rm sph}$ is unjustified, as it implies a
minimum apparent mass $M_{\rm sph}$ greater than the true minimum value by a
factor $1/b\gg 1$. An observer viewing a system of this type from outside the
beam (i.e. by far the more common case) would assign an apparent luminosity
$\sim L_{\rm Edd}$ because of the much lower specific intensity. So even
though even though the system’s apparent luminosity is very high when viewed
along the beam, this substantial near–isotropic luminosity seen in all other
directions must in all cases characterize the source’s effect on its
surroundings, as for example seen in emission line data, which are therefore
the same whether one views the beamed emission or not.

Very recently, X–ray polarimetric observations by Veledina et al. (2023) have
effectively settled two decades of controversy over the cause of ULX
behaviour, concluding ‘unambiguously’ that much of the accretion energy
produced by the stellar–mass X–ray binary Cyg X-3 does follow two narrow
channels (in that particular system pointing away from the observer).

This agrees closely with the picture suggested by King et al. (2001) (see
Lasota & King 2023 for a discussion) which implies that ULX behaviour is
simply a phase of standard X–ray binary evolution characterized by very high
mass transfer rates. This phase occurs whenever a donor star which is more
massive than the accretor fills its Roche lobe. Then mass transfer decreases
the binary separation and amplifies the mass transfer, until this reaches the
(very high) limiting rate set by the donor star’s thermal timescale. I will
show later in the paper that high–redshift supermassive black holes may at
some point have to deal with similarly super–Eddington mass supply rates,
producing effects like those seen in ULXs, and consequent problems for
deducing masses.

For ULXs, King (2009) showed that quite generally,

$b\simeq\frac{73}{{\dot{m}}^{2}}.$ (3)

This formula agrees with simple physical reasoning about the open funnel
geometry (King, 2009), since conditions at large disc radii are set by
$\dot{M}_{\rm tr}$, while if $\dot{m}>1$, conditions at small disc radii
converge to values set by the Eddington rate $\dot{M}_{\rm Edd}$. Assuming
normal stellar–mass black hole and neutron star accretors, equation (3) gives
acceptable fits to all observations of ULXs (see King et al., 2023 for a
review).

##  3 Mass Estimates and Outflows

The results summarized in the last Section, and the scale–free nature of black
hole accretion, suggest that it may be interesting to consider what
observations of high–$z$ quasars would tell us if the sources are assumed
beamed like ULXs, as a result of being fed mass at super–Eddington rates,
rather than radiating isotropically. Since black–hole accretion is scale–free
we can assume that (3) holds for such ULX–like quasars. Note that the argument
that Haiman & Cen (2002) made against quasar beaming does not apply to the
picture I consider here – they considered a quasar model beamed towards the
observer, but which they explicitly assumed had a low total luminosity. In
constrast we saw in the last Section that ULX–type emission is always
accompanied by a substantial near–isotropic component, which would evidently
produce the same emission–line properties as an unbeamed quasar.

The first consequence is that for ULX–like quasars the luminosity alone no
longer requires unusually high black hole masses in order to avoid conflict
with the Eddington limit, just as one had the freedom to consider stellar–mass
ULX models rather than the original IMBH proposal for explaining the original
ULX systems.

But the deduction of large masses in high–redshift quasars generally comes
from the use of virial mass indicators rather than the quasar luminosities (cf
Willott, McLure & Jarvis, 2003). Virial indicators use the substantial
observed velocity widths $\Delta V$ of emission lines such as H$\beta$ and Mg
II to deduce the SMBH mass by assuming the the gas is in bound orbits in the
AGN’s broad–line region (BLR). At low redshift, reverberation mapping of AGN
gives the response time of the BLR to central AGN luminosity changes. This
fixes the BLR size $R$, and allows a black–hole mass estimate

$M\sim\frac{R(\Delta V)^{2}}{G}.$ (4)

At high redshift the much larger expected size of the BLR (several light
years) makes reverberation mapping difficult, so the SMBH mass estimates
proceed by extrapolating the correlation between continuum luminosity and BLR
size observed for low–redshift systems. But a serious problem in using this
virial method to estimate masses for the hypothesised ULX–like quasars is that
the spectra of stellar–mass ULXs are dominated by strong and fast outflows.
The velocity widths then do not correspond to motions of line–emitting
material in bound orbits about the SMBH, but radial winds with the escape
speed from the black hole’s vicinity.

In stellar–mass ULXs these winds remove most of the super–Eddington mass
supply rate $\dot{M}$, and have observed velocity widths $>0.2c$ (e.g. Walton
et al., 2016). This agrees with the values expected for momentum–driven flows,
where the gas velocity $v$ is given by equating outflow and radiation momenta
as

$\dot{M}v\simeq\frac{L_{\rm Edd}}{c},$ (5)

which implies a radial outflow speed

$v\simeq\frac{\eta}{\dot{m}}c,$ (6)

where $\eta\sim 0.1$ is the accretion efficiency. As expected, this gives no
direct information about the black hole mass, as the problem is scale–free.

Momentum–driven ‘UFO’ winds like this, characterized by P Cygni profiles in
their X–ray spectra, are very widely observed in bright AGN (e.g. Tombesi et
al., 2010a, b and many subsequent papers). Since the velocity width is no
longer given by orbital motion of bound gas in the BLR we cannot use it to
deduce the mass of the central supermassive black hole if the AGN are
‘ULX–like’.444 Physically, these UFO winds correspond to cases where the
reverse shock produced where the outflow collides with the host galaxy
interstellar gas is efficiently cooled, typically by the inverse Compton
effect of the AGN radiation field (see King & Pounds, 2015 and King, 2023 for
reviews). This prevents the thermal expansion of the shocked gas driving much
stronger and faster energy–driven outflows. But for a sufficiently large SMBH
mass, Compton shock cooling becomes inefficient, the outflows become
energy–driven, and now drive off all of the accreting gas. This fixes the SMBH
mass, and is the likely origin of the $M-\sigma$ relation.

##  4 Masses and Distances

The quadratic dependence of the beaming factor $b$ (see eqn 3) on
$\dot{m}=\dot{M}/\dot{M}_{\rm Edd}$ means that deducing limits on accretor
masses in disc–wind beamed systems is more complex than for isotropic
emission.

Both $\dot{M}_{\rm Edd}$ and $L_{\rm Edd}$ are proportional to the accretor
mass $M$, so for a given mass supply rate $\dot{M}$ we have

$b\propto\left(\frac{M}{\dot{M}}\right)^{2}$ (7)

giving the apparent (assumed isotropic) luminosity

$L_{\rm sph}\simeq\frac{L_{\rm Edd}}{b}\propto\frac{\dot{M}^{2}}{M},$ (8)

i.e. the tighter beaming for a lower–mass accretor outweighs its lower
Eddington luminosity when we predict its apparent luminosity.

Then the standard procedure of deducing an accretor mass $M_{\rm sph}$ by
arguing that the apparent luminosity $L_{\rm sph}$ cannot exceed the Eddington
limit for the corresponding mass $M_{\rm sph}$ results in an estimate

$M_{\rm sph}\propto\frac{\dot{M}^{2}}{M}$ (9)

for the minimum accretor mass. Thus in an ensemble of beamed sources with
similar super–Eddington mass supply rates $\dot{M}$, the sources with the
highest apparent luminosities actually have the lowest masses because their
beaming is quadratically tighter.

But beaming means that many of these sources will not be observable. In the
simple case where redshifts are modest, to lie in the beam of one such source
one has to search a space volume $\sim 1/Nb$, where $N$ is the spatial number
density of the sources. Then the nearest detectable beamed source is at a
distance

$D\sim\left(\frac{3}{4\pi
Nb}\right)^{1/3}\propto\left(\frac{\dot{M}}{NM}\right)^{2/3}.$ (10)

We see that in observations of a collection of accreting sources with similar
mass supply rates $\dot{M}$, beamed according to (3),

the systems with the highest apparent luminosities and apparent masses are in
reality the most distant, i.e. the most redshifted, of those sources with the
lowest black hole masses555In the stellar–mass case, King (2009) used these
arguments to point out that some extragalactic sources usually assumed to be
AGN might actually be ULXs (‘pseudoblazars’), distinguishable from real AGN
only by the fact that they are not constrained to lie in the dynamical centres
of their host galaxies..

If these same beamed sources were instead assumed to radiate isotropically, we
would deduce that the highest redshift sources were intrinsically the
brightest, and because of the Eddington limit, must have very high black hole
masses. This is exactly the result one gets using the straightforward
interpretation of the data on distant quasars.

These very simple arguments already suggest that ULX–type beaming may have
interesting implications for the discussion of the masses of high–redshift
QSOs. But it is clear that taking these ideas further requires careful
consideration of the evolution of number density $N$ with redshift.

##  5 Discussion

The arguments given in this paper suggest that high–redshift QSO central black
holes may have relatively modest masses if in reality they are fed at
potentially strongly super–Eddington rates.

This is already known to hold for stellar–mass ULX systems. Extending this
picture to high–$z$ QSOs is reasonable, since at the relevant redshifts it is
likely that some black hole masses may not have reached large values, and for
example may lie below the $M-\sigma$ value

$M_{\sigma}\simeq 3.2\times 10^{8}\rm M_{\odot}\sigma_{200}^{4}$ (11)

(e.g. Kormendy & Ho, 2013) corresponding to the velocity dispersion
$\sigma=200\sigma_{200}\,{\rm km\,s^{-1}}$ of their host galaxy bulges. If the
galaxy is strongly disturbed (e.g. by a merger), mass would be supplied at a
dynamical rate

$\dot{M}_{\rm dyn}\lesssim\frac{f_{g}\sigma^{3}}{G},$ (12)

where $f_{g}\simeq 0.16$ is the cosmological gas fraction, giving an Eddington
factor

$\dot{m}\sim\frac{\dot{M}_{\rm dyn}}{\dot{M}_{\rm
Edd}}=\frac{54}{M_{8}^{1/4}}\frac{M_{\sigma}}{M}$ (13)

(cf King & Pounds, 2015; King, 2023). From (3) this implies typical beaming
fractions $b\lesssim 10^{-2}$, with corresponding reductions in the deduced
true SMBH mass $\sim bM_{\rm sph}$.

Evidently the condition for ULX–type disc–wind beaming of the accretion
luminosity is readily satisfied in many high–redshift galaxies. Given the
scale–free property of black hole accretion it would be surprising if Nature
made use of this effect for stellar–mass black holes only.

Future observations should provide direct tests of these ideas. For
stellar–mass ULXs X–ray polarimetry revealed the beaming. For AGN systems at
nonzero redshift the relevant photon energy is softer, and it is unclear
whether a test like this would be possible, and if so, conclusive. The recent
discovery of the very bright QSO system J0529–4351 (Wolf et al., 2024) may
allow direct measurement of its large expected BLR size, and so test the
extrapolations used in current methods for estimating black hole masses in
high–$z$ QSOs. Theoretically, calculations allowing one to discriminate
between emission lines from outflows and from orbiting gas would be very
useful.

## DATA AVAILABILITY

No new data were generated or analysed in support of this research.

## ACKNOWLEDGMENTS

I thank Jean–Pierre Lasota for many discussions of ULX behaviour, and the
referee for thoughtful comments which improved the paper.

## References

  * [1] Begelman, M.C., 2010, MNRAS, 402, 673
  * [2] Begelman, M.C., King, A.R., Pringle, J.E., 2006, MNRAS, 370, 399
  * [3] Begelman, M.C., Volonteri, M., Rees, M.J., 2006, MNRAS, 370, 289
  * [4] Colbert, E.J.M., Mushotzky, R.F., 1999, ApJ 519, 89
  * [5] aiman, Z., Cen, R., 2002, ApJ 578, 702
  * [6] Johnson, J.J.L., Whalen, D.J., Li, H., Holz, D.E, 2013, ApJ, 771,116
  * [7] Kajava. J.J.E., & Poutanen, J., 2009, MNRAS, 398,1450
  * [8] King, A.R., 2009, MNRAS 393, L41
  * [9] King, A.R., & Pounds, 2015, Annual Reviews of Astronomy and Astrophysics, 53, 155
  * [10] King, A.R., 2023, Supermassive Black Holes, Cambridge University Press
  * [11] King, A.R., Lasota, J.-P., Middleton, M., 2023, New Astronomy Reviews 2023, Volume 96, article id. 101672
  * [12] King, A.R., & Muldrew, S.I., 2016, MNRAS, 455, 1211
  * [13] King, A.R., Pringle, J.E., Hofmann, J.A., 2008, MNRAS, 385, 1621
  * [14] Kormendy, J., Ho, L.C., 2013, ARAA, 51, 511
  * [15] Lasota, J.-P., King, A.R., 2023, MNRAS, 526, 2506
  * [16] Lommen, D., Yungelson, L., van den Heuvel, E., Nelemans, G., Portegies Zwart, S. 2005, A&A, 443, 231
  * [17] Miller, M.C., Hamilton, D.P., 2002, MNRAS 330, 232
  * [18] Portegies Zwart, S.F., McMillan, S.L.W., 2002, ApJ 576, 899
  * [19] Regan, J. A., Haehnelt, M.G., 2009, MNRAS, 396, 343
  * [20] Shakura, N.I., Sunyaev, R.A., 1973, A&A 24, 337
  * [21] Tombesi, F., Cappi, M., Reeves, J. N., Palumbo, G. G. C., Yaqoob, T., Braito, V., Dadina, M., 2010a, A&A, 521
  * [22] Tombesi, F., Sambruna, R.M., Reeves, J.N., Braito, V., Ballo, L., Gofford, J., Cappi, M., Mushotzky, R. F., 2010b ApJ 719, 700
  * [23] Veledina, A., Muleri, F., Poutanen, J., and 127 more, arXiv:2303.01174
  * [24] Walton, D.J., Middleton, M.J., Pinto, C., Fabian, A.C., Bachetti, M., Barret, D., Brightman, M., Huerst, F., Harrison, F.A., Miller, J.M., Stern, D. 2016, ApJ 826, L26b
  * [25] Willott, C.J., McLure, R.J., Jarvis M.J., 2003, ApJ, 587, L15
  * [26] Wolf., C., et al, 2024, arXiv:2402.15101


# Skewness of local logarithmic exports

Sung-Gook Choi  Department of Physics, Inha University, Incheon 22212, Korea
Korea Housing & Urban Guarantee Corporation, Busan 48400, Korea  Deok-Sun Lee
School of Computational Sciences, Korea Institute for Advanced Study, Seoul
02455, Korea<EMAIL_ADDRESS>
###### Abstract

The distributions of trade values and relationships among countries and
product categories reflect how countries select their trade partners and
design export portfolios. Here we consider the exporter-importer network and
the exporter-product network with directed links weighted by the logarithm of
the corresponding export values each year from 1962 to 2018, and study how the
weights of the outgoing links from each country are distributed. Such local
logarithmic export distributions by destinations and products are found to
follow approximately the Gaussian distribution across exporters and time,
implying random assignment of export values on logarithmic scale. However, a
non-zero skewness is identified, changing from positive to negative as
exporters have more partner importers and more product categories in their
portfolios. Seeking the origin, we analyze how local exports depend on the
out-degree of exporter and the in-degrees of destinations/products and
formulate their quantitative and measurable relation incorporating randomness,
which uncovers the fundamental nature of the export strategies of individual
countries.

##  I Introduction

Individual countries, though equipped with different resources for production
and put in different environments, commonly fine tune the strategies on how to
allocate resources, invest capital, and make international trade partnership
towards the economic development Hidalgo  _et al._ (2007); Tacchella  _et al._
(2018); Chenery and Taylor (1968); Kuznets (1971); Lei  _et al._ (2015); Choi
and Lee (2019). Such efforts are universal and thus can be imprinted in the
patterns of domestic production, consumption, and international trade.
Especially in international trade are reflected the decisions on what
products, to what countries, and how much each country exports Hummels and
Klenow (2005). The empirical and analytical study of the trade value
distributions may therefore illuminate the underlying principles of the global
economic systems.

The network approach has been instrumental to quantify various nature of
international trade, revealing its far-from-random connection among countries
Serrano and Boguñá (2003); Garlaschelli and Loffredo (2004); Fagiolo  _et al._
(2009), and the broadness of the distributions of trade values by countries
Cha  _et al._ (2010); Almog  _et al._ (2015); Oh and Lee (2017) or products
Hidalgo and Hausmann (2009); Tacchella  _et al._ (2012); Barbier and Lee
(2017). These studies mostly aimed at identifying typical representative
characteristics, consequently neglecting local variations which may be
essential for understanding the principles governing the trade strategies of
individual countries. In this light, here we investigate how the export values
of a country are distributed by its destinations (partner importers) and
product categories. We compile the international trade data-sets from two
sources Feenstra  _et al._ (2005); wto , covering different time periods, to
obtain the annual export value for each pair of an exporter and an importer
country, and for each pair of an exporter country and a product category.
Using the result, we construct the exporter-importer (EI) network and the
exporter-product (EP) network each year, where links are directed, from
exporter countries to their partner importers or to their exported product
categories, and are weighted by the logarithm of the corresponding annual
export values. Then the distribution of the weights of the outgoing links from
each country represents how the values of the exports from the country are
distributed by destinations or products, which we analyze in the present
study.

We find that such local logarithmic export distributions measured for
individual exporters follow approximately the Gaussian distribution, implying
the possibility of the first two moments to fully characterize the whole
distributions. However, the distributions exhibit significant asymmetry and
significant non-zero skewness. Moreover, the skewness shows a systematic
variation - decreases with the out-degree of exporter, the number of partner
importers or of the exported products. In seeking the origin, we first find
that the average export of a country per destination/product increases with
its out-degree, meaning that the whole distribution is pushed relatively to
the right (left) for a hub (non-hub) exporter and thus may give rise to a
thicker left (right) tail. From the network perspective, this finding
characterizes the influence of a source node on the average weight of its
outgoing links. Individual link weights are found to be dependent also on the
in-degrees of the target nodes, the number of connected exporters of
destination/product in the EI/EP network. Yet the dependence is accompanied by
large fluctuations. We derive the formula decomposing a local export into a
deterministic term dependent on the in-degree of the target and a random
noise, which is well defined empirically and leads us to disentangle the
skewness of local exports in terms of the third order joint cumulants of the
in-degrees of the targets and random noises. These results therefore uncover
the mechanisms generating the skewed distribution of local exports but also
provide a quantitative framework to analyze the architecture of weighted
complex networks.

The paper is organized as follows. In Sec. II, we introduce briefly the data-
sets and present the empirical findings on the skewness of the local
logarithmic export distributions. We study in Sec. III the distributions of
the in-degrees of destinations and products receiving links from each exporter
and compare the skewness of such local in-degrees and local exports. In Sec.
IV, we show how the in-degrees of destinations/products and random noises
together generate the skewness of local exports. We discuss our findings and
future works in Sec. V.

##  II Local export distributions

###  II.1 Data-sets and two trade networks

We compile the NBER-UN data Feenstra  _et al._ (2005) and the World Trade
Organization (WTO) data wto  to obtain the annual export value
$V_{cpc^{\prime}}(t)$ in US dollars of a product category $p$ exported by a
country $c$ to another country $c^{\prime}$ in the period $1962\leq t\leq
2018$. The number of countries exporting or importing at least one product
ranges between 154 and 190 and the total number of exported products
fluctuates between 686 and 1458 in the studied period [Appendix A1].

Each country exports selected products to selected destination countries in
different amounts, which constitutes its export portfolio. To represent the
whole collection of such trade relationships, here we consider two kinds of
directed and weighted networks, the exporter-importer(EI) network and the
exporter-product(EP) network each year as follows [Fig. A1]. Nodes in the EI
network are countries exporting or importing at least one product in given
year $t$. A link is assigned from a country $c$ to another country
$c^{\prime}$ if the aggregated export from $c$ to $c^{\prime}$,
$V_{cc^{\prime}}(t)=\sum_{p}V_{cpc^{\prime}}(t)$, is positive. The (binary)
adjacency matrix ${\bf A}^{\rm(EI)}=(A_{cc^{\prime}})$ is given by
$A_{cc^{\prime}}=\theta(V_{cc^{\prime}})$ with $\theta(x)=1$ for $x>0$ and $0$
for $x\leq 0$. The EP network is a bipartite network consisting of two types
of nodes, countries and product categories, and they are connected if
$V_{cp}(t)=\sum_{c^{\prime}}V_{cpc^{\prime}}(t)>0$. Therefore the adjacency
matrix ${\bf A}^{\rm(EP)}=(A_{cp})$ is given by $A_{cp}=\theta(V_{cp})$. The
links in the EI and the EP network are then weighted by the base-10 logarithm
of the corresponding export values as

$\displaystyle v_{cc^{\prime}}(t)$
$\displaystyle=\log_{10}\left(V_{cc^{\prime}}(t)\right),$ $\displaystyle
v_{cp}(t)$ $\displaystyle=\log_{10}\left(V_{cp}(t)\right).$ (1)

The link-weight matrix, or the weighted adjacency matrix, is given by ${\bf
v}^{\rm(EI)}=(v_{cc^{\prime}})$ and ${\bf v}^{\rm(EP)}=(v_{cp})$,
respectively. We take logarithm given that export values span several orders
of magnitude Bhattacharya  _et al._ (2008); Fagiolo  _et al._ (2009); Cha  _et
al._ (2010); Oh and Lee (2017). As we will show, the distributions of $v$’s
are close to the Gaussian distribution and we will investigate their
characteristics. On the other hand, the distribution of the raw values of
export $V$’s are close to power-laws Cha  _et al._ (2010); Oh and Lee (2017).

###  II.2 Distribution of local logarithmic exports

How a country $c$ divides its export among different destinations
($c^{\prime}$) or product categories ($p$) can be represented by the
distributions of its logarithmic export values

$\displaystyle P_{c}^{\rm(EI)}(v)$ $\displaystyle={1\over
q_{c}^{\rm(EI)}}\sum_{c^{\prime}\in\mathtt{I}_{c}^{\rm(EI)}}\delta(v_{cc^{\prime}}-v),$
$\displaystyle P_{c}^{\rm(EP)}(v)$ $\displaystyle={1\over
q_{c}^{\rm(EP)}}\sum_{p\in\mathtt{I}_{c}^{\rm(EP)}}\delta(v_{cp}-v),$ (2)

where $\mathtt{I}_{c}^{\rm(EI)}=\\{c^{\prime}|A_{cc^{\prime}}>0\\}$ and
$\mathtt{I}_{c}^{\rm(EP)}=\\{p|A_{cp}>0\\}$ are the set of the partner
importers and of the exported products of $c$. The out-degree
$q_{c}^{\rm(EI)}$ and $q_{c}^{\rm(EP)}$ are the size of those sets equal to
the number of outgoing links from $c$,

$\displaystyle q_{c}^{\rm(EI)}$
$\displaystyle=\sum_{c^{\prime}}A_{cc^{\prime}}=\sum_{c^{\prime}\in\mathtt{I}_{c}^{\rm(EI)}}1,$
$\displaystyle q_{c}^{\rm(EP)}$
$\displaystyle=\sum_{p}A_{cp}=\sum_{p\in\mathtt{I}_{c}^{\rm(EP)}}1.$ (3)

$q_{c}^{\rm(EI)}$ and $q_{c}^{\rm(EP)}$ can be considered as a measure of the
diversity of country $c$ in the trade partnership and inventory.
$P_{c}^{\rm(EI)}(v)$ and $P_{c}^{\rm(EP)}(v)$ correspond to ego-centric link-
weight distributions in the EI and the EP network and we call them local
(logarithmic) export distributions in that they are measured only for the
outgoing links of $c$.

Figure 1:  Statistics of the logarithm of local exports. (a, b) Distribution
of the logarithm of the exports from each country $c$ (a) to different
destinations and (b) of different product categories in year 1980. Color
varies with the out-degree of $c$. (c, d) The mean and standard deviation of
the logarithmic exports from country $c$ (c) to different destinations and (d)
of different products versus the out-degree of $c$. (e, f) Distribution of the
standardized logarithmic exports from country $c$ (e) to different
destinations and (f) of different products. The dashed lines in (e) and (f)
are $\tilde{P}^{\rm(Gaussian)}(\tilde{v})={e^{-{\tilde{v}^{2}\over
2}}\over\sqrt{2\pi}}$. (g, h) The skewness of the logarithmic exports from $c$
(g) to different destinations and (h) of different products versus the out-
degree of $c$. Shades in (g) and (f) represent the middle 90% range of the
random-version skewness.

The local export distributions $P_{c}(v)$ commonly take a bell shape while the
mean $\langle v\rangle_{c}^{\rm(EI)}$ or $\langle v\rangle_{c}^{\rm(EP)}$ is
different from exporter to exporter [Fig. 1 (a-d)]. We present the results
obtained for year 1980 as an example in most figure plots, but all the
statements of ours are valid throughout the studied period unless stated
otherwise. The standard deviation $\sigma_{v;c}^{\rm(EI)}$ or
$\sigma_{v;c}^{\rm(EP)}$ does not vary with $c$ as significantly as the mean
does. Let us introduce the standardized variable $\tilde{v}$ such that

$v_{c\nu}=\langle
v\rangle_{c}^{\rm(EG)}+\sigma_{v;c}^{\rm(EG)}\tilde{v}_{c\nu}\ {\rm with}\
(\nu,{\rm G})=\left\\{\begin{array}[]{l}(c^{\prime},{\rm I})\ {\rm or}\\\
(p,{\rm P})\end{array}\right.$ (4)

where we used $\nu$ and $G$ to avoid repeating almost the same expressions
such that $\nu$ represents an importer $c^{\prime}$ or a product $p$, and G
represents I or P. Then the standardized distributions
$\tilde{P}_{c}(\tilde{v})$ for different $c$’s look collapsing onto the
Gaussian function of zero mean and unit variance
$\tilde{P}^{\rm(Gaussian)}(\tilde{v})={e^{-{\tilde{v}^{2}\over
2}}\over\sqrt{2\pi}}$ [Fig. 1 (e) and (f)]. If the data collapse were perfect,
$\\{\tilde{v}_{c\nu}|\nu\in\mathtt{I}_{c}\\}$ would be considered as Gaussian
random variables satisfying $\langle\tilde{v}\rangle_{c}=0,\
\langle\tilde{v}^{2}\rangle_{c}=1,\ \langle\tilde{v}^{n}\rangle_{c}=0$ for
$n\geq 3$ with $\langle\tilde{v}^{n}\rangle_{c}={1\over
q_{c}}\sum_{\nu\in\mathtt{I}_{c}}\tilde{v}^{n}_{c\nu}$. We would think that
every country divides its export among destinations and products randomly on
logarithmic scale universally and that the mean and the standard deviation
fully characterize the statistics of the local logarithmic export values.

###  II.3 Skewness

The data collapse is not perfect, however. To quantify the deviation, we
measure the skewness of local logarithmic exports for each exporter $c$ as

$\gamma_{v;c}={{1\over q_{c}}\sum_{\nu\in\mathtt{I}_{c}}(v_{c\nu}-\langle
v\rangle_{c})^{3}\over{\sigma_{v;c}}^{3}}=\langle\tilde{v}^{3}\rangle_{c},$
(5)

where we dropped the superscript (EI) and (EP) but they are assumed
implicitly. The results are in Fig. 1 (g) and (h). Skewness represents the
normalized asymmetry of a distribution such that the skewness is positive
(negative) if the right(left) tail is thicker. The observed
$\gamma_{v;c}^{\rm(EP)}$’s are significant for most exporters $c$, and
$\gamma_{v;c}^{\rm(EI)}$’s seem partly so. To check the significance, we
shuffle weights - logarithmic export values - among links to obtain a thousand
of realizations of weight-shuffled EI and EP network and measure their
skewness. In case of year 1980 presented in Fig. 1 (g) and (h), 140 among 161
$(=87\%)$ exporters with out-degrees not smaller than 10 find their true
skewness $\gamma_{v;c}^{\rm(EP)}$ outside the middle 90% range of the skewness
from randomized networks [Fig. 1 (h)]. 42 ($=26\%$) exporters have true
skewness $\gamma_{v;c}^{\rm(EI)}$ outside the middle 90% range of the random-
version skewness [Fig. 1 (g)]. See Appendix A2 and Fig. A2 for more details.

The skewness of local logarithmic exports tends to decrease, from positive to
negative, with the out-degree of the exporter $c$, more strongly for
$\gamma_{v;c}^{\rm(EP)}$ than for $\gamma_{v;c}^{\rm(EI)}$ [Fig. 1 (g) and
(h)]. Such correlation is robustly observed in the considered period, except
for the weak correlation between $\gamma_{v;c}^{\rm(EI)}$ and
$q_{c}^{\rm(EI)}$ in the middle eighties and nineties [Fig. A3].

When the distribution of a finite-ranged variable is much closer to its lower
(upper) bound, a thicker right (left) tail can be developed, yielding a
positive (negative) skewness. This can explain the observed negative
correlation between $\gamma_{v;c}$ and $q_{c}$; $P_{c}(v)$ moves to the right
on the $v$ axis if $c$ has a large out-degree, as captured by the growth of
the mean logarithmic export $\langle v\rangle_{c}$ with $q_{c}$ [Fig. 1
(a-d)]. Empirically $\gamma_{v;c}=\langle\tilde{v}^{3}\rangle_{c}$ is
negatively correlated with $\langle v\rangle_{c}$ as expected [Fig. A4]. It is
noteworthy that the export value of a country per destination/product grows
with the out-degree of the exporter, as it implies a super-linear growth of
the whole export of a country $V_{c}=\sum_{c^{\prime},p}V_{cc^{\prime}p}$ with
the out-degree $q_{c}^{\rm(EI)}$ and $q_{c}^{\rm(EP)}$. More developed or
richer countries may have many partners and export products Bhattacharya  _et
al._ (2008); Cha  _et al._ (2010); Almog  _et al._ (2015); Garlaschelli and
Loffredo (2004) and therefore our results point out different skewness of
export portfolio between developed and developing countries.

We have just argued that the skewness of local exports $\gamma_{v;c}$ are
influenced by the exporter $c$ via the $c$-dependent mean $\langle
v\rangle_{c}$. Individual local exports $v_{cc^{\prime}}$’s and $v_{cp}$’s
must be dependent also on the importer country $c^{\prime}$ and the product
category $p$ Lei  _et al._ (2015). Exporters have different sets of partners
and export products, possibly depending on the exporters’ out-degrees, which
might be influencing the skewness of local export distributions. In the next
sections, we study how the standardized logarithmic exports $\tilde{v}_{c\nu}$
vary with $\nu$, which will deepen our understanding of the statistics of
local exports and of the export strategies of countries in international
trade.

##  III Local in-degree distributions

What products and to what countries a country exports and how large its value
is may depend on the domestic industrial landscape of the exporter itself but
also on the characteristics of the partner importers and product categories
such as the economic development of and the social and political relationship
with the partner countries, and the technology level of products, and so on.
Among them, how many countries export to a country and how many countries
export a product category can be used as important information when an
exporter country constructs its product portfolio, as they represent the
popularity in international trade, reflecting the collective evaluation by
countries Lei  _et al._ (2015). They correspond to the in-degrees of a country
$c^{\prime}$ and of product $p$ in the EI and the EP network, defined as

$k_{\nu}=\sum_{c}A_{c\nu}=\sum_{c\in\mathtt{O}_{\nu}}1,$ (6)

with $\nu$ being $c^{\prime}$ or $p$, and $\mathtt{O}_{c^{\prime}/p}$ the set
of nodes (countries) sending outgoing links to $c^{\prime}$ or $p$. We take
the in-degree as the representative characteristic of an importer or a product
and investigate how the local logarithmic export value to/of a
destination/product relates to the in-degree of the latter. To this end, we
examine in this section the distributions of the in-degrees of the partner
importers and the exported products of given exporter, which we call local in-
degree distributions, and their correlations, if any, with local export
values.

Figure 2: Statistics of local in-degrees in the EI and the EP network. The in-
degree of a product or a country node means the number of exporter countries
exporting the product or to the country. (a, b) Distribution of the in-degrees
of (a) the country nodes in the EI network and (b) the product nodes in the EP
network receiving links from each exporter country $c$ in year 1980. Color
varies with the out-degree of $c$. Inset: The mean and standard deviation of
local in-degrees versus the out-degree of exporter. (c, d) The skewness of the
in-degrees of (c) the country in the EI network and (d) the product nodes in
the EP network receiving links from $c$ versus the out-degree of the exporter
$c$.

###  III.1 Local in-degree distribution and its moments

Like the local export distributions, the local in-degree distributions
$P_{c}^{\rm(EI)}(k)={1\over
q_{c}^{\rm(EI)}}\sum_{c^{\prime}\in\mathtt{I}_{c}^{\rm(EI)}}\delta(k_{c^{\prime}}^{\rm(EI)}-k)$
and $P_{c}^{\rm(EP)}(k)={1\over
q_{c}^{\rm(EP)}}\sum_{p\in\mathtt{I}_{c}^{\rm(EP)}}\delta(k_{p}^{\rm(EP)}-k)$
show a large variation with exporters $c$, related to the variation of the
mean local in-degree $\langle k\rangle_{c}$ [Fig. 2 (a) and (b)]. Differently
from the mean local logarithmic export, $\langle k\rangle_{c}$ decreases with
the out-degree of exporter. The standard deviations $\sigma_{k;c}^{\rm(EI)}$
and $\sigma_{k;c}^{\rm(EP)}$ increase with $q_{c}$. Such behaviors of the
first two moments and the whole distributions imply that hub exporters select
unpopular destinations/products as well as popular ones while non-hub
exporters tend to select popular partners and products only. In the context of
complex networks, hub exporters and popular destinations/products have links
to and from one another, forming a core, while non-hub exporters and unpopular
destinations/products tend to have links only to and from the nodes in the
core, viewed as periphery, which is characteristic of the core-periphery
structure Lee  _et al._ (2014).

The skewness of the local in-degree distribution $\gamma_{k;c}={1\over
q_{c}}\sum_{\nu\in\mathtt{I}_{c}}(k_{\nu}-\langle
k\rangle_{c})^{3}/\sigma_{k;c}^{3}$ increases with the out-degree of exporter
as shown in Fig. 2 (c) and (d). See Fig. A5 as well for its robustness against
time. As we discussed for the correlation of the skewness of local export and
the out-degree of exporter, the growth of $\gamma_{k;c}$, with $q_{c}$ may be
caused by the shift of local in-degree distributions in the $k$ axis to the
left as the out-degree of exporter increases. Also the distributions are
narrower for non-hub exporters, explaining the larger absolute values of
$\gamma_{k;c}$ for non-hub exporters than for hub ones given that cubic
fluctuations are normalized by the standard deviations in computing the
skewness.

###  III.2 Correlation of local exports and in-degrees at the aggregate level

Figure 3: Correlation between local in-degrees and local logarithmic exports
and their skewness. (a) The joint probability density function (PDF)
$P^{\rm(EI)}(\tilde{k},\tilde{v})=\sum_{c,c^{\prime}}A_{cc^{\prime}}\delta(\tilde{k}_{cc^{\prime}}-\tilde{k})\delta(\tilde{v}_{c\nu}-\tilde{v})/\sum_{c,c^{\prime}}A_{cc^{\prime}}$
of standardized local in-degrees and local logarithmic exports in the EI
network. The orientation of the red (high density) region to the upper right
indicates a positive correlation. (b) The joint PDF
$P^{\rm(EP)}(\tilde{k},\tilde{v})=\sum_{c,p}A_{cp}\delta(\tilde{k}_{cp}-\tilde{k})\delta(\tilde{v}_{cp}-\tilde{v})/\sum_{c,p}A_{cp}$.
(c) Scatter plot of the skewness of local logarithmic exports to different
destinations $\gamma_{v;c}^{\rm(EI)}$ versus the skewness
$\gamma_{k;c}^{\rm(EI)}$ of local in-degrees in the EI network. Color varies
with the out-degree of $c$ as in Fig. 2. (d) Scatter plot of the skewness of
local logarithmic exports of different products $\gamma_{v;c}^{\rm(EP)}$
versus the skewness $\gamma_{k;c}^{\rm(EP)}$ of local in-degrees in the EP
network.

At the aggregate level, the local export values of destinations/products and
their in-degrees are correlated. The Pearson correlation coefficient between
the aggregated standardized values
$\\{\tilde{v}_{c\nu}|c,\nu\in\mathtt{I}_{c}\\}$ and
$\\{\tilde{k}_{c\nu}|c,\nu\in\mathtt{I}_{c}\\}$ is 0.51 and 0.29 for EI and
EP, respectively, in year 1980. Here we introduced the standardized local in-
degrees $\tilde{k}$ for exporter $c$ defined by

$k_{\nu}=\langle k\rangle_{c}+\sigma_{k;c}\tilde{k}_{c\nu}.$ (7)

The joint distributions of the standardized local in-degrees and export in
Fig. 3 (a) and (b) also support their positive correlations.

If the local export value to/of a destination/product is strongly correlated
with the in-degree of the latter, then one can expect

$\tilde{v}_{c\nu}\simeq\tilde{k}_{c\nu}.$ (8)

Then the two skewness should be identical or at least similar
$\gamma_{v;c}=\langle\tilde{v}^{3}\rangle_{c}\simeq\langle\tilde{k}^{3}\rangle_{c}=\gamma_{k;c}$.
In Fig. 3 (c) and (d), however, one can see that they are quite different and
some of them have even different signs. Moreover, $\gamma_{v;c}^{\rm(EP)}$
shows a significant negative correlation with $\gamma_{k;c}^{\rm(EP)}$. Such
discrepancy between the two skewness shows that Eq. (8) needs correction. In
the next section, we present an improved version of Eq. (8) which relates the
two skewness correctly and leads us to understand at a deeper level the
skewness of local logarithmic exports.

##  IV Roles of in-degrees and fluctuations in local exports

###  IV.1 Decomposing local exports

Figure 4: Formulating local exports as functions of local in-degrees for given
exporter. (a) Standardized logarithmic exports
$\\{\tilde{v}_{cc^{\prime}}|c^{\prime}\in\mathtt{I}_{c}\\}$ to different
destinations $c^{\prime}$’s from a selected country $c$ versus their
standardized in-degrees
$\\{\tilde{k}_{cc^{\prime}}|c^{\prime}\in\mathtt{I}_{c}\\}$ in 1980. The slope
$\rho_{c}$ of the fitting line is equal to their Pearson correlation
coefficient. The deviation of $\tilde{v}_{cc^{\prime}}$ from the prediction by
the in-degree, $\rho_{c}\tilde{k}_{cc^{\prime}}$, is set to
$\sqrt{1-\rho^{2}}\xi_{cc^{\prime}}$, by which the random noise
$\xi_{cc^{\prime}}$ is defined. (b) $\\{\tilde{v}_{cp}\\}$ versus
$\\{\tilde{k}_{p}\\}$ for a selected country $c$ in 1980. (c) Pearson
correlation coefficients $\rho_{c}^{\rm(EI)}$ between the logarithmic exports
to different destinations from an exporter $c$ and their in-degrees versus the
out-degree of $c$. (d) Pearson correlation coefficients $\rho_{c}^{\rm(EP)}$
between the logarithmic exports of different products from an exporter $c$ and
their in-degrees versus the out-degree of $c$.

The standardized local exports to/of and in-degrees of destinations/products
are shown as examples for two selected countries in Fig. 4 (a) and (b). The
slope of the best fitting line of the data points
$\\{(\tilde{k}_{c\nu},\tilde{v}_{c\nu})|\nu\in\mathtt{I}_{c}\\}$ corresponds
to their Pearson correlation coefficient, which differs from exporter to
exporter. Also we find that the data points fluctuate considerably deviating
from the fitting line. Such fluctuation is different in magnitude and sign
depending on the target - destination/product - for given exporter. We take
the sum of all the factors contributing to local exports except for the in-
degrees as a random noise and represent the local export to/of a
destination/product as

$\tilde{v}_{c\nu}=\rho_{c}\tilde{k}_{c\nu}+\sqrt{1-{\rho_{c}}^{2}}\xi_{c\nu},$
(9)

where $\rho_{c}^{\rm(EI)}$ and $\rho_{c}^{\rm(EP)}$ are the Pearson
correlation coefficients between local exports and in-degrees given by

$\rho_{c}=\langle\tilde{k}\tilde{v}\rangle_{c}={1\over
q_{c}}\sum_{\nu\in\mathtt{I}_{c}}\tilde{k}_{\nu}\tilde{v}_{c\nu},$ (10)

and $\xi_{cc^{\prime}}$ and $\xi_{cp}$ are random noises of zero mean and unit
strength satisfying

$\displaystyle\langle\xi\rangle_{c}$ $\displaystyle={1\over
q_{c}}\sum_{\nu\in\mathtt{I}_{c}}\xi_{c\nu}=0,$
$\displaystyle\langle\tilde{k}\xi\rangle_{c}$ $\displaystyle={1\over
q_{c}}\sum_{\nu\in\mathtt{I}_{c}}\tilde{k}_{\nu}\xi_{c\nu}=0,$
$\displaystyle\langle\xi^{2}\rangle_{c}$ $\displaystyle={1\over
q_{c}}\sum_{\nu\in\mathtt{I}_{c}}\xi_{c\nu}^{2}=1.$ (11)

The relations in Eq. (11) are derived by using Eq. (10),
$\langle\tilde{k}\rangle_{c}=\langle\tilde{v}\rangle_{c}=0$, and
$\langle\tilde{k}^{2}\rangle_{c}=\langle\tilde{v}^{2}\rangle_{c}=1$. The
graphical meaning and measurment of $\rho$’s and $\xi$’s are presented in Fig.
4 (a). All the quantities in Eq. (9) can be measured empirically, which
advances greatly our study as will be presented in the remaining part. We
expect that the statistics of random noise $\xi_{c\nu}$ encodes the relevant
factors mainly stemming from the exporter $c$ while the standardized in-degree
$\tilde{k}$ reflect those from the destination/product. They are independent
of each other in that their covariance $\langle\tilde{k}\xi\rangle_{c}$
vanishes, but as we will show below, they have non-zero higher-order joint
cumulants.

From the dependence of $\rho_{c}$ on the out-degree $q_{c}$, more significant
in the EP than the EI network [Fig. 4 (c) and (d)], we see that non-hub
exporters determine their export values to/of destinations/products mainly by
referring to their own circumstances encoded by the noise term while the local
exports from hub exporters depend strongly on the popularity of
destinations/products.

###  IV.2 Joint third-order cumulants of in-degree and random noise

Using Eq. (9), one can evaluate the skewness of local exports as

$\displaystyle\gamma_{v;c}$
$\displaystyle=\rho_{c}^{3}\langle\tilde{k}\rangle_{c}+3{\rho_{c}}^{2}\sqrt{1-{\rho_{c}}^{2}}\langle\tilde{k}^{2}\xi\rangle_{c}$
$\displaystyle+3(1-{\rho_{c}}^{2})\rho_{c}\langle\tilde{k}\xi^{2}\rangle_{c}+(1-{\rho_{c}}^{2})^{3/2}\langle\xi^{3}\rangle_{c},$
(12)

where four third-order cumulants,
$\langle\tilde{k}^{3}\rangle_{c},\langle\tilde{k}^{2}\xi\rangle_{c},\langle\tilde{k}\xi^{2}\rangle_{c}$,
and $\langle\xi^{3}\rangle_{c}$ can be measured empirically and capture
asymmetry in the distribution of $\tilde{k}$’s and $\xi$’s. The four cumulants
therefore represent four different modes to generate asymmetry in the local
export distribution. We have seen that
$\gamma_{k;c}=\langle\tilde{k}^{3}\rangle_{c}$ depends on $q_{c}$ in the
opposite way to how $\langle\tilde{v}^{3}\rangle_{c}$ does, which we now
understand can happen if $\rho_{c}$ is sufficiently small as
$\langle\tilde{k}^{3}\rangle_{c}$ is weighted by $\rho_{c}^{3}$ in Eq. (12).
$\tilde{k}$ and $\xi$ have zero covariance $\langle\tilde{k}\xi\rangle_{c}=0$
but we find empirically that they have non-zero third-order joint cumulants
$\langle\tilde{k}^{2}\xi\rangle_{c}$ and $\langle\tilde{k}\xi^{2}\rangle_{c}$
[Fig. 5 and A6]. In principle, $\langle\tilde{k}^{2}\xi\rangle_{c}$ is
positive (negative) when the data points of $(\tilde{k},\tilde{v})$ tend to be
located above (below) the fitting line $\tilde{v}^{\rm(fit)}=\rho\tilde{k}$
when $\tilde{k}^{2}$ is large. If $(\tilde{k},\tilde{v})$ tends to spread more
(less) as $\tilde{k}$ increases, the joint cumulant
$\langle\tilde{k}\xi^{2}\rangle$ will be positive (negative).

Figure 5: Third-order cumulants of local in-degrees and random noises and
their relative magnitudes. (a-f) are obtained from the data of year 1980. (a)
The skewness of random noises $\langle\xi^{3}\rangle_{c}^{\rm(EI)}$ in the
local exports from each exporter $c$ to different destinations versus the out-
degree $q_{c}^{\rm(EI)}$. (b) The skewness of random noises
$\langle\xi^{3}\rangle_{c}^{\rm(EP)}$ in the local exports from exporter $c$
of different products versus the out-degree $q_{c}^{\rm(EP)}$. (c) Third-order
joint cumulants $\langle\tilde{k}^{2}\xi\rangle_{c}^{\rm(EI)}$ and
$\langle\tilde{k}\xi^{2}\rangle_{c}^{\rm(EI)}$ versus $q_{c}^{\rm(EI)}$. (d)
$\langle\tilde{k}^{2}\xi\rangle_{c}^{\rm(EP)}$ and
$\langle\tilde{k}\xi^{2}\rangle_{c}^{\rm(EP)}$ versus $q_{c}^{\rm(EP)}$. (e)
Relative magnitudes of the four third-order cumulants of local in-degrees and
random noises in the EI network. (f) The same plot as (e) for the EP network.
(g) The averages of the relative magnitudes of the third-order cumulants in
the EI networks as functions of time. Errorbars represents the standard
deviations. (h) The same plot as (g) for the EP networks.

The skewness of the random noise $\xi$ is expected to make a considerable
contribution to the whole skewness of local export, for the magnitudes and the
$q_{c}$-dependence of $\langle\xi^{3}\rangle_{c}$ are empirically close to
those of $\langle\tilde{v}^{3}\rangle$ [Fig. 5 (a) and (b)]. It suggests that
the statistics of the random noise $\xi$, whose nature mainly depends on the
exporter itself, may govern the statistics of local exports. However, the
joint cumulants $\langle\tilde{k}^{2}\xi\rangle_{c}$ and
$\langle\tilde{k}\xi^{2}\rangle_{c}$ are not negligible [Fig. 5 (c) and (d)]
but even make bigger contributions to the whole skewness
$\langle\tilde{v}^{3}\rangle_{c}$ for some exporters $c$ [Fig. 5 (e) and (f)].
We measure the contributions of the four third-order cumulants to the skewness
of local exports by the ratio of the absolute values of the four terms in Eq.
(12),
$|\rho^{3}\langle\tilde{k}^{3}\rangle_{c}|,|3{\rho_{c}}^{2}\sqrt{1-{\rho_{c}}^{2}}\langle\tilde{k}^{2}\xi\rangle_{c}|,|3(1-{\rho_{c}}^{2})\rho_{c}\langle\tilde{k}\xi^{2}\rangle_{c}|,|(1-{\rho_{c}}^{2})^{3/2}\langle\xi^{3}\rangle_{c}|$
to their sum
$|\rho_{c}^{3}\langle\tilde{k}\rangle_{c}|+|3{\rho_{c}}^{2}\sqrt{1-{\rho_{c}}^{2}}\langle\tilde{k}^{2}\xi\rangle_{c}|+|3(1-{\rho_{c}}^{2})\rho_{c}\langle\tilde{k}\xi^{2}\rangle_{c}|+|(1-{\rho_{c}}^{2})^{3/2}\langle\xi^{3}\rangle_{c}|$.
Therefore the sum of those four contributions is one ($100\%$) for each
exporter. As expected, the contribution of the skewness of in-degrees is not
significant, $9.0\pm 9.0\%$ to $\gamma_{v;c}^{\rm(EI)}$ and $1.0\pm 1.7\%$ to
$\gamma_{v;c}^{\rm(EP)}$ in year 1980 [Fig. 5 (e) and (f)]. On the other hand,
the contribution of the skewness of the random noise
$\langle\xi^{3}\rangle_{c}$ is $31\pm 25\%$ and $73\pm 26\%$ to
$\gamma_{v;c}^{\rm(EI)}$ and $\gamma_{v;c}^{\rm(EP)}$, respectively.
Interestingly, the contribution of a joint third-order cumulant
$\langle\tilde{k}\xi^{2}\rangle_{c}$ makes a large contribution, $42\pm 20\%$
and $20\pm 19\%$ to $\gamma_{v;c}^{\rm(EI)}$ and $\gamma_{v;c}^{\rm(EP)}$
respectively in year 1980. Across time [Fig. 5 (g) and (h)],
$\langle\xi^{3}\rangle_{c}$ contributes about $32\%$ and
$\langle\tilde{k}\xi^{2}\rangle$ $34\%$ to $\gamma_{v;c}^{\rm(EI)}$. For
$\gamma_{v;c}^{\rm(EP)}$, time variations are stronger than for
$\gamma_{v;c}^{\rm(EI)}$; $\langle\xi^{3}\rangle_{c}$ makes the leading
contribution $40\sim 84\%$ and the next-leading contribution is made by
$\langle\tilde{k}\xi^{2}\rangle_{c}$, ranging $10\%$ to $48\%$.

The empirically measured third-order cumulants allows us to better understand
the skewness of local exports. $\langle\xi^{3}\rangle_{c}$ and
$\langle\tilde{k}\xi^{2}\rangle_{c}$ together contribute $66\pm 4\%$ and
$88\pm 5\%$ to $\gamma_{v;c}^{\rm(EI)}$ and $\gamma_{v;c}^{\rm(EP)}$,
respectively. They tend to decrease with the out-degree of exporter and are
thus responsible for the out-degree-dependence of $\gamma_{v;c}$. As such, it
is mainly via the two cumulants $\langle\tilde{k}\xi^{2}\rangle_{c}$ and
$\langle\xi^{3}\rangle_{c}$ that an exporter and its partners/products affect
the skewness of their local logarithmic exports.

Moreover, the analysis based on Eq. (12) elucidates the nature of the export
strategies of individual countries. As shown in Fig. 5 (c) and (d),
$\langle\tilde{k}\xi^{2}\rangle_{c}$ is negative for most exporters. In 1980,
126 among 161 countries have $\langle\tilde{k}\xi^{2}\rangle_{c}^{\rm(EI)}\leq
0$ and 152 countries have $\langle\tilde{k}\xi^{2}\rangle_{c}^{\rm(EP)}\leq
0$. Such negative $\langle\tilde{k}\xi^{2}\rangle_{c}$ is observed for many
countries throughout the considered period [Fig. A6], and means that the
spread of local exports is weaker for larger $\tilde{k}$. As such, a country
determines its export value to a popular destination or of a popular product
by more referring to its in-degree than to/of an unpopular
destination/product, for the latter of which the circumstances of the exporter
itself can be important. Also we find that
$\langle\tilde{k}^{2}\xi\rangle_{c}$ is positive for many countries, which
implies that the local export to/of quite a popular or quite an unpopular
destination/product tends to be larger than the average. See Fig. A7 for
graphical examples. All these results illuminate how the export values of an
exporter country are allocated to destinations/products, and combined with the
different connection patterns between hub and non-hub exporters, illuminate
the universal and the non-universal aspects of the strategies of countries in
international trade.

##  V Summary and Discussion

We have investigated how a country divides its export value over destinations
and products by analyzing the local logarithmic export distributions. They are
close to the Gaussian distribution, but exhibit significant skewness. We
studied the origin and the implications of such skewness for the export
strategies of countries in international trade. The local export per
destination/product is not constant but large for hub exporters and small for
non-hub exporters, which contributes to the variation of the skewness of local
exports with exporters. We presented the formula relating the local export
value to/of a destination/product to the in-degree of the latter and a random
noise, which are measurable empirically and thereby allow us to decipher the
origin of the skewness of local exports and to understand quantitatively the
export strategies of countries. We showed that the statistics of the random
noise and its third-order correlation with in-degrees are responsible for the
skewness of local exports and that the local exports to/of more popular
destination/products are governed more by their in-degrees than the export
to/of less popular ones are.

Our study aimed at understanding the connection topology and the local link-
weight distributions in weighted networks, here the world trade networks, and
therefore our methods and results can be applied to other weighted complex
networks including natural, information, social, and artificial networks. The
triangular relationships involving exporters, importers, and products can be
represented by hypergraphs Ghoshal  _et al._ (2009), and generalizing our
analysis for them is in progress.

###### Acknowledgements.

We thank Su Do Yi for valuable comments. This work was supported by the
National Research Foundation of Korea (NRF) grants funded by the Korean
Government (No. 2019R1A2C1003486) and a KIAS Individual Grant (No. CG079901)
at Korea Institute for Advanced Study.

##  Appendix A1 Data-sets

In this study, the NBER-UN data-set in the period $1962\leq t\leq 1999$
Feenstra  _et al._ (2005) and the World Trade Organization (WTO) data-set for
$2000\leq t\leq 2018$ wto  are compiled to obtain the set of export values
$\\{V_{cc^{\prime}p}(t)|c,c^{\prime}\in\mathtt{C}(t),p\in\mathtt{P}(t)\\}$ for
$1962\leq t\leq 2018$. Here $\mathtt{C}(t)$ is the set of countries exporting
or importing at least one product category and $\mathtt{P}(t)$ is the set of
trade product categories exported by at least one country in year $t$. To
classify trade products, the 4-digit standard international trade
classification (SITC) sit (2006) is used in Ref. Feenstra  _et al._ (2005).
The WTO data-sets wto  are based on the Harmonized System (HS) hs (2) as the
product classification scheme, and we also select its 4-digit version to make
the total number of product categories similar to that of the 4-digit SITC.
These two schemes of product classification can be mapped to each other uns .

We remark that the data-sets in the period 1984 to 1999 contain some virtual
product categories inserted to fill some gaps in the total values of trade
across different levels of classification Feenstra  _et al._ (2005). The NBER-
UN data-sets record trade over 1000 US dollars. Yet these differences
depending on data source and time period do not affect the main results of the
present study.

Figure A1: Trade networks in 1970 and their global degree and link-weight
distributions.  (a) A subnetwork of the EI network in 1970, consisting of 30
countries having the largest (out$+$in)-degrees connected by 141 links of
standardized local weight (logarithmic export value) larger than 1, i.e.,
$\tilde{v}>1$. (b) A subnetwork of the EP network of 15 countries having the
largest out-degrees and 15 products having the largest in-degrees connected by
133 links of $\tilde{v}>1$. (c) Global distribution of the logarithm of
exports
$P^{\rm(EI)}(v)={\sum_{c,c^{\prime}}A_{cc^{\prime}}\delta(v_{cc^{\prime}}-v)\over\sum_{c,c^{\prime}}A_{cc^{\prime}}}$
in the EI network with $A_{cc^{\prime}}$ the adjacency matrix. (d) Global
distribution of the logarithm of exports $P^{\rm(EP)}(v)$ in the EP network.
(e) Global out-degree distribution
$P^{\rm(EI)}(q)={\sum_{c}\delta(q_{c}^{\rm(EI)}-q)\over\sum_{c}1}$ and the
plot of the rank $R_{c}^{\rm(EI)}$ versus the out-degree $q_{c}^{\rm(EI)}$ of
each exporter $c$ in the EI network. (f) Global out-degree distribution
$P^{\rm(EP)}(q)$ and the plot of the rank $R_{c}^{\rm(EP)}$ versus the out-
degree $q_{c}^{\rm(EP)}$ of exporter $c$ in the EP network. (g) Global in-
degree distribution
$P^{\rm(EI)}(k)={\sum_{c}\delta(k_{c}^{\rm(EI)}-k)\over\sum_{c}1}$ and the
plot of the rank $R_{c}^{\rm(EI)}$ versus the in-degree $k_{c}^{\rm(EI)}$ of
importer $c$ in the EI network. (f) Global in-degree distribution
$P^{\rm(EP)}(k)={\sum_{p}\delta(k_{p}^{\rm(EP)}-k)\over\sum_{p}1}$ and the
plot of the rank $R_{p}^{\rm(EP)}$ versus the in-degree $k_{p}^{\rm(EP)}$ of
product $p$ in the EP network.

##  Appendix A2 Significance of skewness and its correlation with out-degree

Figure A2: Significance of skewness of local logarithmic exports. (a) The
fraction of exporters whose true skewness remain in the middle $x(\%)$ range
of the random-version skewness in 1980. Squares and circles are for the EI and
the EP network. (b) The same plot as (a) in year 2010.

To check the significance of the skewness of local logarithmic exports, we
compare the true skewness with the values of the skewness from
$\mathcal{N}=1000$ realizations of weight-shuffled EI and EP networks. The
out-degrees of each exporter $q_{c}^{\rm(EI)}$ and $q_{c}^{\rm(EP)}$ do not
change by this randomization. Then, for each exporter $c$ we are given
$\mathcal{N}$ values of the random-version skewness
$\\{\gamma_{v;c}^{{\rm(EI;random)}(\ell)}|\ell=1,2,\ldots,\mathcal{N}\\}$ in
the EI network and
$\\{\gamma_{v;c}^{{\rm(EP;random)}(\ell)}|\ell=1,2,\ldots,\mathcal{N}\\}$ in
the EP network, and we check whether the true skewness
$\gamma_{v;c}^{\rm(EI)}$ and $\gamma_{v;c}^{\rm(EP)}$ belong to the middle
$x(\%)$ range of those random-version skewness values. Repeating this for all
exporters, we obtain the fraction $y(\%)$ of exporters whose true skewness are
included in their respective middle $x(\%)$ ranges of the random-version
skewness. If $y$ as a function of $x$ remains close to $x$, then one can say
that true skewness is not distinguished from the random-version ones and
therefore the true skewness is not a significant one. On the other hand, if
$y$ is far smaller than $x$, then it means that the true skewness tends to be
outside the range expected in the randomized networks, demonstrating the
significance of the observed skewness in a sense that it is hardly expected in
randomized networks.

We found that the latter is the case. In Fig. A2, we present the plot of $y$
versus $x$ for the EI and the EP network in year 1980 and 2010. $y$ is far
below the line $y=x$ in case of the EP network. Also for the EI network, $y$
is smaller than $x$ though the difference is not so impressive as for the EP
network. These results suggest that the empirically-observed skewness of local
logarithmic exports are significant.

Figure A3: Time-dependent correlation of the skewness with the average of
local logarithmic exports and with the out-degree of exporter. The Pearson
correlation coefficient between the skewness of local logarithmic exports
$\gamma_{v;c}^{\rm(EI,EP)}$ and the out-degree $q_{c}^{\rm(EI,EP)}$ of
exporter or the average local logarithmic export $\langle
v\rangle_{c}^{\rm(EI,EP)}$ as a function of time $t$.  Figure A4: Plot of the
skewness versus the average of local logarithmic exports. (a) Plot of the
skewness of local logarithmic exports $\gamma_{v;c}^{\rm(EI)}$ versus the
average of local logarithmic exports $\langle v\rangle_{c}^{\rm(EI)}$ in the
EI network in year 1980. (b) The same plot as (a) in the EP network.  Figure
A5: Correlation of skewness of local in-degrees with the out-degree of
exporter and the average local in-degree. The Pearson correlation coefficient
between the skewness of local in-degrees $\gamma_{k;c}^{\rm(EI,EP)}$ and the
out-degree $q_{c}^{\rm(EI,EP)}$ of exporter or the average local in-degree
$\langle k\rangle_{c}^{\rm(EI,EP)}$ as a function of time $t$. Correlations
are observed robustly in the studied period, except for the correlations of
$\gamma_{k;c}^{\rm(EP)}$ with $q_{c}^{\rm(EP)}$ and with $\langle
k\rangle_{c}^{\rm(EP)}$ which change sign in the period 1984-1999 presumably
owing to the virtual product categories added artificially [See Appendix A1].
Figure A6: Robustness against time of the sign of the correlation coefficient
and third-order cumulants.  (a) The portions of exporters showing a non-
negative correlation coefficient between local logarithmic exports and local
in-degrees ($\rho_{c}^{\rm(EI)}\geq 0$), a non-negative cumulant
$\langle\tilde{k}^{2}\xi\rangle_{c}^{\rm(EI)}\geq 0$, and another non-positive
cumulant $\langle\tilde{k}\xi^{2}\rangle_{c}^{\rm(EI)}\neq 0$ of local in-
degrees and random noises are shown respectively for the EI network as
functions of time. (b) The same plots as (a) for the EP network.  Figure A7:
Schematic distribution of $\\{(\tilde{k},\tilde{v})\\}$ corresponding to
positive or negative third-order cumulants. Each panel shows an example of the
distribution of the data points of standardized local exports and in-degrees
$\\{(\tilde{k},\tilde{v})\\}$ exhibiting a given positive or negative
cumulant. The solid line represents $\tilde{v}=\rho\tilde{k}$ with $\rho=1$
for simplicity. Therefore the deviation of the data points from the line
represents the random nooise $\xi$. Note that
$\langle\tilde{k}\rangle=\langle\xi\rangle=0$ and
$\langle\tilde{k}^{2}\rangle=\langle\xi^{2}\rangle=1$.

## References

  * Hidalgo  _et al._ (2007) C. A. Hidalgo, B. Klinger, A.-L. Barabási, and R. Hausmann, Science 317, 482 (2007).
  * Tacchella  _et al._ (2018) A. Tacchella, D. Mazzilli, and L. Pietronero, Nat. Phys. 14, 861 (2018).
  * Chenery and Taylor (1968) H. B. Chenery and L. Taylor, Rev. Econ. Stat. 50, 391 (1968).
  * Kuznets (1971) S. Kuznets,  _Economic Growth of Nations: Total Output and Production Structure_ (Belknap Press of Harvard University Press, 1971).
  * Lei  _et al._ (2015) H. Lei, Y. Chen, R. Li, D. He, and J. Zhang, PLOS ONE 10, 1 (2015).
  * Choi and Lee (2019) S.-G. Choi and D.-S. Lee, Sci. Rep. 9, 15871 (2019).
  * Hummels and Klenow (2005) D. Hummels and P. J. Klenow, Am. Econ. Rev. 95, 704 (2005).
  * Serrano and Boguñá (2003) M. A. Serrano and M. Boguñá, Phys. Rev. E 68, 015101 (2003).
  * Garlaschelli and Loffredo (2004) D. Garlaschelli and M. I. Loffredo, Phys. Rev. Lett. 93, 188701 (2004).
  * Fagiolo  _et al._ (2009) G. Fagiolo, J. Reyes, and S. Schiavo, Phys. Rev. E 79, 036115 (2009).
  * Cha  _et al._ (2010) M.-Y. Cha, J. W. Lee, and D.-S. Lee, J. Korean Phys. Soc. 56, 998 (2010).
  * Almog  _et al._ (2015) A. Almog, T. Squartini, and D. Garlaschelli, New J. Phys. 17, 013009 (2015).
  * Oh and Lee (2017) C.-Y. Oh and D.-S. Lee, Phys. Rev. E 95, 052319 (2017).
  * Hidalgo and Hausmann (2009) C. A. Hidalgo and R. Hausmann, Proc. Natl. Acad. Sci. U.S.A. 106, 10570 (2009).
  * Tacchella  _et al._ (2012) A. Tacchella, M. Cristelli, G. Caldarelli, A. Gabrielli, and L. Pietronero, Sci. Rep. 2, 723 (2012).
  * Barbier and Lee (2017) M. Barbier and D.-S. Lee, J. Stat. Mech. , 123403 (2017).
  * Feenstra  _et al._ (2005) R. C. Feenstra, R. E. Lipsey, H. Deng, A. C. Ma, and H. Mo,  _World Trade Flows: 1962-2000_, Working Paper 11040 (National Bureau of Economic Research, 2005).
  * (18) Http://data.wto.org (2020).
  * Bhattacharya  _et al._ (2008) K. Bhattacharya, G. Mukherjee, J. Saramäki, K. Kaski, and S. S. Manna, J. Stat. Mech. 2008, P02002 (2008).
  * Lee  _et al._ (2014) S. H. Lee, M. Cucuringu, and M. A. Porter, Phys. Rev. E 89, 032810 (2014).
  * Ghoshal  _et al._ (2009) G. Ghoshal, V. Zlatić, G. Caldarelli, and M. E. J. Newman, Phys. Rev. E 79, 066118 (2009).
  * sit (2006) _Standard International Trade Classification Revision 4_ (United Nation Publication, 2006).
  * hs (2) Http://www.wcoomd.org/en/topics/nomenclature/instrument-and-tools/hs-nomenclature-2017-edition.aspx.
  * (24) Https://unstats.un.org/unsd/tradekb/.


information was obtained from private communication with Oliver Gutsche.) The
0.9 (1.4) US dollar cents assumption puts the annual cost of simulation for
CMS in the range of 3.5-6.2 (5.5-10) million US dollars, half of it spent on
executing the Geant4 module. A corollary to this discussion is that
improvements of 1%, 10%, and 35% in the time performance of the Geant4 toolkit
would render 50-80k, 500-800k, 1.8-2.8M US dollars per year of savings to CMS.
Improvements on the order of a 2-5 speed-up factor, as targeted by current R&D
efforts (GeantV [78]), would yield savings on the order of 2-4 (3-6) million
US dollars per year. An important question, rarely addressed even by modern
experiments at the time of detector design and technology selection, is
related to the added costs to the detector construction, commissioning, and
operations that comes from detector choices that maximize physics output in
exchange for expensive and time-consuming simulation and reconstruction
operations.

It is important to mention that the LHC experiments expect their computing
needs to increase by a factor of 10 to 100 in the High-Luminosity LHC (HL-LHC)
era, depending on the solutions developed to face simulation, pile-up, and
reconstruction challenges arising from the high-luminosity environment. In
principle, reconstruction would take a larger fraction of the computing
resources during the HL-LHC era, since the CPU time consumption is predicted
to increase exponentially with the number of pile-up events. However, while
simulation code is highly optimized and offers few non-revolutionary time
performance improvement opportunities, the reconstruction code under
development for the upgraded or new HL-LHC sub-detector systems still offers
low hanging fruit to exploit, at least in the case of CMS. Consequently, a
significant improvement in simulation computing performance is a need in
present times of flat budgets, and so are the research efforts with that goal
in mind.

The Geant4 Collaboration has gone to great lengths to improve the toolkit
computing performance during the last few years, as code was reviewed and
optimized. In 2013, the introduction of event-level multithreading
capabilities in Geant4 brought significant memory savings, as illustrated in
Fig. 55, which shows the CPU time (top) and memory (bottom) consumption ratios
of a CMS standalone simulation application (outside of the CMS software
framework), based on Geant4 version 10.1.p02, executed for 5 GeV electrons in
multithreaded and sequential modes. While time performance does not improve,
deviating from perfect scaling by approximately 10% when executed on 30 cores,
memory consumption improves significantly with 170 MB used in the first event,
and only 30 MB per event used by each additional thread [79].

Figure 55: CPU time (top) and memory (bottom) ratios for a CMS standalone
simulation application (outside of the CMS software framework) based on Geant4
version 10.1.p02 executed for 5 GeV electrons in multithreaded and sequential
modes [79].

Figs. 56, 57 show the percentage change in CPU time performance taking Geant4
version 10.0 as a reference, starting with Geant4 version 9.4.p02 (2010) and
ending with Geant4 version 10.2 (2015) for the standalone CMS application and
a simple calorimeter configuration made of Cu-Scintillator in a 4 Tesla
magnetic field [79]. The study is performed for 50 GeV $e^{-}$, $\pi^{-}$, and
protons, as well as for $H\rightarrow ZZ$ events. In average, the time
performance improvement through the life of the LHC experiments (2010-2015) is
of the order of 35%. All tests were performed on the same hardware (AMC
Opteron 6128 HE @ 2 GHz) using the same operating system. Remarkably, the
percentage time performance improvement during the period of time shown in the
plots is in the double digits, even as the physics models were improved
significantly for accuracy, something that typically comes associated with a
time performance penalty.

Figure 56: Percentage change in CPU time performance with respect to Geant4
version 10.0, starting with Geant4 version 9.4.p02 (2010) and ending with
Geant4 version 10.2 (2015), for a standalone CMS application (outside of the
CMS software framework) and a simple calorimeter configuration made of Cu-
Scintillator in a 4 Tesla magnetic field. The study is performed for 50 GeV
$e^{-}$ and $\pi^{-}$ [79].

Figure 57: Percentage change in CPU time performance with respect to Geant4
version 10.0, starting with Geant4 version 9.4.p02 (2010) and ending with
Geant4 version 10.2 (2015), for a standalone CMS application (outside of the
CMS software framework) and a simple calorimeter configuration made of Cu-
Scintillator in a 4 Tesla magnetic field. The study is performed for protons
and $H\rightarrow ZZ$ events [79].

The cost of simulation presented before refers only to the purchase, life-
cycle, operations, and maintenance of the computing resources allocated to the
task. The values do not include the design, development, validation,
operation, and support of the simulation tool-kits, such as Geant4, or of the
simulation software in the experiments. As a reference, during its 22 years of
existence, the person-power investment in Geant4 has totaled more than 500
person-years, equivalent to about 100 M US dollars, including fringe benefits
and overhead. Additional investments of the same order have been made by the
major $21^{st}$ century HEP experiments on detector specific simulation and
reconstruction software. An interesting corollary is that the cost of the
physics software amounts to a significant fraction of the cost of the
detectors.

It would be an interesting exercise to estimate the cost of running a modern
HEP experiment with and without efficient simulation tool-kits and full
simulation software. The truth is that the experiments as we know them today
would simply not exist without these tools. How much physics would be lost to
a deficient detector design or poor optimization? How would the design and
operation of systems such as data acquisition, distribution, storage, and
analysis workflows be affected? How accurate, efficient, and fast would
reconstruction algorithms, calibration and analysis methods be? How much
person-power and how many years of delay in delivering scientific publications
would it cost to reproduce, without good quality simulation, the level of
accuracy in physics measurements achieved by modern experiments? Is it
possible at all to deliver physics of the quality and accuracy we produce
today without simulation?

##  8 The Future

The accuracy of the simulation software developed for the current generation
of high-energy physics detectors, coupled with the speed of contemporary
computers, has enabled the experiments to perform tasks that scientists could
only have dreamed of before. Simulation helps physicists to design and
optimize detectors for best physics performance, stress-test the computing
infrastructure, program data reconstruction algorithms that perform almost as
in design specifications at the begining of the experiment run, develop data-
driven techniques for calibration and physics analysis, and produce data
samples with the properties predicted by many candidate theories to describe
currently unexplained physical phenomena.

Modern HEP experiments generate and handle an enormous amount of real and
simulated data. For its size and complexity, these data has earned a place in
the world of what is known as Big Data. The experiments at the CERN Large
Hadron Collider (LHC) have produced, reconstructed, stored, transferred, and
analyzed tens of billion of simulated events during the first two runs.
According to Ref. [80], the amount of data collected and stored by the LHC
experiments through the end of 2013 was of the order of 15 PB/year, not so far
from the 180 PB/year uploaded to Facebook, the 98 PB of data in the Google
search index, or the 15 PB/year in videos uploaded to YouTube. Integrated on
time over the last two decades, the cost of simulation and reconstruction in
large modern HEP experiments exceeded the one hundred million dollars mark.

The high instantaneous luminosity required at the LHC experiments, needed to
reach the 3000 fb-1 integrated luminosity milestone associated with the high-
luminosity LHC physics program, will tax heavily the performance of the
reconstruction algorithms. Through the end of the 2030’s, the experiments
expect to collect 150 times more data than in Run 1. The 50 PB of raw data
produced in 2016 will grow to approximately 600 PB in 2026 while the CPU needs
will increase by a factor of about 60. The exact numbers will depend on the
approximations and the loss of information that the experiments are willing to
tolerate to keep computing performance within the limits established by the
available resources. Thus, the effort to improve the computing performance of
the simulation and reconstruction software requires immediate attention, in
order to restrain the increasing demand of computing power within the limits
of flat budgets. Although transistor density growth is more or less keeping up
with Moore’s law, doubling every couple of years, clock speed has been flat
since approximately 2003. Consequently, solutions must be found elsewhere,
leveraging the core count growth in multicore machines, using new generation
coprocessors, and re-engineering code under new programming paradigms based on
concurrency and parallel programming. Coprocessors, or accelerators,
specialize in operations such as floating point math or graphics. A hybrid
computing model would allow to share work across a mixture of computers with
different architectures. Each processor type could be used to perform
different tasks depending on its nature.

In parallel, experiments are transforming their software frameworks to support
event multithreading and task-level parallelization. In the specific case of
Geant4, the release of the first version with multithreading capability in
2013 allowed significant savings in memory, although not in time performance.
For the latter, expert teams are invested in R&D programs to explore the
potential of multithread track-level (particle-level) parallelization,
improved instruction pipelining, data locality, and vectorization for single
instruction multiple data. One example is the GeantV [78] project to develop
the next generation detector simulation toolkit, with a goal set to achieve a
speedup factor of 2 to 5 with respect to Geant4, while enhancing the physics
accuracy of the code and offering fast simulation options that include machine
learning techniques for fast and precise tuning.

Breakthroughs in the design of simulation and reconstruction code, exploiting
the benefits of fine granularity parallelism in applications running in modern
computer architectures, will be essential to address the software and
computing challenges faced by the HEP experiments of the 21st century.

## Acknowledgements

Funding: Fermilab is operated by Fermi Research Alliance, LLC under Contract
No. DE-AC02-07CH11359 with the U.S. Department of Energy.

I would like to express my gratitude to the many people I have consulted on
the material to include in this review. First and foremost, I am eternally
grateful to John Harvey, former leader of the software group at CERN, who
initially asked me a few questions about the impact of simulation in
experimental particle physics, eventually encouraged me to write this article,
and was the first person to read and comment on a draft. I am also grateful to
my Geant4 co-collaborators who built on previous Geant experience and
developed the most widespread and successful detector simulation toolkit the
HEP field has ever known. Their work has yet to receive the recognition and
support it deserves. I want to thank Federico Carminati, Guenther Dissertori,
and Paris Sphicas for reading and commenting on a draft. Thanks to my ATLAS
colleagues John Chapman, Andrea Dotti, Zach Marshall, and Ariel Schwartzman
for pointing me to outstanding ATLAS material on simulation, test beams, jets,
and missing transverse energy. In order of topic appearance in the paper, I
would like to thank Julia Yarba and Alberto Ribon for pointing me to specific
Geant4 physics validation plots and associated information, Sunanda Banerjee
for useful discussions and the figures with the CMS test beam and single
particle response comparisons, Mike Tartaglia and Harrison Prosper for trying
hard to recall details of the D0 test beam experiments, my many D0 colleagues
with whom we navigated the difficult waters of producing high quality physics
measurements with neither a solenoidal field nor fast-enough high-quality
simulation software during Run 1, Liz Sexton-Kennedy and Robert Harris for
discussions on the CDF simulation software and first inclusive jet cross
section results, Soon Young Jun for providing precise information and
documentation on the CDF Monte Carlo tuning effort and for the Geant4
computing performance plots, Marjorie Shapiro for private communication on the
CDF QFL fast simulation software, my CMS colleagues with whom we experienced
the thrill of developing outstanding simulation code and producing the most
precise physics measurements ever in a hadron collider, Kevin Burkett for
pointing me to CDF publications data statistics, Oliver Gutsche for providing
most of the information for the cost evaluation of the simulation operation in
CMS. Last but not least, I want to thank Krzysztof Genser and the rest of the
members of my Fermilab Physics and Detector Simulation group (PDS) for their
hard work in the area of simulation software research, development and
support, as well as the Fermilab Scientific Computing Division and the US
Department of Energy for their continued support of the Fermilab Geant
operations and research programs.

## References

  * [1] S. Agostineli et al. (Geant4 Collaboration), Gean4-a simulation toolkit, Nucl. Instrum. Meth. A 506 (2003) 250–303.  doi:10.1016/S0168-9002(03)01368-8.
  * [2] J. Allison et al. (Geant4 Collaboration), Recent developments in Geant4, Nucl. Instrum. Meth. A 835 (2016) 186–225.  doi:10.1016/j.nima.2016.06.125.
  * [3] R. Ford, W. Nelson, The EGS Code System - Version 3, Stanford Linear Accelerator Center Report SLAC-210.
  * [4] R. Brun, F. Bruyant, M. Maire, A. McPherson, P. Zanarini, Geant3 CERN-DD-EE-84-1.
  * [5] The ATLAS Collaboration, Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC, Phys. Lett. B 716 (2012) 1–29.  doi:10.1016/j.physletb.2012.08.020.
  * [6] The CMS Collaboration, Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC, Phys. Lett. B 716 (2012) 30–61.  doi:10.1016/j.physletb.2012.08.021.
  * [7] D.J. Lange et al. for the CMS Collaboration, Upgrades for the CMS full simulation, J. Phys. Conf. Ser. 608 (2015) 012056.
  * [8] S. Banerjee et al. for the CMS Collaboration, Validation and tuning of the CMS simulation, J. Phys. Conf. Ser. 331 (2011) 032015.
  * [9] V. Daniel Elvira for the CMS Collaboration, Readiness of the CMS detector simulation, NSS/MIC 2007 C07-10-28 (2007) 2081–2085.  doi:10.1109/NSSMIC.2007.4436563.
  * [10] P. J. Clark for the ATLAS Collaboration, The ATLAS Detector Simulation, Nucl. Phys. B - Proceedings Supplements 215 (2011) 85–88.  doi:10.1016/j.nuclphysbps.2011.03.142.
  * [11] Z. Marshall for the ATLAS Collaboration, Validation and performance studies for the ATLAS simulation, J. Phys. Conf. Ser. 219 (2010) 032016.
  * [12] A. Rimoldi for the ATLAS Collaboration, The ATLAS detector simulation application, Nucl. Phys. B - Proceedings Supplements 172 (2007) 49–52.  doi:10.1016/j.nuclphysbps.2007.07.005.
  * [13] G. Grindhammer et al., The Fast Simulation of Electromagnetic and Hadronic Showers, Nucl. Instrum. Meth. A 290 (1990) 469.
  * [14] A. Ferrari, P. Sala, A. Fasso, , J. Ranft, FLUKA: a multi-particle transport code, CERN-2005-10 (2005), INFN/TC_05/11, SLAC-R-773.
  * [15] The MARS Code System http://mars.fnal.gov.
  * [16] C. Adloff et al., Study of the interactions of pions in the CALICE silicon-tungsten calorimeter prototype, JINST 5 (2010) 05007.  doi:10.1088/1748-0221/5/05/P05007.
  * [17] M. Apollonio et al. (HARP Collaboration), Forward production of charged pions with incident $\pi^{+}/-$ on nuclear targets measured at the CERN PS, Nucl. Phys. A 821 (2009) 118–192.  doi:10.1016/j.nuclphysa.2009.01.080.
  * [18] M. Apollonio et al. (HARP Collaboration), Large-angle production of charged pions with incident pion beams on nuclear targets, Phys. Rev. C 80 (2009) 065207.  doi:10.1103/PhysRevC.80.065207.
  * [19] M. G. Catanesi et al. (HARP Collaboration), Large-angle production of charged pions with 3-12 GeV/c incident protons on nuclear targets, Phys. Rev. C 77 (2008) 055207.  doi:10.1103/PhysRevC.77.055207.
  * [20] C. Alt et al. (NA49 Collaboration), Inclusive production of charged pions in p+C collisions at 158 GeV/c beam momentum, Eur. Phys. J. C 49 (2007) 897–917.  doi:10.1140/epjc/s10052-006-0165-7.
  * [21] C. Alt et al. (NA49 Collaboration), Inclusive production of charged pions in p+p collisions at 158 GeV/c beam momentum, Eur. Phys. J. C 45 (2006) 343–381.  doi:10.1140/epjc/s2005-02391-9.
  * [22] T. Anticic et al. (NA49 Collaboration), Inclusive production of protons, anti-protons and neutrons in p+p collisions at 158 GeV/c beam momentum, Eur. Phys. J. C 65 (2010) 9–63.  doi:10.1140/epjc/s10052-009-1172-2.
  * [23] N. Abgrall et al. (NA61/SHINE Collaboration), Measurement of production properties of positively charged kaons in proton-carbon interactions at 31 GeV/c, Phys. Rev. C 85 (2012) 035210.  doi:10.1103/PhysRevC.85.035210.
  * [24] N. Abgrall et al. (NA61/SHINE Collaboration), Measurements of cross sections and charged pion spectra in proton-carbon interactions at 31 GeV/c, Phys. Rev. C 84 (2011) 034604.  doi:10.1103/PhysRevC.84.034604.
  * [25] V.V. Uzhinsky for the Geant4 Hadronic Working Group, The Fritiof (FTF) Model in Geant4, International Conference on Calorimetry for the High Energy Frontier (CHEF 2013) C13-04-22.4 (2013) 260–264.
  * [26] D. Wright, M. Kelsey, The Geant4 Bertini Cascade, Nucl. Instrum. Meth. A 804 (2015) 175–188.  doi:10.1016/j.nima.2015.09.058.
  * [27] Yu. D. Bayukov et al., Angular Dependences Of Inclusive Nucleon Production In Nuclear Reactions At High-energies And Separation Of Contributions From Quasifree And Deep Inelastic Nuclear Processes, Sov. J. Nucl. Phys. 42 (1985) 116.
  * [28] S. Banerjee for the CMS Collaboration, Validation of Physics Models of Geant4 using data from CMS experiment, To be published in Journal of Physics Conference Series (JPCS). Proceedings of the CHEP 2016 Conference, San Francisco https://twiki.cern.ch/twiki/bin/viewauth/CMS/Geant4Validation2016.
  * [29] The ATLAS Collaboration, A measurement of the calorimeter response to single hadrons and determination of the jet energy scale uncertainty using LHC Run-1 pp-collision data with the ATLAS detector, Eur. Phys. J. C (2017) 77:26.  doi:10.1140/epjc/s10052-016-4580-0.
  * [30] ATLAS Public Results https://twiki.cern.ch/twiki/bin/view/AtlasPublic/ ApprovedPlotsTileTestBeamResults.
  * [31] The D0 Collaboration, Beam tests of the D0 uranium liquid argon end calorimeters, Nucl. Instrum. Meth. A 324 (1993) 53–76.  doi:10.1016/0168-9002(93)90965-K.
  * [32] The D0 Collaboration, Jet energy scale determination in the D0 experiment, Nucl. Instrum. Meth. A 763 (2014) 442–475.  doi:10.1016/j.nima.2014.05.044.
  * [33] The D0 Collaboration, Determination of the Absolute Jet Energy Scale in the DZERO Calorimeters, Nucl. Instrum. Meth. A 424 (1999) 352–394.  doi:10.1016/S0168-9002(98)01368-0.
  * [34] A. Bhatti et al., Determination of the jet energy scale at the collider detector at Fermilab, Nucl. Instum. Meth. A 566 (2006) 375–412.  doi:10.1016/j.nima.2006.05.269.
  * [35] The CMS Collaboration, Determination of Jet Energy Calibration and Transverse Momentum Resolution in CMS, JINST 6 (2011) 11002.  doi:10.1088/1748-0221/6/11/P11002.
  * [36] The ATLAS Collaboration, Jet energy measurement with the ATLAS detector in proton–proton collisions at $\sqrt{s}=7$ TeV, Eur. Phys. J. C 73 (2013) 2304.  doi:10.1140/epjc/s10052-013-2304-2.
  * [37] Florian Beaudette for the CMS Collaboration, The CMS Particle Flow Algorithm, International Conference on Calorimetry for the High Energy Frontier (CHEF 2013) C13-04-22.4 (2013) 295–304.
  * [38] The ATLAS Collaboration, Jet energy measurement and its systematic uncertainty in proton–proton collisions at $\sqrt{s}=7$ TeV with the ATLAS detector, Eur. Phys. J. C 75 (2015) 17.  doi:10.1140/epjc/s10052-014-3190-y.
  * [39] The ATLAS Collaboration, Jet energy resolution in proton-proton collisions at $\sqrt{s}=7$ TeV recorded in 2010 with the ATLAS detector, Eur. Phys. J. C 73 (2013) 2306.  doi:10.1140/epjc/s10052-013-2306-0.
  * [40] S. P. Martin, A Supersymmetry Primer arXiv:hep-ph/9709356.
  * [41] The CMS Collaboration, Search for New Physics with Jets and Missing Transverse Momentum in pp collisions at $\sqrt{s}=7$ TeV, JHEP 8 (2011) 155\.  doi:10.1007/JHEP08(2011)155.
  * [42] T. Schum, QCD Processes and Search for Supersymmetry at the LHC, DESY-THESIS-2012-029.
  * [43] The CMS Collaboration, Performance of CMS muon reconstruction in pp collision events at $\sqrt{s}=7$ TeV, JINST 7 (2012) 10002.  doi:10.1088/1748-0221/7/10/P10002.
  * [44] The CMS Collaboration, Performance of electron reconstruction and selection with the CMS detector in proton-proton collisions at $\sqrt{s}=8$ TeV, JINST 10 (2015) 06005.  doi:10.1088/1748-0221/10/06/P06005.
  * [45] The ATLAS Collaboration, Measurement of the muon reconstruction performance of the ATLAS detector using 2011 and 2012 LHC proton-proton collision data, Eur. Phys. J. C 74 (2014) 3130.  doi:10.1140/epjc/s10052-014-3130-x.
  * [46] The ATLAS Collaboration, Electron efficiency measurements with the ATLAS detector using 2012 LHC proton-proton collision data, Submitted to Eur. Phys. J. C arXiv:1612.01456 [hep-ex].
  * [47] The CMS Collaboration, Search for new physics in the multijet and missing transverse momentum final state in proton-proton collisions at $\sqrt{s}=8$ TeV, JHEP 6 (2014) 055.  doi:10.1007/JHEP06(2014)055.
  * [48] The CMS Collaboration, CMS Phase II Upgrade Scope Document, CERN-LHCC-2015-019 https://cds.cern.ch/record/2055167.
  * [49] J.de Favereau et al., DELPHES 3: a modular framework for fast simulation of a generic collider experiment, JHEP 02 (2014) 057.  doi:10.1007/JHEP02(2014)057.
  * [50] The ATLAS Collaboration, ATLAS Phase II Upgrade Scope Document, CERN-LHCC-2015-020 https://cds.cern.ch/record/2055248.
  * [51] WLCG: Worldwide LHC Computing Grid, http://wlcg-public.web.cern.ch/tier-centres.
  * [52] D. Futyan, R. Mankel, Christoph Paus for the CMS Collaboration, The CMS Computing, Software and Analysis Challenge, Journal of Physics: Conference Series 219 (2010) 032008.  doi:10.1088/1742-6596/219/3/032008.
  * [53] The CMS Collaboration, Description and performance of track and primary-vertex reconstruction with the CMS tracker, JINST 9 (2014) 10009.  doi:10.1088/1748-0221/9/10/P10009.
  * [54] The CMS Collaboration, Studies of Tracker Material, CMS Public Note, CMS-PAS-TRK-10-003 https://cds.cern.ch/record/1279138.
  * [55] The ATLAS Collaboration, Photon Conversions at $\sqrt{s}=900$ GeV measured with the ATLAS Detector, ATLAS Conference Note, ATLAS-CONF-2010-007 https://cds.cern.ch/record/1274001.
  * [56] The ATLAS Collaboration, Study of the Material Budget in the ATLAS Inner Detector with $K_{S}^{0}$ decays in collision data at $\sqrt{s}=900$ GeV, ATLAS Conference Note, ATLAS-CONF-2010-019 https://cds.cern.ch/record/1277651.
  * [57] The CMS Collaboration, Identification of b-quark jets with the CMS experiment, JINST 8 (2013) 04013.  doi:10.1088/1748-0221/8/04/P04013.
  * [58] The ATLAS Collaboration, Performance of b-jet identification in the ATLAS experiment, JINST 11 (2016) 04008.  doi:10.1088/1748-0221/11/04/P04008.
  * [59] Torbjorn Sjostrand et al., An Introduction to PYTHIA 8.2, arXiv:1410.3012 [hep-ph]doi:10.1016/j.cpc.2015.01.024.
  * [60] T. Sjostrand, S. Mrenna, P. Skands, PYTHIA 6.4 Physics and Manual, JHEP05 (2006) 026doi:10.1088/1126-6708/2006/05/026.
  * [61] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H.-S. Shao, T. Stelzer, P. Torrielli, M. Zaro, The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP07 (2014) 079doi:10.1007/JHEP07(2014)079.
  * [62] The CMS Collaboration, Performance of the CMS missing transverse energy reconstruction in $pp$ data at $\sqrt{s}=8$ TeV, JINST 10 (2015) 2006\.  doi:10.1088/1748-0221/10/02/P02006.
  * [63] The ATLAS Collaboration, Measurement of $W^{\pm}$ and $Z$-boson production cross sections in pp collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, Phys. Lett. B 759 (2016) 601.  doi:10.1016/j.physletb.2016.06.023.
  * [64] The ATLAS Collaboration, Performance of algorithms that reconstruct missing transverse momentum in $\sqrt{s}=8$ TeV proton–proton collisions in the ATLAS detector, Submitted to Eur. Phys. J. C arXiv:1609.09324 [hep-ex].
  * [65] The ATLAS Collaboration, Measurement of inclusive jet and dijet production in $pp$ collisions at $\sqrt{s}=7$ TeV using the ATLAS detector, Phys. Rev. D 86 (2012) 014022.  doi:10.1103/PhysRevD.86.014022.
  * [66] The CMS Collaboration, Measurement of the Inclusive Jet Cross Section in $pp$ Collisions at $\sqrt{s}=7$ TeV, Phys. Rev. Lett. 107 (2011) 132001.  doi:10.1103/PhysRevLett.107.132001.
  * [67] The CDF Collaboration, Inclusive Jet Cross Section in $\overline{p}p$ Collisions at $\sqrt{s}=1.8$ TeV, Phys. Rev. Lett. 77 (1996) 438.  doi:10.1103/PhysRevLett.77.438.
  * [68] The D0 Collaboration, Inclusive Jet Cross Section in $\overline{p}p$ Collisions at $\sqrt{s}=1.8$ TeV, Phys. Rev. Lett. 82 (1999) 2451.  doi:10.1103/PhysRevLett.82.2451.
  * [69] The ATLAS Collaboration, Measurement of inclusive jet and dijet cross sections in proton-proton collisions at 7 TeV centre-of-mass energy with the ATLAS detector, Eur. Phys. J. C 71 (2011) 1512.  doi:10.1140/epjc/s10052-010-1512-2.
  * [70] The D0 Collaboration, Inclusive Jet Production in $\overline{p}p$ Collisions, Phys. Rev. Lett. 86 (2001) 1707–1712.  doi:10.1103/PhysRevLett.86.1707.
  * [71] The CDF Collaboration, Measurement of the inclusive jet cross section at the Fermilab Tevatron $p\overline{p}$ collider using a cone-based jet algorithm, Phys. Rev. D 78 (2008) 052006.  doi:10.1103/PhysRevD.78.052006.
  * [72] The D0 Collaboration, Measurement of the inclusive jet cross section in $p\overline{p}$ collisions at $\sqrt{s}=1.96$ TeV, Phys. Rev. D 85 (2012) 052006.  doi:10.1103/PhysRevD.85.052006.
  * [73] The CDF Collaboration, Measurement of the Inclusive Jet Cross Section in $\overline{p}p$ Collisions at $\sqrt{s}=1.8$ TeV, Phys. Rev. Lett. 62 (1989) 613.  doi:10.1103/PhysRevLett.62.613.
  * [74] The CDF Collaboration, Inclusive jet cross section in $\overline{p}p$ collisions at $\sqrt{s}=1.8$ TeV, Phys. Rev. Lett. 68 (1992) 1104.  doi:10.1103/PhysRevLett.68.1104.
  * [75] Summary Table https://www-cdf.fnal.gov/physics/physics.html.
  * [76] Year-to-year history of paper submissions https://www-d0.fnal.gov/Run2Physics/WWW/results.htm.
  * [77] CMS Publications versus Time http://cms-results.web.cern.ch/cms-results/public-results/publications-vs-time/.
  * [78] J. Apostolakis et al. (GeantV Development Group), Adaptive track scheduling to optimize concurrency and vectorization in GeantV, J. Phys.: Conf. Ser. 608 (2015) 012003.  doi:10.1088/1742-6596/608/1/012003.
  * [79] Source: Geant4 Testing and Quality Assurance Working Group https://g4cpt.fnal.gov/g4p/prplots/cpu_by_version.html.
  * [80] J. Pearlstein, Information Revolution: Big Data Has Arrived at an Almost Unimaginable Scale, Wired Magazine https://www.wired.com/2013/04/bigdata.