# Where are NANOGrav’s big black holes?

Gabriela Sato-Polito<EMAIL_ADDRESS>School of Natural Sciences,
Institute for Advanced Study, Princeton, NJ 08540, United States  Matias
Zaldarriaga  School of Natural Sciences, Institute for Advanced Study,
Princeton, NJ 08540, United States  Eliot Quataert  Department of
Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

###### Abstract

Multiple pulsar timing array (PTA) collaborations have recently reported the
first detection of gravitational waves (GWs) of nanohertz frequencies. The
signal is expected to be primarily sourced by inspiralling supermassive black
hole binaries (SMBHBs) and these first results are broadly consistent with the
expected GW spectrum from such a population. Curiously, the measured amplitude
of the GW background in all announced results is a bit larger than theoretical
predictions. In this work, we show that the amplitude of the stochastic
gravitational wave background (SGWB) predicted from the present-day abundance
of SMBHs derived from local scaling relations is significantly smaller than
that measured by the PTAs. We demonstrate that this difference cannot be
accounted for through changes in the merger history of SMBHs and that there is
an upper limit to the boost to the characteristic strain from multiple merger
events, due to the fact that they involve black holes of decreasing masses. If
we require the current estimate of the black hole mass density — equal to the
integrated quasar luminosity function through the classic Sołtan argument — to
be preserved, then the currently measured PTA result would imply that the
typical total mass of SMBHs contributing to the background should be at least
$\sim 3\times 10^{10}M_{\odot}$, a factor of $\sim 10$ larger than previously
predicted. The required space density of such massive black holes corresponds
to order 10 $3\times 10^{10}M_{\odot}$ SMBHs within the volume accessible by
stellar and gas dynamical SMBH measurements. By virtue of the GW signal being
dominated by the massive end of the SMBH distribution, PTA measurements offer
a unique window into such rare objects and complement existing electromagnetic
observations.

Evidence for the first detection of low-frequency ($\sim~{}1-100$nHz)
gravitational waves (GWs) has recently been reported by multiple pulsar timing
array (PTA) collaborations. The signature of GWs on the arrival time of pulses
is a correlated red noise with a particular angular dependence given by the
Helling-Downs HD  curve. This signature has been detected with varying levels
of significance by the European PTA and Indian PTA (EPTA:2023fyk , hereafter
EPTA+InPTA), by the North American Nanohertz Observatory for Gravitational
waves (NANOGrav:2023gor , hereafter NANOGrav), Parkes PTA (Reardon:2023gzh ,
hereafter PPTA), and the Chinese Pulsar Timing Array Xu:2023wog , that span
between $2-4\sigma$. Their results have also been shown to be consist with
each other, with discrepancies smaller than $1\sigma$, despite the differences
in both observations and modelling choices
InternationalPulsarTimingArray:2023mzf .

Although the origin of the signal is unknown, the dominant source of the
stochastic gravitational wave background (SGWB) is expected to be a population
of inspiralling supermassive black-hole binaries (SMBHBs) (see, e.g., Ref.
Burke-Spolaor:2018bvk ). Strong observational evidence suggests that most, if
not all, galaxies host SMBHs at their centers 1995ARA&A..33..581K ;
1998Natur.395A..14R . Following the merger of host galaxies, SMBHs form
binaries that may also merge 1980Natur.287..307B , which would produce the
observed PTA signal, albeit through a highly uncertain process.

The most precise and direct measurements of the SMBH population are derived
from observations in our local Universe, where direct kinematic estimates of
their masses are possible kormendy_ho ; mcconnell_ma  for several tens of
black holes, along with various properties of the host galaxies. The tight
correlations between the SMBH mass and host properties likely suggest a shared
evolution and can be extrapolated to the total SMBH population by combining
such scaling relations with galaxy catalogs 2004MNRAS.351..169M ; vika2009 ;
shankar2009 ; shankar2013 . At higher redshifts, inferences of the SMBH
population rely on the active black holes powering AGNs. The total mass
accumulated in SMBHs from AGN relics can be directly related to the total
luminosity they emit through the classic Sołtan argument soltan1982 , while
the cosmic evolution of the SMBH mass function can be theoretically modelled
through the continuity equation 1971ApJ…170..223C ; 1992MNRAS.259..725S ;
2002MNRAS.335..965Y ; 2013MNRAS.428..421S ; Tucci:2016tyc .

The GW spectrum measured by recent PTA analyses is broadly consistent with
expectations from an inspiral driven by GW emission Phinney:2001di , with the
largest deviation in the slope of the spectrum found by NANOGrav. All
announced results, however, found amplitudes larger than expected from
theoretical predictions NANOGrav:2023hfp . In this work, we use the formalism
laid out in Ref. Phinney:2001di  to connect the binaries contributing to the
SGWB to the present-day population of remnant (single) black holes in order to
assess the significance of this discrepancy. This approach yields a simple but
robust model for the SGWB that sheds light on the uncertainties associated
with the merger history of SMBHBs, which enables a clear comparison between
PTAs and other observations of SMBHs.

We show that the difference between predicted and measured values of the SGWB
cannot be accounted for by an increased black hole merger rate for a SMBH mass
function consistent with local observations. We demonstrate that there is a
mathematical upper limit to the boost of the characteristic strain from
multiple mergers, which is only a factor of $\sim 1.64$ due to the fact that
they involve black holes of decreasing masses. Since the predicted amplitude
of the SGWB falls short of the measured value even in the most ideal scenario,
we explore modifications to the mass function prediction that could lead to
this difference.

If the current estimate of the mass density in black holes $\rho_{\rm BH}$,
consistent with the Sołtan argument, is to be preserved, we argue that all
modifications to the SMBH mass function are equivalent to a change in the peak
mass contributing to the SGWB. The PTA measurements would therefore generally
imply that the SGWB is dominated by SMBHs roughly 10 times larger than
previously estimated, of around $M^{\rm GW}_{\rm peak}\sim 3\times
10^{10}M_{\odot}$. We begin by investigating whether a higher intrinsic
scatter $\epsilon_{0}$ in the $M-\sigma$ relation or a higher value for the
turnover $\sigma_{*}$ in the velocity dispersion function might be responsible
for this misestimation, but find that the required values are in significant
disagreement with current observations.

Another reason why the peak BH mass contributing to the SGWB may be larger is
if the $M-\sigma$ relation becomes steeper for high velocity dispersion
galaxies. Indeed, although there are too few objects to derive definitive
conclusions, the $M-\sigma$ does appear to steepen for $\sigma\gtrsim 275$km
s-1kormendy_ho ; mcconnell_ma , a consequence of the galaxy population
“saturating” at high $\sigma$ 2013ApJ…769L…5K ; Thomas:2013bga ;
2012MNRAS.424..224H . The currently favored picture is that these massive
galaxies that host the most massive SMBHs have cores that are formed in the
aftermath of gas-poor major mergers by SMBHBs, which scour the stars at the
center of the merged galaxy 1991Natur.354..212E ; 1997AJ….114.1771F ;
2001ApJ…563…34M ; 2019ApJ…887…10S . This scenario suggests that the high-mass
end of the SMBH mass function would be underpredicted by the $M-\sigma$
relation, but would be more accurately estimated from the luminosity or bulge
mass relations. We show, however, that the amplitude of the SGWB predicted
from the relation between black hole mass and bulge mass also leads to a
smaller value than observed, thus suggesting that this effect is not
sufficient to account for the discrepancy. If a larger characteristic strain
is to be explained by a steeper $M-\sigma$ relation, it must therefore
necessarily imply either a larger abundance of high stellar mass galaxies or a
corresponding change to the $M-M_{\rm bulge}$ relation. Finally, we also
consider a more general modification to the SMBH mass function, in which we
simply boost the mass of the black hole associated with galaxies of a given
velocity dispersion. We show what the required black hole mass is for each
value of $\sigma$, suggesting that the difference may be explained by a small
number of extremely massive black holes or a large number of lighter ones.

This paper is organized as follows. In Sec. I we discuss the SMBH population
inferred from local scaling relations coupled to galaxy catalogs, as well as
the information offered by the quasar luminosity function about this
population. We discuss the model for the characteristic strain in Sec. II and
how it depends on the present-day SMBH population presented in the previous
section. In Sec. III we show the discrepancy between the two sets of
observations and discuss potential sources of this difference, and conclude in
Sec. IV.

##  I Black hole population from scaling relations and AGNs

For nearby sources, black hole masses can be directly estimated from gas or
stellar kinematics, along with host galaxy properties such as velocity
dispersion, bulge stellar mass, and luminosities. The current sample contains
several tens of galaxies kormendy_ho ; mcconnell_ma  and establishes tight
relations between supermassive black holes and their hosts, which are given by
power laws between host properties ($X$) and the black hole mass ($M_{\rm
BH}$)

$\log_{10}M_{\rm BH}=a_{\bullet}+b_{\bullet}\log_{10}X,$ (1)

with a log-normal intrinsic scatter of $\epsilon_{0}$. In order to infer the
total SMBH population, these scaling relations must then be combined with
galaxy catalogs.

Our fiducial model relies on the galaxy stellar velocity dispersion ($\sigma$)
as a proxy for black hole mass. We combine the $M-\sigma$ scaling relation
above (with $X=\sigma/200$km s-1) and combine with the galaxy velocity
dispersion function (VDF), typically parametrized as

$\phi(\sigma)d\sigma=\phi_{*}\left(\frac{\sigma}{\sigma_{*}}\right)^{\alpha}\frac{e^{-(\sigma/\sigma_{*})^{\beta}}}{\Gamma(\alpha/\beta)}\beta\frac{d\sigma}{\sigma},$
(2)

which corresponds to a generalization of the Schechter function. While many
host galaxy properties have been shown to correlate with the mass of the
central SMBH, recent work suggests that velocity dispersion offers a more
reliable, and perhaps more fundamental, prediction of black hole mass
kormendy_ho ; 2016ApJ…831..134V ; Marsden:2020yyd . However, this may not
necessarily hold for the most massive SMBHs — the velocity dispersion of
galaxies appears to saturate above $\sigma\gtrsim 275$km s-1, which might lead
to an underprediction of the massive end of the SMBH mass function relative to
bulge mass ($M_{\rm bulge}$) or luminosity scaling relations 2007ApJ…662..808L
; 2019ApJ…887…10S . Hence, while we focus on the velocity dispersion in this
discussion, we also estimate the SMBH mass function using the bulge mass
scaling relation.

In the absence of scatter in the $M-\sigma$ relation, the SMBH mass function
can be computed by directly combining Eqs. 1 and 2. Substituting velocity
dispersion for the black hole mass leads to

$\phi(M_{\rm BH})dM_{\rm BH}=\phi_{*}\left(\frac{M_{\rm
BH}}{M_{*}}\right)^{\alpha/b_{\bullet}}\frac{e^{-(\frac{M_{\rm
BH}}{M_{*}})^{\beta/b_{\bullet}}}}{\Gamma(\alpha/\beta)}\frac{\beta}{b_{\bullet}}\frac{dM_{\rm
BH}}{M_{\rm BH}},$ (3)

where we have defined
$\log_{10}M_{*}=a_{\bullet}+b_{\bullet}\log_{10}\sigma_{*}/200$km s-1. In the
presence of scatter, we now have that the mass function is given by a
convolution between the velocity dispersion function and the probability
distribution function $p(\log_{10}M_{\rm BH}|\log_{10}\sigma)$

$\phi(M_{\rm BH})=\int d\sigma\frac{p(\log_{10}M_{\rm
BH}|\log_{10}\sigma)}{M_{\rm BH}\log(10)}\phi(\sigma),$ (4)

where we assume $p$ to be log-normal

$\begin{split}p(\log_{10}&M_{\rm
BH}|\log_{10}\sigma)=\frac{1}{\sqrt{2\pi}\epsilon_{0}}\\\
&\times\exp\left[-\frac{1}{2}\left(\frac{\log_{10}M_{\rm
BH}-a_{\bullet}-b_{\bullet}\log_{10}\sigma}{\epsilon_{0}}\right)^{2}\right].\end{split}$
(5)

Figure 1: SMBH mass function computed from a combination of VDFs and GSMF with
local scaling relations reported in the literature. The VDFs correspond to
Refs. 2010MNRAS.404.2087B  and 2011ApJ…737L..31B (B10 and B11, respectively),
the GSMF is given in 2015MNRAS.454.4027D  (DS15), while the adopted scaling
relations are presented in Refs. kormendy_ho  and mcconnell_ma  (KH13 and
MM13, respectively). For the sake of visual clarity, we have omitted the error
bands of most SMBH mass functions in the plot, and show only for the B10 VDF
result since it corresponds to our fiducial choice. The grey dashdotted line
shows an example of a SMBH mass function that matches the characteristic
strain observed by PTAs, computed by increasing the intrinsic scatter to
$\epsilon_{0}=0.62$.

We predict the SMBH mass function from various combinations of velocity
dispersion functions and $M-\sigma$ relations, which we describe next. Ref.
2010MNRAS.404.2087B  (B10) uses spectroscopic data from SDSS of $\sim 250000$
galaxies across $4681$deg2 of the sky and redshifts $z\lesssim 0.3$. We choose
the VDF fit to all galaxies and including only $\sigma>125$km s-1, since the
low-velocity dispersion end is subject to greater uncertainties and systematic
errors (e.g. sample completeness), but are not expected to host SMBHs that
significantly contribute to the SGWB. Ref. 2011ApJ…737L..31B  (B11) uses
photometric data from the UKIDSS Ultra-Deep Survey (UDS) and the NEWFIRM
Medium Band Survey of the COSMOS field (NMBS COSMOS), covering 0.77deg2 and
0.21deg2, respectively. Stellar velocity dispersions are estimated from the
radius, Sérsic index, and stellar mass using the virial theorem. We use the
results measured between redshift $0.3\leq z\leq 0.6$ and fit the generalized
Schechter function given in Eq. 2. Finally, Ref. 2022ApJ…939…90T  (T20)
reports spectroscopic measurements of the VDF from the Large Early Galaxy
Astrophysics Census (LEGA-C) survey, which targets massive galaxies in the
COSMOS field, from which we use the measurements between $0.6<z\leq 0.7$ to
fit Eq. 2. Despite the differences in redshift, Refs. 2011ApJ…737L..31B  and
2022ApJ…939…90T , which explore broader redshift ranges, do not find evidence
for substantial evolution in the VDF, and we therefore compare the three
different measurements as representatives of the present-day VDF. We use two
different fits to the $M-\sigma$ relation from Refs. kormendy_ho  (KH13) and
mcconnell_ma  (MM13, fit to all galaxies).

We also estimate the SMBH mass function using the bulge mass scaling relation.
We use the galaxy stellar mass function (GSMF) from Ref. 2015MNRAS.454.4027D
(DS15), which was measured from a complete sample of half a million galaxies
from the SDSS survey. A significant source of systematic error in measurements
of the stellar mass function arises in the determination of the total flux in
a given bandpass and the impact of sky background subtraction. Ref.
2015MNRAS.454.4027D  focuses on estimating flux corrections by stacking
samples of similar galaxies selected by various properties, while accounting
for a variety of potential sources of systematic biases. The measured stellar
mass function lies roughly in between previously published results
2009MNRAS.398.2177L ; Bernardi:2013mqa . We combine this GSMF with the scaling
relations from MM13 estimated from bulges with dynamical masses, and from
KH13, where we assume a bulge fraction of 1.

In Fig. 1, we show multiple estimates of the SMBH mass function using a
variety of different types of observations. All choices lead to similar black-
hole abundances for the relevant mass range probed by the SGWB (which, as we
show in the next section, is around $\sim 3\times 10^{9}M_{\odot}$). For the
estimates based on velocity dispersion, this can be attributed to the fact
that direct mass measurements for such objects are relatively robust and
correspond to galaxies with $\sigma\sim 300$km s-1. Low velocity dispersion
galaxies $\sigma\lesssim 150$km s-1 are impacted by sample completeness, while
galaxies with high velocity dispersions $\sigma\gtrsim 400$km s-1 are rare and
few objects can be found in observed catalogs. The typical galaxies
contributing to the SGWB are not impacted by either effects or by
morphological selection criteria, being dominated by early-type galaxies.

We choose as our fiducial model the combination of the VDF from B10 and the
$M-\sigma$ relation from MM13. In summary, $X=\sigma/200$km s-1 in Eq. 1 and
we have chosen fiducial parameters: $\phi_{*}=(2.61\pm 0.16)\times
10^{-2}$Mpc-3, $\sigma_{*}=159.6\pm 1.5$ km s-1, $\alpha=0.41\pm 0.02$, and
$\beta=2.59\pm 0.04$ for the VDF, and $a_{\bullet}=8.32\pm 0.05$,
$b_{\bullet}=5.64\pm 0.32$, and $\epsilon_{0}=0.38$, for the $M-\sigma$
relation.

Beyond our local Universe, only active SMBHs powering active galactic nuclei
(AGN) can be observed. Assuming that the present-day population of SMBHs are
relics of AGNs, the total mass density accumulated in black holes can be
related to the energy density emitted by AGNs through the Sołtan argument
soltan1982 . If all black holes accrete with a constant radiative efficiency
$\epsilon_{r}$, the quasar bolometric luminosity ($L_{\rm bol}$) can be
related to the SMBH accretion rate

$L_{\rm bol}=\epsilon_{r}\dot{M}c^{2},$ (6)

and the integrated quasar luminosity function therefore corresponds to the
integrated mass density in SMBHs

$\displaystyle\rho_{\rm BH}$ $\displaystyle=\int dM_{\rm BH}\phi(M_{\rm
BH})M_{\rm BH}$ (7)
$\displaystyle=\frac{1-\epsilon_{r}}{\epsilon_{r}c^{2}}\int_{0}^{\infty}dz\frac{dt}{dz}\int
d\log L_{\rm bol}\phi(L_{\rm bol},z)L_{\rm bol},$ (8)

with a relatively consistent value found throughout the literature
2002MNRAS.335..965Y ; 2004MNRAS.351..169M ; shankar2009 ; 2007ApJ…654..731H .
For instance, Ref. 2007ApJ…654..731H  found a value of

$\rho_{\rm BH}=4.81^{+1.24}_{-0.99}\left(\frac{0.1}{\epsilon_{r}}\right)\times
10^{5}M_{\odot}\text{Mpc}^{-3},$ (9)

where the radiative efficiency is typically assumed to be a free parameter,
fit to explain the relic black hole population.

##  II Black hole population from PTAs

The characteristic strain of the gravitational-wave background produced by
inspiralling SMBHBs can be written as the integral of the emission from
binaries belonging to each mass and redshift interval, summed over the entire
population of sources. That is,

$h^{2}_{c}(f)=\frac{4G}{\pi c^{2}f^{2}}\int
d\mathcal{M}\int\frac{dz}{1+z}\frac{dn}{dzd\mathcal{M}}\left.\left(\frac{dE_{\rm
gw}}{d\log f_{r}}\right)\right|_{f_{r}=(1+z)f}$ (10)

where $f_{r}$ and $f$ are the emitted and observed GW frequencies, and
$dn/dzd\mathcal{M}$ is the number of black hole binaries from redshifts
between $[z,z+dz]$ and chirp masses $[\mathcal{M},\mathcal{M}+d\mathcal{M}]$.
For a circular orbit in the Newtonian regime, the energy spectrum depends only
on the chirp mass of the binary, and is given by

$\frac{dE_{\rm gw}}{d\log f_{r}}=\frac{1}{3G}(G\mathcal{M})^{5/3}(\pi
f_{r})^{2/3}.$ (11)

Switching from chirp mass to total mass $M_{\rm BH}$, and defining

$\epsilon({M_{\rm
BH}})\equiv\frac{\eta}{3}\left(\frac{v}{c}\right)^{2}=\frac{\eta}{3}\frac{(GM_{\rm
BH}\pi f)^{2/3}}{c^{2}},$ (12)

where $\eta=q/(1+q)^{2}$ is the symmetric mass ratio and $q=m_{1}/m_{2}$.
Using this definition, the characteristic strain can be rewritten as

$\begin{split}h^{2}_{c}(f)=&\frac{4G}{\pi f^{2}}\int dM_{\rm
BH}\int\frac{dz}{(1+z)^{1/3}}\int dq\frac{dn}{dzdM_{\rm BH}dq}\\\ &\times
M_{\rm BH}\epsilon({M_{\rm BH}},q)\end{split}$ (13)

We can interpret $\epsilon$ as an average efficiency of GW emission.

Figure 2: Contribution to the mass density in black holes (red) and to the
SGWB (dark blue) per logarithmic mass bin, assuming the fiducial set of
parameters (which includes scatter in the $M-\sigma$ relation). The peak of
the each kernel is marked by the vertical dashed lines and the exact values
are quoted above.

A crucial step, pointed out in Ref. Phinney:2001di , is that the remnant black
hole population — resultant from the merger history described by
$\frac{dn}{dzdM_{\rm BH}dq}$ — must correspond to the present-day single SMBH
population. Assuming all black holes have undergone a single merger event in
their lifetimes, the mass function of mergers integrated over all redshifts
and mass ratios must correspond to the single BH mass function ($\phi$) today

$\phi(M_{\rm BH})=\frac{dn}{dM_{\rm BH}}=\int dz\int dq\frac{dn}{dzdM_{\rm
BH}dq}.$ (14)

We will demonstrate in Sec. III.2 that the assumption of a single merger event
provides a reasonable estimate of the SGWB, since repeated mergers do not
significantly increase GW emission. This is due to the fact that each earlier
merger must occur between black holes of decreasing masses and therefore
contribute progressively less to the characteristic strain.

Assuming that the mass and redshift dependence are separable, such that
$\frac{dn}{dzdM_{\rm BH}dq}=p_{z}(z)p_{q}(q)\frac{dn}{dM_{\rm BH}}$, we can
write the characteristic strain of the SGWB in terms of the local SMBH mass
function

$h^{2}_{c}(f)=\frac{4G}{\pi f^{2}}\langle(1+z)^{-1/3}\rangle\int dM_{\rm
BH}\phi(M_{\rm BH})M_{\rm BH}\langle\epsilon\rangle({M_{\rm BH}}),$ (15)

where we have defined the redshift and mass ratio averages

$\langle(1+z)^{-1/3}\rangle=\int dz\frac{p_{z}(z)}{(1+z)^{1/3}}$ (16)

and

$\langle\epsilon\rangle({M_{\rm BH}})=\int dq\ p_{q}(q)\epsilon({M_{\rm
BH}},q).$ (17)

We choose the following parametrizations for the redshift and mass-ratio
distribution functions

$p_{z}(z)=z^{\gamma}e^{-z/z_{*}},\quad\text{and}\quad p_{q}(q)=q^{\delta},$
(18)

with fiducial values of $\gamma=0.5$, $z_{*}=0.3$, and $\delta=-1$, an assumed
minimum value of $q_{\rm min}=0.1$, and where both functions are then
normalized to unity. These values are chosen to roughly reproduce the
simulation-based results of Ref. Kelley:2016gse , such that the median
redshift of the SGWB signal is sourced at $z\sim 0.3$ and over $90\%$ of the
contribution comes from $z\lesssim 1$. Notice that the redshift integral is
very insensitive to variations in the probability distribution of mergers. Our
fiducial set of parameters leads to $\langle(1+z)^{-1/3}\rangle=0.9$, while
its upper limit is 1 and, for instance, a much broader redshift distribution
with $z_{*}=2$ results in $\langle(1+z)^{-1/3}\rangle=0.68$. The maximum
efficiency of $\eta=0.25$ corresponds to equal mass ratio mergers, while our
fiducial log-uniform distribution leads to $\langle\eta\rangle=0.18$.

Using the black hole mass function predicted from Eq. 4 and our fiducial set
of parameters, we show in Fig. 2 the contribution to the BH mass density and
the characteristic strain per logarithmic mass. This illustrates that there is
a characteristic mass that primarily dictates the observed signals and that
the SGWB receives the majority of its contribution from significantly higher
masses than the mass density, since it is weighted by a factor of $M_{\rm
BH}^{5/3}$.

Figure 3: Characteristic strain, black hole mass density, and their ratios as
a function of the SMBH mass that most contributes to the SGWB ($M^{\rm
GW}_{\rm peak}$), and the corresponding intrinsic scatter $\epsilon_{0}$ in
the upper axis. In the bottom panel, the value of the ratio between
characteristic strain and the fiducial estimate are shown in the left axis,
while the black hole mass density is shown on the right.

The black hole mass density and characteristic strain can now be computed
using Eqs. 7 and 15 by directly integrating over the BH mass function,

$\displaystyle\rho_{{\rm BH},0}=$
$\displaystyle\phi_{*}M_{*}\frac{\Gamma\left(\nicefrac{{(\alpha+b_{\bullet})}}{{\beta}}\right)}{\Gamma\left(\nicefrac{{\alpha}}{{\beta}}\right)},$
(19) $\displaystyle h^{2}_{c,0}(f)=$
$\displaystyle\frac{4G\langle(1+z)^{-1/3}\rangle}{\pi c^{2}f^{2}}\rho_{{\rm
BH},0}\epsilon(M_{*})\frac{\Gamma\left(\nicefrac{{(\alpha+5b_{\bullet}/3)}}{{\beta}}\right)}{\Gamma\left(\nicefrac{{(\alpha+b_{\bullet})}}{{\beta}}\right)},$
(20)

where the gamma function $\Gamma$ results from the mass integral and the
assumed functional form of the velocity dispersion function, and the subscript
$0$ denotes the value without scatter in the $M-\sigma$ relation.

By combining Eqs. 7 and 15 with 4 we can show that scatter boosts the mass
density by a factor of $\rho_{\rm BH}=\rho_{\rm
BH,0}e^{\epsilon_{0}^{2}\log^{2}(10)/2}$ and the characteristic strain by
$h^{2}_{c}=h^{2}_{c,0}e^{25/18\epsilon_{0}^{2}\log^{2}(10)}$. We find the
following values for the black hole mass density and characteristic strain

$\begin{split}\rho_{\rm BH}=&4.5\times
10^{5}M_{\odot}\text{Mpc}^{-3}\left(\frac{M_{*}}{5.8\times
10^{7}M_{\odot}}\right)\\\
&\times\left(\frac{1.5}{e^{\frac{1}{2}\epsilon_{0}^{2}\log^{2}(10)}}\right)\\\
h^{2}_{c}(f)=&1.2\times
10^{-30}\left(\frac{f}{\text{yr}^{-1}}\right)^{-4/3}\left(\frac{M_{*}}{5.8\times
10^{7}M_{\odot}}\right)^{2/3}\\\ &\times\left(\frac{\rho_{\rm BH}}{4.5\times
10^{5}M_{\odot}\text{Mpc}^{-3}}\right)\left(\frac{2}{e^{\frac{8}{9}\epsilon_{0}^{2}\log^{2}(10)}}\right).\end{split}$
(21)

Notice that the mass scale $M_{*}$ corresponds to the turnover in the velocity
dispersion function, but it is not the mass that most contributes to the SGWB
or the mass density. These were shown in Fig.2 and can be derived analytically
in absence of scatter from Eqs. 3, 19, and 20, by computing the peak of the
mass integral kernel. The mass that most contributes to the BH mass density
and SGWB, per logarithmic mass bin, are given by

$\begin{split}M^{\rho}_{\rm
peak}=&M_{*}\left[\frac{(\alpha+b_{\bullet})}{\beta}\right]^{b_{\bullet}/\beta}\\\
M^{\rm GW}_{\rm
peak}=&M_{*}\left[\frac{(\alpha+\nicefrac{{5b_{\bullet}}}{{3}})}{\beta}\right]^{b_{\bullet}/\beta},\end{split}$
(22)

which then scales as $\log M_{\rm peak}\propto\epsilon^{2}_{0}$ in the
presence of scatter.

One interesting quantity to consider is the ratio between the characteristic
strain squared and the black hole mass density. Although the mass density can
only be inferred from the quasar luminosity function up to a constant factor
given by the radiative efficiency, the ratio $h^{2}_{c}/\rho_{\rm BH}$ cancels
the amplitude of the mass function and therefore any dependence on
$\epsilon_{r}$. If we therefore consider a black hole mass function that leads
to a $\rho_{\rm BH}$ consistent with the Sołtan argument for some value of
$\epsilon_{r}$, and a value of $h_{c}^{2}$, then this ratio can only be
altered through changes to the mass function that lead to a different $M_{\rm
peak}^{\rm GW}$. Increasing the characteristic strain while maintaining the
current estimate of $\rho_{\rm BH}$ therefore implies modifications to the
high-mass tail of the SMBH mass function, that can be generally described as a
change to the peak mass contributing to the SGWB. The ratio as a function of
the peak mass is shown on the top panel of Fig. 3, while $h^{2}_{c}$ and
$\rho_{\rm BH}$ are shown on the bottom.

Figure 4: Comparison between the characteristic strain predicted using galaxy
scaling relations and galaxy catalogs (open symbols), and the values reported
by PTA collaborations (filled dots). The left vertical axis shows the
characteristic strain value at the reference frequency of $f=$yr-1 for a fixed
power law of $\gamma=13/3$, while the right vertical axis shows the ratio with
respect to the NANOGrav measurement. The discrepancy between the predicted and
the three measured values ranges from 2–4.5$\sigma$. Red open symbols
correspond to the scaling relations from MM13, while blue open symbols
correspond to KH13. The open circles, crosses, and diamonds use the VDF from
B10, B11, and T20, respectively, while the open stars were computed using the
GSMF from DS15.

##  III Can direct SMBH observations and GWs be reconciled?

The predicted characteristic strain at $f=$yr-1 for all of the combinations of
velocity dispersion and stellar mass functions with scaling relations
described in Sec. I, and the values measured by PTAs are shown in Fig. 4. The
discrepancy between the different predictions and measurements varies between
$2-4.5\sigma$, assuming Gaussian errors. Based on the discussion presented in
Sec.II, we explore two aspects of the SGWB prediction that may be subject to
change: the SMBH mass function and the merger history. We find that an
increase in the merger rate of black holes is insufficient to boost the
predicted characteristic strain to the value measured by PTAs.

Motivated by the general agreement between the black hole mass density
inferred from local scaling relations and from the quasar luminosity function
through Sołtan’s argument for a radiative efficienty around $\epsilon_{r}\sim
0.1$, we begin by considering a scenario in which we require the local mass
density of BHs to be unchanged. We also consider more general variations to
the BH abundance, such as the existence of a new population of SMBHs, and
infer their required mass in order to match the PTA measurement.

Figure 5: Cumulative contribution of residuals of the $M-\sigma$ relation. The
dark blue solid line shows the residuals of the scaling relation fit from the
data presented in Ref. mcconnell_ma , the dotted line shows a Gaussian fit to
the residuals, and the light blue solid line shows the cumulative contribution
to the characteristic strain from the same fit. Note that the scatter in the
observed data includes both an intrinsic component and a contribution from
measurement errors. The results shown above include both and therefore
correspond to a slightly larger scatter (0.43) than the intrinsic (0.38).
Figure 6: Parameter values required to match both the measured characteristic
strain from PTAs and the inferred black hole mass density. The panel on the
left shows the curves of constant $h_{c}$ (in dark blue) and $\rho_{\rm BH}$
(in red). The best-fit parameters from direct observations of SMBHs and host-
galaxy properties correspond to the point $x_{\rm fid}=(159.6,0.38)$. The
points $x_{1}=(123,0.84)$ and $x_{2}=(150,0.62)$ correspond to sets of
parameters where the predicted $h^{2}_{c}$ and $\rho_{\rm BH}$ exactly match
the central fiducial value, and where they match within a 90% confidence
interval, respectively. The points in the panel on the right are the BH data
used in Ref. mcconnell_ma . The inner (outer) contours show the regions where
50% (90%) of the contribution to the SGWB come from for the set of parameters
$x_{1}$ in dashed and $x_{2}$ in solid lines. The red curve corresponds to the
fiducial VDF and $M-\sigma$ relation, but with an additional population of BHs
associated with galaxies of a given $\sigma$. Each point along the red curve
therefore shows the required BH mass if a single BH population were to account
for the $h_{c}$ value measured by NANOGrav. The light blue dotted line and
shaded band shows Eq. 23 for $d_{\bullet}=10.5^{+1.5}_{-3.6}$.

###  III.1 Changes to galaxy scaling relations

As noted in Sec. II, if we require the characteristic strain to be fixed to
its measured value and the BH mass density to be consistent with the current
estimates, then the relevant outcome of any variation of the VDF or the
$M-\sigma$ parameters is to change the mass that most contributes to the SGWB.
Fig.3 therefore allows us to infer that the characteristic strain must be
dominated by $\sim 3\times 10^{10}M_{\odot}$ black holes, roughly 10 times
larger than previously predicted.

We explore variations in parameters that lead to such enhancements in the
abundance or the typical masses of the most massive SMBHs, such as the
intrinsic scatter $\epsilon_{0}$ in the $M-\sigma$ relation and the
characteristic velocity dispersion $\sigma_{*}$. A larger $\sigma_{*}$ boosts
the abundance of high velocity dispersion galaxies, while a higher intrinsic
scatter in the $M-\sigma$ relation allows for more massive black holes to be
associated with galaxies of a given $\sigma$. It is worth noting that, since
the characteristic strain is weighted by $M^{5/3}$, it will always receive
contributions from more massive outliers of the $M-\sigma$ relation, even for
the best-fit value of $\epsilon_{0}$ from local observations. This feature is
shown in Fig. 5, using the black hole data from Ref. mcconnell_ma . Assuming
the best-fit value of scatter results in 50% of the contribution to $h_{c}$
being sourced from black holes that lie in the 5% tail of the $M-\sigma$
relation.

We show in Fig. 6 the curves of constant black hole mass density and
characteristic strain in the $\sigma_{*}$ and $\epsilon_{0}$ parameter space.
Requiring both values to be consistent with current measurements therefore
requires particular values of such parameters, which correspond to the
intersection of the curves shown. We highlight two points of intersection:
$x_{1}$, which results in an exact match to the central value of both the
measured $h^{2}_{c}$ and inferred $\rho_{\rm BH}$, while $x_{2}$ corresponds
to a set of parameters that lead to predicted values of $h^{2}_{c}$ and
$\rho_{\rm BH}$ that are consistent with observations at a 90% confidence
level.

The point $x_{2}$ has a very similar value of $\sigma_{*}$ to the fiducial,
but a significantly larger intrinsic scatter in the $M-\sigma$ relation. Since
the effect of scatter on the final SMBH mass function is to increase the
abundance of the most massive black holes, this appears to be a
straightforward mechanism to boost the characteristic strain, while minimizing
the change to the BH mass density. However, the required value of scatter is
in significant disagreement with local observations. This can be seen in the
residuals of the $M-\sigma$ relation and, equivalently, in the right panel of
Fig. 6. The points with error bars show the data set used in Ref. mcconnell_ma
, while the contours show the regions that contribute 50% and 90% of the SGWB.
This illustrates that changing the parameters of the VDF or the $M-\sigma$
relation in order to increase the GW prediction leads to the conclusion that
the majority of the contribution to the SGWB would be sourced by host galaxies
with associated SMBHs that have never been observed.

Motivated by the apparent upwards trend of the $M-\sigma$ relation for high
velocity dispersion galaxies, we implement an $M-\sigma$ relation given by the
double power law

$M_{\rm BH}=10^{a_{\bullet}}\left(\frac{\sigma}{200\text{km
s}^{-1}}\right)^{b_{\bullet}}+10^{c_{\bullet}}\left(\frac{\sigma}{200\text{km
s}^{-1}}\right)^{d_{\bullet}},$ (23)

where $a_{\bullet}$, $b_{\bullet}$, and a single value of scatter
$\epsilon_{0}$ are fixed to the values fit in MM13, and $c_{\bullet}$ is fixed
such that the two power laws intercept at $\sigma=250$km s-1. In order to
match $h_{c}(f=\text{yr}^{-1})=2.4^{+0.7}_{-0.6}\times 10^{-15}$, we find
$d_{\bullet}=10.5^{+1.5}_{-3.6}$. It has been argued in the literature that an
upwards trend in the $M-\sigma$ relation may be driven by the velocity
dispersion saturating at high values (which is not expected or observed for
the other host galaxy properties), and therefore that bulge mass or luminosity
scaling relations offer a more reliable prediction of the massive end of the
black hole mass function. However, it is worth emphasizing that since we find
a consistent value of the characteristic strain for predictions based on
$\sigma$ and $M_{\rm bulge}$, then the steepening required to match the PTA
measurements must necessarily be larger than this effect, as it would also
require a boost in the $M_{\rm bulge}$-based prediction.

Although the BH mass density has been a useful point of comparison between
SMBH mass functions estimated from local scaling relations and inferences
based on the quasar luminosity function, it is not a directly measured
quantity. If we instead allow $\rho_{\rm BH}$ to assume any value,
$\sigma_{*}$ and $\epsilon_{0}$ may lie in any point along the dark blue curve
in Fig. 6. We may also suppose that the measured VDF and $M-\sigma$ relation
are fixed, but that there is a new population of SMBHs with a typical mass
$M_{\rm new}$ hosted by galaxies with $\sigma_{\rm new}$, in addition to the
standard population. The red curve on the right panel in Fig. 6 shows the
required value of $M_{\rm new}$ as a function of $\sigma=\sigma_{\rm
new}$111We assume a Gaussian bump centered at each value of $\sigma$ with
width of 10km s-1 on top of the standard $M-\sigma$ relation, boosting SMBH
masses by the amount required to match $h_{c}$ measured by NANOGrav. The red
curve shows the required SMBH mass when galaxies with each value of $\sigma$
are boosted.. As expected from our previous discussion, boosting the masses of
SMBHs near the peak of the mass kernel of the SGWB leads to the smallest
change relative to its fiducial value, while high velocity dispersion hosts
require SMBHs of increasing masses to overcome their decreasing abundances and
the required SMBH masses associated with low velocity dispersion hosts
plateaus at $\sim 10^{9}M_{\odot}$ as a consequence of the plateau in the VDF,
leading to values that are in significant disagreement with observations.

In Appendix B we compare the population of SMBHs required to explain the PTA
results with those measured locally and inferred using local scaling
relations. Fig. 9 (also Fig. 1) shows that the best-fit PTA results can be
explained by a SMBH space-density in the local universe corresponding to a
yet-undetected population of order 10 $3\times 10^{10}M_{\odot}$ SMBHs within
100 Mpc (a volume inspired by the MASSIVE survey of massive galaxies and SMBHs
2014ApJ…795..158M ). This required space-density of very massive BHs is
independent of whether the population is due to non-Gaussian scatter in BH-
galaxy correlations or the steepening in such correlations highlighted above.

###  III.2 Changes to the merger history

Figure 7: Impact of multiple merger events on the predicted characteristic
strain. The panel on the left shows the redshift distribution of the first 3
merger events, assuming a fixed timescale of $\tau=4$ Gyr between subsequent
mergers. The panel on the right shows the cumulative increase of the predicted
characteristic strain as a function of number of mergers. The orange points
show the increase for instantaneous mergers, the red points correspond to a 4
Gyr time-delay, and the grey arrow shows the upper limit computed in Eq. 24.

In the previous section, we assumed that all black holes in the present-day
Universe have undergone one merger event in order to connect the single SMBH
mass function to the binary total mass function. Here we investigate whether
changes to the merger history of SMBHs can sufficiently increase the amplitude
of the SGWB.

We begin with simple and unphysical, but informative example. Ignoring the
time scale for the black holes to merge and assuming that all mergers occur
between equal mass black holes, we can estimate the boost to the
characteristic strain from multiple mergers by noting that at each earlier
generation $N$, the number of black holes will double and their masses will be
halved. That is, we assume that the present-day single black holes of masses
$M$ were formed from from equal-mass parents of masses $M/2$ that were twice
as abundant, and so on. This leads to the following convergent series

$\left.h^{2}_{c}\right|_{N\rightarrow\infty}=\left.h^{2}_{c}\right|_{N=1}\left(\sum_{N=1}^{\infty}2^{-2(N-1)/3}\right)\approx
2.7\left.h^{2}_{c}\right|_{N=1}.$ (24)

Hence, even in the extreme scenario of an infinite number of instantaneous
mergers, the characteristic strain $h_{c}$ can only be a factor of 1.64 larger
than our previously predicted value, which is still lower than the best-fit
value of any of the published PTA results.

Instead of an upper limit, we will now compute a fiducial estimate of this
boost, assuming that they occur with a certain distribution of mass ratios and
redshifts. We consider two compounding effects at each previous merger: the
shift of the redshift distribution and the reduced efficiency of GW emission
relative to an equal mass merger. We assume that the redshift distribution for
the latest merger $p_{N=1}$ is given by Eq. 18, and compute the distribution
of the previous merger $p_{N+1}$ by considering a fixed timescale $\tau$
between subsequent mergers. Every redshift $z$ can therefore be mapped onto
the redshift of the previous merger $z_{N+1}=z_{N}+\Delta z(\tau)$. Since some
fraction of objects will never merge (i.e. $t(z)+\tau$ may be larger than the
age of the Universe), the new distribution function is normalized to the
fraction that merges, which is not necessarily 1. We show in Fig. 7 a case
with $\tau=4$ Gyr and the resulting redshift distribution of each previous
merger.

The effect of unequal mass mergers can be computed by considering that the
progenitors of each black hole with a final mass $M_{\rm BH}$ will have masses
$M_{\rm BH}/(1+q)$ and $qM_{\rm BH}/(1+q)$, where $q$ is drawn from the
distribution function defined in Eq. 18. We show in App. A that the boost to
the characteristic strain after $N$ mergers predicted through a Monte Carlo
sampling of the mass ratio distribution is nearly identical to the prediction
with a fixed mass ratio given by the mean $\langle q\rangle$. Furthermore, we
derive an equivalent series to Eq. 24 for the unequal mass case. The highest
boost is obtained for $q=1$ given by 24 but it does not change much with $q$,
in the limit $q\rightarrow 0$ the boost factor is $12/5=2.4$.

The combined effect of the finite merger timescale and the unequal mass
mergers is a smaller boost to the characteristic strain as a function of
number of mergers and a convergence to a lower maximum value, since the finite
merger time implies a cutoff to the contributions. The cumulative contribution
of mergers is shown in the right panel of Fig. 7. The timescale of 4 Gyr is
inspired by the state-of-the-art simulations presented in Ref. Kelley:2016gse
, which is the average lifetime of a $\sim 10^{9}M_{\odot}$ black hole. Note
that this is still optimistic in several ways. For instance, it neglects the
timescale of the galaxy/halo merger and that black holes grow through
accretion, which would result in a significantly lower mass for the parent
SMBHs at earlier times. We therefore conclude that, while an increase in the
merger rate may alleviate the tension between the predicted amplitude of the
SGWB and the values measured by PTAs, it is insufficient to fully account for
the discrepancy.

##  IV Conclusion

The detection of nanohertz gravitational waves by several PTA collaborations
has recently opened a new window into the supermassive black hole population,
that is complementary to the picture offered by existing direct local
measurements and by quasars in the distant Universe. The GW measurement is
particularly sensitive to the most massive SMBHs, which are rare and therefore
challenging to find through direct observations.

By comparing GW and electromagnetic measurements of SMBHs (following
Phinney:2001di ), we find an inconsistency in the inferred SMBH population,
leading to a lower predicted value of the SGWB for mass functions consistent
with local observations. We show that this difference cannot be accounted for
by boosting the merger rate of SMBHs and is therefore likely a consequence of
an underestimated mass function of SMBHs. We show that if the current estimate
of the black hole mass density, consistent with the integrated luminosity of
quasars through the Sołtan argument, is to be preserved, then the PTA
measurement generally implies that the SMBH mass that most contributes to the
SGWB must be $\sim 3\times 10^{10}M_{\odot}$, about ten times higher than
would be naively predicted using local measurements of the SMBH mass function.
This would suggest that the most massive SMBHs found in our local Universe
comprise the typical population contributing to the SGWB.

We explore a few potential avenues that could lead to a higher prediction of
the high mass tail of the SMBH mass function. We begin by modifying the
intrinsic scatter in the $M-\sigma$ relation and the turnover scale in the
VDF, but find that the required parameters to account for the PTA measurement
are inconsistent with the residuals of the $M-\sigma$ fit. We then show that a
steeper $M-\sigma$ relation for high velocity dispersion galaxies can account
for the missing high-mass population. Although it has been argued that the
$M-\sigma$ relation should steepen as a consequence of the saturation of the
velocity dispersion at high values in gas-poor mergers, we emphasize that
since both the $\sigma$\- and $M_{\rm bulge}$-based predictions lead to
consistent and underpredicted values of the characteristic strain, then the
required boost to the abundance of the most massive black holes must
necessarily be larger than this effect. This exploration leads to a directly
testable prediction for the local abundance of the most massive SMBHs and
therefore highlights the importance of populating this region of parameter
space with local observations. Quantitatively, the best-fit PTA results
require a SMBH space-density in the local universe that corresponds to a yet-
undetected population of order 10 $3\times 10^{10}M_{\odot}$ or 1
$10^{11}M_{\odot}$ SMBHs within 100 Mpc (a volume in which stellar and gas
dynamical BH measurements are feasible). However, masses significantly above a
few $10^{10}M_{\odot}$ are disfavored since they would likely be detected as
individual sources by PTAs.

Although the GW signal alone offers only an integrated measurement of the SMBH
mass function, we showed that combining GWs with local observations of SMBHs
and the quasar luminosity function leads to a greater insight about this black
hole population. Crucially, GW measurements probe the massive end of the SMBH
mass function, which are rare and therefore elusive in electromagnetic
observations. While some black holes of masses $\sim 10^{10}M_{\odot}$ have
been found locally, the low surface brightness of the cored inner region
hinders direct mass measurements. Furthermore, since these extremely massive
SMBHs are most likely hosted by galaxies that appear to have been formed from
repeated gas-poor mergers, they are also unlikely to be seen as the brightest
quasars. GWs may therefore be offering a unique perspective on the most
massive black holes in the Universe.

###### Acknowledgements.

We would like to thank Jenny Greene, Scott Tremaine, Luke Zoltan Kelley, and
Nick Kokron, for helpful discussions. EQ and MZ thank Roger Blandford for a
cryptic comment about the NANOGrav results at the Simons Symposium on
Multiscale Physics, which motivated this work. GSP gratefully acknowledges
support from the Friends of the Institute for Advanced Study Fund. MZ is
supported by NSF 2209991 and NSF-BSF 2207583. EQ is supported in part by a
Simons Investigator Award from the Simons Foundation.

## References

  * (1) R. W. Hellings and G. S. Downs, “Upper limits on the isotropic gravitational radiation background from pulsar timing analysis.,” Astrophys. J. Lett. 265 (Feb., 1983) L39–L42.
  * (2) EPTA Collaboration, J. Antoniadis et al., “The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals,” arXiv:2306.16214 [astro-ph.HE].
  * (3) NANOGrav Collaboration, G. Agazie et al., “The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background,” Astrophys. J. Lett. 951 no. 1, (2023) L8, arXiv:2306.16213 [astro-ph.HE].
  * (4) D. J. Reardon et al., “Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array,” Astrophys. J. Lett. 951 no. 1, (2023) L6, arXiv:2306.16215 [astro-ph.HE].
  * (5) H. Xu et al., “Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I,” Res. Astron. Astrophys. 23 no. 7, (2023) 075024, arXiv:2306.16216 [astro-ph.HE].
  * (6) International Pulsar Timing Array Collaboration, G. Agazie et al., “Comparing recent PTA results on the nanohertz stochastic gravitational wave background,” arXiv:2309.00693 [astro-ph.HE].
  * (7) S. Burke-Spolaor et al., “The Astrophysics of Nanohertz Gravitational Waves,” Astron. Astrophys. Rev. 27 no. 1, (2019) 5, arXiv:1811.08826 [astro-ph.HE].
  * (8) J. Kormendy and D. Richstone, “Inward Bound—The Search For Supermassive Black Holes In Galactic Nuclei,” Annu. Rev. Astron. Astrophys. 33 (Jan., 1995) 581.
  * (9) D. Richstone, E. A. Ajhar, R. Bender, G. Bower, A. Dressler, S. M. Faber, et al., “Supermassive black holes and the evolution of galaxies.,” Nature (London) 385 no. 6701, (Oct., 1998) A14, arXiv:astro-ph/9810378 [astro-ph].
  * (10) M. C. Begelman, R. D. Blandford, and M. J. Rees, “Massive black hole binaries in active galactic nuclei,” Nature (London) 287 no. 5780, (Sept., 1980) 307–309.
  * (11) J. Kormendy and L. C. Ho, “Coevolution (Or Not) of Supermassive Black Holes and Host Galaxies,” Annu. Rev. Astron. Astrophys. 51 no. 1, (Aug., 2013) 511–653, arXiv:1304.7762 [astro-ph.CO].
  * (12) N. J. McConnell and C.-P. Ma, “Revisiting the Scaling Relations of Black Hole Masses and Host Galaxy Properties,” Astrophys. J.  764 no. 2, (Feb., 2013) 184, arXiv:1211.2816 [astro-ph.CO].
  * (13) A. Marconi, G. Risaliti, R. Gilli, L. K. Hunt, R. Maiolino, and M. Salvati, “Local supermassive black holes, relics of active galactic nuclei and the X-ray background,” Mon. Not. R. Astron. Soc. 351 no. 1, (June, 2004) 169–185, arXiv:astro-ph/0311619 [astro-ph].
  * (14) M. Vika, S. P. Driver, A. W. Graham, and J. Liske, “The Millennium Galaxy Catalogue: the Mbh-Lspheroid derived supermassive black hole mass function,” Mon. Not. R. Astron. Soc. 400 no. 3, (Dec., 2009) 1451–1460, arXiv:0908.2102 [astro-ph.CO].
  * (15) F. Shankar, D. H. Weinberg, and J. Miralda-Escudé, “Self-Consistent Models of the AGN and Black Hole Populations: Duty Cycles, Accretion Rates, and the Mean Radiative Efficiency,” Astrophys. J.  690 no. 1, (Jan., 2009) 20–41, arXiv:0710.4488 [astro-ph].
  * (16) F. Shankar, “Black hole demography: from scaling relations to models,” Classical and Quantum Gravity 30 no. 24, (Dec., 2013) 244001, arXiv:1307.3289 [astro-ph.CO].
  * (17) A. Soltan, “Masses of quasars.,” Mon. Not. R. Astron. Soc. 200 (July, 1982) 115–122.
  * (18) A. Cavaliere, P. Morrison, and K. Wood, “On Quasar Evolution,” Astrophys. J.  170 (Dec., 1971) 223.
  * (19) T. A. Small and R. D. Blandford, “Quasar evolution and the growth of black holes.,” Mon. Not. R. Astron. Soc. 259 (Dec., 1992) 725–737.
  * (20) Q. Yu and S. Tremaine, “Observational constraints on growth of massive black holes,” Mon. Not. R. Astron. Soc. 335 no. 4, (Oct., 2002) 965–976, arXiv:astro-ph/0203082 [astro-ph].
  * (21) F. Shankar, D. H. Weinberg, and J. Miralda-Escudé, “Accretion-driven evolution of black holes: Eddington ratios, duty cycles and active galaxy fractions,” Mon. Not. R. Astron. Soc. 428 no. 1, (Jan., 2013) 421–446, arXiv:1111.3574 [astro-ph.CO].
  * (22) M. Tucci and M. Volonteri, “Constraining supermassive black hole evolution through the continuity equation,” Astron. Astrophys. 600 (2017) A64, arXiv:1603.00823 [astro-ph.GA].
  * (23) E. S. Phinney, “A Practical theorem on gravitational wave backgrounds,” arXiv:astro-ph/0108028.
  * (24) NANOGrav Collaboration, G. Agazie et al., “The NANOGrav 15 yr Data Set: Constraints on Supermassive Black Hole Binaries from the Gravitational-wave Background,” Astrophys. J. Lett. 952 no. 2, (2023) L37, arXiv:2306.16220 [astro-ph.HE].
  * (25) J. Kormendy and R. Bender, “The Lvprop$\sigma$8 Correlation for Elliptical Galaxies with Cores: Relation with Black Hole Mass,” Astrophys. J. Lett. 769 no. 1, (May, 2013) L5.
  * (26) J. Thomas, R. P. Saglia, R. Bender, P. Erwin, and M. Fabricius, “The Dynamical Fingerprint of Core Scouring in Massive Elliptical Galaxies,” Astrophys. J. 782 no. 1, (2014) 39, arXiv:1311.3783 [astro-ph.GA].
  * (27) J. Hlavacek-Larrondo, A. C. Fabian, A. C. Edge, and M. T. Hogan, “On the hunt for ultramassive black holes in brightest cluster galaxies,” Mon. Not. R. Astron. Soc. 424 no. 1, (July, 2012) 224–231, arXiv:1204.5759 [astro-ph.CO].
  * (28) T. Ebisuzaki, J. Makino, and S. K. Okumura, “Merging of two galaxies with central black holes,” Nature (London) 354 no. 6350, (Nov., 1991) 212–214.
  * (29) S. M. Faber, S. Tremaine, E. A. Ajhar, Y.-I. Byun, A. Dressler, K. Gebhardt, et al., “The Centers of Early-Type Galaxies with HST. IV. Central Parameter Relations.,” Astron. J. 114 (Nov., 1997) 1771, arXiv:astro-ph/9610055 [astro-ph].
  * (30) M. Milosavljević and D. Merritt, “Formation of Galactic Nuclei,” Astrophys. J.  563 no. 1, (Dec., 2001) 34–62, arXiv:astro-ph/0103350 [astro-ph].
  * (31) N. Sahu, A. W. Graham, and B. L. Davis, “Revealing Hidden Substructures in the M BH-$\sigma$ Diagram, and Refining the Bend in the L-$\sigma$ Relation,” Astrophys. J.  887 no. 1, (Dec., 2019) 10, arXiv:1908.06838 [astro-ph.GA].
  * (32) R. C. E. van den Bosch, “Unification of the fundamental plane and Super Massive Black Hole Masses,” Astrophys. J.  831 no. 2, (Nov., 2016) 134, arXiv:1606.01246 [astro-ph.GA].
  * (33) C. Marsden, F. Shankar, M. Ginolfi, and K. Zubovas, “The case for the fundamental $M_{\rm BH}$-$\sigma$ relation,” arXiv:2004.00098 [astro-ph.GA].
  * (34) T. R. Lauer, S. M. Faber, D. Richstone, K. Gebhardt, S. Tremaine, M. Postman, et al., “The Masses of Nuclear Black Holes in Luminous Elliptical Galaxies and Implications for the Space Density of the Most Massive Black Holes,” Astrophys. J.  662 no. 2, (June, 2007) 808–834, arXiv:astro-ph/0606739 [astro-ph].
  * (35) M. Bernardi, F. Shankar, J. B. Hyde, S. Mei, F. Marulli, and R. K. Sheth, “Galaxy luminosities, stellar masses, sizes, velocity dispersions as a function of morphological type,” Mon. Not. R. Astron. Soc. 404 no. 4, (June, 2010) 2087–2122, arXiv:0910.1093 [astro-ph.CO].
  * (36) R. Bezanson, P. G. van Dokkum, M. Franx, G. B. Brammer, J. Brinchmann, M. Kriek, et al., “Redshift Evolution of the Galaxy Velocity Dispersion Function,” Astrophys. J. Lett. 737 no. 2, (Aug., 2011) L31, arXiv:1107.0972 [astro-ph.CO].
  * (37) R. D’Souza, S. Vegetti, and G. Kauffmann, “The massive end of the stellar mass function,” Mon. Not. R. Astron. Soc. 454 no. 4, (Dec., 2015) 4027–4036, arXiv:1509.07418 [astro-ph.GA].
  * (38) L. Taylor, R. Bezanson, A. van der Wel, A. Pearl, E. F. Bell, F. D’Eugenio, et al., “The Velocity Dispersion Function for Massive Quiescent and Star-forming Galaxies at 0.6 ¡ z $\leq$ 1.0,” Astrophys. J.  939 no. 2, (Nov., 2022) 90, arXiv:2210.00144 [astro-ph.GA].
  * (39) C. Li and S. D. M. White, “The distribution of stellar mass in the low-redshift Universe,” Mon. Not. R. Astron. Soc. 398 no. 4, (Oct., 2009) 2177–2187, arXiv:0901.0706 [astro-ph.CO].
  * (40) M. Bernardi, A. Meert, R. K. Sheth, V. Vikram, M. Huertas-Company, S. Mei, and F. Shankar, “The massive end of the luminosity and stellar mass functions: Dependence on the fit to the light profile,” Mon. Not. Roy. Astron. Soc. 436 (2013) 697, arXiv:1304.7778 [astro-ph.CO].
  * (41) P. F. Hopkins, G. T. Richards, and L. Hernquist, “An Observational Determination of the Bolometric Quasar Luminosity Function,” Astrophys. J.  654 no. 2, (Jan., 2007) 731–753, arXiv:astro-ph/0605678 [astro-ph].
  * (42) L. Z. Kelley, L. Blecha, and L. Hernquist, “Massive Black Hole Binary Mergers in Dynamical Galactic Environments,” Mon. Not. Roy. Astron. Soc. 464 no. 3, (2017) 3131–3157, arXiv:1606.01900 [astro-ph.HE].
  * (43) C.-P. Ma, J. E. Greene, N. McConnell, R. Janish, J. P. Blakeslee, J. Thomas, and J. D. Murphy, “The MASSIVE Survey. I. A Volume-limited Integral-field Spectroscopic Study of the Most Massive Early-type Galaxies within 108 Mpc,” Astrophys. J.  795 no. 2, (Nov., 2014) 158, arXiv:1407.1054 [astro-ph.GA].
  * (44) NANOGrav Collaboration, G. Agazie et al., “The NANOGrav 15 yr Data Set: Bayesian Limits on Gravitational Waves from Individual Supermassive Black Hole Binaries,” Astrophys. J. Lett. 951 no. 2, (2023) L50, arXiv:2306.16222 [astro-ph.HE].
  * (45) J. Thomas, C.-P. Ma, N. J. McConnell, J. E. Greene, J. P. Blakeslee, and R. Janish, “A 17-billion-solar-mass black hole in a group galaxy with a diffuse core,” Nature (London) 532 no. 7599, (Apr., 2016) 340–342, arXiv:1604.01400 [astro-ph.GA].

Figure 8: Efficiency of GW emission relative to the equal mass case. Symbols
corresponds to the results of merger tree simulations and the lines the series
of equation 27.

##  Appendix A Efficiency of unequal mass mergers

Assume all mergers have a fixed mass ratio $q$ and corresponding $\eta$. Let
us define $S(N,q)$ the efficiency ratio relative to a single equal mass merger
of a tree of mergers that goes down $n$ levels. All merger have an efficiency
ratio relative to the equal mass case of $4\eta$. After one merger the two
progenitor black holes have masses $M_{\rm BH}/(1+q)$ and $M_{\rm BH}q/(1+q)$
so:

$S(0,q)=4\eta\\\ \\\ ;\\\ \\\
S(1,q)=4\eta\left(1+\frac{1}{(1+q)^{5/3}}+\frac{q^{5/3}}{(1+q)^{5/3}}\right).$
(25)

After a second level, the black holes split again and there is one black hole
with mass $M_{\rm BH}1/(1+q)^{2}$, two with masses $M_{\rm BH}q/(1+q)^{2}$,
and one with $M_{\rm BH}(q/(1+q))^{2}$, so that:

$S(2,q)=4\eta\left(1+\frac{1+q^{5/3}}{(1+q)^{5/3}}+(\frac{1+q^{5/3}}{(1+q)^{5/3}})^{2}\right).$
(26)

To go down each one level, each of the black holes split into pieces by
factors $1/(1+q)$ and $q/(1+q)$, so there are many possible masses. The number
of black holes of the different masses however are given by the binomial
coefficients so that in general:

$S(N,q)=4\eta\sum_{n=0}^{N}\left(\frac{1+q^{5/3}}{(1+q)^{5/3}}\right)^{n}=4\eta\frac{1-\left(\frac{1+q^{5/3}}{(1+q)^{5/3}}\right)^{N+1}}{1-\left(\frac{1+q^{5/3}}{(1+q)^{5/3}}\right)}$
(27)

For a fixed $N$, $S(N,q)$ is maximum for $q=1$.

Rather than label the tree by the number of levels, one could label it by the
mass of the largest black hole in that level, $M_{\rm BH}/(1+q)^{N+1}$. The
analytic prediction in the figure simply plots $(1+q)^{N+1}$ vs $S(N,q)$ where
$q$ is the average value for each distribution, where $\langle
q\rangle=[0.3909,0.55,1]$ for the log, linear and equal mass cases (assumes
$q_{min}=0.1$). In these formulas $N$ does not need to be an integer. The
discrete points show the results of averages of Monte Carlo realization where
the mass ratio was sampled form the appropriate distribution. Each subsequent
point corresponds to one additional level of the tree.

##  Appendix B Implication for local observations

Figure 9: Total number of SMBHs above a given minimum mass $M_{\rm min}$.
Solid lines with error bands correspond to predictions using the fiducial VDF
and $M-\sigma$ relation (grey) and the fiducial VDF and an $M-\sigma$ relation
given by Eq. 23 with the best-fit value that matches the characteristic strain
measured by PTAs (light blue). The dashed line with dots correspond to the
fiducial mass function, with a Gaussian bump above $M_{\rm min}$ such that the
amplitude is chosen to match the measured $h_{c}$. The dark blue curve shows
the SMBHs observed so far.

The higher amplitude of the SGWB found by PTA observations may imply a
“missing” populations of SMBHs in the current estimates of the mass function.
While the PTA signal alone cannot distinguish their properties, requiring
consistency with the current estimate of the black hole mass density would
suggest an increase in the high-mass end of the SMBH mass function. If this is
the case, we should expect an increase in the number of very massive SMBHs in
nearby galaxies and we provide a rough estimate of this population.

Inspired by the ongoing observations of the MASSIVE Survey 2014ApJ…795..158M ,
we show in Fig. 9 the number of SMBHs above a certain minimum mass $M_{\rm
min}$ within a survey with a spherical volume of 100 Mpc radius. The solid
lines correspond to predictions for a single mass function, the lower computed
with out fiducial VDF and $M-\sigma$ relation, and the higher with the
fiducial VDF and a steeper $M-\sigma$ relation that reproduced the central
value of $h_{c}$ observed by NANOGrav, and error bands which correspond to
Poisson errors. The dashed line corresponds to a SMBH mass function identical
to the fiducial, but with a Gaussian bump added above $M_{\rm min}$ with width
of 0.1 dex and with a height chosen so as to match $h_{c}$. Each point in the
dashed line therefore corresponds to a different mass function, and should be
interpreted as a statement that the characteristic strain observed by NANOGrav
may be explained by a mass function that predicts $\sim 60$ SMBHs with masses
above $\gtrsim 10^{10}M_{\odot}$ or $\sim 2$ SMBHs with masses above $\gtrsim
10^{11}M_{\odot}$ within a distance of 100Mpc. However, masses significantly
above $10^{10}M_{\odot}$ are disfavored by upper limits in individual SMBHB
searches, as they would likely be detected as point sources in PTA
observations NANOGrav:2023pdq . The line of observed black holes includes the
data from Ref. mcconnell_ma , with the addition of the SMBH reported in Ref.
2016Natur.532..340T . Note that this does not correspond to a complete sample
within the assumed volume, unlike the sample measured by the MASSIVE Survey.


# On $p$-adic solutions to KZ equations, ordinary crystals, and
$p^{s}$-hypergeometric solutions

Alexander Varchenko⋆ and  Vadim Vologodsky∘

###### Abstract.

In this paper we consider the KZ connection associated with a family of
hyperelliptic curves of genus $g$ over the ring of $p$-adic integers
$\mathbb{Z}_{p}$. Then the dual connection is the Gauss-Manin connection of
that family. We observe that the Gauss-Manin connection has an ordinary
$F$-crystal structure and its unit root subcrystal is of rank $g$. We prove
that all local flat sections of the KZ connection annihilate the unite root
subcrystal, and the space of all local flat sections of the KZ connection is a
free $\mathbb{Z}_{p}$-module of rank $g$.

We also consider the reduction modulo $p^{s}$ of the unit root subcrystal for
any $s\geq 1$. We prove that its annihilator is generated by the so-called
$p^{s}$-hypergeometric flat sections of the KZ connection. In particular, that
means that the reduction modulo $p^{s}$ of an arbitrary local flat section of
the KZ connection over $\mathbb{Z}_{p}$ is a linear combination of the
$p^{s}$-hypergeometric flat sections.

⋆Department of Mathematics, University of North Carolina at Chapel Hill
Chapel Hill, NC 27599-3250, USA

∘ Institute for Advanced Study, Princeton, NJ 08540, USA

††footnotetext: ⋆E-mail<EMAIL_ADDRESS>
∘E-mail<EMAIL_ADDRESS>supported by a Simons Foundation Investigator
Grant No. 622511 through Bhargav Bhatt.

###### Contents

  1. 1 Introduction
    1. 1.1 KZ equations
    2. 1.2 Gauss-Manin connection
    3. 1.3 Solutions modulo $p^{s}$
    4. 1.4 Formal $p$-adic schemes over $\mathbb{Z}_{p}$
    5. 1.5 $F$-crystal
    6. 1.6 $F$-crystal and $p^{s}$-hypergeometric sections
    7. 1.7 Pairing on $\mathcal{U}^{\prime}$
    8. 1.8 Local flat sections
  2. 2 Proof of Theorem 1.3
    1. 2.1 Some cohomology classes
    2. 2.2 Map $C_{s}$ and $p^{s}$-hypergeometric solutions
    3. 2.3 Proof of Theorem 1.3
  3. 3 Proof of Theorem 1.6

##  1\. Introduction

The Knizhnik-Zamolodchikov (KZ) equations were discovered in conformal field
theory, see [KZ, EFK]. Versions of these equations have been observed in
representation theory, algebraic geometry, combinatorics, field theory, … In
[CF, DF, DJMM, SV1] the KZ equations were identified with equations for flat
sections of suitable Gauss-Manin connections, and solutions to the KZ
equations were constructed in the form of multidimensional hypergeometric
integrals.

In [SV2, V3], the KZ equations were considered modulo $p^{s}$ where $p$ is a
prime number, and polynomial solutions to the KZ equations modulo $p^{s}$ were
constructed by an elementary procedure as suitable $p^{s}$-approximations of
the hypergeometric integrals.

In [VV], the question was addressed of whether all solutions to the KZ
equations in characteristic $p$ are generated by the $p$-hypergeometric
functions. The affirmative answer was demonstrated in a non-trivial example.

In [V3, VZ], the KZ equations associated with a family of hyperelliptic curves
of genus $g$ were studied. It was demonstrated that the module of
$p^{s}$-hypergeometric solutions has a $p$-adic limit as $s$ tends to
infinity. This limit is a KZ-invariant rank $g$ subbundle $\mathcal{W}$ of the
associated KZ connection.

In this paper we describe how the $p^{s}$-hypergeometric solutions and
subbundle $\mathcal{W}$ are related to the ordinary $F$-crystal of the Gauss-
Manin connection of that family of hyperelliptic curves. More precisely, we
consider the KZ connection associated with the family of hyperelliptic curves
over the ring of $p$-adic integers $\mathbb{Z}_{p}$. Then the dual connection
is the Gauss-Manin connection of that family. We observe that the Gauss-Manin
connection has an ordinary $F$-crystal structure and its unit root subcrystal
is of rank $g$. We prove that all local flat sections of the KZ connection
annihilate the unite root subcrystal, the space of all local flat sections of
the KZ connection is a free $\mathbb{Z}_{p}$-module of rank $g$, the space of
all local flat sections generate the subbundle $\mathcal{W}$.

We also consider the reduction modulo $p^{s}$ of the unit root subcrystal for
any $s\geq 1$. We prove that its annihilator is generated by the
$p^{s}$-hypergeometric flat sections of the KZ connection. In particular, that
means that the reduction modulo $p^{s}$ of an arbitrary local flat section of
the KZ connection over $\mathbb{Z}_{p}$ is a linear combination of the
$p^{s}$-hypergeometric solutions.

The general project is to analyze solutions of arbitrary KZ equations over
$\mathbb{Z}_{p}$ and relate them to $p^{s}$-hypergeometric solutions and
associated $F$-crystals. At this stage, we consider first non-trivial examples
to establish a map of the interrelations between these concepts.

###  1.1. KZ equations

Let $g$ be a positive integer. Denote $n=2g+1$.

In this paper we consider the following example of the KZ equations. Let
$\Omega_{ij}$ be the $n\times n$-matrix $(c_{kl})$ with
$c_{ij}=c_{ji}=-c_{ii}=-c_{jj}=1$ and all other entries equal to $0$. Let
$z=(z_{1},\dots,z_{n})$. Define

$\displaystyle H_{i}(z)=\sum_{j\neq i}\frac{\Omega_{ij}}{z_{i}-z_{j}},\quad
i=1,\dots,n,$

the Gaudin Hamiltonians. Consider the system of differential and algebraic
equations for column vectors $I(z)=(I_{1}(z)$, …, $I_{n}(z))^{\intercal}$:

(1.1) $\displaystyle\phantom{aaaaaa}2\partial_{i}I\ -\
\,H_{i}(z)\,I\,=\,0,\quad i=1,\dots,n,\quad I_{1}(z)+\dots+I_{n}(z)=0,$

where $\partial_{i}$ denotes $\frac{\partial}{\partial z_{i}}$;
$M^{\intercal}$ denotes the transpose matrix of a matrix $M$.

Let $\mathbb{A}^{n}$ be affine space with coordinates $z_{1},\dots,z_{n}$. Let
$S$ be the complement to the union of all diagonal hyperplanes defined by
equations $z_{i}=z_{j}$, $i\neq j$. The KZ equations define a flat KZ
connection $\nabla^{\operatorname{KZ}}$ on the trivial bundle $S\times V\to S$
whose fiber is the vector space $V$ of $n$-vectors with the zero sum of
coordinates. Denote this bundle by $\mathcal{V}$ and call the pair
$(\mathcal{V},\nabla^{\operatorname{KZ}})$ the KZ bundle.

######  Remark.

The system of equations (1.1) is the system of the original KZ equations
associated with the Lie algebra $\mathfrak{sl}_{2}$ and the subspace of
singular vectors of weight $n-2$ of the tensor power
$({\mathbb{C}}^{2})^{{\otimes}n}$ of the two-dimensional irreducible
$\mathfrak{sl}_{2}$-module, up to a gauge transformation, see this example in
[V1, Section 1.1]. The original KZ equations depend on a parameter $\kappa$.
In our paper, $\kappa=2$, see the coefficient of $\partial_{i}$ in (1.1).

###  1.2. Gauss-Manin connection

Consider the family $X$ of hyperelliptic curves $X(z)$ defined by the affine
equation

(1.2) $\displaystyle y^{2}=(x-z_{1})\dots(x-z_{n}).$

For $z\in S$, the curve $X(z)$ is a smooth projective curve of genus $g$. The
KZ connection $\nabla^{\operatorname{KZ}}$ is dual to the Gauss-Manin
connection $\nabla^{\operatorname{GM}}$ over $S$ with fiber
$H^{1}_{dR}(X(z))$.

For example, the flat sections of the KZ equations over the field
${\mathbb{C}}$ are labeled by flat sections $\gamma(z)$ of the Gauss-Manin
connection over $S$ with fiber $H_{1}(X(z),{\mathbb{C}})$ and given by the
following construction. Define the master function

(1.3) $\displaystyle\Phi(x,z)=\big{(}(x-z_{1})\dots(x-z_{n})\big{)}^{-1/2}$

and the ${n}$-vector of hyperelliptic integrals

(1.4) $\displaystyle
I^{(\gamma)}(z)=(I_{1}(z),\dots,I_{n}(z))^{\intercal},\qquad
I_{j}=\int_{\gamma(z)}\frac{dx}{(x-z_{j})y}=\int_{\gamma(z)}\frac{\Phi(x,z)dx}{x-z_{j}}\,.$

######  Theorem 1.1.

The vector $I^{(\gamma)}(z)$ satisfies the KZ equations (1.1). All complex
solutions to the KZ equations have this form. The complex vector space of
solutions of the form (1.4) is $n-1$-dimensional.

Theorem 1.1 is a classical statement discussed as an example in [V1, Section
1.1], also see [V1, Formula (1.3)].

The solutions in (1.4) are called the hypergeometric solutions to the KZ
equations.

###  1.3. Solutions modulo $p^{s}$

Let $p$ be an odd prime number. For an integer $s\geq 1$ define the master
polynomial

$\displaystyle\Phi_{s}(x,z)=\big{(}(x-z_{1})\dots(x-z_{n})\big{)}^{(p^{s}-1)/2}$

which can be considered as a $p^{s}$-approximation of the master function
$\Psi(x,z)$. For $\ell=1,\dots,g$, define the column $n$-vector

$\displaystyle
Q^{s,\ell}(z)=(Q^{s,\ell}_{1},\dots,Q^{s,\ell}_{n})^{\intercal}$

as the coefficient of $x^{\ell p^{s}-1}$ in the column $n$-vector of
polynomials
$\big{(}\frac{\Phi_{s}(x,z)}{x-z_{1}},\dots,\frac{\Phi_{s}(x,z)}{x-z_{n}}\big{)}^{\intercal}$.
If $\ell>g$, then the polynomial $\frac{\Phi_{s}(x,z)}{x-z_{i}}$ does not have
the monomial $x^{\ell p^{s}-1}$.

######  Theorem 1.2 ([V3, Theorems 4.4], [VZ, Theorem 6.9]).

The column $n$-vectors $Q^{s,\ell}(z),$ $\ell=1,\dots,g$, of polynomials in
$z$ are solutions of the KZ equations (1.1) modulo $p^{s}$. They generate a
free $(\mathbb{Z}/p^{s})[z]$-module of rank $g$.

The $n$-vectors $Q^{s,\ell}(z)$, $\ell=1,\dots,g$, are called the
$p^{s}$-hypergeometric solutions to the KZ equations (1.1).

###  1.4. Formal $p$-adic schemes over $\mathbb{Z}_{p}$

Recall some definitions. A formal $p$-adic scheme $Y/\mathbb{Z}_{p}$ is a
collection of schemes $Y_{s}/(\mathbb{Z}/p^{s})$ together with isomorphisms

$\displaystyle
Y_{s+1}\times_{\operatorname{Spec}\mathbb{Z}/p^{s+1}}\operatorname{Spec}\mathbb{Z}/p^{s}\simeq
Y_{s}\,.$

If $f_{1},\dots,f_{k}\in\mathbb{Z}_{p}[z_{1},\dots,z_{n}]$, then

$\displaystyle\mathbb{Z}_{p}[z_{1},\dots,z_{n},f_{1}^{-1},\dots,f_{k}^{-1}]^{\wedge}_{p}\
:\,=\
\lim_{\leftarrow}(\mathbb{Z}/p^{s})[z_{1},\dots,z_{n},f_{1}^{-1},\dots,f_{k}^{-1}]\,$

is the $p$-adic completion of
$\mathbb{Z}_{p}[z_{1},\dots,z_{n},f_{1}^{-1},\dots,f_{k}^{-1}]$. The formal
$p$-adic scheme associated with
$\mathbb{Z}_{p}[z_{1},\dots,z_{n},f_{1}^{-1},\dots,f_{k}^{-1}]^{\wedge}_{p}$
is denoted by

$\displaystyle\operatorname{Spf}\,\mathbb{Z}_{p}[z_{1},\dots,z_{n},f_{1}^{-1},\dots,f_{k}^{-1}]^{\wedge}_{p}.$

Let

$\mathfrak{S}={\operatorname{Spf}}\,{\mathbb{Z}}_{p}[z_{1},\dots,z_{n},(z_{i}-z_{j})^{-1},\;1\leq
i\neq j\leq n]^{\wedge}_{p}$

be the formal $p$-adic scheme of the complement to the union of diagonal
hyperplanes $z_{i}-z_{j}=0$, $1\leq i,j\leq n$.

In what follows we consider the KZ bundle $\mathcal{V}$ over the formal
$p$-adic scheme $\mathfrak{S}$ and its subschemes and reductions modulo
$p^{s}$.

In [V3, Section 5.1], a $g\times g$-matrix $A(z)$ was defined with entries in
$\mathbb{Z}[z]$ and determinant having a nonzero projection to
${\mathbb{F}}_{p}[z]$. The matrix $A(z)$ projected to ${\mathbb{F}}_{p}[z]$ is
the Hasse–Witt matrix of the curve $X(z)$ considered over ${\mathbb{F}}_{p}$.
The fact that $\det A(z)$ is nonzero over ${\mathbb{F}}_{p}$ means that the
curve $X(z)$ is ordinary in the terminology of [DI].

Let

$\mathfrak{D}={\operatorname{Spf}}\,{\mathbb{Z}}_{p}[z_{1},\dots,z_{n},(\det
A)^{-1},(z_{i}-z_{j})^{-1},\;1\leq i\neq j\leq n]^{\wedge}_{p}$

be the open formal $p$-adic subscheme of $\mathfrak{S}$ defined by $\det
A(z)\neq 0$.

Denote

(1.5) $\displaystyle\mathcal{O}(\mathfrak{S})$ $\displaystyle=$
$\displaystyle{\mathbb{Z}}_{p}[z_{1},\dots,z_{n},(z_{i}-z_{j})^{-1},\;1\leq
i\neq j\leq n]^{\wedge}_{p}\,.$ $\displaystyle\mathcal{O}(\mathfrak{D})$
$\displaystyle=$ $\displaystyle{\mathbb{Z}}_{p}[z_{1},\dots,z_{n},(\det
A)^{-1},(z_{i}-z_{j})^{-1},\;1\leq i\neq j\leq n]^{\wedge}_{p}\,.$

According to [VZ, Section 6.5], for any $s\geq 1$, the column vectors
$Q^{s,\ell}(z)$, $\ell=1,\dots,g$, generate a rank $g$ flat subbundle of the
KZ bundle $\mathcal{V}_{s}$ over $\mathfrak{D}_{s}$. The subbundle is denoted
by $\mathcal{W}_{s}$ . According to [VZ, Theorem 6.9], we have
$\mathcal{W}_{s+1}|_{\mathfrak{D}_{s}}=\mathcal{W}_{s}$ as subbundles of
$\mathcal{V}_{s}$. Thus, the sequence of bundles $\mathcal{W}_{s}$ has a
$p$-adic limit as $s\to\infty$, denoted by $\mathcal{W}$, where
$\mathcal{W}\subset\mathcal{V}|_{\mathfrak{D}}$.

###  1.5. $F$-crystal

Consider our family $X$ of hyperelliptic curves over $\mathfrak{S}$. Let
$H^{1}_{dR}(X/\mathfrak{S})$ be the relative de Rham cohomology of the family.
This is a bundle of rank $2g$ equipped with the Gauss-Manin connection
$\nabla^{\operatorname{GM}}$. We denote it by $\mathcal{E}$.

We identify $\mathcal{V}$ with $\mathcal{E}^{*}$ dual to $\mathcal{E}$, cf.
formula (1.4). This identification allows us to identify the flat sections of
the KZ connection with the $\nabla^{\operatorname{GM}}$-flat linear
functionals $I:\mathcal{E}\to\mathcal{O}(\mathfrak{S})$.

The map $z_{i}\mapsto z_{i}^{p}$, $i=1,\dots,n$, extends to a well-defined
algebra endomorphism
$F:\mathcal{O}(\mathfrak{S})\to\mathcal{O}(\mathfrak{S})$. Denote

$\displaystyle
F^{*}\mathcal{E}=\mathcal{O}(\mathfrak{S}){\otimes}_{\mathcal{O}(\mathfrak{S}),F}\mathcal{E}.$

Thus, $F^{*}\mathcal{E}$ is the quotient of
$\mathcal{O}(\mathfrak{S})\otimes_{\mathbb{Z}}\mathcal{E}$ by the subgroup
generated by $r\otimes r^{\prime}v-F(r^{\prime})r\otimes v$,
$r,r^{\prime}\in\mathcal{O}(\mathfrak{S})$, $v\in\mathcal{E}$. Let

$\displaystyle\phi:F^{*}\mathcal{E}\to\mathcal{E}$

be the crystalline Frobenius, see [BO]. The morphism $\phi$ respects the
connections. The pair $(\mathcal{E},\phi)$ is called an $F$-crystal.

The $F$-crystal $(\mathcal{E},\phi)$ restricted to $\mathfrak{D}$ is ordinary
in the sense of [DI] since $\det A(z)\neq 0$ on $\mathfrak{D}$. In particular,
by [DI], there is a unique exact sequence of bundles with connections over
$\mathfrak{D}$,

(1.6) $\displaystyle
0\to\mathcal{U}\hookrightarrow\mathcal{E}\to\mathcal{E}/\mathcal{U}\to 0,$

such that

  1. (i)

$\operatorname{rank}\,\mathcal{U}=g$,

  2. (ii)

$\phi\,|_{\mathcal{F}^{*}\mathcal{U}}:F^{*}\mathcal{U}\to\mathcal{U}$ is an
isomorphism,

  3. (iii)

$\phi\,|_{F^{*}(\mathcal{E}/\mathcal{U})}:F^{*}(\mathcal{E}/\mathcal{U})\to\mathcal{E}/\mathcal{U}$
has the form $p\cdot\phi_{1}$, where
$\phi_{1}:F^{*}(\mathcal{E}/\mathcal{U})\to\mathcal{E}/\mathcal{U}$ is an
isomorphism, and $p\,\cdot$ is multiplication by $p$.

  4. (iv)

$\mathcal{U}$ equals
$\bigcap_{i}\phi^{i}((F^{i})^{*}\mathcal{E})\subset\mathcal{E}$, where
$\phi^{i}$ is viewed as $(F^{i})^{*}\mathcal{E}\to\mathcal{E}$.

Note that (iv) is implied by (ii) and (iii).

The pair $(\mathcal{U},\phi)$ is called the unit root subcrystal of
$(\mathcal{E},\phi)$. We refer to sequence (1.6) as the slope sequence.

The KZ bundle $\mathcal{V}=\mathcal{E}^{*}$ has the dual $F$-crystal
structure,

(1.7) $\displaystyle
p\cdot(\phi^{*})^{-1}:F^{*}(\mathcal{E}^{*})\to\mathcal{E}^{*},$

denoted by $(\mathcal{V},\psi)$. Here $\phi^{*}$ is the morphism dual to
$\phi$. The restriction of the dual $F$-crystal to $\mathfrak{D}$ is ordinary.
Let

(1.8) $\displaystyle
0\to\mathcal{U}^{\prime}\hookrightarrow\mathcal{V}\to\mathcal{V}/\mathcal{U}^{\prime}\to
0$

be the corresponding unique exact sequence of bundles with connections over
$\mathfrak{D}$ with properties (i)-(iv). Here $\mathcal{U}^{\prime}$ is the
annihilator of $\mathcal{U}$.

###  1.6. $F$-crystal and $p^{s}$-hypergeometric sections

Recall the $p^{s}$-hypergeometric solutions $Q^{s,\ell},$ $\ell=1,\dots,g$,
and the bundle $\mathcal{W}_{s}$ generated by the them.

Let
$\mathcal{U}^{\prime}_{s}:=\mathcal{U}^{\prime}{\otimes}(\mathbb{Z}/p^{s})$ be
the reduction of $\mathcal{U}^{\prime}$ modulo $p^{s}$.

######  Theorem 1.3.

For $s\geq 1$, the $p^{s}$-hypergeometric solutions $Q^{s,\ell},$
$\ell=1,\dots,g$, form a flat basis of the
$\Gamma(\mathfrak{D}_{s},\mathcal{O}_{\mathfrak{D}_{s}})$-module
$\mathcal{U}^{\prime}_{s}$, that is,
$\mathcal{W}_{s}=\mathcal{U}^{\prime}_{s}$.

Theorem 1.3 is proved in Section 2.

######  Corollary 1.4.

The KZ bundle $\mathcal{U}^{\prime}$ is the $p$-adic limit of the sequence of
KZ bundles $\mathcal{W}_{s}$ as $s\to\infty$. ∎

Theorem 1.3 and Corollary 1.4 give a new proof of [VZ, Theorem 6.9] that the
sequence of bundles $\mathcal{W}_{s}$ has a $p$-adic limit as $s\to\infty$,
and the limit is a KZ-invariant subbundle of $\mathcal{V}|_{\mathfrak{D}}$.
Notice that the proof in [VZ] is based on matrix Dwork–type congruences, while
the proof of Theorem 1.3 is based on properties of ordinary $F$-crystals.

###  1.7. Pairing on $\mathcal{U}^{\prime}$

We have the Poincare pairing on $\mathcal{E}=H^{1}_{dR}(X/\mathfrak{S})$ and
the dual Poincare pairing $(\cdot,\cdot)$ on $\mathcal{V}=\mathcal{E}^{*}$.
Using properties (ii), (iii), and the fact that
$(\phi(a),\phi(b))=p\cdot(a,b)$ for all $a,b\in\mathcal{V}$, we conclude that
$(\cdot,\cdot)|_{\mathcal{U}^{\prime}}=0$ and obtain the following corollary
of Theorem 1.3.

######  Corollary 1.5.

For $s\geq 1$, the $p^{s}$-hypergeometric solutions $Q^{s,\ell},$
$\ell=1,\dots,g$, generate a Lagrangian subbundle of $\mathcal{V}_{s}$ with
respect to the Poincare pairing $(\cdot,\cdot)$. ∎

One can calculate the pairing $(\cdot,\cdot)$ explicitly, cf. [VV, Formula
(4.30)]. Then equations $(Q^{s,\ell},Q^{s,m})=0$, $\ell,m=1,\dots,g$, become a
system of explicit congruences modulo $p^{s}$ for rational functions in $z$,
cf. orthogonal relations [VV, Formula (1.14)] for $p$-hypergeometric
solutions.

###  1.8. Local flat sections

Here we discuss the KZ bundle locally. Let ${\mathbb{K}}$ be a field of
characteristic $p$, $W({\mathbb{K}})$ its ring of Witt vectors. The reader may
assume that ${\mathbb{K}}={\mathbb{F}}_{p}$ and then
$W({\mathbb{F}}_{p})={\mathbb{Z}}_{p}$.

Choose a point $a=(a_{1},\dots,a_{n})\in\mathfrak{D}(W({\mathbb{K}}))$. Denote
$t_{i}=z_{i}-a_{i}$, $t=(t_{1},\dots,t_{n})$,

$\displaystyle
R=W({\mathbb{K}})[[t_{1},\dots,t_{n}]],\qquad\mathcal{V}_{\operatorname{loc}}=\\{(v_{1},\dots,v_{n})\in
R^{\oplus n}\mid v_{1}+\dots+v_{n}=0\\}.$

Formulas (1.1) well-define a KZ connection $\nabla^{\operatorname{KZ}}$ on
$\mathcal{V}_{\operatorname{loc}}$. We study the flat sections of this
connection.

Let $H^{1}_{dR}(X/R)$ be the relative de Rham cohomology of our family of
hyperelliptic curves. This is a free $R$-module of rank $2g$ equipped with the
Gauss-Manin connection $\nabla^{\operatorname{GM}}$. We denote it by
$\mathcal{E}_{\operatorname{loc}}$.

We identify the $R$-module $\mathcal{V}_{\operatorname{loc}}$ with the
$R$-module $\mathcal{E}_{\operatorname{loc}}^{*}$ dual to
$\mathcal{E}_{\operatorname{loc}}$ and identify the KZ connection
$\nabla^{\operatorname{KZ}}$ with the connection on
$\mathcal{E}_{\operatorname{loc}}^{*}$ dual to the Gauss-Manin connection
$\nabla^{\operatorname{GM}}$. This identification allows us to identify the
flat sections of the KZ connection with the $\nabla^{\operatorname{GM}}$-flat
linear functionals $I:\mathcal{E}_{\operatorname{loc}}\to R$.

Consider the algebra endomorphism $F:R\to R$ defined by the formula
$t_{i}\mapsto t_{i}^{p}$, $i=1,\dots,n$. Let
$(\mathcal{E}_{\operatorname{loc}},\phi)$ be the associated $F$-crystal. The
$F$-crystal $(\mathcal{E}_{\operatorname{loc}},\phi)$ is ordinary since
$a\in\mathfrak{D}(W({\mathbb{K}}))$. According to [DI], there is a unique
exact sequence of $R$-modules with connections,

(1.9) $\displaystyle
0\to\mathcal{U}_{\operatorname{loc}}\hookrightarrow\mathcal{E}_{\operatorname{loc}}\to\mathcal{E}_{\operatorname{loc}}/\mathcal{U}_{\operatorname{loc}}\to
0.$

such that

  1. (i)

$\operatorname{rank}\,\mathcal{U}_{\operatorname{loc}}=g$,

  2. (ii)

$\phi\,|_{\mathcal{F}^{*}\mathcal{U}_{\operatorname{loc}}}:F^{*}\mathcal{U}_{\operatorname{loc}}\to\mathcal{U}_{\operatorname{loc}}$
is an isomorphism,

  3. (iii)

$\phi\,|_{F^{*}(\mathcal{E}_{\operatorname{loc}}/\mathcal{U}_{\operatorname{loc}})}:F^{*}(\mathcal{E}_{\operatorname{loc}}/\mathcal{U}_{\operatorname{loc}})\to\mathcal{E}_{\operatorname{loc}}/\mathcal{U}_{\operatorname{loc}}$
has the form $p\cdot\phi_{1}$, where
$\phi_{1}:F^{*}(\mathcal{E}_{\operatorname{loc}}/\mathcal{U}_{\operatorname{loc}})\to\mathcal{E}_{\operatorname{loc}}/\mathcal{U}_{\operatorname{loc}}$
is an isomorphism.

  4. (iv)

$\mathcal{U}_{\operatorname{loc}}$ equals
$\bigcap_{i}\phi^{i}((F^{i})^{*}\mathcal{E}_{\operatorname{loc}})\subset\mathcal{E}_{\operatorname{loc}}$,
where $\phi^{i}$ is viewed as
$(F^{i})^{*}\mathcal{E}_{\operatorname{loc}}\to\mathcal{E}_{\operatorname{loc}}$
,

  5. (v)

$\mathcal{U}_{\operatorname{loc}}$ has a $\nabla^{\operatorname{GM}}$-flat
basis.

Compare this list with properties (i)-(iv) in Section 1.5 and notice the new
property (v).

The second main result of this paper is the following theorem.

######  Theorem 1.6.

  1. $\operatorname{(i)}$

If $I:\mathcal{E}_{\operatorname{loc}}\to R$ is a flat section of the KZ
bundle, then $I|_{\mathcal{U}_{\operatorname{loc}}}=0$.

  2. $\operatorname{(ii)}$

The space of flat sections of the KZ connection on
$\mathcal{V}_{\operatorname{loc}}$ is a free $W({\mathbb{K}})$-module of rank
$g$.

Theorem 1.6 is deduced in Section 3 from paper [VV] where it is shown that the
space of solutions of the KZ equations (1.1) in characteristic $p$ is
$g$-dimensional.

######  Corollary of Theorems 1.3 and 1.5.

Let $I^{1},\dots,I^{g}$ be a basis of the module of flat sections of
$\mathcal{V}_{\operatorname{loc}}$ . Then there exists a unique matrix
$(b_{s,j}^{i})\in\operatorname{Mat}_{g}((W({\mathbb{K}})/p^{s})[[t]])$, such
that $\det(b_{s,j}^{i})\not\equiv 0$ modulo $p$, and

(1.10) $\displaystyle I^{i}\ \equiv\
\sum_{j=1}^{g}\,b_{s,j}^{i}\,Q^{s,j}\pmod{p^{s}},\qquad i=1,\dots,g.$

Moreover, each $b_{s,j}^{i}$ is a quasi-constant, that is $db_{s,j}^{i}=0$. ∎

######  Remark.

One can show that for any $s>1$, there exists a flat map
$\mathcal{E}_{\operatorname{loc}}\to R/p^{s}$ such that
$I|_{\mathcal{U}_{\operatorname{loc}}}\neq 0$, cf. part (i) of Theorem 1.6.
That is, for any $s>1$, not every solution to the KZ equation modulo $p^{s}$
is a linear combination of $p^{s}$-hypergeometric solutions.

It is shown in [VV], that every flat map $\mathcal{E}_{\operatorname{loc}}\to
R/p$ is a linear combination of $p$-hypergeometric solutions.

######  Remark.

It is an interesting problem to choose a natural basis $I^{1},\dots,I^{g}$ and
describe the matrix $(b^{i}_{s,j})$ explicitly. Cf. [V2, Corollary 7.4], [SlV,
Theorem 9.8], [SmV, Theorem 1.2].

##  2\. Proof of Theorem 1.3

###  2.1. Some cohomology classes

Let
$\nabla^{\operatorname{GM}}_{i}:=\nabla^{\operatorname{GM}}_{\frac{\partial}{\partial
z_{i}}}$ denote the differentiation with respect to the Gauss-Manin connection
along the vector field $\frac{\partial}{\partial z_{i}}$. Define the
differential forms

$\displaystyle\omega_{i}=\frac{dx}{(x-z_{i})y}\,,\quad i=1,\dots,n.$

The following facts are well-known.

######  Lemma 2.1.

  1. $\operatorname{(i)}$

The elements $[\omega_{i}]\in H^{1}_{dR}(X/\mathfrak{S})$, $i=1,\dots,n-1$,
form a basis of $H^{1}_{dR}(X/\mathfrak{S})$ and

(2.1) $\displaystyle[\omega_{1}]+\dots+[\omega_{n}]=0.$

  2. $\operatorname{(ii)}$

We have

(2.2) $\displaystyle\nabla^{\operatorname{GM}}_{i}\,[\omega_{j}]$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\,\frac{[\omega_{i}]-[\omega_{j}]}{z_{i}-z_{j}}\,,\quad
i\neq
j,\qquad\nabla^{\operatorname{GM}}_{i}\,[\omega_{i}]=\frac{1}{2}\,\sum_{j\neq
i}\,\frac{[\omega_{i}]-[\omega_{j}]}{z_{i}-z_{j}}\,.$

######  Corollary 2.2.

Formulas (2.2) imply the duality of the KZ connection
$\nabla^{\operatorname{KZ}}$ and the connection dual to the Gauss-Manin
connection $\nabla^{\operatorname{GM}}$. ∎

###  2.2. Map $C_{s}$ and $p^{s}$-hypergeometric solutions

For $s\geq 1$, define a generalized Cartier map
$C_{s}:\mathcal{E}\to((\mathbb{Z}/p^{s})[z])^{g}$ first on generators
$[\omega_{i}]$, $i=1,\dots,n$, and then extend it to $\mathcal{E}$ by
linearity.

We have

$\displaystyle\omega_{i}=\frac{dx}{(x-z)y}=\frac{1}{y^{p^{s}}}\,\frac{y^{p^{s}-1}dx}{x-z_{i}}=\frac{1}{y^{p^{s}}}\,{\sum}_{k}c_{k,i}(z)x^{k}dx,$

for suitable $c_{k,i}(z)\in\mathbb{Z}[z]$. Define

(2.3) $\displaystyle
C_{s}:[\omega_{i}]\mapsto(c_{p^{s}-1,i}(z),c_{2p^{s}-1,i}(z),\dots,c_{gp^{s}-1,i}(z)).$

Recall the $p^{s}$-hypergeometric solutions
$Q^{s,\ell}(z)=(Q^{s,\ell}_{1},\dots,Q^{s,\ell}_{n})^{\intercal}$,
$\ell=1,\dots,g$. Clearly

(2.4) $\displaystyle Q^{s,\ell}(z)=(c_{\ell p^{s}-1,1},c_{\ell
p^{s}-1,2},\dots,c_{\ell p^{s}-1,n}).$

######  Lemma 2.3.

Let $\omega=\sum_{i}g_{i}(z)\omega_{i}$ be such that $\omega=\frac{\partial
f}{\partial x}(x,z)dx$ for some function $f(x,z)$, then $C_{s}(\omega)=0$.

###### Proof.

By assumption we have

$\displaystyle\frac{\partial f}{\partial
x}(x,z)dx=\frac{1}{y^{p^{s}}}{\sum}_{i,k}g_{i}(z)c_{k,i}(z)x^{k}dx.$

Hence

$\displaystyle\frac{\partial}{\partial
x}\big{(}y^{p^{s}}f(x,z))dx\equiv{\sum}_{i,k}g_{i}(z)c_{k,i}(z)x^{k}dx\pmod{p^{s}}.$

Thus ${\sum}_{i}g_{i}(z)c_{\ell p^{s}-1,i}(z)\equiv 0\pmod{p^{s}}$ for
$\ell=1,\dots,g,$ and $C_{s}(\omega)=0$. ∎

###  2.3. Proof of Theorem 1.3

Let $U\subset X$ be the affine curve defined by equation (1.2). Set

$\displaystyle\mathcal{G}_{s}=\ker\big{(}H^{1}(X_{s}/\mathfrak{D}_{s})\to
H^{1}(U_{s}/\mathfrak{D}_{s})\big{)},$ [1][1][1]For any nonempty affine
$W\subset X$, we have
$\mathcal{G}_{s}=\ker\big{(}H^{1}(X_{s}/\mathfrak{D}_{s})\to
H^{1}(W_{s}/\mathfrak{D}_{s})\big{)}$ and $\mathcal{G}_{s}$ does not depend on
the choice of $W$.

. By Lemma 2.3, we have $C_{s}(\mathcal{G}_{s})=0$. Hence the generalized
Cartier map $C_{s}$ descends to a map

(2.5) $\displaystyle
C_{s}:H^{1}(X_{s}/\mathfrak{D}_{s})/\mathcal{G}_{s}\to\mathcal{O}(\mathfrak{D})^{g}.$

Since the $p^{s}$-hypergeometric solutions $Q^{s,\ell}$, $\ell=1,\dots,g$, are
linear independent, this map is onto. Hence,
$\operatorname{rank}\,\mathcal{G}_{s}\leq g$. [2][2][2]Alternatively, one
checks that $(\cdot,\cdot)|_{\mathcal{G}_{s}}=0$, where $(\cdot,\cdot)$ is the
Poincare form, and hence, $\operatorname{rank}\,\mathcal{G}_{s}\leq g$.

######  Lemma 2.4.

We have $\mathcal{U}_{s}\subset\mathcal{G}_{s}$.

###### Proof.

Every class in $H^{1}(U_{s}/\mathfrak{D}_{s})$ is represented by a global
1-form. This implies that the Frobenius on $H^{1}(U_{s}/\mathfrak{D}_{s})$ is
divisible by $p$. The lemma follows from (iv) in Section 1.5. ∎

By (i) in Section 1.5, we have $\operatorname{rank}\,\mathcal{U}_{s}=g$. Hence
$\mathcal{G}_{s}=\mathcal{U}_{s}$. Theorem 1.3 is proved.

##  3\. Proof of Theorem 1.6

The theorem is deduced below from paper [VV] where it is shown that the space
of solutions of the KZ equations (1.1) in characteristic $p$ is
$g$-dimensional.

Without loss of generality we may assume that $\mathbb{K}$ is algebraically
closed. Denote $\Lambda=W(K)$. Recall the $F$-crystal
$(\mathcal{V}_{\operatorname{loc}},\psi)$ of the $R$-module
$\mathcal{V}_{\operatorname{loc}}=\mathcal{E}_{\operatorname{loc}}^{*}$. We
have the corresponding slope exact sequence of $R$-modules with connections,

(3.1) $\displaystyle
0\to\mathcal{U}^{\prime}\hookrightarrow\mathcal{V}_{\operatorname{loc}}\to\mathcal{V}_{\operatorname{loc}}/\mathcal{U}^{\prime}\to
0.$

We need to show that every flat section of $\mathcal{V}_{\operatorname{loc}}$
lies in $\mathcal{U}^{\prime}$.

Recall from [DI] that there exists a basis $u_{1},\dots,u_{g}$,
$e_{1},\dots,e_{g}$ of $\mathcal{V}_{\operatorname{loc}}$ such that

  1. (i)

the collection $u_{1},\dots,u_{g}$ is a flat basis of $\mathcal{U}^{\prime}$,

  2. (ii)

for $i=1,\dots,g$, we have

$\displaystyle\nabla^{\operatorname{KZ}}(e_{i})=\sum_{j=1}^{g}d\log(q_{ij}){\otimes}u_{j}\,,$

where $q_{ij}\in R^{\times}$.

Let $f=\sum_{i}(a_{i}(t)u_{i}+b_{i}(t)e_{i})$ be a flat section of
$\mathcal{V}_{\operatorname{loc}}$. Our goal is to show that all $b_{i}(t)$
are equal to zero. That statement implies the theorem.

Using $\nabla^{\operatorname{KZ}}(f)=0$, we observe that $db_{i}(t)=0$ and
hence $b_{i}\in\Lambda$.

We assume that
$(b_{1},b_{2},\dots,b_{g})=p^{s}(b^{\prime}_{1},\dots,b^{\prime}_{g})$, where
the vector $(b^{\prime}_{1},\dots,b^{\prime}_{g})$ is nonzero modulo $p$.

Writing $\nabla^{\operatorname{KZ}}(f)=0$ we obtain

(3.2) $\displaystyle
da_{i}+p^{s}{\sum}_{j=1}^{g}\,b^{\prime}_{j}\,d\log(q_{ij})=0$

for all $i$. Reduce all terms of this identity modulo $p^{s+1}$. Let

(3.3) $\displaystyle
d\tilde{a}_{i}+p^{s}{\sum}_{j=1}^{g}\,\tilde{b}_{j}\,d\log(\tilde{q}_{ij})=0$

be the resulting identity. Here
$\tilde{b}_{1},\dots,\tilde{b}_{g}\in{\mathbb{K}}$. Denote
$\eta_{i}=\sum_{j=1}^{g}\tilde{b}_{j}\,d\log(\tilde{q}_{ij})$.

We use the following well-known lemma (see e.g. [I, Corollary 2.3.14 and
Equation (2.2.5)]).

######  Lemma 3.1.

Let $s$ be a positive integer and $\eta\in\Omega^{1}({\mathbb{K}}[[t]])$ a
closed 1-form. Identifying $\Omega^{1}({\mathbb{K}}[[t]])$ with
$p^{s}\Omega^{1}(W_{s+1}({\mathbb{K}})[[t]])$ we denote by $p^{s}\eta$ the
corresponding form in $p^{s}\Omega^{1}(W_{s+1}({\mathbb{K}})[[t]])$. Assume
that $p^{s}\eta$ is exact. Let
$C:\Omega^{1}_{\text{closed}}({\mathbb{K}}[[t]])\to\Omega^{1}({\mathbb{K}}[[t]])$
be the Cartier map. Then, for each positive $i$, the 1-form

$\displaystyle C^{i}(\eta):=\underbrace{C\circ\dots\circ C}_{\text{i}}(\eta)$

is well-defined, closed, and $C^{s+1}(\eta)=0$. ∎

Return to the exact form $p^{s}\eta_{i}$ . By the lemma, we have
$C^{s+1}(\eta_{i})=0$. Using the fact that
$C(d\log(\tilde{q}_{ij}))=d\log(\tilde{q}_{ij})$ we obtain that

$\displaystyle
C\left({\sum}_{j=1}^{g}\,(\tilde{b}_{j})^{1/p^{s}}\,d\log(\tilde{q}_{ij})\right)=0.$

This implies that
${\sum}_{j=1}^{g}\,(\tilde{b}_{j})^{1/p^{s}}\,d\log(\tilde{q}_{ij})$ is exact.
Therefore,
${\sum}_{j=1}^{g}\,(\tilde{b}_{j})^{1/p^{s}}\,d\log(\tilde{q}_{ij})+d\tilde{a}_{i}=0$
for some $\tilde{a}_{i}\in{\mathbb{K}}[[t]]$. Then
$\tilde{f}=\sum_{i}(\tilde{a}_{i}u_{i}+(\tilde{b}_{i})^{1/p^{s}}e_{i})$ is a
flat section of $\mathcal{V}_{\operatorname{loc}}{\otimes}{\mathbb{F}}_{p}$
with at least one nonzero $e_{i}$-component.

This contradicts to [VV] where it is shown that the space of flat sections of
$\mathcal{V}_{\operatorname{loc}}{\otimes}{\mathbb{F}}_{p}$ is of rank $g$ and
hence is generated by linear combinations of $u_{1},\dots,u_{g}$ with constant
coefficients. Theorem 1.6 is proved.

######  Remark.

Assume that ${\mathbb{K}}=\overline{{\mathbb{F}}_{p}}$. Then, for some $r$, we
have $C^{r}(\eta_{i})=\eta_{i}$. In particular, the sequence
$\eta_{i},C(\eta_{i}),C^{2}(\eta_{i}),\cdots$ is periodic. On the other hand,
by the lemma, we have $C^{s+1}(\eta_{i})=0$. We infer $\eta_{i}=0$ for all
$i$. Then $\tilde{f}=\sum_{i}\tilde{b}_{i}e_{i}$ is a flat section of
$\mathcal{V}_{\operatorname{loc}}{\otimes}{\mathbb{F}}_{p}$ .

## References

  * [1]
  * [BO] P. Berthelot, A. Ogus, Notes on Crystalline Cohomology, Series: Mathematical Notes, Princeton University Press, v. 256, 1–256, https://www.jstor.org/stable/j.ctt130hk6f
  * [CF] P. Christe, R. Flume, The four-point correlations of all primary operators of the $d=2$ conformally invariant $SU(2)$ $\sigma$-model with Wess-Zumino term, Nucl. Phys. B282, (1987) 466–494
  * [DJMM] E. Date, M. Jimbo, A. Matsuo, T. Miwa, Hypergeometric-type integrals and the $\mathfrak{sl}(2,{\mathbb{C}})$ Knizhnik-Zamolodchikov equation, International Journal of Modern Physics B, Vol. 04, No. 05, (1990) 1049–1057
  * [DI] P. Deligne, L. Illusie, Cristaux ordinaires et coordonnées canoniques, Séminaire de Géométrie Algébrique d’Orsay 1976?78; Lecture Notes in Mathematics Volume: 868 (1981), 80–127
  * [DF] V. Dotsenko, V. Fateev, Conformal algebra and multipoint correlation functions in $2D$ statistical models, Nuclear Physics B, Vol. 240, Issue 3, 15 October 1984, 312–348
  * [EFK] P. Etingof, I. Frenkel, A. Kirillov, Lectures on representation theory and Knizhnik–Zamolodchikov equations, Mathematical Surveys and Monographs, 58, AMS, Providence, RI, 1998. 198 pp. ISBN: 0-8218-0496-0
  * [I] L. Illusie, Complexe de de Rham-Witt et cohomologie cristalline, Ann. Sci. École Norm. Sup., Serie 4, vol. 12, (1979), no. 4, 501–661
  * [KZ] V. Knizhnik and A. Zamolodchikov, Current algebra and the Wess-Zumino model in two dimensions, Nucl. Phys. B247 (1984), 83–103
  * [SV1] V. Schechtman, A. Varchenko, Arrangements of hyperplanes and Lie algebra homology, Invent. Math. 106 (1991), 139–194
  * [SV2] V. Schechtman, A. Varchenko, Solutions of KZ differential equations modulo $p$, Ramanujan J. 48 (2019), no. 3, 655–683; https://doi.org/10.1007/s11139-018-0068-x
  * [SlV] A. Slinkin, A. Varchenko, Hypergeometric Integrals Modulo $p$ and Hasse–Witt Matrices, Arnold Math J., 7 (2021), 267–311, https://doi.org/10.1007/s40598-020-00168-2
  * [SmV] A. Smirnov, A. Varchenko, The $p$-adic approximations of vertex functions via $3D$-mirror symmetry, arXiv:2302.03092, 1–22
  * [V1] A. Varchenko, Special functions, KZ type equations, and Representation theory, CBMS, Regional Conference Series in Math., n. 98, AMS (2003)
  * [V2] A. Varchenko, Hyperelliptic integrals modulo $p$ and Cartier-Manin matrices, Pure and Applied Math Quarterly, 16 (2020), n. 3, 315–336
  * [V3] A. Varchenko, Notes on solutions of KZ equations modulo $p^{s}$ and $p$-adic limit $s\to\infty$, with Appendix written jointly with S. Sperber, Contemporary Mathematics, 780 (2022), 309–347, https://doi.org/10.1090/conm/780/15695
  * [VV] A. Varchenko, V. Vologodsky, Finding all solutions to the KZ equations in characteristic $p$, arXiv:2405.05159, 1–40
  * [VZ] A. Varchenko, W. Zudilin, Congruences for Hasse–Witt matrices and solutions of $p$-adic KZ equations, Pure and Applied Mathematics Quarterly, Volume 20 (2024), Number 1, Special Issue dedicated to Corrado De Concini, 565 – 597, DOI: https://dx.doi.org/10.4310/PAMQ.2024.v20.n1.a13