# Formation of Magnetically Truncated Accretion Disks in 3D Radiation-
Transport Two-Temperature GRMHD Simulations

M.T.P. Liska  Institute for Theory and Computation, Harvard University, 60
Garden Street, Cambridge, MA 02138, USA  G. Musoke  Anton Pannekoek Institute
for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam,
The Netherlands  A. Tchekhovskoy  Center for Interdisciplinary Exploration &
Research in Astrophysics (CIERA), Physics & Astronomy, Northwestern
University, Evanston, IL 60202, USA  O. Porth  Anton Pannekoek Institute for
Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The
Netherlands  A. M. Beloborodov  Physics Department and Columbia Astrophysics
Laboratory, Columbia University, 538 West 120th Street, New York, NY 10027,
USA  Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1,
D-85741, Garching, Germany

###### Abstract

Multi-wavelength observations suggest that the accretion disk in the hard and
intermediate states of X-ray binaries (XRBs) and active galactic nuclei (AGN)
transitions from a cold, thin disk at large distances into a hot, thick flow
close to the black hole. However, the formation, structure and dynamics of
such truncated disks are poorly constrained due to the complexity of the
thermodynamic, magnetic, and radiative processes involved. We present the
first radiation-transport two-temperature general relativistic
magnetohydrodynamic (GRMHD) simulations of truncated disks radiating at $\sim
35\%$ of the Eddington luminosity with and without large-scale poloidal
magnetic flux. We demonstrate that when a geometrically-thin accretion disk is
threaded by large-scale net poloidal magnetic flux, it self-consistently
transitions at small radii into a two-phase medium of cold gas clumps floating
through a hot, magnetically dominated corona. This transition occurs at a
well-defined truncation radius determined by the distance out to which the
disk is saturated with magnetic flux. The average ion and electron
temperatures in the semi-opaque corona reach, respectively, $T_{i}\gtrsim
10^{10}K$ and $T_{e}\gtrsim 5\times 10^{8}K$. The system produces radiation,
powerful collimated jets and broader winds at the total energy efficiency
exceeding $90\%$, the highest ever energy extraction efficiency from a
spinning black hole by a radiatively efficient flow in a GRMHD simulation.
This is consistent with jetted ejections observed during XRB outbursts. The
two-phase medium may naturally lead to broadened iron line emission observed
in the hard state.

##  1 Introduction

Observations spanning the radio to gamma-ray frequency range show that black
hole (BH) X-ray binaries (XRBs) cycle through different spectral states of
accretion over the course of months to years (e.g Esin1997; Remillard2006). BH
XRBs spend most of their time in a dim, quiescent state. However, from time to
time they go into outburst during which they transition into orders of
magnitude more luminous spectral states some of which are associated with
powerful collimated outflows, or jets. During such spectral state transitions,
the spectrum of the source changes from a power-law–like ‘hard’ spectrum to a
black-body–like ‘soft’ spectrum. During the outbursts, the luminosity $L$ can
reach a substantial fraction of the Eddington limit, $L_{Edd}$, at which the
radiation forces become comparable with gravity. In some extreme cases, such
as in extra-galactic ultra-luminous X-Ray sources (ULXs) and in the galactic
XRB SS433, the inferred bolometric luminosity can even exceed the Eddington
limit (e.g. Fabrika2004; Fuchs2006; Begelman2006; Middleton20201), which can
potentially be a sign of super-Eddington accretion (Sadowski2015_Jet). In
fact, it is conceivable that most black hole growth occurs during outbursts
when the disk is accreting at or above the Eddington limit (e.g. volonteri07).
In particular, the ‘luminous-hard state’ is very interesting since it is
associated with the most powerful transient jets in XRBs (e.g., Fender2004)
and Fanaroff-Riley type 2 (FRII, fr74) jets in luminous AGN (e.g.
Tchekhovskoy2015; 2020ARA&A..58..407D).

Accretion disks are generally described by either a geometrically-thin,
optically-thick disk (ss73; Novikov1973) or a geometrically-thick, optically-
thin radiatively inefficient accretion flow (RIAF) such as an advection
dominated accretion flow (ADAF, see Narayan1994). In a geometrically-thin disk
the density and optical depth are high. The balance between the radiative
cooling and turbulent dissipation rates determines the geometric thickness, or
the scaleheight, of the disk. Such geometrically-thin disk models can explain
thermal emission in the high-soft state. In a RIAF, the density is so low that
the plasma decouples into a two-temperature fluid, in which the Coulomb
collisions between the ions and electrons become too rare to equilibrate the
temperatures of the two species. This thermal decoupling, along with the
preferential heating of the ions, prevents the disk from radiating a
dynamically significant amount of dissipated energy. RIAFs can explain
emission in the quiescent and hard spectral states (e.g.
AbramowiczFragile2013, and references therein).

When a quiescent black hole accretion system goes into an outburst, it
transitions into a hybrid state that features both the thermal ‘soft’ and non-
thermal ‘hard’ components in its spectrum. This suggests the presence of both
a standard accretion disk and a RIAF (e.g. Remillard2006). In the literature,
‘truncated’ disk models form a popular paradigm describing such hybrid states
(e.g. Esin1997; Ferreira2006; Done2007; Marcel2018A), including their quasi-
periodic variability (e.g. Ingram2009). In these truncated disk models, the
outer part of the disk is geometrically-thin while the inner part is
geometrically-thick, hot, and described by an RIAF. We refer to the inner, hot
RIAF–like part of the system as the corona. The nature of the corona is
actively debated, and its origin and geometry – often represented as a hot
‘lamppost’ – are not definitively known. The corona can Comptonize its own
cyclo-synchrotron emission in addition to thermal emission from a thin
accretion disk (e.g. Wardzinsky2000; DelSanto2013). It is important to note
that the power-law emission attributed to the corona can make up a significant
fraction of the total X-ray luminosity in BHXRBs and AGN and hence in such
cases the hot coronal gas is likely to be also dynamically important (e.g.
Reynolds_Fabian1997; Fukumura_Kazanas2007; Wilkins_Fabian2012, and references
therein).

A major point of contention concerns the behavior of the truncation radius,
$r_{\rm tr}$, which marks the transition from the outer geometrically-thin,
radiatively efficient disk to the inner geometrically-thick, radiatively
inefficient corona. It is not clear if the softening of the spectrum during a
state transition is primarily caused by the truncation radius moving inwards
(Fender2004) or the corona shrinking in size (e.g. Kara2019), and what roles
Compton-upscattered X-ray emission from the base of the jet (e.g. Markoff2005)
and/or relativistic winds (Beloborodov1999) play. In addition, the presence of
broadened iron lines (e.g. Reis2010), which are typically associated with a
reservoir of cold optically-thick gas near the innermost stable circular orbit
(ISCO), in the low-hard state of several BHXRB systems is inconsistent with
the standard picture of a truncated disk, which does not predict the presence
of any cold gas within the truncation radius (e.g. Done2006; Done2007).

Simulating truncated accretion disks remains extremely challenging because
both radiation and thermal decoupling between ions and electrons is expected
to play a prominent role. Current MHD simulations only partly addressed disk
truncation by simplifying the (radiative) physics with an ad-hoc cooling
function which led to a smooth transition (Hogg2017; Hogg2018). Here, we use
the state-of-the-art radiation-transport two-temperature GRMHD simulations to
demonstrate that truncated accretion disks naturally form in the presence of
large-scale net poloidal magnetic flux. In Sections 2 and 3 we describe our
GRMHD code H-AMR and the initial conditions, before presenting our results and
concluding in Sections 4 and 5.

##  2 Numerical Setup

We use for our simulations the graphical processing unit (GPU) accelerated
GRMHD code H-AMR (Liska2018A; Liska2020; see Porth2019 for code comparison).
H-AMR is triple level parallelized using a hybrid CUDA-OpenMP-MPI framework
that solves the equations of GRMHD in the Kerr-Schild foliation. Magnetic
fields are evolved using a staggered grid approached based on Gardiner2005.
The spatial reconstruction of primitive variables from cell centers to cell
faces is done using a third order accurate PPM method (Collela1984), which
leads to second order overall spatial accuracy (cross-dimensional terms are
not taken into account). Integration in time is performed using the second
order IMEX2A scheme described in Mckinney2013. H-AMR features a flexible
adaptive mesh refinement (AMR) framework that allows for the refinement both
in space and time. Namely, our local adaptive time-stepping approach allows
different blocks to have different time steps even at the same refinement
level (Liska2020).

Figure 1: A transverse slice through density with magnetic field lines shown
in black for our model with poloidal magnetic flux (RADPOL, left) and model
with toroidal magnetic flux (RADTOR, right) right before enabling the
radiation transport. While the RADTOR simulation forms a thin accretion disk
with a constant scaleheight, the inner disk in the RADPOL simulation enters
the MAD state and accretion proceeds through non-axisymmetric RTI modes.

We have recently implemented on-the-fly radiation transport into H-AMR using
an M1 (two-moment) closure, following the implementations in Koral
(Scadowski2013; Sadowski2017) and HARMRAD (Mckinney2013). We include the
energy-averaged bound-free, free-free, and cyclo-Synchrotron absorption
($\kappa_{abs}$) and emission ($\kappa_{em}$) opacities in addition to the
electrons scattering ($\kappa_{es}$) opacity as given in Mckinney2017. We use
grey-opacities energy-averaged over a diluted blackbody spectrum
(Mckinney2017) to provide as accurate as possible opacities within a single
energy M1-framework. In addition, we assume that only electrons contribute and
set $\kappa_{ion}=0$. This is valid since in high-density regions where bound-
free opacities ($\kappa_{bf}\propto\rho$) become dominant, Coulomb collisions
typically equilibrate the plasma to a single-temperature.

Figure 2: A transverse slice through RADPOL. The right hemisphere shows a
zoom-in of the left hemisphere. The mesh-block boundaries are illustrated
using white lines. Each mesh-block has $N_{r}\times N_{\theta}\times
N_{\phi}=36\times 36\times 102$ cells.

The M1 closure works well in the optically-thick disk body and outside of the
thin disk where the radiation moves predominantly in one direction
(perpendicular to the disk surface). Since this is not always valid in the
corona we further discuss the physical implications of the M1 approximation on
our results in section 5. To approximate Compton cooling of the corona we
implemented both blackbody and photon number conserving Comptonization in
H-AMR (Sadowski2015). By evolving the photon number ($n_{rad}$), a more
accurate radiation temperature can be calculated in the optically thin corona
($T_{rad}\sim(E_{rad}/(2.7k_{b}n_{rad})$ instead of
$T_{rad}=(E_{rad}/a)^{1/4}$). This increases the accuracy of the temperature
dependent emission, absorption and Compton scattering coefficients in the
plasma, and might play an important effect when the radiation spectrum
deviates from a blackbody such as in an optically thin corona.

Figure 3: The presence of large-scale poloidal magnetic flux leads to the
development of a two-phase medium: a low-density, thick, hot corona-like
accretion flow with patches of cold gas floating through it (top row, model
RADPOL) within $r\lesssim 20r_{g}$. In contrast, in the absence of the
poloidal magnetic flux the cold, thin flow extends down to the black hole
(bottom row, model RADTOR). The three columns, from left to right, show
vertical slices through the density, ion temperature, and electron
temperature; the right hemispheres show zoom-ins on the left hemispheres.
Magnetic field lines are shown in black, jet boundary in white
(${p_{mag}}={\rho}c^{2}$), and the last scattering surface in pink
($\tau_{es}=1$, integrated in vertical direction). In RADPOL the accretion
disk transitions into a magnetic pressure supported corona within $r\lesssim
20r_{g}$, which decouples into a 2-temperature plasma of hot, radiating,
electrons and very hot, non-radiating, ions. Optically thick patches of cold
gas, visible as high density regions with $\tau_{es}\gtrsim 1$ within
$r\lesssim 25r_{g}$, cover a significant fraction of the surface area in the
corona near the midplane. In RADTOR the plasma in the disk remains strongly
coupled and optically thick. Even above the photosphere, the plasma is much
colder than in RADPOL.

To allow for the development of a two-phase medium we have implemented two-
temperature thermodynamics in H-AMR, and combined it with our M1 radiation
scheme. In H-AMR, we evolve the ion and electron entropy tracer,
$\kappa={p_{e,i}}/{\rho^{\Gamma-1}}$, where $p_{e,i}$ is the electron and ion
pressure, respectively, $\rho$ is the fluid-frame rest-mass density. In this
work, we adopt $\Gamma=5/3$ for both electrons and ions: we have verified that
this is a safe choice by carrying out simulations with variable polytropic
index for the electrons. The total dissipation is calculated as the difference
between the internal energies computed from the total energy equation and from
the entropy tracers (Ressler2015). For the purposes of partitioning the
dissipated energy between the ions and electrons, we assume here that the
dissipation ultimately occurs via magnetic reconnection in unresolved current
sheets. The dissipated energy is divided between the ions and electrons using
an analytical prescription derived from particle-in-cell simulations
(Rowan2017). We account for the energy exchange via Coulomb collisions between
ions and electrons through an implicit source term (Sadowski2017). Due to the
low radiative efficiency of ions as compared to electrons, when a two-
temperature fluid forms, we assume that radiative cooling reduces the electron
entropy and keeps the ion entropy constant.

H-AMR allows broad flexibility in the choice of a coordinate system. Here, we
adopt spherical polar coordinates, $r,\theta,\phi$, and choose a uniform grid
in the $\log r$, $\theta$, and $\phi$ variables. We use a base resolution of
$1020\times 432\times 288$ cells in the three dimensions, respectively. We
then use AMR to increase this resolution in the regions of interest, as we
describe below. We place the inner radial boundary inside of the event horizon
(with $\geq 5$ radial cells inside of the event horizon at the base level:
this ensures that the inner radial grid boundary is properly causally-
disconnected from the black hole exterior). We place the outer radial boundary
at $R_{out}=10^{4}r_{g}$. To resolve the structure of the thin equatorial
disk, we quadruple the resolution for $5r_{g}\lesssim r\lesssim 120r_{g}$
within $7.5^{\circ}$ of the midplane to an effective resolution of $4080\times
1728\times 1152$ cells by adding 2 layers of static mesh refinement. In
addition, we use 4 levels of local adaptive time-stepping, which adaptively
sets the timestep based on the local Courant condition (Liska2020). This
increases the timestep in the outer disk, which speeds up and improves the
accuracy of the simulations (Chatterjee2019). To prevent cell squeezing near
the poles ($\theta=0,\pi$), we progressively reduce the $\phi$-resolution to
$N_{\phi}=18$ within $30^{\circ}$ from each pole (we leave the resolutions in
the other two dimensions unaffected). We use outflow boundary conditions in
$r$-, transmissive boundary conditions in $\theta$\- and periodic boundary
conditions in $\phi$-directions. The grid features in total $\sim 0.45\times
10^{9}$ cells and achieves a computational speed of effectively $1.5\times
10^{7}$ zone-cycles/s/GPU (taking into account the speedup from the local
adaptive timestepping). This is an order of magnitude less compared to the
non-radiative version of H-AMR. The grid is visualized in Fig. 2.

##  3 Physical Setup

To study the effect of large-scale magnetic flux on the structure of the inner
disk, we consider two models. We initialize the model RADPOL with a purely
poloidal magnetic field and model RADTOR with a purely toroidal magnetic
field. Both models feature a rapidly spinning black hole of dimensionless spin
$a=0.9375$ and use a $\Gamma=5/3$ for ions and electrons. Since ideal GRMHD is
unable to describe the physical processes responsible for injecting gas in the
jet funnel we use drift-frame density floors (Ressler2017) to ensure that
${\rho c^{2}}\geq{p_{mag}}/12.5$ in the jets.

Model RADPOL starts with a thick equatorial torus in hydrostatic equilibrium
with the inner radius at $r_{in}=12.5r_{g}$, density maximum at
$r_{max}=25r_{g}$, and outer radius $r_{out}=200r_{g}$ (Fishbone1976). The
magnetic field has the vector potential $A_{\phi}\propto(\rho-0.05)^{2}r^{2}$.
We normalize this magnetic field such that
$\beta=\max{p_{gas}}/\max{p_{mag}}=30$ where $p_{gas}$ and $p_{mag}$ are the
gas and magnetic pressures, and ‘$\max$’ denotes the maximum taken over the
torus. This setup is intended to represent the hard intermediate state, during
which a hard-state thick torus collapses into a thin accretion disk due to
radiative cooling (e.g. Esin1997). This radiative collapse is modelled by
cooling the torus on the orbital timescale using a cooling function (see Sec.
4). This cooling process is more thoroughly described in Liska2018C, who used
the same setup but ran the simulation at a lower resolution and for a shorter
duration.

Model RADTOR starts with a thin equatorial accretion disk on Keplerian orbits
in approximate hydrostatic equilibrium whose density distribution,
$\rho(r,z)\propto r^{-1}\exp{(-{z^{2}}/{2h^{2}})}$, extends from an inner
radius, $r_{in}=6r_{g}$, to the outer radius, $r_{out}=76r_{g}$ at the
constant $h/r=0.03$. Magnetic field is described by a vector potential
$A_{\theta}\propto(\rho-0.0005)r^{2}$, which is normalized to give an
approximately uniform $\beta\sim 7$. This setup is intended to represent the
high soft state, where neither large-scale poloidal magnetic flux nor jets are
present, presumably because all net poloidal magnetic flux has diffused out
after the disk has settled into a geometrically thin disk (e.g. Begelman2014).

Figure 4:  The absorption heating rate ($\lambda_{abs}$, left column),
emission cooling rate ($\lambda_{em}$, middle column), and Compton cooling
rate ($\lambda_{sc}$, right column) for models RADPOL (top row) and RADTOR
(bottom row). The right hemispheres show zoom-ins of the left hemisphere.
Magnetic field lines are shown in black, jet boundary in white
(${p_{mag}}={\rho}c^{2}$), and the last scattering surface in pink
($\tau_{es}=1$, integrated in vertical direction). For clarity, we have
removed data from the jet, whose thermodynamics cannot be modeled in GRMHD as
explained in Sec. 4. Cooling in the corona of RADPOL ($r\lesssim 25r_{g}$) is
concentrated near the equatorial plane, where the plasma has electron
temperatures of $T_{e}\gtrsim 5\times 10^{8}K$. Here Compton cooling is
expected to give rise to hard non-thermal emission that dominates over thermal
emission. Since the plasma is optically thick (and self-absorbed) in model
RADTOR, a more detailed ray-tracing analysis will be necessary to calculate
the hardness of the spectrum.

##  4 Results

We evolve both models in two stages. During the first stage, we use neither
the radiation transport nor two-temperature thermodynamics. Instead, we use a
cooling function (Noble2009) to cool the disk to a target thermal scaleheight
of $({h}/{r})_{init}=0.03$: this cooling process takes $\sim 10^{3}r_{g}/c$
for RADPOL and mimics the catastrophic cooling of a BHXRB disk undergoing an
outburst as described in Liska2018C. We initialize RADTOR at the target
temperature profile corresponding to the same disk scaleheight,
$({h}/{r})_{init}=0.03$. We evolve RADPOL for $t\sim 188,840r_{g}/c$ and
RADTOR for $t\sim 136,075r_{g}/c$, which gives both models sufficient time to
reach an approximately constant accretion rate on the BH. This allows us to
select a well-defined accretion rate (with respect to the Eddington rate)
before we turn on radiation.

The left panel in Figure 1 shows that RADPOL accumulates large-scale poloidal
magnetic flux in the inner disk: this is similar to the model described in
Liska2018C, which is essentially the same as RADPOL, apart from a small disk
tilt, numerical grid, and duration. Model RADTOR does not generate nor advect
any significant large scale poloidal magnetic flux as was suggested in
previous work considering geometrically thick disks (Liska2018B). The reason
for this is unclear, but it might involve the inability of geometrically thin
disks to generate and/or advect magnetic flux inwards (e.g. Lubow1994). It is
also possible that the simulation was not run long enough to capture this
effect for geometrically thin disks, which have a 2 orders of magnitude longer
viscous time than the disk presented in Liska2018B. The flux accumulation in
RADPOL is accompanied by a decrease in density at $r\lesssim 20r_{g}$ and the
development of a magnetically arrested disk (MAD, Tchekhovskoy2011): the
accretion in the inner disk proceeds through magnetic Rayleigh-Taylor
instabilities (RTI), as is expected in a MAD (e.g. Tchekhovskoy2011;
Mckinney2012). Neither of these two effects are observed for RADTOR (right
panel of Fig. 1) and highlights the importance of the net poloidal magnetic
flux for driving the disk truncation.

In the second stage, we set $M_{BH}=10M_{\odot}$, renormalize the density such
that $\langle\dot{M}\rangle\sim 0.35\dot{M}_{Edd}$, and set
$T_{rad}=T_{e}=T_{i}$, where we define $\dot{M}_{Edd}=L_{Edd}/\eta$ with
$\eta=\eta_{NT}=0.178$ and $L_{Edd}={4\pi GM_{BH}}/{\kappa_{es}}$. We assume
solar abundances of hydrogen, helium and metals ($X=0.7,Y=0.28,Z=0.02$). We
then run both simulations for another $t\sim 1.35\times 10^{4}r_{g}/c$ in full
radiation-transport two-temperature GRMHD: for time interval
$t\in[188840,203360]r_{g}/c$ in RADPOL and $t\in[136076,148997]r_{g}/c$ in
RADTOR. We use blackbody Comptonization for both models, except during
$t\in[191840,193650]r_{g}/c$ in RADPOL, where we use photon-number conserving
Comptonization. In addition, we use for RADPOL during
$t\in[193650,196550]r_{g}/c$ a lower density floor of ${\rho
c^{2}}\geq{p_{mag}}/25$. Additionally, for the time interval of
$t\in[188840,192630]r_{g}/c$ in RADPOL we limit the ratio of the electron to
ion temperatures to $T_{e}/T_{i}\gtrsim 10^{-2}$ instead of the usual
$T_{e}/T_{i}\gtrsim 10^{-4}$. These tests validate that the results presented
in this article are not sensitive to the Comptonization routine employed or
the value of the density floors as can be witnessed by the animations on our
YouTube channel and in Figure 10. However, setting the temperature floor to
$T_{e}/T_{i}\gtrsim 0.01$ artificially increased the temperature of the
electrons in current sheets near the midplane within $r\lesssim 5r_{g}$ of the
black hole.

The top row in Figure 3 (see also our YouTube channel for animations of the
figure) shows that in RADPOL the inner disk decouples into a hot, optically
thin, two-temperature plasma and can thus be associated with a corona. Since
the physical processes such as pair creation (e.g. Levinson2018) responsible
for mass loading of the jet cannot be modelled in GRMHD simulations, which
artificially inject gas in the jet for the scheme to remain numerically
stable, the jet thermodynamics fall beyond the scope of this work. Electrons
in the corona reach temperatures of $T_{e}\gtrsim 5\times 10^{8}K$ while the
ions reach temperatures of $T_{i}\gtrsim 10^{10}K$ irrespective of the
Comptonization routine employed. Dissipation, most likely caused by magnetic
reconnection in current sheets near the midplane (e.g. RipperdaLiskaEtAl2021),
leads to localised heating of the plasma to temperatures of $T_{i}\sim
10^{12}K$ and $T_{e}\sim 10^{9}K$, which suggest that electrons radiate their
heat locally (see also Beloborodov2017). A significant part of this heated
plasma flows out as a wind along the magnetic field lines (Fig. 7). On the
other hand, pockets of cold gas survive and fall into the black hole
predominantly along the equatorial plane. As explained in Sec. 5 this has
interesting consequences for reflection modelling.

Figure 4 shows the cooling rate of the plasma due to emission
($\lambda_{em}=\rho\kappa_{em}a_{rad}T_{e}^{4}$) and inverse Compton
scattering
($\lambda_{sc}=\rho\kappa_{es}E_{rad}\,{4k_{b}T_{e}}/{m_{e}c^{2}}$), in
addition to the heating rate through absorption
($\lambda_{abs}=\rho\kappa_{abs}E_{rad}$). As can be seen, the cooling is
predominantly concentrated near the equatorial plane which has an optial depth
of order unity ($\tau_{es}\sim 1$). Since the Compton scattering rate exceeds
the emission rate, it is conceivable that most radiation actually gets Compton
up-scattered before leaving the corona. We estimate that the Compton-$y$
parameter ($y=\tau_{es}\,{4k_{b}T_{e}}/{m_{e}c^{2}}$), which measures the
average fractional energy gain of photons due to Comptonization (e.g.
Rybicki1986), can easily exceed $y\gtrsim 1$ in the $T_{e}\gtrsim 5\times
10^{8}K$ plasma near the midplane of the corona. Pinpointing the exact value
of the Compton y-parameter is non-trivial since absorption will play an
important role in some optically thick patches near the midplane. A more
detailed calculation that takes into account emission, absorption and
scattering in the plasma along photon geodesics will be required to further
constrain the Compton-$y$ parameter. Nevertheless, this work suggests that the
corona will up-scatter and significantly harden thermal radiation from the
disk and cyclo-Synchrotron emission from the corona (see also Dexter2021;
Scepi2021). In contrast, the bottom row in Figure 3 shows that in RADTOR the
ions and electrons are strongly coupled in the cold disk, which remains
optically thick to both scattering and absorption. Since this disk extends all
the way to the event horizon, we expect a strong thermal emission component in
the emergent spectrum. Compton up-scattering of disk photons by the winds
sandwiching the inner disk to temperatures necessary to significantly harden
the spectrum is not expected since these winds are predominantly cold
($T_{e}\lesssim 5\times 10^{7}K$). The hot pockets of gas with temperatures
$T_{e}\gtrsim 10^{8}K$ within $r\lesssim 5r_{g}$ in figure 3 are relatively
rare and seem to be more prominent at later times when the disk becomes
numerically under-resolved due to thermal collapse (Sec. 7.2).

Figure 5: The first demonstration of highly efficient energy extraction from a
rapidly spinning black hole by a radiatively efficient accretion flow in a
radiation GRMHD simulation. From top to bottom, the panels show the accretion
rate in Eddington units, normalized magnetic flux on black hole and outflow
efficiencies (total in blue, jet in black, wind in cyan, and radiative in
purple; see the legend) for the radiative stage of the simulations. Both
RADPOL and RADTOR accrete at $\dot{M}/\dot{M}_{Edd}\sim 0.35$. The red lines
demarcate the time intervals when we use either photon-conserving
Comptonization ($t\in[191840,193650]r_{g}/c$) or a higher density floor
($t\in[193650,196550]r_{g}/c$). While RADPOL floods the black hole with the
large-scale poloidal magnetic flux and reaches the MAD state, RADTOR does not
develop any net large-scale poloidal magnetic flux. In RADPOL, the total
efficiency exceeds $\eta_{tot}\sim 90\%$, for the first time demonstrating
highly efficient energy extraction from a radiatively efficient accretion flow
out of spinning black hole in a radiation GRMHD simulation. The efficiencies
of the jet and winds reach respectively $\eta_{jet}\sim 50\%$ and
$\eta_{wind}\sim 15\%$. In contrast, in the RADTOR simulation, the total
efficiency is much smaller, $\eta_{tot}\sim 40\%$, with a much smaller
$\eta_{jet}\lesssim 1\%$, and and approximately the same $\eta_{wind}\sim
20\%$. Both models achieve similar radiative outflow efficiencies,
$\eta_{rad}\sim 20\%$.

We show the accretion rate $\dot{M}=-\int_{0}^{2\pi}\int_{0}^{\pi}\rho
u^{r}\sqrt{-g}d\theta d\phi$, with the integral evaluated at $r=5r_{g}$ to
avoid the contamination by the density floors, and dimensionless black hole
poloidal magnetic flux $\phi_{\rm
BH}=0.5{\int_{0}^{2\pi}\int_{0}^{\pi}|B^{r}|\sqrt{-g}d\theta
d\phi}/{\langle\dot{M}\rangle}^{1/2}$, with the integral evaluated at the
event horizon, $r=r_{H}\simeq 1.3r_{g}$, in the upper and middle panels of
Figure 5. In the RADPOL model, the accretion disk transports so much of the
initial poloidal magnetic flux to the black hole that it becomes as strong as
gravity, obstructs the accretion, reaches a saturated value, and leads to a
MAD. The normalized magnetic flux $\phi_{BH}\sim 33$ is about two-thirds of
that in thick accretion disks (Tchekhovskoy2011) and consistent with the MAD
saturation value in thin accretion disks (Cheng, Liska et al 2022, in prep).
The lower panel of Figure 5 shows that RADPOL launches powerful jets and
winds. Namely, we show the jet ($\eta_{jet}$) and wind ($\eta_{wind}$) energy
efficiencies, which we obtain by normalizing the electromagnetic plus fluid
energy fluxes at $r=5r_{g}$ in the jets and winds by the mass accretion rate
(see Liska2018C). We define the jet-wind boundary by the
${p_{mag}}={\rho}c^{2}$ criterion. In RADPOL, we find $\eta_{wind}\sim 15\%$
and $\eta_{jet}\sim 50\%$. In RADTOR, we find $\eta_{wind}\sim 20\%$ and
$\eta_{jet}\lesssim 1\%$. In both models the associated winds are
significantly hotter than the disk. We compute the radiative efficiency from
the normalized radiation energy flux measured at $r=100r_{g}$ and find
$\eta_{rad}\sim 20\%$ in both models. This is remarkably close to the
Novikov1973 value $\eta_{NT}\sim 17.8\%$.

As Figure 3 shows, in RADPOL the excess magnetic flux remains in the inner
disk, $r\lesssim r_{\rm tr}\approx 20r_{g}$. The left column in Figure 6 shows
that this leads to the formation of a magnetic pressure supported low-density
inner disk coupled to a radiation pressure supported high-density outer disk.
The transition between these two regimes occurs at a well-defined ‘magnetic
truncation’ radius, $r_{\rm tr}$, which is set by the distance out to which
the inner disk is flooded by the poloidal magnetic flux. Similar to the
single-temperature simulation stage (Fig. 1), here the accretion in the inner
disk also proceeds through magnetic RTI and is the radiative analog of the MAD
(see also Teixeira2017). On the other hand, the disk in RADTOR remains
radiation pressure dominated and does not contain any large-scale poloidal
magnetic flux. We define the effective viscosity as
$\alpha_{eff}=-\langle{v_{r}v_{\phi}\rangle_{\rho}}/{\langle({h}/{r})_{init}^{2}v_{\phi}^{2}\rangle_{\rho}}$,
where $\langle\dots\rangle_{\rho}$ is the density-weighted angular average,
and $v_{r}$ and $v_{\phi}$ are the physical radial and azimuthal velocity
components, respectively. Figure 6(g,h) shows that $\alpha_{eff}\sim 10^{1}$
in RADPOL and $\alpha_{eff}\sim 10^{-2}$ in RADTOR. This illustrates that the
presence of large scale magnetic flux, as in RADPOL, leads to more rapid
accretion of gas (e.g. Ferreira2006). Please note that even when
$\alpha_{eff}\sim 100$ the ratio between the radial and azimuthal velocities
becomes $v_{r}/v_{\phi}=(h/r)^{2}\alpha\lesssim 0.1$, suggesting the infalling
gas has a significant toroidal velocity component capable of producing a
broadened iron line.

In both models the scaleheight of the radiation pressure supported part of the
disk shrinks to a density averaged scaleheight
${h}/{r}=\langle|\theta-\pi/2|\rangle_{\rho}\sim 0.015-0.02$ (Fig. 6 panels g
and h and YouTube channel) at which point (when the disk scaleheight or MRI
wavelength become resolved by less than $\sim 8$ cells beyond $r\gtrsim
5r_{g}$) we terminate the simulation. This runaway cooling effect suggests
that both disks are subject to viscous and/or thermal instabilities
(Lightman1974; Shakura1976). However, in RADPOL magnetic pressure stabilizes
the inner disk against thermal collapse (see also Sadowski2016_mag;
Jiang2019), which settles into a scaleheight of ${h}/{r}\sim 0.05-0.2$. Please
note that the outer disk in both RADPOL and RADTOR is radiation pressure
dominated and thus remains thermally unstable. Future simulations with higher
resolutions and longer runtimes will investigate whether magnetic pressure or
other physical processes eventually stabilize the entire accretion disk
against runaway collapse (e.g. Begelman2007). This was supported by radiative
GRMHD simulations of thin accretion disks seeded with equipartition strength
magnetic fields (e.g. Sadowski2016_mag; Lancova2019).

Figure 6: Low-density, magnetically-dominated, rapidly accreting, thicker
corona–like accretion flow develops near the black hole only in the presence
of large-scale poloidal magnetic flux. From top to bottom, the panels show
density-averaged midplane density $\bar{\rho}=\int\rho^{2}dV/\int\rho dV$
(panels a and b); the ratios of density-averaged pressures, $\bar{p}=\int\rho
pdV/\int\rho dV$ (panels c and d); effective $\alpha$-viscosity (panels e and
f; dashed lines indicate negative values); and scaleheights (panels g and h)
time-averaged over the radiative evolution stage of the simulations. The left
panels show RADPOL and right panels RADTOR results. The vertical magenta lines
show the location of the ISCO. While in both models the outer disk ($r\gtrsim
20r_{g}$) is radiation pressure dominated, the inner disk in RADPOL becomes
magnetic pressure dominated and forms a two-temperature plasma (with
$\bar{p}_{e}\ll\bar{p}_{i}$). This transition to a magnetic pressure supported
corona in RADPOL is accompanied by a sharp drop in density and an increase in
the effective scaleheight. The presence of large scale poloidal magnetic flux
leads to rapid accretion in RADPOL as evidenced by the extremely high value of
$\alpha_{eff}$.

##  5 Discussion and Conclusion

In this work we have presented the first radiation-transport two-temperature
GRMHD simulations of luminous sub-Eddington accretion disks that accrete close
to the Eddington limit ($L/L_{edd}\gtrsim 0.01$). Previous work in GRMHD only
addressed highly sub-Eddington disks with
$\lambda=\langle\dot{M}\rangle/\dot{M}_{Edd}\ll 0.01$ (e.g. Sadowski2017;
Chael2017; Chael2019; Dexter2021). More specifically, the simulations in this
article accrete at the Eddington ratio
$\lambda=\langle\dot{M}\rangle/\dot{M}_{Edd}\sim 0.35$. One of the
simulations, model RADTOR, had no net poloidal magnetic flux, while the other
one, model RADPOL, became saturated with poloidal magnetic flux and entered a
MAD state. These models were chosen to demonstrate the importance of magnetic
flux saturation on the (thermo-)dynamics of the plasma and do not attempt to
describe a possible transition between RADPOL and RADTOR and/or to rule out
other physical mechanisms that can lead to spectral hardening or disk
truncation (e.g. Meyer1994). We demonstrated that the saturation of magnetic
flux in the inner disk of model RADPOL creates a sharp transition between a
radiation pressure supported outer disk and magnetic pressure supported inner
corona. The sharp transition between the two occurs at a magnetic truncation
radius of $r_{\rm tr}\sim 20r_{g}$. Powerful jets, winds, and radiative
outflows with a combined efficiency of $\eta_{tot}\gtrsim 90\%$ emerge in the
presence of this poloidal magnetic flux. All of these characteristics are
consistent with accretion in the hard (-intermediate state) of BHXRBs.

The corona is best described by a MAD where radial magnetic pressure gradients
prevent ordered quasi-axisymmetric accretion of gas (e.g. Narayan2003;
Igumenschev2003; Igumenshchev2008; Tchekhovskoy2011; Mckinney2012;
Begelman2021). Accretion of gas in the corona proceeds through magnetic RTI
instabilities originating at the truncation radius. While a significant
fraction of the infalling gas gets heated by magnetic reconnection to
$T_{i}\gtrsim 10^{10}K$ and $T_{e}\gtrsim 5\times 10^{8}K$, we also observe
cold, optically thick, clumps of gas falling into the black hole (Fig.3). This
could potentially explain broad iron line emission observed in the hard
spectral states of XRBs and AGN (e.g. Reis2010) without the need for the
corona to recondense into a thin accretion disk near the black hole (e.g.
Meyer-Hofmeister2011).

Our model differs substantially from most truncated accretion disk models in
the literature (e.g. Esin1997; Ferreira2006; Begelman2014). In such models
heat conducted from the corona into the upper layers of the disk leads to the
disk’s gradual evaporation into a corona (e.g. Meyer1994; liu98; Qian2007).
However, in magnetically truncated disks the transition from a disk into a
corona occurs at a well-defined radius which is determined by the location of
the magnetospheric radius (e.g. Igumenshchev2009). Since accretion within the
magnetospheric radius proceeds through RTI modes, we expect that the power
spectrum of the Comptonized emission from the corona will differ substantially
from the outer disk’s blackbody emission. In future work we will address how
the cold-phase gas transitions into hot-phase gas in the corona. At the moment
of writing, we suspect that magnetic reconnection in equatorial current sheets
heats a small fraction of the cold-phase gas to extremely high temperatures,
while leaving most of the cold-phase gas unaffected. The hot-phase case
subsequently flows out along magnetic field lines (Sec. 7.1). This is a
conceivable scenario in a magnetic pressure supported plasma, but will require
further analysis (Liska et al 2022C in prep).

Spectral state transitions (e.g., Fender2004) can potentially be explained by
the outwards diffusion of magnetic flux (e.g. Lubow1994; Begelman2014). From
an observational perspective, the truncation radius needs to move inwards to
explain the softening of the spectrum during a hard-to-soft state transition
(e.g. Esin1997) while the magnetic flux needs to diffuse out to explain the
lack of any jets in the high-soft state (e.g. Fender2004). This has never been
observed in GRMHD simulations since magnetic flux diffusion is expected to
occur on much longer timescales than can be captured by state-of-the-art GRMHD
models. Nevertheless, it is conceivable that outwards diffusion of the
magnetic flux in our model will progressively reduce the magnetic truncation
radius, which will cause the spectrum to become softer. However, even when the
entire disk leaves the MAD state, residual large scale poloidal magnetic flux
might still be present and launch powerful winds from the inner disk. Such
winds could be significantly hotter compared to the winds in RADTOR and could
potentially explain the hard emission tail observed in the soft
(-intermediate) states of BHXRBs. In addition, other processes such as shocks
(Musoke2022, Liska et al 2022B in prep) and emission from jets (e.g.
Markoff2005) can lead to significant hardening of the spectrum.

Interestingly, Kinch2021 used an iterative post-processing approach
(Kinch2019) to generate spectra of thin accretion disks, which were simulated
in GRMHD using a cooling function (Noble2009; Noble2011). In these simulations
the complex radiative processes responsible for cooling the disk and corona
(defined as gas outside of the $\tau_{es}\lesssim 1$ surface) were
approximated with cooling functions (Noble2009; Kinch2020). The temperature of
the radiation emitting electrons in the corona was determined by matching the
cooling rate in the corona to the Compton emission rate of the electrons.
Making a direct comparison between Kinch2021 and RADPOL/RADTOR simulations is
difficult since both the initial configuration and numerical methods differ
substantially. The thin disk in Kinch2021 is, for example, much thicker
($h/r\sim 0.05-0.06$) and contains a significant amount of poloidal magnetic
flux, but stays below the MAD saturation value, suggesting neither RADPOL or
RADTOR can describe it. In addition, our radiative GRMHD simulations model the
radiative processes in the disk and corona such as emission, absorption, and
scattering self-consistently. We also account for thermal decoupling between
ions and electrons in the corona, and include feedback between radiation and
gas. For example, in our radiative models the local dissipation rate is not
guaranteed to match the local emission rate which can lead to runaway thermal
collapse of the accretion disk (Sec. 4).

By ray-tracing our simulation data (e.g Krawczynski_Henric2021) we will be
able place unique constraints on the relative contributions of the cold-phase
and hot-phase gas. However, generating self-consistent spectra will face
several challenges. For example, magnetic reconnection in collisionless plasma
leads to the formation of non-thermal particle populations (e.g Sironi_2014;
Beloborodov2017; Sironi2020; Navin2021). Since ideal GRMHD is unable to model
the full energy distribution of such particles, deviation from the thermal
spectrum is a free parameter that should be quantified by future particle-in-
cell (PIC) simulations. In addition, it is likely that there are some
inaccuracies in the electron temperatures presented in this work. These stem
from the intrinsic challenges associated with quantifying the dissipation rate
in a highly magnetized plasma, and from the inherent limitations related to
our two-moment (M1) radiation scheme which treats radiation as a highly
collisional fluid. This works best in regions of high optical depth where
emission and scattering of photons is localized. However, in the optically
thin corona there are multiple radiation fields which M1 cannot capture. These
include thermal radiation from the cold-phase gas and disk, Compton-scattered
radiation from the hot-phase gas in the corona, and reflected coronal
radiation from the cold-phase gas and disk. This can lead to irradiation of
the dense gas clumps in the corona by each other and the hot-phase gas, which
can potentially disperse the plasma to a less clumpy state than predicted by
M1. In other words, future work will be necessary to determine if the clumpy
structure in the corona survives when properly accounting for all radiation
fields, and, if the irradiated clumps will indeed be able to produce a
relativistic broadened iron line (which will require an elaborate calculation
to determine the ionisation state of the plasma).

##  6 Acknowledgements

We thank Michiel van der Klis, Adam Ingram, Sera Markoff, Ramesh Narayan, and
Eliot Quataert for insightful discussions. An award of computer time was
provided by the Innovative and Novel Computational Impact on Theory and
Experiment (INCITE) program under award PHY129. This research used resources
of the Oak Ridge Leadership Computing Facility, which is a DOE Office of
Science User Facility supported under Contract DE-AC05-00OR22725. This
research also used HPC and visualization resources provided by the Texas
Advanced Computing Center (TACC) at The University of Texas at Austin, which
contributed to our results via the LRAC allocation AST20011 and through
Texascale Days (http://www.tacc. utexas.edu). We acknowledge PRACE for
awarding us access to JUWELS Booster at GCS@JSC, Germany. ML was supported by
the John Harvard Distinguished Science Fellowship, GM is supported by a
Netherlands Research School for Astronomy (NOVA), Virtual Institute of
Accretion (VIA) postdoctoral fellowship and AT by the National Science
Foundation grants AST-2009884, AST-1815304 and AST-1911080.

##  7 Appendix

###  7.1 Formation of outflows

Both RADPOL and RADTOR form radiative and mechanically driven outflows. In
figure 7 we use streamlines to illustrate the direction of the radiation
velocity (red) and gas velocity (black). Both of these outflows move out
radially or vertically from the disk’s midplane. In future work we plan to
address the mass and energy ejection rates of these outflows and if these
outflows indeed become gravitational unbound. Note that radiation is treated
as a single fluid in the M1 approximation and thus only the net radiative
energy flow is captured. In reality, the radiation field would consist out of
a blackbody-component originating in the disk and Compton-scattered radiation
originating from the corona which can be reflected of the disk.

Figure 7: A transverse slice through density in RADPOL and RADTOR. The right
hemisphere is a zoom-in of the left hemisphere. Velocity streamlines are
illustrated in black and radiation streamlines are illustrated in red. In
model RADPOL both the radiation and gas flow out radially. In RADTOR the gas
velocity shows a chaotic structure, while the radiation escapes mostly
vertical.

###  7.2 Convergence

To quantify the convergence of our models we plot in Figure 9 the radial
dependence of the mass accretion rate and MRI quality factors. In model RADPOL
the mass accretion rate reaches a steady value within $r\lesssim 25r_{g}$.
This suggests RADPOL achieved inflow equilibrium in the MAD region of the disk
and corona. In RADTOR no inflow equilibrium can be achieved since the disk
collapses into an infinitely thin slab due to radiative cooling. The MRI
quality factors ($Q_{r}$, $Q_{\theta}$, $Q_{\phi}$) are defined as the number
of cells per MRI-wavelength ($\lambda_{MRI}$) weighted by $\sqrt{b^{2}}\rho$
($Q_{r,\theta.\phi}=|\lambda_{MRI}|_{\sqrt{b^{2}\rho}}/N_{r,\theta,\phi}$).
Resolving the MRI wavelength by at least 10 cells, and ideally by more than 20
cells(e.g. Shiokawa2012; Porth2019), is necessary to capture the the growth of
the longest-wavelength MRI modes in the $\theta-$ and $\phi-$ direction. Both
RADPOL and RADTOR have $Q_{\phi}\gtrsim 20$ throughout the domain. However, in
RADTOR $Q_{r}$ and $Q_{\theta}$ reach a marginal $Q_{r}\sim Q_{\theta}\sim
6-10$ for $r\lesssim 10$. Since the initial setup only contains a toroidal
magnetic field it is not unexpected that the poloidal components of the
magnetic field become challenging to resolve, especially after the inner disk
has undergone thermal collapse. We suspect that this can potentially heat up
the plasma above the photosphere in RADTOR (Fig. 8).

Figure 8: A transverse slice through density (left panels), ion temperature
(middle panels) and electrons temperature (right panels) in RADTOR at early
times (upper panels) and late times (lower panels). The right hemisphere is a
zoom-in of the left hemisphere. Magnetic field lines are visualised in black
and the last scattering surface is pink. At late times the inner disk is
subject to runaway thermal collapse, which causes the temperature of the gas
above the photosphere (cyan lines) to modestly increase. This could be caused
by insufficient resolution in the disk. An animation of this figure is
available. The top panels of the animation correspond to the large scale
portions of the figure panels while the bottom panels in the animation are the
small scale portions on the right side of each figure panel. The animation
shows the RADTOR model in the time interval [136076, 149037]. The real-time
duration of the animation is 52 seconds. Note that the full resolution
versions are available at our YouTube channel.  Figure 9: Accretion rate
$\dot{M}$ (upper panels) and MRI quality factors $Q_{r,\theta,\phi}$ (lower
panels) as function of radius averaged over the final $1000r_{g}/c$ runtime of
the respective simulation. Dashed lines show negative values. While RADPOL
achieves inflow equilibrium within $r\lesssim 25r_{g}$, RADTOR is thermally
unstable and never achieves inflow equilibrium. The MRI quality factors give
the number of cells per MRI wavelength. RADPOL is sufficiently resolved
throughout the domain while RADTOR becomes underresolved ($Q\lesssim 10$)
within $r\sim 10r_{g}$ as the disk undergoes runaway thermal collapse.

###  7.3 Effect of Comptonization and density floors

To quantify the effects of different treatments for the Comptonization and a
lower density floor, we plot in figure 10 a snapshot of each simulation, which
was restarted at $t=188839r_{g}/c$. From these figures and the animations in
the supplementary materials we conclude there is no strong dependence on
either of these two numerical choices for the temperature of the disk and/or
corona supporting the conclusions presented in this paper. However, the
temperature of the jet strongly depends on the value of the density floors.
This is an inherent limitation of GRMHD and will need to be addressed by
particle-in-cell simulations (e.g. Levinson2018).

Figure 10: A transverse slice through density (left panels), ion temperature
(middle panels) and electron temperature (right panel) for model RADPOL with
photon non-conserving comptonization (RADPOL, top panels), with photon-
conserving comptonization (RADPOL-COMP, middle panels) and with lower density
floors (RADPOL-FLOOR, lower panels). The right hemisphere is a zoom-in of the
left hemisphere. Magnetic field lines are black while the location of the jet-
disk boundary is white, and the $\tau_{es}=1$ surface is pink. An animation of
this figure is available. The top panels of the animation correspond to the
large scale portions of the figure panels while the bottom panels in the
animation are the small scale portions on the right side of each figure panel.
The animation shows the RADPOL model in the time interval [188839, 203360],
except for the time intervals [191840, 193650] and [193650, 196550] which show
the RADPOL_COMP and RADPOL_FLOOR models, respectively. The real-time duration
of the animation is 58 seconds. Note that the full resolution versions are
available at our YouTube channel. Based on the animation, we can conclude that
the temperature of the disk and corona remains relatively insensitive to these
choices in the numerics.


# Lower bounds for the directional discrepancy with respect to an interval of
rotations

Dmitriy Bilyk  School of Mathematics, University of Minnesota, Minneapolis, MN
55455<EMAIL_ADDRESS><EMAIL_ADDRESS>and  Michelle Mastrianni

###### Abstract.

We show that the lower bound for the optimal directional discrepancy with
respect to the class of rectangles in $\mathbb{R}^{2}$ rotated in a restricted
interval of directions $[-\theta,\theta]$ with $\theta<\frac{\pi}{4}$ is of
the order at least $N^{1/5}$ with a constant depending on $\theta$.

###### Key words and phrases:

discrepancy, average decay of the Fourier transform

###### 2000 Mathematics Subject Classification:

Primary 11K38, 42B05; Secondary 52A10

Both authors have been supported by the National Science Foundation (grant DMS
2054606).

##  1\. introduction

In the present paper we discuss the directional discrepancy in the plane.
Consider a set $\Omega\subset\left[-\frac{\pi}{4},\frac{\pi}{4}\right]$, which
we shall call the _allowed rotation set_. Let the class of sets
$\mathcal{R}_{\Omega}$ contain all rectangles in $\mathbb{R}^{2}$ which make
an angle of $\omega$ with the $x$-axis, where $\omega\in\Omega$:

$\mathcal{R}_{\Omega}=\\{\text{rectangles }R:\text{ one side of }R\text{ makes
angle }\omega\in\Omega\text{ with the x-axis}.\\}$ (1.1)

We then define the directional discrepancy of an $N$-point set
$P\subset[0,1]^{2}$ in the directions given by $\Omega$ as the extremal
discrepancy of $P$ with respect to elements of $\mathcal{R}_{\Omega}$. In
other words,

$D(P,\mathcal{R}_{\Omega})=\sup_{R\in\mathcal{R}_{\Omega}}\big{|}|P\cap
R|-N\cdot\textup{vol}(R\cap[0,1]^{2})\big{|},$ (1.2)

and the optimal directional discrepancy is

$D(N,\mathcal{R}_{\Omega})=\inf_{|P|=N}D(P,\mathcal{R}_{\Omega}).$ (1.3)

The two extreme cases of the directional discrepancy are well studied and
yield very different behavior of optimal asymptotic estimates.

  * •

No rotations. Axis-parallel rectangles. When $\Omega=\\{0\\}$ is a singleton
(i.e. all rectangles point in a fixed direction), one recovers the classical
case of axis-parallel rectangles, which results in very small, logarithmic
discrepancy [Sch, Le]:

$D(N,\\{0\\})\approx\log N.$ (1.4)

However more complex classes of test set yield much larger discrepancy
estimates:

  * •

All directions. Arbitrarily rotated rectangles. On the other hand, the case of
all directions, i.e. $\Omega=\left[-\frac{\pi}{4},\frac{\pi}{4}\right]$ (that
is, when $\mathcal{R}_{\Omega}$ consists of all arbitrarily rotated
rectangles) is also fairly well-understood: we have polynomial discrepancy
with bounds

$N^{1/4}\lesssim\mathcal{D}(N,\mathcal{R}_{[-\frac{\pi}{4},\frac{\pi}{4}]})\lesssim
N^{1/4}\sqrt{\log N}.$ (1.5)

Both the lower and upper bound are due to Beck [Beck87]; in §2.1 we outline
Beck’s proof of the lower bound in the case of all rotations. Perhaps
unsurprisingly, discrepancy estimates are almost identical in the case of
discs in $\mathbb{R}^{2}$ [Mo].

Our main goal in this paper is to understand what happens “in between” these
two extreme cases (a single direction and all possible directions). The
broader question is about the interplay between the geometry of the set
$\Omega$ of allowed rotations and the sharp discrepancy bounds with respect to
$\mathcal{R}_{\Omega}$ (or, even more generally, the influence of the
geometric properties of the class of test sets on discrepancy).

Various different choices of $\Omega$, as well as their geometric
characteristics, are of interest, such as Cantor-type sets, different infinite
sequences of directions, and sets with given (Hausdorff or Minkowski)
dimension. Upper bounds for some of these cases have been studied in [BMPS11,
BMPS16], where the arguments rely heavily on the Diophantine properties of the
allowed rotation set $\Omega$.

However, we are not aware of any results about lower bounds for these
intermediate cases in the literature. Already the simplest case when $\Omega$
is an _interval_ of directions is not easy to understand and produces many
open questions. In this paper we prove the first result for this case: namely,
a lower bound of the order $N^{1/5}$. Our main result, proven in §3.2, is
stated more precisely below.

######  Theorem 1.1.

Let $\Omega=[-\theta,\theta]$ with $\theta<\frac{\pi}{4}$, i.e. $\Omega$ is a
restricted interval of rotations. Then there is a constant $\gamma>0$ such
that for all $N\geq\gamma\theta^{-5}$,

$D(N,\mathcal{R}_{\Omega})\geq cN^{1/5}$

for some constant $c$ independent of $\theta$.

One can easily restate this result as a bound for all $N\in\mathbb{N}$, but
with a constant that depends on the length of the interval.

######  Corollary 1.2.

Let $\Omega=[-\theta,\theta]$ with $\theta<\frac{\pi}{4}$. There exists a
constant $c^{\prime}>0$ such that for all $N\in\mathbb{N}$,

$D(N,\mathcal{R}_{\Omega})\geq c^{\prime}\theta N^{1/5}.$

It is natural, in accordance with (1.4), that the constant decays as $\theta$
goes to zero.

Part of the proof of these results is based on and inspired by the arguments
recently developed by Brandolini and Travaglini in [BT] for discrepancy
estimates with respect to dilations and translations of a convex body with
given local convexity and smoothness properties.

It remains unclear whether or not it is possible to improve this lower bound
for directional discrepancy on restricted intervals. In particular, one might
expect that it should have order $N^{\frac{1}{4}}$, just as in the case of all
rotations, but this question is wide open.

We proceed in §2 with an exposition of the Fourier transform decay estimates,
which are used, in particular, in the proof of the $\Omega(N^{1/4})$ lower
bound (1.5) for the discrepancy in the extreme case of all rotations, see
§2.1. We adapt the Fourier bounds to our problem, see Proposition 2.1, and
explain why the proof technique for all rotations does not generalize easily
to the case of a restricted interval. In §3.1 we outline some relevant results
and ideas from [BT] and then give a proof of Theorem 1.1 in §3.2.

Throughout the paper, we shall use the notation $A\lesssim B$ meaning that
there exists an absolute constant $C$ such that $A\leq CB$ (the implicit
constant is assumed to be independent of $N$ and $\theta$ and other relevant
parameters; in particular, we shall carefully trace the dependence on
$\theta$). The relation $A\approx B$ means that $A\lesssim B$ and $B\lesssim
A$. The constants $c$, $c^{\prime}$, $c_{i}$ etc. are not necessarily the same
from line to line.

######  Remark 1.3.

The proof of Theorem 1.1 presented in §3 employs Fourier series, and
therefore, treats the unit square $[0,1]^{2}$ as the torus $\mathbb{T}^{2}$,
and from that viewpoint, rectangles in the class $\mathcal{R}_{\Omega}$ are
“periodic” rectangles. This gives rise to a slightly different notion of
discrepancy than the one defined in (1.2), where one looks at the intersection
of rectangles with the unit square rather than periodic extensions. It is easy
to see, however, that the lower bound for the discrepancy with respect to
periodic rectangles immediately yields a lower bound for the former
discrepancy (1.2). This can be seen from Figure 1. Indeed, if a periodic
rectangle has large discrepancy, then at least one of at most four
intersections with the unit square also has large discrepancy, possibly with
different constant.

Figure 1. On the left, a rotated rectangle intersected with the unit square;
on the right, the same rectangle viewed as periodic.

Finally, we remark that the lower bounds we obtain are for the $L^{2}$
discrepancy, i.e., appropriate quadratic averages of the discrepancy over
rectangles in $\mathcal{R}_{\Omega}$, which naturally implies lower bounds for
the supremum (extremal) discrepancy as defined in (1.2)–(1.3).

##  2\. Fourier transform techniques

Our proof will rely on the average decay estimates for the Fourier transform
of the test sets. This technique is often employed in various geometric
problems (including discrepancy) in which rotational invariance or curvature
is present, see e.g. [Beck87, BCIT, BT, IoLi]. In particular, the Fourier
transform method is used in Beck’s proof [Beck87] of the lower bound for the
discrepancy with respect to rectangles rotated in arbitrary directions (1.5),
which is beautifully presented in [Ma, Sec. 7.1] and which we briefly explain
in §2.1.

We now present the average Fourier decay bounds relevant to our problem.
Instead of rectangles $\mathcal{R}_{[-\theta,\theta]}$, we shall restrict our
attention just to the class of rotated squares
$\mathcal{S}_{[-\theta,\theta]}$. Let $\mathds{1}_{r,\nu}$ denote the
indicator function of the square centered at the origin with “radius” $r$ and
making angle $\nu$ with the $x$-axis, i.e. the square $[-r,r]^{2}$ rotated by
the angle $\nu$ counterclockwise. We shall be interested in the behavior of
the following quantity:

$\varphi_{R,\theta}(\xi)=\frac{1}{R{\theta}}\int_{R/2}^{R}\int_{-\theta}^{\theta}|\widehat{\mathds{1}}_{r,\nu}(\xi)|^{2}\hskip
2.84544ptd\nu\hskip
2.84544ptdr,\,\,\,\,\,\,\,\xi=(\xi_{1},\xi_{2})\in\mathbb{R}^{2},$ (2.1)

where $R>0$ and $\theta\in(0,\pi/4]$, i.e. this is the squared Fourier
transform of the indicator of the square averaged over rotations in the
interval $[-\theta,\theta]$ and dilations between $R/2$ and $R$. The relevance
of this quantity to discrepancy estimates will become apparent in, e.g., §2.1
and §3.2.

For a fixed direction $\nu=0$, it is easy to compute

$\widehat{\mathds{1}}_{r,0}(\xi)=\Big{(}\int_{-r}^{r}e^{-2\pi
ix_{1}\xi_{1}}dx_{1}\Big{)}\Big{(}\int_{-r}^{r}e^{-2\pi
ix_{2}\xi_{2}}dx_{2}\Big{)}=\frac{\sin(2\pi\xi_{1}r)}{\pi\xi_{1}}\cdot\frac{\sin(2\pi\xi_{1}r)}{\pi\xi_{1}}.$
(2.2)

Zeros of this function explain the need for averaging over dilations: it
“smears” the function, thus eliminating the zeros and allowing one to obtain
uniform lower bounds for $\varphi_{R,\theta}(\xi)$.

We have the following lower bounds for the average Fourier decay
$\varphi_{R,\theta}(\xi)$:

######  Proposition 2.1.

Let $\theta\in(0,\pi/4)$ and $R>0$.

  1. (i)

There exists a constant $c^{\prime}>0$ such that

$\varphi_{R,\theta}(\xi)\gtrsim\frac{R}{\theta|\xi|^{3}}$ (2.3)

whenever $\operatorname{arg}(\xi)\in(-\theta/2,\theta/2)$ and
$|\xi|\geq\frac{c^{\prime}}{\theta R}$.

  2. (ii)

There exists a constant $c^{\prime\prime}>0$ such that

$\varphi_{R,\theta}(\xi)\gtrsim\frac{1}{|\xi|^{4}}$ (2.4)

for all $\xi\in\mathbb{R}^{2}$ with $|\xi|\geq\frac{c^{\prime\prime}}{R}$.

Before we proceed to the proof of this proposition we would like to make a
couple remarks:

######  Remark 2.2.

In the case when $\theta=\pi/4$, i.e. if we consider squares rotated in all
directions, the decay estimate (2.3) of part (i) holds independently of
$\operatorname{arg}(\xi)$, i.e. for all $\xi$ large enough. Precisely this
estimate allows one to obtain the lower bound (1.5) for the discrepancy of
arbitrarily rotated rectangles, see the discussion in §2.1. When
$\theta<\pi/4$, however, i.e. in the case of restricted intervals, this
estimate holds only on a sector with aperture slightly smaller than the
interval of rotations, which prevents one from obtaining the same discrepancy
bound and leads to much more delicate arguments.

######  Remark 2.3.

Estimate (2.4) of part (ii) of the proposition actually holds for each square,
without the need to average over rotations. In fact, even more generally, such
an estimate holds for any planar convex set, see e.g. Theorem 24 in [BT] as
well as the discussion in §3.1. Nevertheless, we shall present a simple proof
for the square below.

###### Proof of Proposition 2.1.

We start with part (i). The beginning of the argument, as well as the
notation, closely follows Lemma 7.5 in [Ma]. As we already computed in (2.2),

$|\widehat{\mathds{1}}_{r,0}(\xi)|^{2}=\frac{\sin^{2}(2\pi\xi_{1}r)\sin^{2}(2\pi\xi_{2}r)}{{\pi^{4}}\xi_{1}^{2}\xi_{2}^{2}}.$
(2.5)

Note that
$|\widehat{\mathds{1}}_{r,\nu}(\xi)|^{2}=|\widehat{\mathds{1}}_{r,0}(\widetilde{\xi})|^{2}$
where
$\widetilde{\xi}=(\widetilde{\xi}_{1},\widetilde{\xi}_{2})=(\xi_{1}\cos\nu+\xi_{2}\sin\nu,-\xi_{1}\sin\nu+\xi_{2}\cos\nu).$
For the moment, we will assume $\xi=(|\xi|,0)$. In this case,

$\displaystyle|\widehat{\mathds{1}}_{r,\nu}(\xi)|^{2}=\left(\frac{\sin(2\pi|\xi|r\cos\nu)}{\pi|\xi|\cos\nu}\right)^{2}\,\left(\frac{\sin(2\pi|\xi|r\sin\nu)}{\pi|\xi|\sin\nu}\right)^{2}.$

Note that for $\nu\in[-\theta,\theta]$, $|\sin\nu|\approx|\nu|$ and
$\cos\nu\approx 1$. Hence, restricting to the interval $(0,\theta)$ due to
symmetry, we can write

$\varphi_{R,\theta}(\xi)\approx\frac{2}{R}\int_{R/2}^{R}\frac{1}{\theta}\int_{0}^{\theta}\frac{\sin^{2}(2\pi|\xi|r\cos\nu)\sin^{2}(2\pi|\xi|r\sin\nu)}{|\xi|^{4}\nu^{2}}\hskip
2.84544ptd\nu\hskip 2.84544ptdr$ (2.6)

Assume that $\theta>\frac{1}{R|\xi|}$, i.e. $|\xi|>\frac{1}{R\theta}$.
Consider the integral in (2.6) with the integration in $\nu$ restricted to the
interval from $\frac{1}{R|\xi|}$ to $\theta$:

$\displaystyle\begin{split}\frac{2}{R}\int_{R/2}^{R}\frac{1}{\theta}\int_{\frac{1}{R|\xi|}}^{\theta}&\frac{\sin^{2}(2\pi|\xi|r\cos\nu)\sin^{2}(2\pi|\xi|r\sin\nu)}{|\xi|^{4}\nu^{2}}d\nu
dr\lesssim\frac{2}{R}\int_{R/2}^{R}\frac{1}{\theta}\int_{\frac{1}{R|\xi|}}^{\theta}\frac{1}{|\xi|^{4}\nu^{2}}d\nu
dr\\\
&=\frac{2}{R}\int_{R/2}^{R}\frac{1}{\theta}\Big{[}-\frac{1}{\theta|\xi|^{4}}+\frac{R}{|\xi|^{3}}\Big{]}dr\lesssim\frac{R}{\theta|\xi|^{3}}.\end{split}$
(2.7)

Now let us consider the integration interval $(0,\frac{1}{R|\xi|})$. On this
interval, we have $\sin(2\pi r|\xi|\sin\nu)\approx R|\xi|\nu$. In addition,
for any $\nu$, we have

$\displaystyle\frac{2}{R}\int_{R/2}^{R}\sin^{2}(2\pi r|\xi|\cos\nu)dr\approx
1.$ (2.8)

The upper bound is clear, while the lower bound follows from the averaging
over $r$. Indeed, upon rescaling and setting $C=\pi|\xi|R\cos\nu$, one obtains
the integral

$\int_{1}^{2}\sin^{2}(Ct)\,dt=\frac{1}{2}-\frac{1}{2}\int_{1}^{2}\cos(2Ct)\,dt\geq\frac{1}{2}-\frac{1}{C}\geq\frac{1}{4},$
(2.9)

if $|\xi|R$ is large enough.

We thus have

$\displaystyle\frac{2}{R\theta}$
$\displaystyle\int_{R/2}^{R}\int_{0}^{\frac{1}{R|\xi|}}\frac{\sin^{2}(2\pi|\xi|r\cos\nu)\sin^{2}(2\pi|\xi|r\sin\nu)}{|\xi|^{4}\nu^{2}}d\nu
dr\approx\frac{2}{R\theta}\int_{R/2}^{R}\int_{0}^{\frac{1}{R|\xi|}}\frac{\sin^{2}(2\pi|\xi|r\cos\nu)R^{2}|\xi|^{2}\nu^{2}}{|\xi|^{4}\nu^{2}}d\nu
dr$
$\displaystyle=\frac{1}{\theta}\int_{0}^{\frac{1}{R|\xi|}}\frac{R^{2}|\xi|^{2}\nu^{2}}{|\xi|^{4}\nu^{2}}\Big{(}\frac{2}{R}\int_{R/2}^{R}\sin^{2}(|\xi|r\cos\nu)dr\Big{)}d\nu\approx\frac{1}{\theta}\int_{0}^{\frac{1}{R|\xi|}}\frac{R^{2}}{|\xi|^{2}}d\nu\approx\frac{R}{\theta|\xi|^{3}}.$

Hence, putting this together with (2.7), we have

$\varphi_{R,\theta}(\xi)\approx\frac{R}{\theta|\xi|^{3}}.$ (2.10)

This proves estimate (2.3) in the specific case $\xi=(|\xi|,0)$. We now show
how to remove this assumption.

Let $\xi=(\xi_{1},\xi_{2})$ with
$\theta_{0}=\operatorname{arg}(\xi)=\tan^{-1}(\frac{\xi_{2}}{\xi_{1}})=\alpha\in(0,\frac{\theta}{2})$,
and denote $\xi^{\prime}=(|\xi|,0)$ where $|\xi|$ is the length of $\xi$.
Then, obviously,
$\widehat{\mathds{1}}_{r,\nu}(\xi)=\widehat{\mathds{1}}_{r,\nu-\alpha}(\xi^{\prime})$.
Therefore,

$\displaystyle\frac{1}{\theta}\int_{0}^{\theta}|\widehat{\mathds{1}}_{r,\nu}(\xi)|^{2}d\nu=\frac{1}{\theta}\int_{0}^{\theta}|\widehat{\mathds{1}}_{r,\nu-\alpha}(\xi^{\prime})|^{2}d\nu=\frac{1}{\theta}\int_{-\alpha}^{\theta-\alpha}|\widehat{\mathds{1}}_{r,\nu}(\xi^{\prime})|^{2}d\nu\geq\frac{1}{2}\cdot\frac{2}{\theta}\int_{0}^{\theta/2}|\widehat{\mathds{1}}_{r,\nu}(\xi^{\prime})|^{2}d\nu.$

Hence,

$\varphi_{R,\theta}(\xi)\gtrsim\varphi_{R,\theta/2}(\xi^{\prime})\gtrsim\frac{R}{\theta|\xi|^{3}}$

when $|\xi|\gtrsim\frac{1}{R\theta}$, which finishes the proof of part (i) of
Proposition 2.1.

We now turn to the proof of part (ii). Observe that it is enough to prove the
bound (2.4) for one direction, e.g. $\nu=0$, without restrictions on
$\operatorname{arg}(\xi)$. Recall that $|\widehat{\mathds{1}}_{r,0}(\xi)|^{2}$
is given by (2.5). If both $|\xi_{1}|$, $|\xi_{2}|\gtrsim\frac{1}{R}$, then

$\frac{2}{R}\int_{R/2}^{R}{\sin^{2}(2\pi\xi_{1}r)\sin^{2}(2\pi\xi_{2}r)}dr\approx
1,$

which can be shown similarly to estimate (2.8) in the proof of part (i).
Therefore, in this case,

$\displaystyle\frac{2}{R}\int_{R/2}^{R}|\widehat{\mathds{1}}_{r,0}(\xi)|^{2}\,dr$
$\displaystyle=\frac{2}{R}\int_{R/2}^{R}\frac{\sin^{2}(2\pi\xi_{1}r)\sin^{2}(2\pi\xi_{2}r)}{{\pi^{4}}\xi_{1}^{2}\xi_{2}^{2}}dr\gtrsim\frac{1}{\xi_{1}^{2}\xi_{2}^{2}}\geq\frac{2}{(\xi_{1}^{2}+\xi_{2}^{2})^{2}}=\frac{2}{|\xi|^{4}}.$

On the other hand, if $|\xi_{1}|\lesssim\frac{1}{R}$, but
$|\xi|\geq|\xi_{2}|\gtrsim\frac{1}{R}$, then

$\frac{\sin^{2}(2\pi\xi_{1}r)}{\pi^{2}\xi_{1}^{2}}\approx r^{2}\approx R^{2},$

and once again invoking (2.8) we get

$\displaystyle\frac{2}{R}\int_{R/2}^{R}|\widehat{\mathds{1}}_{r,0}(\xi)|^{2}\,dr$
$\displaystyle\approx\frac{R^{2}}{\xi_{2}^{2}}\cdot\frac{2}{R}\int_{R/2}^{R}{\sin^{2}(2\pi\xi_{2}r)}dr\approx\frac{R^{2}}{\xi_{2}^{2}}\gtrsim\frac{1}{\xi_{2}^{4}}\approx\frac{1}{|\xi|^{4}}.$

∎

Since in practice we are going to use the decay estimates of Proposition 2.1
when $R\approx 1$ is fixed, we state this case as a separate corollary.

######  Corollary 2.4.

Let $\theta\in(0,\pi/4)$ and assume that $R\approx 1$ is fixed, e.g.
$R=\frac{1}{16}$. Then

$\varphi_{R,\theta}(\xi)\gtrsim\begin{cases}\frac{1}{\theta|\xi|^{3}}&\textup{
if }-\theta/2\leq\operatorname{arg}(\xi)<\theta/2\textup{ and
}|\xi|\gtrsim\frac{1}{\theta},\\\ \frac{1}{|\xi|^{4}}&\textup{ if
}|\xi|\gtrsim 1.\end{cases}$ (2.11)

###  2.1. The case of all rotations

We would like to explain briefly how Fourier estimates akin to (2.3) are used
to prove the classical lower bound for the discrepancy with respect to
$\mathcal{R}_{[-\frac{\pi}{4},\frac{\pi}{4}]}$, i.e. all rotated rectangles.
This would explain the use of the average Fourier decay and demonstrate the
difficulties of the restricted case. The argument is originally due to Beck in
[Beck87].

######  Theorem 2.5.

For the class of rectangles rotated in all directions
$\mathcal{R}_{[-\frac{\pi}{4},\frac{\pi}{4}]}$, we have

$D(N,\mathcal{R}_{[-\frac{\pi}{4},\frac{\pi}{4}]})\gtrsim N^{1/4}.$ (2.12)

###### Proof.

To simplify matters, we prove the bound for the class
$\mathcal{S}_{[-\frac{\pi}{4},\frac{\pi}{4}]}$ of rotated squares; since
$\mathcal{S}_{[-\frac{\pi}{4},\frac{\pi}{4}]}\subseteq\mathcal{R}_{[-\frac{\pi}{4},\frac{\pi}{4}]}$,
the lower bound also follows for rotated rectangles. The proof makes use of an
amplification method which exploits the fact that the $L_{2}$-discrepancy is a
convolution of two functions.

Let $P_{N}$ be any $N$-point set in $[0,1]^{2}$. If
$S(q,r,\nu)\subset[0,1]^{2}$ is an axis-parallel square centered at $q$ of
side length $2r$ and making angle $\nu$ with the $x$-axis, it is easy to see
that the local discrepancy of $P_{N}$ with respect to $S(q,r,\nu)$ can be
expressed as a convolution

$D(P_{N},S(q,r,\nu))=\int_{\mathbb{R}^{2}}\mathds{1}_{r,\nu}(q-p)d\mu(p)=(\mathds{1}_{r,\nu}\ast\mu)(q),$
(2.13)

where $\mu$ is the measure on $R^{2}$ given by
$\mu=\sum_{i=1}^{N}\delta_{p_{i}}-N\lambda$, with $\lambda$ denoting the
restriction of the usual Lebesgue measure on $\mathbb{R}^{2}$ to $[0,1]^{2}$.
Then by Plancherel’s theorem, the $L^{2}$ average of the local discrepancies
with respect to shifts $q$ satisfies

$\displaystyle\|D(P_{N},R(q,r,\nu)\|_{L^{2}(dq)}^{2}$
$\displaystyle=\int_{\mathbb{R}^{2}}D(P_{N},S(q,r,\nu))^{2}dq$ (2.14)
$\displaystyle=\int_{\mathbb{R}^{2}}\widehat{(\mathds{1}_{r}\ast\mu)}(\xi)^{2}d\xi=\int_{\mathbb{R}^{2}}|\widehat{\mathds{1}}_{r,\nu}(\xi)|^{2}|\widehat{\mu}(\xi)|^{2}d\xi,$

i.e. on the Fourier side the “point component” $\widehat{\mu}$ and the “shape
component” $\widehat{\mathds{1}}_{r,\nu}$ are completely separated.

Now, if $r=\frac{1}{2\sqrt{N}}$, then the local discrepancy
$|D(P_{N},S(q,r,\nu))|\geq\frac{1}{4}$, since “counting part” is always
integer, while the “area part” is $\frac{1}{4}$. So in this case we have the
trivial lower bound

$\int_{\mathbb{R}^{2}}D(P_{N},S(q,r,\nu))^{2}dq\gtrsim 1.$ (2.15)

The amplification idea is as follows. If
$|\widehat{\mathds{1}}_{r,\nu}(\xi)|^{2}$ would grow linearly in $r$ for any
value of $\xi$, then by changing $r=\frac{1}{2\sqrt{N}}$ to $r\approx 1$ we
could amplify the trivial lower bound for the square of $L_{2}$-discrepancy to
$\Omega(\sqrt{N})$, thus giving a discrepancy lower bound of
$\Omega(N^{1/4})$. This linear growth phenomenon does not hold for fixed $\nu$
and $r$. Indeed, expression (2.5) for $|\widehat{\mathds{1}_{r,0}}(\xi)|^{2}$
shows that zeros and concentration near the axes create a problem. But this
can be fixed by averaging with respect to dilations and translations, i.e. by
considering the quantity $\varphi_{R,\pi/4}(\xi)$.

As mentioned earlier, for all rotations (2.3) holds without any restrictions
on $\operatorname{arg}(\xi)$, and even more generally one gets

$\varphi_{R,\pi/4}(\xi)\approx\min\Big{(}R^{4},\frac{R}{|\xi|^{3}}\Big{)}\text{
for all }R,$ (2.16)

which does satisfy the linear growth alluded to earlier

$\varphi_{aR,\theta}(\xi)\geq ka\cdot\varphi_{R,\theta}(\xi).$ (2.17)

This easily yields the discrepancy lower bound (2.12) for all rotations. Thus,
letting $R_{0}=\frac{1}{2\sqrt{N}}$ and $R\in[R_{0},\frac{1}{2}),$ we have

$\displaystyle\frac{2}{R}\int_{R/2}^{R}\frac{2}{\pi}\int_{-\pi/4}^{\pi/4}\|D(P_{N},$
$\displaystyle R(q,r,\nu)\|_{L^{2}(dq)}^{2}\,d\nu
dr=\frac{2}{R}\int_{R/2}^{R}\frac{2}{\pi}\int_{-\pi/4}^{\pi/4}\int_{\mathbb{R}^{2}}|\widehat{\mathds{1}}_{r,\nu}(\xi)|^{2}|\widehat{\mu}(\xi)|^{2}d\xi
d\nu dr$
$\displaystyle=\int_{\mathbb{R}^{2}}|\widehat{\mu}(\xi)|^{2}\varphi_{R,2\pi}(\xi)d\xi\geq
k\frac{R}{R_{0}}\int_{\mathbb{R}^{2}}|\widehat{\mu}(\xi)|^{2}\varphi_{R_{0},\pi}(\xi)d\xi\gtrsim
R\sqrt{N},$

where the last inequality follows from the trivial discrepancy bound (2.15). ∎

The fact that (2.16) or, more preceisely, (2.3) holds only on a sector for
restricted intervals of directions prevents the amplification argument above
from working in this setting. In the following section we shall at least
partially overcome this difficulty.

##  3\. Restricted intervals of rotations

We turn to the case of a restricted interval of rotations. We prove Theorem
1.1 inspired by techniques from a recent paper of Brandolini and Travaglini in
[BT, 2022].

As is shown in the previous section, the favorable decay estimates for
$\varphi_{R,\theta}$ only hold on a sector of Fourier space. Thus, if we
wanted to prove an analogue of Theorem 2.5 for the restricted interval of
rotations $[-\theta,\theta]$ using the same amplification technique, we would
need to integrate the analogue of (2.14) on the sector
$S_{\theta}=\\{\xi=(\xi_{1},\xi_{2}):\tan^{-1}(\frac{\xi_{2}}{\xi_{1}})\in(-\frac{\theta}{2},\frac{\theta}{2})\\}$
rather than all of $\mathbb{R}^{2}$, i.e. consider the expression

$\displaystyle\int_{S_{\theta}}|\widehat{\mathds{1}}_{r,\nu}(\xi)|^{2}|\widehat{\mu}(\xi)|^{2}d\xi$
(3.1)

But there is a problem: in order to show that the $L_{2}$ discrepancy is of
roughly the same order as the above expression, we would need to show that the
point component $|\widehat{\mu}(\xi)|^{2}$ carries at least a constant
fraction of its mass on that sector. (Note that in the proof for all
rotations, it was not necessary to know anything about the behavior of
$|\widehat{\mu}(\xi)|^{2}$.) Since $\mu$ is the sum of $N$ dirac masses at
each point in the point set $P_{N}$, the behavior of $\widehat{\mu}$ varies
wildly depending on $P_{N}$. It is therefore difficult to conclude anything
about the relationship between the $L_{2}$ discrepancy and (3.1).

######  Remark 3.1.

The decay estimates for $\varphi_{R,\theta}$ in the proof of the amplification
lemma match those for a disc, for $\xi$ large enough. Classes of discs and
classes of all rotated rectangles tend to behave very similarly in discrepancy
settings because of their shared rotational invariance and Fourier transform
properties. In the restricted setting, the decay estimates agree with those
for a disc only on a restricted sector of Fourier space.

######  Remark 3.2.

Understanding the behavior of $|\widehat{\mu}(\xi)|^{2}$ is of interest in
other problems. For example, in [St] it is shown that one of the first $4N$
Fourier frequencies of the measure $\mu$ on the torus is not $0$: in fact,
this holds with the same constant $4$ holds on all manifolds, not just the
torus. On the torus, one can prove this directly using the following result of
Montgomery [Mo]: let $p_{1},\cdots,p_{N}$ be $N$ points in $\mathbb{T}^{2}$.
Then

$\sum_{\begin{subarray}{c}|k_{1}|\leq X_{1}\\\ |k_{2}|\leq X_{2}\\\
(k_{1},k_{2})\neq(0,0)\end{subarray}}\Big{|}\sum_{j=1}^{N}e^{2\pi im\cdot
p_{j}}\Big{|}^{2}\geq NX_{1}X_{2}-N^{2}$

for any positive numbers $X_{1},X_{2}$. The proof is quick and uses the Fejér
kernel to lower-bound the left hand side. One then wonders if other kernels
could be employed to obtain information about $\widehat{\mu}$ on sectors or
annuli. Not much progress has yet been made in this direction.

Fortunately, using techniques based [BT], we can avoid using the amplification
technique entirely and obtain a lower bound of $\Omega(N^{1/5})$ for the
directional discrepancy on any restricted interval of rotations. First let us
briefly summarize the results of [BT] which deal with a different, but
somewhat related notion of discrepancy.

###  3.1. Discrepancy with respect to rotations and dilations of a convex body

The paper [BT] studies the problem of discrepancy with respect to translations
and dilations (but not rotations!) of a given planar convex body with minimal
smoothness and curvature assumptions (or in absence thereof). The authors use
the following regularity measure:

######  Definition 3.3.

Let $C\subset\mathbb{T}^{2}$ be a convex body. For every unit vector
$\Theta=(\cos\theta,\sin\theta)$ and $\delta>0$, let

$\gamma_{\Theta}(\delta)=\Big{\\{}x\in C:x\cdot\Theta=\inf_{y\in
C}(y\cdot\Theta)+\delta\Big{\\}},$

in other words $\gamma_{\Theta}(\delta)$ is the chord in $C$ perpendicular to
the vector $\Theta$ and a distance $\delta$ from $\partial C$. Denote by
$|\gamma_{\Theta}(\delta)|$ its length. This quantity measures smoothness and
convexity of the boundary $\partial C$ in the direction $\Theta$.

It is shown in [BT] via a simple geometric argument that for any convex body
$C$ there exists a $\delta_{0}$ and $c>0$ such that for $0<\delta<\delta_{0}$
and every direction $\Theta$, we have

$\displaystyle|\gamma_{\Theta}(\delta)|\geq c\delta.$ (3.2)

This observation together with Theorem 24 in [BT] provides an alternative way
to prove and generalize part (ii) of Proposition 2.1.

We also note that if $C$ has $C^{2}$ boundary (e.g., a disc), then again there
is a $\delta_{0}$ and $c>0$ such that for $0<\delta<\delta_{0}$ and each
$\Theta$,

$\displaystyle|\gamma_{\Theta}(\delta)|\geq c\delta^{1/2}.$ (3.3)

The main result of [BT] on the discrepancy bound is stated below. Essentially,
it states that if the boundary of a convex body $C$ behaves like a disc on
some interval of directions (in the sense that $\gamma_{\theta}(\delta)$ grows
at least like $\delta^{1/2}$), and if we know the behavior of
$\gamma_{\theta}(\delta)$ elsewhere on the boundary, we can obtain a lower
bound on the discrepancy with respect to the class of all translations and
dilations of $C$.

######  Theorem 3.4 ([BT], Theorem 7).

Let $C$ be a convex body and let $\Theta=(\cos\omega,\sin\omega)$. Assume
there are constants $\delta_{0},c_{1},c_{2}>0$, $1/2\leq\sigma\leq 1$ and an
interval $\Omega$ in $(-\pi,\pi)$ such that for every $0<\delta\leq\delta_{0}$
we have

$\begin{cases}|\gamma_{-\Theta}(\delta)|+|\gamma_{\Theta}(\delta)|>c_{1}\delta^{1/2},&\omega\in\Omega\\\
|\gamma_{-\Theta}(\delta)|+|\gamma_{\Theta}(\delta)|>c_{2}\delta^{\sigma}&\omega\notin\Omega.\end{cases}$

Then there exists $c>0$ such that for every set $P_{N}$ of $N$ points in
$\mathbb{T}^{2}$ we have

$D_{2}(P_{N},\mathcal{T}(C))\geq cN^{2/(2\sigma+3)},$

where $\mathcal{T}(C)=\\{\lambda C+t:0\leq\lambda\leq
1,t\in\mathbb{T}^{2}\\}$.

Due to (3.2), Theorem 3.4 automatically gives, for a large class of convex
bodies $C$, a discrepancy lower bound of order at least $N^{1/5}$ with respect
to the class $\mathcal{T}(C)$. In this setting of discrepancy with respect to
classes of convex bodies, we can think of the “rotations” component as being
encoded in the smoothness of the boundary $\partial C$ on the interval
$\Omega$. Thus, this setting is at least heuristically similar to the class of
all translations, dilations and rotations (in the interval $\Omega$) of a
rectangle $R$.

###  3.2. Proof of Theorem 1.1

The proof of Theorem 3.4 exploits the decomposition (2.14) of the $L_{2}$
discrepancy that was used in the proof of Theorem 2.5. We also make use of an
analogue this decomposition in our proof of Theorem 1.1.

######  Remark 3.5.

As described in Remark 1.3, in the proof of Theorem 1.1, similarly to Theorem
3.4, we deal with the periodic version of discrepancy, i.e. we transfer the
problem to the torus $\mathbb{T}^{2}$ and the discrete Fourier space
$\mathbb{Z}^{2}$. Obviously, the Fourier coefficients of indicators of squares
are equal to their Fourier transforms evaluated at integer points. Hence,
decay estimates given in Corollary 2.4 continue to hold in this case.

The following fact, based on a classical lemma of Montgomery [Mo], provides
the lower bound for the “point component” in the proof .

######  Lemma 3.6 ([BT], Lemma 25).

Let $B$ be a neighborhood of the origin. Then there is a positive constant $c$
such that for every convex symmetric body $U$ in $\mathbb{R}^{2}$ and every
finite set $\\{p_{j}\\}_{j=1}^{N}\subset\mathbb{T}^{2}$ we have

$\sum_{m\in(U\setminus B)\cap\mathbb{Z}^{2}}\Big{|}\sum_{j=1}^{N}e^{2\pi
im\cdot p_{j}}\Big{|}^{2}\geq N\cdot\emph{\text{area}}(U)/4-cN^{2}.$

We can now proceed with the proof of the main theorem for discrepancy on
restricted intervals.

###### Proof of Theorem 1.1.

Let $P=\\{p_{1},\cdots,p_{N}\\}$ be an $N$-point set in $[0,1]^{2}$, and let
$\Omega=[-\theta,\theta]$ with $\theta<\pi/4$ be the interval of allowed
directions.

Consider the axis-parallel rectangle $R_{0}$ with vertices
$(\pm\frac{X}{2},\pm\frac{Y}{2})$, satisfying the conditions

  1. (1)

$XY=\kappa N$ for some constant $\kappa$,

  2. (2)

$\theta X\gtrsim Y$, and

  3. (3)

$Y\gtrsim\frac{1}{\theta}$.

Furthermore, let

$\psi=\frac{Y}{X}\hskip 14.22636pt\text{ and }\hskip
14.22636ptM=M_{\theta}=\left[\frac{\theta}{2\psi}\right]=\left[\frac{\theta}{2}\frac{X}{Y}\right].$

For every $-M\leq j\leq M$ we consider (in the Fourier space) the rotated
rectangles $R_{j}=r_{j\psi}R_{0}$, where $r_{j\psi}$ is the rotation by angle
$j\psi$ about the origin. The union $\bigcup_{j=-M}^{M}R_{j}$ then
approximates the sector in the plane between angles $-\theta/2$ and
$\theta/2$, as depicted in Figure 2 (recall that this is the sector on which
$\varphi_{R,\theta}(\xi)$ satisfies a better Fourier decay estimate, see
Corollary 2.4).

Figure 2. The approximating rectangles in Fourier space when $M=2$.

For every $m=(m_{1},m_{2})\in\mathbb{Z}^{2}$, let

$\Phi(m)=\sum_{j=-M}^{M}\mathds{1}_{R_{j}}(m),$

so that $\Phi(m)$ keeps track of how many rectangles the lattice point $m$
lies in. A simple counting argument shows that if $|m|\geq Y$ and
$m\in\bigcup_{j=-M}^{M}R_{j}$,

$\displaystyle\Phi(m)\lesssim\frac{Y}{|m|\psi}\approx\frac{X}{|m|},$ (3.4)

and if $|m|\leq Y$, we can crudely bound $\Phi(m)$ by the total number of
rectangles $2M+1$, so

$\displaystyle\Phi(m)\lesssim\frac{\theta X}{Y},$ (3.5)

where we have used the fact that that $2M+1\approx\theta\frac{X}{Y}$ since
$M\approx\theta\frac{X}{Y}\gtrsim 1$ by the assumption that $\theta X\gtrsim
Y$.

Our goal now is to find a constant $\rho$ (depending on $N$ and $\theta$) so
that

$\rho\Phi(m)\leq\varphi_{R,\theta}(m)$ (3.6)

with some fixed $R\approx 1$, say, $R=\frac{1}{16}$ for $|m|\gtrsim 1$.
According to Corollary 2.4, it would suffice to have

$\displaystyle\rho\Phi(m)\lesssim\begin{cases}\theta^{-1}|m|^{-3}&\textup{ if
}-\theta/2\leq\text{arg}(m)<\theta/2\textup{ and }|m|\gtrsim\theta^{-1}\\\
|m|^{-4}&\textup{ otherwise.}\end{cases}$ (3.7)

When $m\in\cup_{j}R_{j}$ satisfies $|m|\geq Y\gtrsim\frac{1}{\theta}$, it
would suffice to take

$\displaystyle\rho\lesssim\frac{1}{\theta X^{3}}.$ (3.8)

Indeed, in this case, since $m\in\cup_{j}R_{j}$ implies that $|M|\leq X$, we
see that

$\rho\cdot
c\frac{X}{|m|}\lesssim\theta^{-1}X^{-2}|m|^{-1}\lesssim\theta^{-1}|m|^{-3},$

which matches (3.4).

When $|m|\leq Y$, one could take $\rho$ to satisfy

$\displaystyle\rho\lesssim\frac{1}{\theta XY^{3}}.$ (3.9)

In this case one obtains

$\rho\cdot c\frac{\theta X}{Y}\lesssim Y^{-4}\lesssim|m|^{-4}.$

in agreement with (3.5)

Combining (3.8) and (3.9), we choose

$\displaystyle\rho=c\min\Big{(}\frac{1}{\theta X^{3}},\frac{1}{\theta
XY^{3}}\Big{)},$ (3.10)

again for $c$ independent of $\theta$.

We are now ready to choose the parameters which would satisfy conditions
(1)–(3). Since $XY=\kappa N$ we obtain $X\approx\kappa^{3/5}N^{3/5}$ and
$Y=\kappa^{2/5}N^{2/5}$ by setting equal the arguments of the min function in
(3.10). Therefore, $\rho\approx\theta^{-1}\kappa^{-9/5}N^{-9/5}$.

Condition (2) that $\theta\frac{X}{Y}\gtrsim 1$ then requires that $\theta
N^{1/5}\gtrsim 1$, so the remainder of the proof holds under the assumption
that $N\gtrsim\theta^{-5}$. Note that then condition (3) holds since $Y\approx
N^{2/5}\gtrsim\frac{1}{\theta^{2}}\gtrsim\frac{1}{\theta}.$

By the above discussion, for any lattice point $m$ with $|m|$ large enough
(greater than an absolute constant), we have

$\varphi_{R,\theta}(m)=\frac{2}{R}\int_{R/2}^{R}\frac{1}{|\Omega|}\int_{\Omega}\int_{\mathbb{R}^{2}}|\widehat{\mathds{1}}_{r,\nu}(m)|^{2}d\nu
dr\geq\rho\Phi(m)=\rho\sum_{j=-M}^{M}\mathds{1}_{R_{j}}(m).$ (3.11)

We proceed in a fashion similar to the proof of Theorem 2.5. Defining the
discrepancy measure $\mu=\sum_{i=1}^{N}\delta_{p_{i}}-N\lambda$ where
$\lambda$ is the Lebesgue measure on the torus $\mathbb{T}^{2}$, considering
an analog of (2.14) for the $L^{2}$ discrepancy with respect to translations,
and averaging over dilations and rotations, one obtains

$\displaystyle\begin{split}\frac{2}{R}\int_{R/2}^{R}\frac{1}{\Omega}\int_{\Omega}\|D(P_{N},R(q,r,\nu)\|_{L^{2}(dq)}^{2}d\nu
dr&=\frac{2}{R}\int_{R/2}^{R}\frac{1}{|\Omega|}\int_{\Omega}\sum_{m\neq
0}|\widehat{\mathds{1}}_{r,\nu}(m)|^{2}|\widehat{\mu}(m)|^{2}d\nu dr\\\
&=\sum_{m\neq 0}\Big{|}\sum_{j=1}^{N}e^{2\pi
imp_{j}}\Big{|}\frac{2}{R}\int_{R/2}^{R}\frac{1}{|\Omega|}\int_{\Omega}|\widehat{\mathds{1}}_{r,\nu}(m)|^{2}d\nu
dr\\\ &=\sum_{m\neq 0}\Big{|}\sum_{j=1}^{N}e^{2\pi
imp_{j}}\Big{|}\cdot\varphi_{R,\theta}(m)\\\ &\geq\sum_{|m|\geq
K}\Big{|}\sum_{j=1}^{N}e^{2\pi
imp_{j}}\Big{|}^{2}\rho\sum_{j=-M}^{M}\mathds{1}_{R_{j}}(m)\hskip
14.22636pt\text{ (by \ref{eqn:sumofrectslower})}\\\
&=\rho\sum_{j=-M}^{M}\sum_{\begin{subarray}{c}|m|\geq K,\\\ m\in
R_{j}\end{subarray}}\Big{|}\sum_{j=1}^{N}e^{2\pi imp_{j}}\Big{|}^{2}.\\\
\end{split}$ (3.12)

In the computation above, $K$ is an absolute constant large enough so that the
decay estimates in Corollary 2.4 apply. Using Lemma 3.6 and taking $\kappa$
large enough (depending only on $c$ from Lemma 3.6, but independent of $N$ and
$\theta$), the expression above can be bounded below by

$\displaystyle\rho\sum_{j=-M_{\theta}}^{M_{\theta}}(N\cdot\text{area}(R_{j})/4-cN^{2})$
$\displaystyle\geq\rho(2M_{\theta}+1)\left(\frac{1}{4}\kappa
N^{2}-cN^{2}\right)$ $\displaystyle\approx\theta^{-1}N^{-9/5}\cdot\frac{\theta
X}{Y}\cdot N^{2}\approx
N^{-\frac{9}{5}+\frac{3}{5}-\frac{2}{5}+2}=N^{\frac{2}{5}}.$

Taking square roots, we conclude that for $N\gtrsim\theta^{-5}$,

$D(N,\mathcal{R}_{\Omega})\gtrsim cN^{1/5},$

where $c$ is independent of $\theta$.

∎

###### Proof of Corollary 1.2.

Theorem 1.1 automatically implies $D(N,\mathcal{R}_{\Omega})\geq c\theta
N^{1/5}$ for $N\geq\gamma\theta^{-5}$. As in (2.15) in the proof of Theorem
2.5, it is easy to see that for small rectangles of area $\frac{1}{4N}$, the
discrepancy is at least $\frac{1}{4}$. Thus, when $N\leq\gamma\theta^{-5}$, we
trivially have $D(N,\mathcal{R}_{\theta})\geq\frac{1}{4}\geq
c\gamma^{-1/5}\theta N^{1/5}.$ ∎

It is still unclear whether this is the best lower bound for the discrepancy
with respect to restricted intervals. Indeed, in the setting of translations
and dilations of a convex set, [BT] prove that the bound is sharp by finding
an explicit point set that exhibits an upper bound of order $N^{1/5}$;
however, the proof techniques do not seem to translate in this setting. It is
natural to expect that one should be able to improve the bound to order
$N^{1/4}$, as in the case of all rotations. Diophantine properties play a
large role in discrepancy (see [BMPS11, BMPS16]), and since in any interval
one can find a number of any Diophantine type, it seems on a heuristic level
that the case of all rotations and the case of a restricted interval should
not differ much with respect to the discrepancy. There remain many open
problems in the area of directional discrepancy, and any further understanding
of the relationship between the allowed rotation set $\Omega$ and the
discrepancy of $\mathcal{R}_{\Omega}$ would be very interesting.

## References

  * [Beck87] J. Beck. _Irregularities of distribution I._ Acta Math., 159, 1–49 (1987).
  * [BMPS11] D. Bilyk, X. Ma, J. Pipher and C. Spencer. _Directional discrepancy in two dimensions_ , Bulletin of the London Mathematical Society, 43 (6), 1151–1166 (2011).
  * [BMPS16] D. Bilyk, X. Ma, J. Pipher, and C. Spencer. _Diophantine Approximations and Directional Discrepancy of Lattices_ , Transactions of the AMS, 68 (6), 3871–3897 (2016).
  * [BCIT] L. Brandolini, L. Colzani, A. Iosevich, G. Travaglini. _The rate of convergence of Fourier expansions in the plane: a geometric viewpoint._ , Math. Z., 242 (4), 709–724 (2002).
  * [BT] L. Brandolini and G. Travaglini. _Irregularities of Distribution and Geometry of Planar Convex Sets_ , Adv. Math. 396, Paper No. 108162, 40 pp., (2022).
  * [Ch] B. Chazelle, _The Discrepancy Method: Randomness and Complexity_ , Cambridge University Press (2000).
  * [IoLi] A. Iosevich, E. Liflyand, _Decay of the Fourier transform. Analytic and geometric aspects._ Birkhäuser/Springer, Basel (2014).
  * [Le] M. Lerch, _Question 1547_ , L’Intermediaire Math. 11, 144–145 (1904).
  * [Ma] Jiri Matousek, _Geometric Discrepancy: An Illustrated Guide_ , Algorithms and Combinatorics, Vol. 18 (2010).
  * [Mo] H. Montgomery, _Ten Lectures on the interface between Analytic Number Theory and Harmonic Analysis_ , American Mathematical Society (1994).
  * [Sch] W. M. Schmidt, _Irregularities of distribution. VII_ , Acta Arith., 21 45– 50 (1972).
  * [St] S. Steinerberger. _Spectral limitations of quadrature rules and generalized spherical designs_ , International Mathematics Research Notices, 2021 (16), 12265–12280 (2021).