KamLAND-Zen Collaboration

# Search for Majorana Neutrinos with the Complete KamLAND-Zen Dataset

S. Abe  Present address: Kamioka Observatory, Institute for Cosmic-Ray
Research, The University of Tokyo, Hida, Gifu 506-1205, Japan  Research Center
for Neutrino Science, Tohoku University, Sendai 980-8578, Japan  T. Araki
Research Center for Neutrino Science, Tohoku University, Sendai 980-8578,
Japan  K. Chiba  Research Center for Neutrino Science, Tohoku University,
Sendai 980-8578, Japan  T. Eda  Research Center for Neutrino Science, Tohoku
University, Sendai 980-8578, Japan  M. Eizuka  Research Center for Neutrino
Science, Tohoku University, Sendai 980-8578, Japan  Y. Funahashi  Research
Center for Neutrino Science, Tohoku University, Sendai 980-8578, Japan  A.
Furuto  Research Center for Neutrino Science, Tohoku University, Sendai
980-8578, Japan  A. Gando  Research Center for Neutrino Science, Tohoku
University, Sendai 980-8578, Japan  Y. Gando  Research Center for Neutrino
Science, Tohoku University, Sendai 980-8578, Japan  Department of Human
Science, Obihiro University of Agriculture and Veterinary Medicine, Obihiro,
Hokkaido 080-8555, Japan  S. Goto  Research Center for Neutrino Science,
Tohoku University, Sendai 980-8578, Japan  T. Hachiya  Research Center for
Neutrino Science, Tohoku University, Sendai 980-8578, Japan  K. Hata  Research
Center for Neutrino Science, Tohoku University, Sendai 980-8578, Japan  K.
Ichimura  Research Center for Neutrino Science, Tohoku University, Sendai
980-8578, Japan  S. Ieki  Research Center for Neutrino Science, Tohoku
University, Sendai 980-8578, Japan  H. Ikeda  Research Center for Neutrino
Science, Tohoku University, Sendai 980-8578, Japan  K. Inoue  Research Center
for Neutrino Science, Tohoku University, Sendai 980-8578, Japan  K.
Ishidoshiro  Research Center for Neutrino Science, Tohoku University, Sendai
980-8578, Japan  Y. Kamei  Present address: Center for Advanced Photonics,
RIKEN, Wako, Saitama, 351-0198, Japan  Research Center for Neutrino Science,
Tohoku University, Sendai 980-8578, Japan  N. Kawada  Research Center for
Neutrino Science, Tohoku University, Sendai 980-8578, Japan  Y. Kishimoto
Research Center for Neutrino Science, Tohoku University, Sendai 980-8578,
Japan  M. Koga  Research Center for Neutrino Science, Tohoku University,
Sendai 980-8578, Japan  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo Institutes for Advanced Study, The
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  A. Marthe  Research
Center for Neutrino Science, Tohoku University, Sendai 980-8578, Japan  Y.
Matsumoto  Research Center for Neutrino Science, Tohoku University, Sendai
980-8578, Japan  T. Mitsui  Present address: Faculty of Health Sciences,
Butsuryo College of Osaka, Sakai, Osaka 593-8328, Japan  Research Center for
Neutrino Science, Tohoku University, Sendai 980-8578, Japan  H. Miyake
Research Center for Neutrino Science, Tohoku University, Sendai 980-8578,
Japan  Graduate Program on Physics for the Universe, Tohoku University, Sendai
980-8578, Japan  D. Morita  Research Center for Neutrino Science, Tohoku
University, Sendai 980-8578, Japan  R. Nakajima  Research Center for Neutrino
Science, Tohoku University, Sendai 980-8578, Japan  K. Nakamura  Present
address: Faculty of Health Sciences, Butsuryo College of Osaka, Sakai, Osaka
593-8328, Japan  Research Center for Neutrino Science, Tohoku University,
Sendai 980-8578, Japan  R. Nakamura  Research Center for Neutrino Science,
Tohoku University, Sendai 980-8578, Japan  R. Nakamura  Research Center for
Neutrino Science, Tohoku University, Sendai 980-8578, Japan  J. Nakane
Research Center for Neutrino Science, Tohoku University, Sendai 980-8578,
Japan  T. Ono  Research Center for Neutrino Science, Tohoku University, Sendai
980-8578, Japan  H. Ozaki  Research Center for Neutrino Science, Tohoku
University, Sendai 980-8578, Japan  K. Saito  Research Center for Neutrino
Science, Tohoku University, Sendai 980-8578, Japan  T. Sakai  Research Center
for Neutrino Science, Tohoku University, Sendai 980-8578, Japan  I. Shimizu
Research Center for Neutrino Science, Tohoku University, Sendai 980-8578,
Japan  J. Shirai  Research Center for Neutrino Science, Tohoku University,
Sendai 980-8578, Japan  K. Shiraishi  Research Center for Neutrino Science,
Tohoku University, Sendai 980-8578, Japan  A. Suzuki  Research Center for
Neutrino Science, Tohoku University, Sendai 980-8578, Japan  K. Tachibana
Research Center for Neutrino Science, Tohoku University, Sendai 980-8578,
Japan  K. Tamae  Research Center for Neutrino Science, Tohoku University,
Sendai 980-8578, Japan  H. Watanabe  Research Center for Neutrino Science,
Tohoku University, Sendai 980-8578, Japan  K. Watanabe  Research Center for
Neutrino Science, Tohoku University, Sendai 980-8578, Japan  S. Yoshida
Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan
S. Umehara  Research Center for Nuclear Physics, Osaka University, Ibaraki,
Osaka 567-0047, Japan  K. Fushimi  Department of Physics, Tokushima
University, Tokushima 770-8506, Japan  K. Kotera  Graduate School of
Integrated Arts and Sciences, Tokushima University, Tokushima 770-8502, Japan
Y. Urano  Present address: Institute for Materials Research, Tohoku
University, 2-1-1 Katahira, Aoba-ku, Sendai, Miyagi, 980-8577, Japan  Graduate
School of Integrated Arts and Sciences, Tokushima University, Tokushima
770-8502, Japan  B.E. Berger  Nuclear Science Division, Lawrence Berkeley
National Laboratory, Berkeley, California 94720, USA  B.K. Fujikawa  Kavli
Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study, The University of Tokyo,
Kashiwa, Chiba 277-8583, Japan  Nuclear Science Division, Lawrence Berkeley
National Laboratory, Berkeley, California 94720, USA  J.G. Learned  Department
of Physics and Astronomy, University of Hawaii at Manoa, Honolulu, Hawaii
96822, USA  J. Maricic  Department of Physics and Astronomy, University of
Hawaii at Manoa, Honolulu, Hawaii 96822, USA  Z. Fu  Massachusetts Institute
of Technology, Cambridge, Massachusetts 02139, USA  S. Ghosh  Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139, USA  J. Smolsky
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
L.A. Winslow  Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA  Y. Efremenko  Kavli Institute for the Physics and Mathematics of
the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Department of Physics and
Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA  H.J.
Karwowski  Triangle Universities Nuclear Laboratory, Durham, North Carolina
27708, USA;
Physics Departments at Duke University, Durham, North Carolina 27708, USA;
North Carolina Central University, Durham, North Carolina 27707, USA;
and The University of North Carolina at Chapel Hill, Chapel Hill, North
Carolina 27599, USA  D.M. Markoff  Triangle Universities Nuclear Laboratory,
Durham, North Carolina 27708, USA;
Physics Departments at Duke University, Durham, North Carolina 27708, USA;
North Carolina Central University, Durham, North Carolina 27707, USA;
and The University of North Carolina at Chapel Hill, Chapel Hill, North
Carolina 27599, USA  W. Tornow  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Triangle Universities Nuclear Laboratory, Durham, North Carolina 27708, USA;
Physics Departments at Duke University, Durham, North Carolina 27708, USA;
North Carolina Central University, Durham, North Carolina 27707, USA;
and The University of North Carolina at Chapel Hill, Chapel Hill, North
Carolina 27599, USA  S. Dell’Oro  Center for Neutrino Physics, Virginia
Polytechnic Institute and State University, Blacksburg, Virginia 24061, USA
T. O’Donnell  Center for Neutrino Physics, Virginia Polytechnic Institute and
State University, Blacksburg, Virginia 24061, USA  J.A. Detwiler  Kavli
Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study, The University of Tokyo,
Kashiwa, Chiba 277-8583, Japan  Center for Experimental Nuclear Physics and
Astrophysics, University of Washington, Seattle, Washington 98195, USA  S.
Enomoto  Kavli Institute for the Physics and Mathematics of the Universe
(WPI), The University of Tokyo Institutes for Advanced Study, The University
of Tokyo, Kashiwa, Chiba 277-8583, Japan  Center for Experimental Nuclear
Physics and Astrophysics, University of Washington, Seattle, Washington 98195,
USA  M.P. Decowski  Nikhef and the University of Amsterdam, Science Park,
Amsterdam, the Netherlands  K.M. Weerman  Nikhef and the University of
Amsterdam, Science Park, Amsterdam, the Netherlands  C. Grant  Boston
University, Boston, Massachusetts 02215, USA  Ö. Penek  Boston University,
Boston, Massachusetts 02215, USA  H. Song  Boston University, Boston,
Massachusetts 02215, USA  A. Li  Halıcıoğlu Data Science Institute, Department
of Physics, University of California San Diego, La Jolla, California, 92093,
USA  S.N. Axani  Department of Physics and Astronomy, University of Delaware,
Newark, Delaware 19716, USA  M. Garcia  Department of Physics and Astronomy,
University of Delaware, Newark, Delaware 19716, USA  M. Sarfraz  Department of
Physics and Astronomy, University of Delaware, Newark, Delaware 19716, USA

###### Abstract

We present a search for neutrinoless double-beta ($0\nu\beta\beta$) decay of
136Xe using the full KamLAND-Zen 800 dataset with 745 kg of enriched xenon,
corresponding to an exposure of $2.097$ ton yr of 136Xe. This updated search
benefits from a more than twofold increase in exposure, recovery of photo-
sensor gain, and reduced background from muon-induced spallation of xenon.
Combining with the search in the previous KamLAND-Zen phase, we obtain a lower
limit for the $0\nu\beta\beta$ decay half-life of $T_{1/2}^{0\nu}>3.8\times
10^{26}$ yr at 90% C.L., a factor of 1.7 improvement over the previous limit.
The corresponding upper limits on the effective Majorana neutrino mass are in
the range 28–122 meV using phenomenological nuclear matrix element
calculations.

The search for neutrinoless double-beta ($0\nu\beta\beta$) decay is an active
field of research in nuclear and particle physics Agostini2023 . Observation
of this decay would prove the Majorana nature of neutrinos and demonstrate the
non-conservation of lepton number, a key feature of Leptogenesis, which can
explain the imbalance of the matter over antimatter in the Universe
Fukugita1986 . In the model of light Majorana neutrino exchange, the
$0\nu\beta\beta$ decay rate is proportional to the square of the effective
Majorana neutrino mass
$\left<m_{\beta\beta}\right>\equiv\left|\Sigma_{i}U_{ei}^{2}m_{\nu_{i}}\right|$.
The most stringent limit to date has been provided by KamLAND-Zen, reaching
the inverted mass ordering (IO) region below 50 meV for the first time
Abe2023a . An improved search covering the unexplored IO region would
stringently test theoretical models predicting $\left<m_{\beta\beta}\right>$
in this range Harigaya2012 ; Asaka2020 ; Asai2020 , with a significantly
increased chance to observe $0\nu\beta\beta$ decay Agostini2017 .

The KamLAND-Zen (KamLAND Zero-Neutrino Double-Beta Decay) Gando2012a ;
Gando2012b ; Gando2013a ; Asakura2016 ; Gando2016 ; Gando2019 ; Abe2023a
experiment utilizes the low-background KamLAND detector to search for
$0\nu\beta\beta$ decay. Enriched xenon is dissolved in liquid scintillator
(Xe-LS) contained in a teardrop-shaped inner balloon (IB). The IB is
constructed from 25-$\mu$m-thick transparent nylon film and is suspended at
the center of the KamLAND detector. The IB is surrounded by 1 kton of LS
(Outer LS) contained in a 13-m-diameter spherical outer balloon. Scintillation
light is viewed by 1,325 17-inch and 554 20-inch photomultiplier tubes (PMTs)
mounted on an 18-m-diameter spherical stainless-steel tank. In the previous
phase of the experiment (KamLAND-Zen 400) with 381 kg of enriched xenon, we
probed the quasidegenerate neutrino mass region Gando2016 . To further improve
the search sensitivity, the KamLAND-Zen collaboration upgraded the experiment
(KamLAND-Zen 800) to 745 kg of enriched xenon, nearly twice the target mass of
the previous phase. To hold the increased xenon, we constructed a larger and
cleaner 3.80-m-diameter IB Gando2021 . The Xe-LS consists of 82% decane and
18% pseudocumene (1,2,4-trimethylbenzene) by volume, 2.4 g/liter of the fluor
PPO (2,5-diphenyloxazole), and $(3.13\pm 0.01)$% by weight of isotopically
enriched xenon gas. The isotopic abundances in the enriched xenon are
$(90.85\pm 0.13)\%$ 136Xe and $(8.82\pm 0.01)\%$ 134Xe.

KamLAND-Zen 800 completed data-taking in January 2024 to start dismantling the
current KamLAND detector in preparation for a detector upgrade, KamLAND2-Zen.
This Letter reports on the $0\nu\beta\beta$ decay analysis with the full
KamLAND-Zen 800 dataset collected between February 5, 2019, and January 12,
2024, and its combination with the earlier KamLAND-Zen 400 data.

The experimental signal is a characteristic peak in the energy distribution,
produced by the two $e^{-}$ emitted from 136Xe $\beta\beta$ decay. For
hypothetical $0\nu\beta\beta$ decays, their sum energy is always 2.458 MeV
(Q-value of 136Xe $\beta\beta$ decay) Redshaw2007 , while for $2\nu\beta\beta$
decays the sum has a continuous spectrum up to the Q-value. Event vertex and
energy are reconstructed based on the timing and charge distributions of
photoelectrons (p.e.) recorded by the PMTs. Energy calibration is performed
using 214Bi $\beta$’s and $\gamma$’s from 222Rn ($\tau=5.5$ days), introduced
during the initial filling of the IB with Xe-LS, and 2.225 MeV $\gamma$’s from
capture of spallation neutrons by protons. The detector Monte Carlo (MC)
simulation is based on Geant4 Agostinelli2003 ; Allison2006  and is tuned to
reproduce the timing and charge distributions observed in the data. The
estimated energy and vertex resolutions in the Xe-LS are 6.7%/$\sqrt{E({\rm
MeV})}$ and 13.7 cm/$\sqrt{E({\rm MeV})}$, respectively.

During the earlier period in KamLAND-Zen 800, we observed an increasing trend
in the number of low photo-detection efficiency PMTs caused by age-related
gain decrease. The resulting degradation in energy resolution could increase
the background from $2\nu\beta\beta$ tails. In March 2020, we started a
continuous campaign modifying the input circuitry of the KamLAND front-end
electronics (FEE) to counter the PMT gain decrease. The 1 p.e. gain and timing
offsets due to the circuit modification are corrected for each PMT based on
photon counting with high-statistics background events. Owing to these
restoration efforts, we have succeeded in preventing the degradation of the
energy resolution by a factor of 1.2, which would have lead to a twofold
increase of the $2\nu\beta\beta$ background.

Figure 1: Vertex distribution of candidate singles data (SD) events (black
points) in the $0\nu\beta\beta$ energy window $2.35<E<2.70\,{\rm MeV}$. The
thick black and red lines indicate the shape of the IB and the 1.57-m-radius
spherical volume, respectively. The dot-dashed line indicates the nylon belt
suspending the IB. The thin dashed lines illustrate the shape of the equal-
volume spherical half-shells, which compose the 2.5-m-radius spherical
fiducial volume. The high-count region at the IB bottom within the yellow
solid line indicates the hot spot and is removed through a 0.7 m radius cut.

We apply the following cuts to select $\beta\beta$ decay events: (i) The
reconstructed vertex must be within 2.5 m of the detector center and 0.7 m
away from the bottom hot spot on the IB, which is outlined in Fig. 1. (ii)
Muons and events within 2 ms after muons are rejected. (iii) Bi-Po decay pairs
are eliminated by a delayed coincidence tag, which requires the time and
distance between the prompt Bi and delayed Po decay candidates to be less than
1.9 ms and 1.7 m, respectively, and by rejecting events with a double pulse
identification within a single event acquisition window. Those cuts remove
$(99.95\pm 0.03)$% of 214Bi-214Po events, and $(97.9\pm 0.5)$% of 212Bi-212Po.
(iv) Reactor $\overline{\nu}_{e}$’s identified by a delayed coincidence of
positrons and neutron-capture $\gamma$’s are rejected. (v) Poorly
reconstructed events are rejected. These events are tagged using a vertex-
time-charge discriminator which measures how well the observed PMT time-charge
distributions match those expected based on the reconstructed vertex. The
overall selection inefficiency is less than 0.1%.

Background sources external to the Xe-LS are dominated by radioactive
impurities on the IB film. The contamination levels of 238U and 232Th are
$(4.3\pm 0.2)\times 10^{-12}$ g/g and $(2.4\pm 0.1)\times 10^{-11}$ g/g,
respectively, which are roughly a factor of 10 smaller than those measured on
the previous IB Gando2016 . The reference calculations for 238U and 232Th
mentioned here assume secular equilibrium for comparison with the previously
reported values. In the earlier period of the dataset, we found an increase in
the background rate at the IB bottom, possibly due to the settling of dust
particles containing radioactive impurities. To avoid this possible
background, we tag and remove this high-background period from the dataset
using machine learning algorithms, as discussed in Ref. Abe2023a ; Li2023 .

Solar neutrinos are an intrinsic background source for KamLAND-Zen. Based on
neutrino flux calculations using the standard solar model Serenelli2011 , the
elastic scattering (ES) of 8B solar neutrinos on electrons in the Xe-LS is
estimated to be $(4.9\pm 0.2)\times 10^{-3}$ (ton day)-1, including the effect
of three-flavor neutrino oscillations. The charged-current (CC) interactions
on 136Xe produce $e^{-}$, 136Cs, and $\gamma$’s from the excited states, and
the subsequent decays of 136Cs ($\tau=19.0$ days, $Q=2.548$ MeV) create a
background peak around 2.0 MeV in visible energy, mostly overlapping with the
resolution tail of $2\nu\beta\beta$ decays. The interaction rate is expected
to be $(0.8\pm 0.1)\times 10^{-3}$ (ton day)-1 based on the cross section
calculated in Ref. Ejiri2014 ; Frekers2013 .

Cosmic-ray muons produce neutrons and radioactive isotopes through nuclear
spallation of carbon and xenon, which decay by emitting $\beta$’s or
$\gamma$’s. To remove these events, we apply cuts based on space and time
correlation with muon and neutron capture $\gamma$ events. In the
$0\nu\beta\beta$ window, decays of 10C ($\tau=27.8$ s, $Q=3.65$ MeV) and 6He
($\tau=1.16\,{\rm s}$, $Q=3.51\,{\rm MeV}$) dominate the muon spallation
backgrounds. To reduce these short-lived backgrounds, we remove events within
150 ms after muons, and events reconstructed within 1.6 m of neutron vertices
for 180 s. In addition, we reject remaining background events by the “shower”
likelihood method based on reconstructed muon shower profiles Li2014 ; Li2015a
; Li2015b ; Zhang2016 ; Abe2023b . The overall rejection efficiencies for 10C
and 6He are $>$$99.7$% and $(97.3\pm 1.5)\%$, respectively. To reduce the
137Xe ($\tau=5.5\,{\rm min}$, $Q=4.17\,{\rm MeV}$) background, we remove
events reconstructed within 1.6 m for 27 min of vertices identified as neutron
captures on 136Xe, which produce high energy $\gamma$’s ($Q=4.03\,{\rm MeV}$).
This cut removes $(74\pm 7)\%$ of 137Xe.

The products of muon spallation with lifetimes greater than $\mathcal{O}$(100
s), denoted as “long-lived products”, are attributed to the decay of heavy
isotopes produced by xenon spallation. Xenon spallation can be tagged by
detecting multiple neutrons. To characterize the long-lived products, we
define a likelihood ratio, $R_{L}=L_{\rm spa}/(L_{\rm spa}+L_{\rm acc})$. Here
$L_{\rm spa}$ and $L_{\rm acc}$ are the probability density functions (PDFs)
for long-lived muon-spallation pairs and accidental pairs, respectively. These
PDFs are constructed as a function of neutron multiplicity, distance to
neutron vertices, and time interval from preceding muons. The cut value on
$R_{L}$ is optimized using an estimate with MC simulation tools. FLUKA
BOHLEN2014211 ; Ferrari:898301  is used to calculate the spallation isotope
yields, and Geant4 to calculate the radioactive decays and energy spectra of
the isotopes, considering their sequential decay chain. The predicted
backgrounds are primarily from 132I, 130I, 124I, 122I, 118Sb, 110In, 88Y, and
from other isotopes in smaller amounts (Table IX in Ref. Abe2023b ). The total
background rate in the $0\nu\beta\beta$ window is $30.0\pm 2.2$ (ton yr)-1 in
total, in units of 136Xe exposure. Events that are not classified as coming
from long-lived backgrounds are referred to as “singles data” (SD), and the
others are referred to as “long-lived data” (LD).

Figure 2: Energy spectra of selected $\beta\beta$ candidates within a
1.57-m-radius spherical volume drawn together with best-fit backgrounds, the
$2\nu\beta\beta$ decay spectrum, and the 90% C.L. upper limit for
$0\nu\beta\beta$ decay of (a) singles data (SD), and (b) long-lived data (LD).
The LD exposure is about 10% of the SD exposure.

We improved the neutron identification for post-muon events to enhance the
tagging efficiency of the long-lived backgrounds. We identify these neutrons
with a second, deadtime-free electronics system, called MoGURA Abe2023b .
Previously, the neutron vertex was reconstructed using only PMTs with normal
gain. In this analysis, we added low gain PMTs to the photon counting,
achieving an increase in the mean number of PMT hits from $181$ to $201$ for
2.225 MeV $\gamma$’s from neutron capture, after subtracting afterpulse noise
contributions. The resulting neutron tagging efficiency is estimated to be
$(74.5\pm 0.4)\%$ based on the neutron capture time distribution,
approximately 2% better than in the previous analysis. Our MC study shows that
$(47.1\pm 8.7)\%$ of long-lived spallation backgrounds are classified as LD,
whereas only $9.0$% of uncorrelated events are mis-classified.

The total livetime for SD and LD in KamLAND-Zen 800 is $1131$ days and $111$
days, respectively. The total exposure of SD, which is sensitive to the
$0\nu\beta\beta$ signal, is $2097$ kg yr of 136Xe, more than double the amount
of the previous search Abe2023a . The $0\nu\beta\beta$ decay rate is estimated
from a simultaneous likelihood fit to the binned energy spectra of SD and LD
between 0.5 and 4.8 MeV. The binned energy spectra are made by dividing the
2.5-m-radius fiducial volume into 40 equal-volume bins (illustrated in Fig.
1), considering the time variation of event rates. The contributions from
major backgrounds in the Xe-LS, such as 85Kr, 40K, 210Bi, the 228Th-208Pb
subchain of the 232Th series, and long-lived spallation products, are free
parameters and are left unconstrained in the spectral fit. The contributions
from the 222Rn-210Pb subchain of the 238U series and short-lived spallation
products are constrained by their independent measurements. The parameters of
the detector energy response model are constrained by 214Bi data. The energy
spectral distortion parameter (a linear rescaling) for the long-lived
spallation products is constrained by its estimated uncertainty (Fig. 10 in
Ref. Abe2023b ).

Figure 3: Allowed region for the 136Xe $0\nu\beta\beta$ rate and the long-
lived spallation background rate in the energy region $2.35<E<2.70\,{\rm MeV}$
($0\nu\beta\beta$ window). The red dot and contour lines correspond,
respectively, to the best-fit and the $1\sigma$, 90%, 95%, 99% C.L. intervals
for two degrees of freedom. The horizontal line indicates the MC-based
prediction.

Figure 2 shows the energy spectra of selected candidate SD and LD events
within a 1.57-m-radius spherical volume together with the best-fit curves. The
exposure of 136Xe for SD in this reduced volume is 1.13 ton yr. The best-fit
background contributions are summarized in Table 1. The $0\nu\beta\beta$
signal best-fit is 0 events, indicating no event excess over the background
expectation. We obtain a 90% confidence level (C.L.) upper limit on the number
of $0\nu\beta\beta$ decays of $<$ $10.0$ events ($<$ $8.0$ events in the range
$2.35<E<2.70\,{\rm MeV}$), which corresponds to a limit of $<$ $8.9$ (ton
yr)-1 in units of 136Xe exposure, or $T_{1/2}^{0\nu\beta\beta}>3.4\times
10^{26}$ yr (90% C.L.). An analysis based on the Feldman-Cousins procedure
PhysRevD.57.3873_FCmethod  gives a slightly stronger limit of $4.3\times
10^{26}$ yr (90% C.L.), indicating a limited impact of the physical boundary
on the $0\nu\beta\beta$ rate at low statistics. An MC simulation of an
ensemble of experiments assuming the best-fit background spectrum rwithout a
$0\nu\beta\beta$ signal indicates a median sensitivity of $2.3\times 10^{26}$
yr. The probability of obtaining a limit stronger than that reported here is
24%. In addition to the frequentist analyses above, we also performed a
statistical analysis within the Bayesian framework, assuming a flat prior for
$1/T_{1/2}^{0\nu\beta\beta}$. The Bayesian limit and median sensitivity are
$3.4\times 10^{26}$ yr and $2.4\times 10^{26}$ yr (90% C.L.), respectively.

Figure 4: Effective Majorana neutrino mass $\left<m_{\beta\beta}\right>$ as a
function of the lightest neutrino mass $m_{\rm lightest}$. The dark shaded
regions are predictions based on best-fit values of neutrino oscillation
parameters for (a) the normal ordering (NO) and (b) the inverted ordering
(IO), and the light shaded regions indicate the $3\sigma$ ranges calculated
from oscillation parameter uncertainties DellOro2014 ; NuFIT2020 . The
horizontal lines indicate 90% C.L. upper limits on
$\left<m_{\beta\beta}\right>$ with 136Xe from KamLAND-Zen (this work),
considering an improved phase space factor calculation Kotila2012 ; Stoica2013
and phenomenological NME calculations: shell model Menendez2018 ;
PhysRevC.93.024308horoi ; PhysRevC.101.044315.coraggio2020 ;
PhysRevC.105.034312.coraggio2022  (dot-dashed lines), quasiparticle
randomphase approximation (QRPA) PhysRevC.87.064302.mustonen ;
PhysRevC.91.024613.Hyv ; PhysRevC.98.064325.simko2018 ; PhysRevC.97.045503fang
; PhysRevC.102.044303terasaki (dotted lines), energy-density functional (EDF)
theory PhysRevLett.105.252503.rodriguez ; PhysRevLett.111.142501lopez ;
PhysRevC.95.024305song  (solid lines), interacting boson model (IBM)
PhysRevC.91.034304barea ; PhysRevD.102.095016deppisch  (dashed lines). (c) The
corresponding limits for 136Xe, 76Ge Agostini2020 , and 130Te Adams2022 .
Three theoretical predictions in the IO region are also shown: (A)Harigaya2012
, (B)Asaka2020 , (C)Asai2020 . Table 1: Summary of the estimated and best-fit
background contributions in the energy region $2.35<E<2.70\,{\rm MeV}$ within
the 1.57-m-radius spherical volume. In total, $59$ events were observed.

[h]  Background Estimated Best-fit 136Xe $2\nu\beta\beta$ - $27.60$ Residual
radioactivity in Xe-LS 238U series $0.08\pm 0.01$ $0.08$ 232Th series - $1.35$
External (Radioactivity in IB) 238U series - $6.61$ 232Th series - $0.04$
Neutrino interactions 8B solar $\nu$ $e^{-}$ ES $3.57\pm 0.09$ $3.57$
Spallation products Long-lived $17.82\pm 1.30$ $\dagger$ $19.39$ 10C $0.00\pm
0.05$ $0.00$ 6He $0.54\pm 0.27$ $0.56$ 137Xe $0.72\pm 0.60$ $0.72$

  * $\dagger$

Estimation based on the spallation MC study. This event rate constraint is not
applied in the spectral fit.

Figure 3 shows the allowed region of the 136Xe $0\nu\beta\beta$ rate and the
long-lived spallation background rate from a combined fit of the KamLAND-Zen
400 and 800 datasets, giving a limit of

$T_{1/2}^{0\nu\beta\beta}>3.8\times 10^{26}\,\textrm{yr~{}~{}~{}(90\% C.L.)}.$
(1)

The best-fit long-lived spallation background rate is $30.2\pm 4.5$ (ton
yr)-1indicating good consistency between the MC-based prediction and the LD
analysis. This combined analysis has a sensitivity of $2.6\times 10^{26}$ yr,
and the probability of obtaining a stronger limit is $29$%. The lower limit on
$T_{1/2}^{0\nu\beta\beta}$ can be converted to an upper limit on $\langle
m_{\beta\beta}\rangle$ assuming that the $0\nu\beta\beta$ decay mechanism is
dominated by exchange of light Majorana neutrinos. Assuming the axial coupling
constant $g_{A}\simeq 1.27$, we obtain an upper limit of

$\langle m_{\beta\beta}\rangle<(28-122)\,\textrm{meV~{}~{}~{}(90\% C.L.)},$
(2)

using an improved phase space factor calculation of $14.54\times
10^{-15}\,(\textrm{yr})^{-1}$ Kotila2012 ; Stoica2013 , and phenomenological
nuclear matrix element (NME) calculations, summarized in Table 2. In this
letter, we follow the current standard to use those phenomenological models,
although recent theoretical work has produced new NME calculations
belley2023ab ; Jokiniemi2023  incorporating both quenching physics
Faessler_2008 ; PhysRevC.87.014315  and short-range contact
PhysRevLett.120.202001  contributions, which are missing in the calculations
in Table 2. Considering the calculations from Ref. belley2023ab ;
Jokiniemi2023  gives ranges of upper limits on $\langle m_{\beta\beta}\rangle$
of $71-125$ meV (ab initio belley2023ab ), $25-60$ meV (pnQRPA Jokiniemi2023
), and $44-138$ meV (shell model Jokiniemi2023 ). Figure 4 illustrates the
allowed range of $\langle m_{\beta\beta}\rangle$ as a function of the lightest
neutrino mass for the NME values in Table 2. The shaded regions include the
uncertainties in mixing elements $U_{ei}$ and the neutrino mass splitting for
each hierarchy. Also drawn are experimental limits from $0\nu\beta\beta$ decay
searches in other nuclei Agostini2020 ; Adams2022 . The limit using 136Xe by
KamLAND-Zen has provided the most stringent tests of the effective Majorana
mass in the IO region and several theoretical models Harigaya2012 ; Asaka2020
; Asai2020 . It also provides the strongest constraint on $m_{\rm lightest}$
considering extreme cases of the combination of CP phases and the
uncertainties from neutrino oscillation parameters DellOro2014 ; NuFIT2020 .
We obtain a 90% C.L. upper limit of $m_{\rm lightest}<(84-353)\,{\rm meV}$.

In conclusion, we presented a $0\nu\beta\beta$ decay search with the complete
KamLAND-Zen dataset including improvements in the rejection of muon-induced
long-lived spallation backgrounds. This search places the most stringent
constraints to date on $\left<m_{\beta\beta}\right>$, and reaches the inverted
mass ordering region below 50 meV with half of the phenomenological nuclear
matrix element calculations. The sensitivity is limited by background from
xenon spallation products and $2\nu\beta\beta$ decay. In the near future, we
plan to upgrade the KamLAND detector with more 136Xe mass, enhanced energy
resolution to minimize the leakage of the $2\nu\beta\beta$ tail into the
$0\nu\beta\beta$ analysis window, and enhanced rejection of spallation
products through improved neutron tagging and particle identification methods.

Table 2: Summary of commonly used phenomenological nuclear matrix elements ($M^{0\nu}$) for 136Xe $0\nu\beta\beta$ decay and the corresponding 90% C.L. upper limits on the effective Majorana neutrino mass $\left<m_{\beta\beta}\right>$ from KamLAND-Zen data. The limits shown in Eq. (2) use the minimum and maximum. |  Ref. |  $M^{0\nu}$ |  $\langle m_{\beta\beta}\rangle$ (meV)
---|---|---|---
Shell model | Menendez2018  | 2.28, 2.45 | 59.4, 55.3
PhysRevC.93.024308horoi  | 1.63, 1.76 | 83.1, 77.0
PhysRevC.101.044315.coraggio2020 ; PhysRevC.105.034312.coraggio2022  | 2.39 | 56.7
QRPA | PhysRevC.87.064302.mustonen  | 1.55 | 87.4
PhysRevC.91.024613.Hyv  | 2.91 | 46.6
PhysRevC.98.064325.simko2018  | 2.71 | 50.0
PhysRevC.97.045503fang  | 1.11, 1.18 | 122, 115
PhysRevC.102.044303terasaki  | 3.38 | 40.1
EDF theory | PhysRevLett.105.252503.rodriguez  | 4.20 | 32.3
PhysRevLett.111.142501lopez  | 4.77 | 28.4
PhysRevC.95.024305song  | 4.24 | 32.0
IBM | PhysRevC.91.034304barea  | 3.25 | 41.7
PhysRevD.102.095016deppisch  | 3.40 | 39.9

###### Acknowledgements.

The KamLAND-Zen experiment is supported by JSPS KAKENHI Grants No. 21000001,
No. 26104002, and No. 19H05803; the U.S. National Science Foundation awards
no. 2110720 and no. 2012964; the Heising-Simons Foundation; the Dutch Research
Council (NWO); and under the U.S. Department of Energy (DOE) Grant No. DE-
AC02-05CH11231, as well as other DOE and NSF grants to individual
institutions. The Kamioka Mining and Smelting Company has provided service for
activities in the mine. We acknowledge the support of NII for SINET4.

## References

  * (1) M. Agostini, G. Benato, J. A. Detwiler, J. Menéndez, and F. Vissani,  Rev. Mod. Phys. 95, 025002 (2023).
  * (2) M. Fukugita and T. Yanagida,  Phys. Lett. B 174, 45 (1986).
  * (3) S. Abe et al., (KamLAND-Zen Collaboration),  Phys. Rev. Lett. 130, 051801 (2023).
  * (4) K. Harigaya, M. Ibe, and T. T. Yanagida,  Phys. Rev. D 86, 013002 (2012).
  * (5) T. Asaka, Y. Heo, and T. Yoshida,  Phys. Lett. B 811, 135956 (2020).
  * (6) K. Asai,  The European Physical Journal C 80, 76 (2020).
  * (7) M. Agostini, G. Benato, and J. A. Detwiler,  Phys. Rev. D 96, 053001 (2017).
  * (8) A. Gando et al., (KamLAND-Zen Collaboration),  Phys. Rev. C 85, 045504 (2012).
  * (9) A. Gando et al., (KamLAND-Zen Collaboration),  Phys. Rev. C 86, 021601 (2012).
  * (10) A. Gando et al., (KamLAND-Zen Collaboration),  Phys. Rev. Lett. 110, 062502 (2013).
  * (11) K. Asakura et al., (KamLAND-Zen Collaboration),  Nucl. Phys. A 946, 171 (2016).
  * (12) A. Gando et al., (KamLAND-Zen Collaboration),  Phys. Rev. Lett. 117, 082503 (2016).
  * (13) A. Gando et al., (KamLAND-Zen Collaboration),  Phys. Rev. Lett. 122, 192501 (2019).
  * (14) Y. Gando et al., (KamLAND-Zen Collaboration),  JINST 16, P08023 (2021).
  * (15) M. Redshaw, E. Wingfield, J. McDaniel, and E. G. Myers,  Phys. Rev. Lett. 98, 053003 (2007).
  * (16) S. Agostinelli et al.,  Nucl. Instr. Meth. A 506, 250 (2003).
  * (17) J. Allison et al.,  IEEE Trans. Nucl. Sci. 53, 270 (2006).
  * (18) A. Li et al.,  Phys. Rev. C 107, 014323 (2023).
  * (19) A. M. Serenelli, W. C. Haxton, and C. Peña-Garay,  Astrophys. J. 743, 24 (2011).
  * (20) H. Ejiri and S. R. Elliott,  Phys. Rev. C 89, 055501 (2014).
  * (21) D. Frekers, P. Puppe, J. H. Thies, and H. Ejiri,  Nucl. Phys. A 916, 219 (2013).
  * (22) S. W. Li and J. F. Beacom,  Phys. Rev. C 89, 045801 (2014).
  * (23) S. W. Li and J. F. Beacom,  Phys. Rev. D 91, 105005 (2015).
  * (24) S. W. Li and J. F. Beacom,  Phys. Rev. D 92, 105033 (2015).
  * (25) Y. Zhang et al., (Super-Kamiokande Collaboration),  Phys. Rev. D 93, 012004 (2016).
  * (26) S. Abe et al., (KamLAND-Zen Collaboration),  Phys. Rev. C 107, 054612 (2023).
  * (27) T. Böhlen et al.,  Nuclear Data Sheets 120, 211 (2014).
  * (28) A. Ferrari and et al.,  FLUKA: A multi-particle transport code (program version 2005) CERN Yellow Reports: Monographs (CERN, Geneva, 2005).
  * (29) G. J. Feldman and R. D. Cousins,  Phys. Rev. D 57, 3873 (1998).
  * (30) S. Dell’Oro, S. Marcocci, and F. Vissani,  Phys. Rev. D 90, 033005 (2014).
  * (31) Nufit 5.0, available at http://www.nu-fit.org (2020).
  * (32) J. Kotila and F. Iachello,  Phys. Rev. C 85, 034316 (2012).
  * (33) S. Stoica and M. Mirea,  Phys. Rev. C 88, 037303 (2013); updated in arXiv:1411.5506v3 .
  * (34) J. Menéndez,  J. of Phys. G 45, 014003 (2018).
  * (35) M. Horoi and A. Neacsu,  Phys. Rev. C 93, 024308 (2016).
  * (36) L. Coraggio, A. Gargano, N. Itaco, R. Mancino, and F. Nowacki,  Phys. Rev. C 101, 044315 (2020).
  * (37) L. Coraggio et al.,  Phys. Rev. C 105, 034312 (2022).
  * (38) M. T. Mustonen and J. Engel,  Phys. Rev. C 87, 064302 (2013).
  * (39) J. Hyvärinen and J. Suhonen,  Phys. Rev. C 91, 024613 (2015).
  * (40) F. S̆imkovic, A. Smetana, and P. Vogel,  Phys. Rev. C 98, 064325 (2018).
  * (41) D.-L. Fang, A. Faessler, and F. S̆imkovic,  Phys. Rev. C 97, 045503 (2018).
  * (42) J. Terasaki,  Phys. Rev. C 102, 044303 (2020).
  * (43) T. R. Rodríguez and G. Martínez-Pinedo,  Phys. Rev. Lett. 105, 252503 (2010).
  * (44) N. L. Vaquero, T. R. Rodríguez, and J. L. Egido,  Phys. Rev. Lett. 111, 142501 (2013).
  * (45) L. S. Song, J. M. Yao, P. Ring, and J. Meng,  Phys. Rev. C 95, 024305 (2017).
  * (46) J. Barea, J. Kotila, and F. Iachello,  Phys. Rev. C 91, 034304 (2015).
  * (47) F. F. Deppisch, L. Graf, F. Iachello, and J. Kotila,  Phys. Rev. D 102, 095016 (2020).
  * (48) M. Agostini et al., (GERDA Collaboration),  Phys. Rev. Lett. 125, 252502 (2020).
  * (49) D. Q. Adams et al., (CUORE Collaboration),  Nature 604, 53 (2022).
  * (50) A. Belley, T. Miyagi, S. R. Stroberg, and J. D. Holt,  arXiv:2307.15156v1.
  * (51) L. Jokiniemi, B. Romeo, P. Soriano, and J. Menéndez,  Phys. Rev. C 107, 044305 (2023).
  * (52) A. Faessler et al.,  J. Phys. G: Nucl. Part. Phys. 35, 075104 (2008).
  * (53) J. Barea, J. Kotila, and F. Iachello,  Phys. Rev. C 87, 014315 (2013).
  * (54) V. Cirigliano et al.,  Phys. Rev. Lett. 120, 202001 (2018).


# Optimal incentives in a limit order book: a SPDE control approach111This
work benefits from the financial support of the Chaire Deep Finance and
Statistics. The authors would like to thank René Aïd (Université Paris-
Dauphine) for the insightful discussions they had on the subject.

Bastien Baldacci222Quantitative Advisory Solutions,
<EMAIL_ADDRESS>Philippe Bergault333École Polytechnique,
CMAP, 91128, Palaiseau, France<EMAIL_ADDRESS>(Corresponding author.)

###### Abstract

With the fragmentation of electronic markets, exchanges are now competing in
order to attract trading activity on their platform. Consequently, they
developed several regulatory tools to control liquidity provision /
consumption on their liquidity pool. In this paper, we study the problem of an
exchange using incentives in order to increase market liquidity. We model the
limit order book as the solution of a stochastic partial differential equation
(SPDE) as in [12]. The incentives proposed to the market participants are
functions of the time and the distance of their limit order to the mid-price.
We formulate the control problem of the exchange who wishes to modify the
shape of the order book by increasing the volume at specific limits. Due to
the particular nature of the SPDE control problem, we are able to characterize
the solution with a classic Feynman-Kac representation theorem. Moreover, when
studying the asymptotic behavior of the solution, a specific penalty function
enables the exchange to obtain closed-form incentives at each limit of the
order book. We study numerically the form of the incentives and their impact
on the shape of the order book, and analyze the sensitivity of the incentives
to the market parameters.

Key words: Stochastic control, SPDEs, make-take fees, financial regulation,
limit order book.

##  1 Introduction

At the end of the twentieth century, the electronification of financial
markets started first with stock exchanges, before spreading to every asset
class. On most order-driven markets, computers now handle almost all the
market activity, organized around all-to-all limit order books. The recent
years have also seen the fragmentation of electronic markets, meaning that the
same asset can be traded in several platforms, thus increasing the competition
between exchanges to attract trading activity on their liquidity pools, see
for example [20]. The usual mechanism put in place by exchanges is the make-
take fees system: a fee rebate is associated with executed limit orders while
transaction costs are applied to market orders.444An other way for the
exchange to attract liquidity on its platform is to set an appropriate tick
size, see [7, 14] for instance. This way, exchanges charge in an asymmetric
way liquidity provision and liquidity consumption. One unintended effect of
these incentives mechanism has been the rise of high-frequency trading
activity, see among others [2], which is now responsible for the majority of
liquidity provision on electronic markets. Indeed, traditional market making
companies are now competing with all market participants who can post
liquidity-providing orders. One of the main concerns regarding this activity
is the quality of liquidity provided in times of stress, when the high-
frequency traders tend to leave the market, see [8, 21, 22, 23]. The influence
of the make-take fees mechanism has been studied mostly from an empirical and
data-driven viewpoint, see [1, 10, 11, 19] showing the dependence of depths,
volumes, or price impact of the order book on the make-take fees structure.

The recent work of [15] proposes a make-take fees design, combining the
Avellaneda-Stoikov [3] model and a Principal-Agent framework with drift
control,555There is an extensive literature on this topic, see the pioneering
work of [25] and recent developments such as [16, 17]. to increase market
liquidity on a single asset. In their formulation, the exchange offers a
contract, resulting from a Stackelberg equilibrium between the market maker
and the exchange, to the market maker in order to decrease the bid-ask spread
offered on the asset. The optimal contract proposed is a remuneration indexed
on the number of transactions realized on the asset as well as its dynamics.
This contract is obtained in quasi closed-form, depending on the market
maker’s inventory trajectory and market volatility. It is shown numerically
that the contract effectively reduces the spreads and subsequently increases
liquidity provision. This work has been extended in various ways, taking into
account several market participants, the specificities of options market and
enabling trading on dark liquidity pools, see [4, 5, 6]. However, the very
design of make-take fees through a Principal-Agent has three main practical
drawbacks. First, the form of the incentives is a Markovian function of the
inventory of the market maker. While it is reasonable for the exchange to
track the inventory of the agent on the platform, it neglects the fact that
she may have submitted buy or sell orders on other platforms where the asset
is listed. It induces a possible arbitrage opportunity for the market maker,
who can take advantage of the fact that the exchange can only see the portion
of her inventory on its platform to earn a higher incentive. Second, the
market making model used in these articles is the one of Avellaneda and
Stoikov [3], which is far more suited to Over-The-Counter (OTC) markets than
to limit order books.666See [9] or [18] for more details. In particular, there
is no notion of discrete price limits. Third, the use of agent-based models
for incentives design introduces the notion of risk-aversion for the
participants. While the risk-aversion parameter of the exchange can be chosen
in order to produce a certain type of incentives, the one of the agent is
completely unknown to the exchange. As the values of the incentives are highly
dependent on this parameter, its under or overestimation can lead to very
different contracts.

In this article, we propose an incentives design with a complete different
paradigm. We do not rely on an agent-based model, in the sense that the
incentives do not emerge from the interactions between market participants and
the exchange. We consider a single asset whose order book dynamics is governed
by a system of stochastic partial differential equations (SPDE) on the bid and
ask sides, described in the recent work of Cont and Müller [12]. Two great
advantages of this model are its accuracy with respect to market order-flow
dynamics and its tractability via a parametrization using low-dimensional
diffusion processes. We consider that the exchange is looking for the best
incentives policy to offer to liquidity providers in order to increase the
rate of buy or sell order submissions, so that the shape of the order book is
modified. Ideally, the exchange wishes to have the maximum possible volume on
the bid and ask sides at the first limits. Its willingness to increase the
market liquidity decreases with the depth of the order book. Thus, we avoid
the issues of selecting the risk-aversion parameters of the agents. We
formulate the stochastic control problem of the exchange and prove a
verification theorem characterizing the optimal incentives function. The main
result of this article is that, given the form of the control problem of the
exchange, the Backward Stochastic Partial Differential Equation (BSPDE)
characterizing the problem boils down in fact to a linear parabolic Partial
Differential Equation (PDE). This simplifies drastically the control problem,
as it admits a classical Feynman-Kac representation. Moreover, with the choice
of suitable intensity and penalty functions, the asymptotic solution of the
control problem is obtained in closed-form by solving a simple second order
linear ordinary differential equation (ODE) with constant coefficients. It is
of particular interest for the practical design of incentives as they are now
analytically obtained as functions of the distance from the mid-price and the
market parameters. In particular, the flexibility of the model allows the
exchange to modify specific limits of the order book and to design various
forms of incentives (increase overall liquidity, more emphasis on first and
second limit etc.). We illustrate numerically the form of the incentives
produced by the model and its influence on the shape of the limit order book
(i.e. on market liquidity) by simulating the associated SPDE controlled by the
optimal incentives. We also study numerically the dependence of our model to
the different market parameters.

The article is designed as follows. In Section 2, we describe the order book
dynamics and the optimization problem of the exchange. In particular, we
propose a slightly modified version of the limit order book model introduced
in [12], to take into account the impact of incentives. In Section 3, we
provide a BSPDE representation of the optimal control problem, and show that
it admits a unique solution. We also prove a verification theorem
characterizing the optimal control of the exchange. Section 4 is devoted to
the asymptotic analysis of the control problem and the associated closed-form
solutions, while Section 5 outlines the main numerical results of the
incentives model. The proofs are relegated to the Appendix.

##  2 The model

Throughout the paper, we consider a filtered probability space
$\left(\Omega,\mathcal{F},\mathbb{P};\mathbb{F}=(\mathcal{F}_{t})_{t\geq
0}\right)$ satisfying the usual conditions. We assume this probability space
to be large enough to support all the processes we introduce.

###  2.1 A Cont-Müller [12] limit order book

Let us introduce a market consisting of a single asset, whose price process is
denoted by $(S_{t})_{t\in[0,T]}$ over a time horizon $[0,T]$, with $T>0.$ We
consider a model for the limit order book (LOB) based on the work of Cont and
Müller [12]. In this paper, the authors introduce a LOB modelled by a density
$u:(t,x)\in[0,T]\times[-L,L]\mapsto u_{t}(x)$ representing the volume of
orders available at time $t$ at price $S_{t}+x$ per unit price, at the bid for
$x<0$ and at the ask for $x>0$. For mathematical convenience, we limit the
range of the argument $x$ to a bounded interval $[-L,L]$, with $L>0$, meaning
that one cannot propose a bid (resp. ask) price lower (resp. higher) than
$S_{t}-L$ (resp. $S_{t}+L$).777As $L$ can be taken arbitrarily large, this
assumption is not restrictive for practical applications. The function $u$
satisfies the following stochastic partial differential equation
(SPDE):888$u_{t}$ is the density of the order book at time $t$, i.e. the
volume available at the first limit at the ask (for instance) at time $t$ is
given by $\int_{0}^{\delta}u_{t}(x)dx$, where $\delta$ is the tick size of the
asset. See [12] for more details.

$\displaystyle du_{t}(x)=\left(\vphantom{\beta_{b}\nabla
u_{t}(x)+\alpha_{b}u_{t}(x)-f^{b}(x)}\eta_{a}\Delta u_{t}(x)+\beta_{a}\nabla
u_{t}(x)+\alpha_{a}u_{t}(x)+f^{a}(x)\right)dt+\sigma_{a}u_{t}(x)dW^{a}_{t}\quad\forall
x\in(0,L),$ (1) $\displaystyle du_{t}(x)=\left(\eta_{b}\Delta
u_{t}(x)-\beta_{b}\nabla
u_{t}(x)+\alpha_{b}u_{t}(x)-f^{b}(x)\right)dt+\sigma_{b}u_{t}(x)dW^{b}_{t}\quad\forall
x\in(-L,0),$

with boundary conditions:

$\displaystyle u_{t}(0^{+})=u_{t}(0^{-})=u_{t}(-L)=u_{t}(L)=0\quad\forall
t\in[0,T],$ (2)

with $u_{0}:[-L,L]\rightarrow u_{0}(x)$ fixed, and by convention $u_{t}(x)\leq
0$ for $x<0$ and $u_{t}(x)\geq 0$ for $x>0$, where, as mentioned above, $x$ is
the (signed) distance from the mid-price. The process
$(W^{a}_{t},W^{b}_{t})_{t\in[0,T]}$ is a 2-dimensional Brownian motion
(possibly correlated), and the parameters are taken as follows :
$\eta_{a},\eta_{b},\sigma_{a},\sigma_{b}>0$, $\beta_{a},\beta_{b}\geq 0$ and
$\alpha_{a},\alpha_{b}\leq 0$. $f^{a}:(0,L)\rightarrow\mathbb{[}0,+\infty)$
and $f^{b}:(-L,0)\rightarrow[0,+\infty)$ are two functions. The interpretation
is given below:

  * •

$f^{b}$ (resp. $f^{a}$) represents the rate of buy (resp. sell) order
submissions at a distance $x$ from the mid-price;

  * •

the term $\alpha_{b}u_{t}(x)$ (resp. $\alpha_{b}u_{t}(x)$) represents the
proportional cancellation of limit buy (resp. sell) orders at a distance $x$
from the mid-price;

  * •

the convection term $-\beta_{b}\nabla u_{t}(x)$ (resp. $\beta_{a}\nabla
u_{t}(x)$) represents the replacement of buy (resp. sell) orders by orders
closer to the mid-price;

  * •

the diffusion term $\eta_{b}\Delta u_{t}(x)$ (resp. $\eta_{a}\Delta u_{t}(x)$)
represents the cancellation and symmetric replacement of orders at a distance
$x$ from the mid-price;

  * •

the multiplicative noise term $\sigma_{b}u_{t}(x)dW^{b}_{t}$ (resp.
$\sigma_{a}u_{t}(x)dW^{a}_{t}$) accounts for the high-frequency submissions
and cancellations associated with HFT orders.

We will now considered a slightly modified version of the above model, to take
into account the use of incentives from the exchange to optimize the shape of
the LOB. More precisely, we assume that $u$ is now solution to the following
SPDE:

$\displaystyle du_{t}(x)=\left(\eta_{a}\Delta u_{t}(x)+\beta_{a}\nabla
u_{t}(x)+\alpha_{a}u_{t}(x)+f^{a}(x,Z^{a}_{t}(x))\right)dt+\sigma_{a}u_{t}(x)dW^{a}_{t}\
\forall x\in(0,L),$ (3) $\displaystyle du_{t}(x)=\left(\eta_{b}\Delta
u_{t}(x)-\beta_{b}\nabla
u_{t}(x)+\alpha_{b}u_{t}(x)-f^{b}(x,Z^{b}_{t}(x))\right)dt+\sigma_{b}u_{t}(x)dW^{b}_{t}\
\forall x\in(-L,0),$

where $f^{a},f^{b}:(0,L)\times[0,+\infty)\rightarrow[0,+\infty)$ verify:

  * •

$f^{a}$ and $f^{b}$ are $C^{1}$;

  * •

$f^{a}(x,.)$ and $f^{b}(x,.)$ are increasing for all $x\in(0,L)$;

  * •

$f^{a}(x,.)$ and $f^{b}(x,.)$ are concave for all $x\in(0,L)$.

######  Remark 1.

Although $f^{a}$ and $f^{b}$ are typically decreasing in $x$, our model does
not require it.

In Equation (3), the processes $Z^{a}_{t}(x)$ et $Z^{b}_{t}(x)$ are incentives
that the exchange commits to give to the trader placing a unit order at a
distance $x$ from the mid-price at time $t$, if the order is executed (hence
the incentive is not received by the trader at time $t$ but at the time of the
trade). They belong respectively to the sets of admissible controls
$\mathcal{A}^{a}$ and $\mathcal{A}^{b}$ defined by

$\displaystyle\mathcal{A}^{j}\\!:=\\!\bigg{\\{}Z:\Omega\times[0,T]\times(0,L)\rightarrow[0,+\infty)\Big{|}Z$
$\displaystyle\text{ is
}\mathcal{P}\otimes\mathcal{B}\big{(}(0,L)\big{)}\text{-measurable,}$
$\displaystyle\text{ and
}\mathbb{E}\\!\left[\\!\int_{0}^{T}\\!\\!\\!\int_{0}^{L}f^{j}\left(x,Z_{t}(x)\right)^{2}dxdt\right]\\!<\\!+\infty\\!\bigg{\\}}$

for $j\in\\{a,b\\},$ where $\mathcal{P}$ denotes the $\sigma$-algebra of
$\mathbb{F}$-predictable subsets of $\Omega\times[0,T]$, and
$\mathcal{B}\big{(}(0,L)\big{)}$ denotes the Borelian sets of $(0,L)$.

###  2.2 The optimization problem

The exchange has two symmetrical problems, on the bid side and ask side of the
LOB, respectively. We only consider here the problem on the ask side (the
problem on the bid side is treated similarly). The exchange aims at maximizing

$\displaystyle
J(Z^{a}):=\mathbb{E}\left[\int_{0}^{T}\left(\int_{0}^{L}\left(u^{Z^{a}}_{t}(x)-g(x,Z^{a}_{t}(x))\right)dx\right)dt\right],$

over the set of admissible controls $\mathcal{A}^{a}$, where
$g:(0,L)\times[0,+\infty)\rightarrow[0,+\infty)$ is a penalty function which
verifies the following assumptions:

  * •

$g(x,0)=0$ for all $x\in(0,L)$;

  * •

$g$ is nondecreasing in both its variables;

  * •

$g$ is $C^{1}$;

  * •

$g(x,.)$ is convex for all $x\in(0,L)$.

Hence, the exchange wants to find $\hat{Z}^{a}\in\mathcal{A}^{a}$ such that

$\displaystyle
J(\hat{Z}^{a})=\underset{Z^{a}\in\mathcal{A}^{a}}{\sup}J(Z^{a}).$ (4)

######  Remark 2.

Ideally, the exchange would like the LOB to be as full as possible. However,
providing incentives induces a cost for the exchange. The function $g$ is a
penalty function that allows to take this cost into account. The first
assumption $g(x,0)=0$ is perfectly natural, as providing no incentive induces
no cost. Similarly, $g$ should be nondecreasing in its second variable, as a
higher incentive induces a higher cost. We also assume that $g$ is
nondecreasing in its first variable, to model the fact that the exchange is
mainly interested in having a lot of liquidity at the first limits of the LOB,
and is therefore ready to offer a higher incentive at those limits –
basically, we penalize “less” the cost induced by incentives at the first
limits than the cost induced by incentives at limits further away from the
mid-price.

##  3 Main results

Let us define the Hamiltonian
$H:(0,L)\times\mathbb{R}\times[0,+\infty)\times\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$
by

$H(x,u,z,p,q)=u-g(x,z)+(\alpha_{a}u+f^{a}(x,z))p+\sigma_{a}uq.$

From the properties of $g$ and $f^{a}$, we get the following lemma:

######  Lemma 1.

For all $(x,p,q)\in[0,L]\times[0,+\infty)\times\mathbb{R}$, the function

$(u,z)\in[0,+\infty)\times[0,+\infty)\mapsto H(x,u,z,p,q)$

is concave.

We introduce the differential operator $\mathcal{L}$ given by

$\mathcal{L}u_{t}(x)=\eta_{a}\Delta u_{t}(x)+\beta_{a}\nabla u_{t}(x).$

Let us then consider the following backward stochastic partial differential
equation (BSPDE) $p(t,x)\in\mathbb{R}$ and $q(t,x)\in\mathbb{R}$:

$\displaystyle
dp_{t}(x)=-\left\\{\mathcal{L}^{*}p_{t}(x)+\partial_{u}H\left(x,u_{t}(x),Z^{a}_{t}(x),p_{t}(x),q_{t}(x)\right)\right\\}dt+q_{t}(x)dW^{a}_{t},\quad\forall(t,x)\in(0,T)\times(0,L),$
(5)

where $\mathcal{L}^{*}$ is the adjoint operator of $\mathcal{L}$, with
boundary conditions:

$\displaystyle p_{T}(x)=0,\quad\forall x\in(0,L),$ (6) $\displaystyle
p_{t}(x)=0,\quad\forall(t,x)\in(0,T)\times\\{0,L\\}.$

Note that we can write the BSPDE (5) as

$\displaystyle dp_{t}(x)=-\left\\{\eta_{a}\Delta p_{t}(x)-\beta_{a}\nabla
p_{t}(x)+\left(1+\alpha_{a}p_{t}(x)+\sigma_{a}q_{t}(x)\right)\right\\}dt+q_{t}(x)dW^{a}_{t},\quad\forall(t,x)\in(0,T)\times(0,L).$
(7)

The following result is proved in Appendix A.1:

######  Theorem 1.

Let $\hat{Z}^{a}\in\mathcal{A}^{a}$, and denote by $\hat{u}$ the associated
solution to SPDE (3). Assume there exists a solution $(\hat{p},\hat{q})$ to
the BSPDE (7) with boundary conditions (6). Assume also that:

  1. (i)

For all $t\in[0,T)$ and $x\in(0,L)$,

$H(x,\hat{u}_{t}(x),\hat{Z}^{a}_{t}(x),\hat{p}_{t}(x),\hat{q}_{t}(x))=\underset{z\in[0,+\infty)}{\sup}\
H(x,\hat{u}_{t}(x),z,\hat{p}_{t}(x),\hat{q}_{t}(x));$

  2. (ii)

For all $t\in[0,T]$ and $x\in[0,L]$,

$\hat{p}_{t}(x)\geq 0;$

  3. (iii)

For all $Z^{a}\in\mathcal{A}^{a}$,

$\mathbb{E}\left[\int_{0}^{L}\int_{0}^{T}\left(u^{Z^{a}}_{t}(x)-\hat{u}_{t}(x)\right)^{2}\hat{q}_{t}(x)^{2}dtdx\right]<+\infty,$

where $u^{Z^{a}}$ is the solution to (3) controlled by $Z^{a}$;

  4. (iv)

For all $Z^{a}\in\mathcal{A}^{a}$,

$\mathbb{E}\left[\int_{0}^{L}\int_{0}^{T}\left(u^{Z^{a}}_{t}(x)\right)^{2}\hat{p}_{t}(x)^{2}dtdx\right]<+\infty;$

then $\hat{Z}^{a}$ is an optimal control for the control problem (4).

As the boundary conditions (6) are deterministic, we look for a solution to
the BSPDE (7) such that $q(t,x)=0$ for all $(t,x)\in[0,T]\times[0,L]$. In that
case, $p$ solves the following linear parabolic PDE:999From now on, we write
$p(t,x)$ instead of $p_{t}(x)$ to highlight the fact that $p$ is
deterministic.

$\displaystyle\partial_{t}p(t,x)+\eta_{a}\Delta p(t,x)-\beta_{a}\nabla
p(t,x)+\alpha_{a}p(t,x)+1=0,$ (8)

with boundary conditions:

$\displaystyle p(T,x)=0\ \forall x\in(0,L)$ (9) $\displaystyle p(t,x)=0\
\forall(t,x)\in(0,T)\times\\{0,L\\}.$

This is a classical boundary value problem, and the solution admits a Feynman-
Kac representation of the following form:

$\displaystyle
p(t,x)=\mathbb{E}\left[\int_{t}^{T\wedge\gamma^{t,x}}e^{\alpha_{a}(s-t)}ds\right],$
(10)

where we consider a new Brownian motion $\widehat{W}$ and the following SDE:

$\displaystyle dX_{s}=-\beta_{a}ds+\sqrt{2\eta_{a}}d\widehat{W}_{s},$ (11)

and $\gamma^{t,x}$ denotes the stopping time

$\gamma^{t,x}:=\inf\left\\{s\geq
t\left|X^{t,x}_{s}\not\in(0,L)\right.\right\\},$

where $(X^{t,x}_{s})_{s\geq 0}$ is the solution to SDE (11) starting at time
$t$ and at point $x$.

As $H$ is concave in $z$, the optimal control is then given by the solution
to:

$\displaystyle\partial_{z}H\left(x,u_{t}(x),z,p(t,x),q(t,x)\right)=0$
$\displaystyle\iff$
$\displaystyle-\partial_{z}g(x,z)+p(t,x)\partial_{z}f^{a}(x,z)=0,$

if such a solution exists.

######  Remark 3.

If the above equation has a solution for all $t$ and $x$, it is clear that the
corresponding optimal control satisfies condition $(i)$ in Theorem 1.
Furthermore, notice from Equation (10) that $p(t,x)\geq 0$ for all $t$ and
$x$, hence condition $(ii)$ is also satisfied. Note that we obtain here a
solution such that $q(t,x)=0$ for all $t$ and $x$, and therefore condition
$(iii)$ is also naturally verified. Finally, as $p$ is clearly bounded from
Equation (10), Theorem 6.7 in [13] guarantees that condition $(iv)$ is
satisfied.

##  4 Long-term behaviour and closed-form formula

###  4.1 Asymptotics of $p$

In order to study the asymptotic behaviour of the function $p$, we consider
the stationary equation:

$\displaystyle\eta_{a}\Delta\bar{p}(x)-\beta_{a}\nabla\bar{p}(x)+\alpha_{a}\bar{p}(x)+1=0,$
(12)

with boundary conditions:

$\displaystyle\bar{p}(x)=0\ \forall x\in(0,L).$ (13)

The Feynman-Kac representation of the solution is given by:

$\displaystyle\bar{p}(x)=\mathbb{E}\left[\int_{0}^{\gamma^{0,x}}e^{\alpha_{a}s}ds\right].$
(14)

The following proposition is straightforward by dominated convergence:

######  Proposition 1.

For all $x\in[0,L]$, we have

$p(t,x)\underset{T\rightarrow+\infty}{\longrightarrow}\bar{p}(x).$

Moreover, notice that Equation (12) is a simple second order linear ODE, and
together with conditions (13), its solution is given explicitly by:

$\displaystyle\bar{p}(x)=-\frac{1}{\alpha_{a}}+\mu_{+}e^{\nu_{+}x}+\mu_{-}e^{\nu_{-}x}$
(15)

for all $x\in[0,L]$, with

$\nu_{+}=\frac{\beta_{a}+\sqrt{\beta_{a}^{2}-4\eta_{a}\alpha_{a}}}{2\eta_{a}}\qquad\text{and}\qquad\nu_{-}=\frac{\beta_{a}-\sqrt{\beta_{a}^{2}-4\eta_{a}\alpha_{a}}}{2\eta_{a}},$

and

$\mu_{+}=\frac{1}{\alpha_{a}}-\mu_{-}\qquad\text{and}\qquad\mu_{-}=\frac{1}{\alpha_{a}}\frac{e^{\nu_{+}L}-1}{e^{\nu_{+}L}-e^{\nu_{-}L}}.$

###  4.2 An example with closed-form solution

We propose a penalty function $g$ and an intensity function $f^{a}$ with the
following forms:

$g(x,z)=\bar{A}ze^{\bar{a}x},$

and

$f^{a}(x,z)=\lambda z^{r}e^{-\kappa x}+\lambda_{0}(x)$

for all $(x,z)\in(0,L)\times[0,+\infty)$, where
$\bar{A},\bar{a},\lambda,\kappa>0$, $r\in(0,1)$, and $\lambda_{0}$ is a non-
negative function representing the rate of sell order submissions at a
distance $x$ from the mid-price in the absence of incentives.

######  Remark 4.

These two functions deserve a few comments:

  * •

For a fixed $x\in(0,L)$, the intensity function $f^{a}(x,.)$ has both a
constant part, $\lambda_{0}(x)$, corresponding to the intensity of orders
arrival in the absence of incentives (when
$\left(Z^{a}_{t}(x)\right)_{t\in[0,T]}$ is identically equal to $0$), and a
part that depends on the incentive: $\lambda z^{r}e^{-\kappa x}$.This second
part is exponentially decreasing in $x$ to model the fact that the incentive
has less and less impact as we move away from the mid-price: indeed, the
incentive is only perceived by the liquidity provider if her order is
executed, which is more likely when the order is close to the mid-price. It is
concave in $z$ to model the decrease of the marginal liquidity gain derived
from each additional unit of incentive.

  * •

It seems natural to choose a penalty function $g$ that is linear in $z$, as
$g$ should intuitively be proportional to the cost incurred for the platform.
It is increasing in $x$ to model the fact that the exchange is more inclined
to pay to increase liquidity at the first limits than far away from the mid-
price. Here, we take the example of an exponential increase, but one could
reasonably choose to work with any function $g(x,z)=zh(x)$ for any $C^{1}$
increasing function $h$ (the choice of $h$ depends on how much the platform
wishes to concentrate liquidity on the first limits).

######  Proposition 2.

In this framework, the optimal control $\hat{Z}^{a}\in\mathcal{A}^{a}$ is
given by

$\hat{Z}^{a}_{t}(x)=\left(\frac{p(t,x)\lambda
r}{\bar{A}e^{(\bar{a}+\kappa)x}}\right)^{\frac{1}{1-r}}.$

When $T\rightarrow+\infty$, this control reaches a stationary state given in
closed-form by

$\displaystyle\hat{Z}^{a}(x)=\left(\frac{\bar{p}(x)\lambda
r}{\bar{A}e^{(\bar{a}+\kappa)x}}\right)^{\frac{1}{1-r}}.$ (16)

##  5 Numerical simulations

###  5.1 Model parameters

In this section, we apply our model to the case of a (fictitious) symmetric
limit order book with a tick size given by $\delta=0.01\ \$$, with the
following parameters:

  * •

Diffusion parameters: $\eta_{b}=\eta_{a}=10^{-3}\ \$^{2}\cdot\min^{-1}$.

  * •

Convection parameters: $\beta_{b}=\beta_{a}=2\cdot 10^{-2}\ \$\cdot\min^{-1}$.

  * •

Cancellation parameters: $\alpha_{b}=\alpha_{a}=-0.2\min^{-1}$.

  * •

Noise parameters: $\sigma_{b}=\sigma_{a}=0.3\min^{-\frac{1}{2}}$.

  * •

Correlation between $W^{a}$ and $W^{b}$: $\rho=-0.05$.

  * •

Rate of submissions of the form

$f^{a}(x,z)=f^{b}(-x,z)=\lambda\sqrt{z}e^{-\kappa
x}+\lambda_{0}e^{-\kappa_{0}x},$

with:

    * $\lambda_{0}=50000\ \$^{-1}\cdot\min^{-1},$

    * $\kappa_{0}=50\ \$^{-1},$

    * $\lambda=630000\ \$^{-\frac{3}{2}}\cdot\min^{-1},$

    * $\kappa=100\ \$^{-1}.$

This means that, in the absence of incentives, unit limit orders arrive at a
distance of $0.01\ \$$ from the mid-price at a constant rate of
$\delta\lambda_{0}e^{-0.01\kappa_{0}}\simeq 303\min^{-1},$ and at a distance
of $0.05\ \$$ from the mid-price at a constant rate of
$\delta\lambda_{0}e^{-0.05\kappa_{0}}=41\min^{-1}.$ An incentive of $z=0.01\
\$$ per unit order increases the rate at a distance of $0.01\ \$$ from the
mid-price by $\delta\lambda\sqrt{z}e^{-0.01\kappa}\simeq 232\min^{-1},$ and at
a distance of $0.05\ \$$ from the mid-price by
$\delta\lambda\sqrt{z}e^{-0.05\kappa}=4\min^{-1}.$

Regarding the objective function we consider the following:

  * •

Boundary on the LOB: $L=0.11\ \$$ (this corresponds to $10$ limits on each
side).

  * •

Time horizon: $T=30\min$. This horizon ensures convergence toward the
stationary incentives at time $t=0$.

  * •

Penalty function of the form $g(x,z)=\bar{A}ze^{\bar{a}x},$ with:

    * $\bar{A}=4200\ \$^{-2},$

    * $\bar{a}=50\ \$^{-1}.$

######  Remark 5.

The parameters
$\eta_{b},\eta_{a},\beta_{b},\beta_{a},\alpha_{b},\alpha_{a},\sigma_{b},\sigma_{a},\lambda_{0},$
and $\kappa_{0}$ can easily be calibrated on historical data. The parameters
$\lambda$ and $\kappa$ can be calibrated if the exchange already has
historical data with different levels of incentives; otherwise, some tests and
heuristics should be used. Finally, $\bar{A}$ and $\bar{a}$ should be
calibrated depending on the desired shape for the limit order book and the
amount the exchange is ready to spend; simulations can be carried out in order
to find the appropriate values for these parameters.

###  5.2 Optimal incentives and shape of the order book

We compute the asymptotic optimal incentives using Equation (16), and plot
them in Figure 1. As we observe, the optimal incentives decrease at a very
fast rate as we get away from the mid-price. For more clarity, we report these
incentives in Table 1 (the incentives are symmetric at the bid and at the
ask). This fast decrease is expected, as the exchange wants to concentrate
liquidity at the first limits of the order book.

Figure 1: Optimal incentive for a unit order at each limit of the limit order book. Limit | Incentive ($)
---|---
$1st$ | $1.90\cdot 10^{-3}$
$2nd$ | $8.56\cdot 10^{-6}$
$3rd$ | $2.11\cdot 10^{-8}$
$4th$ | $4.01\cdot 10^{-11}$
$5th$ | $6.41\cdot 10^{-14}$
$6th$ | $8.96\cdot 10^{-17}$
$7th$ | $1.09\cdot 10^{-19}$
$8th$ | $1.12\cdot 10^{-22}$
$9th$ | $8.86\cdot 10^{-26}$
$10th$ | $3.84\cdot 10^{-29}$
Table 1: Optimal incentive for a unit order at each limit of the limit order
book.

We then simulate the limit order book given by equation (3) with and without
incentives, to obtain its average shape at time $T$ and compare the two cases.
The shape is plotted in Figure 2. We see that the use of incentives allows to
increase the liquidity level of the order book, particularly at the first
limits. Note that, in the absence of incentives, the first limit is only the
$5$-th more filled of the order book, whereas with incentives it becomes the
$3$-rd more filled.

Figure 2: Shape of the limit order book at time $T$.

###  5.3 Impact of the different parameters

We now play with the parameters of Equation (3) to study their impact on the
optimal incentives and on the limit order book.

####  5.3.1 Diffusion parameter

Let us start with the diffusion parameter. We now divide it by 2, and set
$\eta_{b}=\eta_{a}=5\cdot 10^{-4}\ \$^{2}\cdot\min^{-1}$. The optimal
incentives are now plotted in Figure 3 and reported in Table 2.

We observe that a decrease in the diffusion parameter leads to an increase in
the incentives. Indeed, reducing $\eta_{b}$ and $\eta_{a}$ improves the impact
of the incentives, as the new orders posted are less likely to be replaced by
orders further away from the mid-price. This is confirmed by the shape of the
order book in Figure 4, where we see that the overall level of the order book
is higher than in Figure 2, particularly at the first limits.

Figure 3: Optimal incentive for a unit order at each limit of the limit order
book.

Limit | Incentive ($)
---|---
$1st$ | $3.51\cdot 10^{-2}$
$2nd$ | $1.73\cdot 10^{-4}$
$3rd$ | $4.77\cdot 10^{-7}$
$4th$ | $1.03\cdot 10^{-9}$
$5th$ | $1.90\cdot 10^{-12}$
$6th$ | $3.15\cdot 10^{-15}$
$7th$ | $4.68\cdot 10^{-18}$
$8th$ | $6.02\cdot 10^{-21}$
$9th$ | $6.14\cdot 10^{-24}$
$10th$ | $3.55\cdot 10^{-27}$
Table 2: Optimal incentive for a unit order at each limit of the limit order
book.
Figure 4: Shape of the limit order book at time $T$.

####  5.3.2 Convection rate

We now turn to the convection rate and multiply it by 5, setting
$\beta_{b}=\beta_{a}=0.1\ \$\cdot\min^{-1}$. The optimal incentives are now
plotted in Figure 5 and reported in Table 3.

Figure 5: Optimal incentive for a unit order at each limit of the limit order book. Limit | Incentive ($)
---|---
$1st$ | $2.21\cdot 10^{-3}$
$2nd$ | $1.18\cdot 10^{-5}$
$3rd$ | $3.53\cdot 10^{-8}$
$4th$ | $8.37\cdot 10^{-11}$
$5th$ | $1.74\cdot 10^{-13}$
$6th$ | $3.29\cdot 10^{-16}$
$7th$ | $5.80\cdot 10^{-19}$
$8th$ | $9.36\cdot 10^{-22}$
$9th$ | $1.28\cdot 10^{-24}$
$10th$ | $1.09\cdot 10^{-27}$
Table 3: Optimal incentive for a unit order at each limit of the limit order
book.

We observe that an increase in the convection rate leads to an increase in the
incentives. Indeed, increasing $\beta_{b}$ and $\beta_{a}$ means that orders
will naturally move towards the first limit, hence increasing the incentives
even at the last limits will have a good impact on liquidity. The overall
level of the order book in Figure 6 is lower than the one in Figure 2, because
orders at the first limit are also more rapidly matched by a market order.

Figure 6: Shape of the limit order book at time $T$.

####  5.3.3 Cancellation rate

Finally, we study the impact of the cancellation rate. We multiply it by 2,
setting $\alpha_{b}=\alpha_{a}=-0.4\min^{-1}$. The optimal incentives are now
plotted in Figure 7 and reported in Table 4.

Figure 7: Optimal incentive for a unit order at each limit of the limit order book. Limit | Incentive ($)
---|---
$1st$ | $1.64\cdot 10^{-2}$
$2nd$ | $7.57\cdot 10^{-5}$
$3rd$ | $1.95\cdot 10^{-7}$
$4th$ | $3.90\cdot 10^{-10}$
$5th$ | $6.69\cdot 10^{-13}$
$6th$ | $1.01\cdot 10^{-15}$
$7th$ | $1.36\cdot 10^{-18}$
$8th$ | $1.57\cdot 10^{-21}$
$9th$ | $1.41\cdot 10^{-24}$
$10th$ | $7.07\cdot 10^{-28}$
Table 4: Optimal incentive for a unit order at each limit of the limit order
book.

We observe that an increase in the convection rate also leads to an increase
in the incentives. Indeed, when $\alpha_{b}$ and $\alpha_{a}$ increase, the
exchange needs to compensate the higher cancellation rate by setting more
attractive incentives. This allows to get roughly the same level and shape for
the order book with incentives in the case with high cancellation rate plotted
in Figure 8 as in the first case in Figure 2.

Figure 8: Shape of the limit order book at time $T$.

##  6 Conclusion

In this paper, we tackled the problem of incentives design in a limit order
book. We represented this problem as an optimal control problem for a SPDE,
using the limit order book model of Cont and Müller [12]. Due to the specific
form of the control problem faced by the exchange, we showed that the solution
can be expressed with a Feynman-Kac representation. Moreover, we gave an
example of a class of intensity and penalty functions under which the
asymptotic behavior of the control problem leads to closed-form incentives at
each limit of the order book. We then simulated a fictitious limit order book
using these incentives to illustrate numerically their influence on its shape.
In particular, we showed that these incentives increase liquidity at the first
limits of the order book. Finally, we performed a sensitivity analysis with
respect to the market parameters.

##  Appendix A Appendix

###  A.1 Proof of Theorem 1

We follow closely the proof of Theorem 10.2 in [24].

Let $Z^{a}\in\mathcal{A}^{a}$, and for the sake of simplicity, denote by $u$
the associated solution to SPDE (3) : $u:=u^{Z^{a}}$. For
$(t,x)\in[0,T)\times(0,L)$, we introduce

$\displaystyle\hat{H}_{t}(x)$
$\displaystyle:=H\left(x,\hat{u}_{t}(x),\hat{Z}^{a}_{t}(x),\hat{p}_{t}(x),\hat{q}_{t}(x)\right)$
$\displaystyle=\hat{u}_{t}(x)-g(x,\hat{Z}^{a}_{t}(x))+\left(\alpha_{a}\hat{u}_{t}(x)+f^{a}(x,\hat{Z}^{a}_{t}(x))\right)\hat{p}_{t}(x)+\sigma_{a}\hat{u}_{t}(x)\hat{q}_{t}(x),$

and

$\displaystyle H_{t}(x)$
$\displaystyle:=H\left(x,u_{t}(x),Z^{a}_{t}(x),\hat{p}_{t}(x),\hat{q}_{t}(x)\right)$
$\displaystyle=u_{t}(x)-g(x,Z^{a}_{t}(x))+\left(\vphantom{\hat{Z}^{a}_{t}(x)}\alpha_{a}u_{t}(x)+f^{a}(x,Z^{a}_{t}(x))\right)\hat{p}_{t}(x)+\sigma_{a}u_{t}(x)\hat{q}_{t}(x).$

Hence we have

$J(\hat{Z}^{a})-J(Z^{a})=I,$

where

$\displaystyle I$
$\displaystyle:=\mathbb{E}\left[\int_{0}^{T}\int_{0}^{L}\left\\{\hat{u}_{t}(x)-u_{t}(x)-\left(g(x,\hat{Z}^{a}_{t}(x))-g(x,Z^{a}_{t}(x))\right)\right\\}dxdt\right]$
$\displaystyle=\mathbb{E}\\!\bigg{[}\int_{0}^{T}\\!\\!\\!\int_{0}^{L}\\!\\!\Big{\\{}\hat{H}_{t}(x)\\!-\\!H_{t}(x)\\!-\\!\big{(}\alpha_{a}(\hat{u}_{t}(x)-u_{t}(x))+f^{a}(x,\hat{Z}^{a}_{t}(x))-f^{a}(x,Z^{a}_{t}(x))\big{)}\hat{p}_{t}(x)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\sigma_{a}(\hat{u}_{t}(x)\\!-\\!u_{t}(x))\hat{q}_{t}(x)\Big{\\}}dxdt\bigg{]}.$

Moreover, denoting by $\tilde{u}_{t}(x)$ the process

$\tilde{u}_{t}(x):=u_{t}(x)-\hat{u}_{t}(x),$

we know from the boundary condition (6) that

$\displaystyle 0$
$\displaystyle=\mathbb{E}\left[\int_{0}^{L}\hat{p}_{T}(x)\tilde{u}_{T}(x)dx\right]$
$\displaystyle=\mathbb{E}\left[\int_{0}^{L}\left(\hat{p}_{0}(x)\tilde{u}_{0}(x)+\int_{0}^{T}\left\\{\tilde{u}_{t}(x)d\hat{p}_{t}(x)+\hat{p}_{t}(x)d\tilde{u}_{t}(x)+\sigma_{a}\tilde{u}_{t}(x)\hat{q}_{t}(x)dt\right\\}\right)dx\right]$
$\displaystyle=\mathbb{E}\Bigg{[}\int_{0}^{L}\int_{0}^{T}\bigg{\\{}\tilde{u}_{t}(x)\left(-\mathcal{L}^{*}\hat{p}_{t}(x)-\widehat{\partial_{u}H}_{t}(x)\right)$
$\displaystyle\qquad+\hat{p}_{t}(x)\left(\mathcal{L}\tilde{u}_{t}(x)+\left(\alpha_{a}\tilde{u}_{t}(x)+f^{a}(x,Z^{a}_{t}(x))-f^{a}(x,\hat{Z}^{a}_{t}(x))\right)\right)+\sigma_{a}\tilde{u}_{t}(x)\hat{q}_{t}(x)\bigg{\\}}dtdx\Bigg{]},$

where

$\widehat{\partial_{u}H}_{t}(x)=\partial_{u}H\left(x,\hat{u}_{t}(x)\hat{Z}^{a}_{t}(x),\hat{p}_{t}(x),\hat{q}_{t}(x)\right)\quad\forall(t,x)\in[0,T)\times(0,L).$

Hence

$\displaystyle J(\hat{Z}^{a})-J(Z^{a})$
$\displaystyle=I-\mathbb{E}\left[\int_{0}^{L}\hat{p}_{T}(x)\tilde{u}_{T}(x)dx\right]$
$\displaystyle=\mathbb{E}\Bigg{[}\int_{0}^{T}\left(\int_{0}^{L}\left\\{\tilde{u}_{t}(x)\mathcal{L}^{*}\hat{p}_{t}(x)-\hat{p}_{t}(x)\mathcal{L}\tilde{u}_{t}(x)\right\\}dx\right)dt$
$\displaystyle\qquad+\int_{0}^{L}\left(\int_{0}^{T}\left\\{\hat{H}_{t}(x)-H_{t}(x)+\widehat{\partial_{u}H}_{t}(x)\tilde{u}_{t}(x)\right\\}dt\right)dx\Bigg{]}.$

As $\tilde{u}_{t}(x)=\hat{p}_{t}(x)=0$ $\forall(t,x)\in[0,T)\times\\{0,L\\}$
a.s, we have a.s.

$\int_{0}^{L}\tilde{u}_{t}(x)\mathcal{L}^{*}\hat{p}_{t}(x)dx=\int_{0}^{L}\hat{p}_{t}(x)\mathcal{L}\tilde{u}_{t}(x)dx,$

and therefore

$J(\hat{Z}^{a})-J(Z^{a})=\mathbb{E}\Bigg{[}\int_{0}^{L}\left(\int_{0}^{T}\left\\{\hat{H}_{t}(x)-H_{t}(x)+\widehat{\partial_{u}H}_{t}(x)\tilde{u}_{t}(x)\right\\}dt\right)dx\Bigg{]}.$

Finally by Lemma 1, we have

$H_{t}(x)-\hat{H}_{t}(x)\leq\widehat{\partial_{u}H}_{t}(x)\tilde{u}_{t}(x)+\widehat{\partial_{z}H}_{t}(x)\left(Z^{a}_{t}(x)-\hat{Z}^{a}_{t}(x)\right),$

where

$\widehat{\partial_{z}H}_{t}(x)=\partial_{z}H\left(x,\hat{u}_{t}(x)\hat{Z}^{a}_{t}(x),\hat{p}_{t}(x),\hat{q}_{t}(x)\right)\quad\forall(t,x)\in[0,T)\times(0,L).$

But by $(i),$

$\widehat{\partial_{z}H}_{t}(x)\left(Z^{a}_{t}(x)-\hat{Z}^{a}_{t}(x)\right)\leq
0,$

and therefore

$\hat{H}_{t}(x)-H_{t}(x)+\widehat{\partial_{u}H}_{t}(x)\tilde{u}_{t}(x)\geq
0,$

and

$J(\hat{Z}^{a})\geq J(Z^{a}).$

## References

  * Angel et al. [2011] J.J. Angel, L.E. Harris, and C.S. Spatt.  Equity trading in the 21st century.  _The Quarterly Journal of Finance_ , 1(01):1–53, 2011.
  * Arnuk and Saluzzi [2010] S. Arnuk and J. Saluzzi.  Concept release on equity market structure.  Technical report, The Federal Register, 2010.
  * Avellaneda and Stoikov [2008] M. Avellaneda and S. Stoikov.  High-frequency trading in a limit order book.  _Quantitative Finance_ , 8(3):217–224, 2008.
  * Baldacci et al. [2019a] B. Baldacci, I. Manziuk, T. Mastrolia, and M. Rosenbaum.  Market making and incentives design in the presence of a dark pool: a deep reinforcement learning approach.  _ArXiv preprint arXiv:1912.01129_ , 2019a.
  * Baldacci et al. [2019b] B. Baldacci, D. Possamaï, and M. Rosenbaum.  Optimal make take fees in a multi market maker environment.  _ArXiv preprint arXiv:1907.11053_ , 2019b.
  * Baldacci et al. [2019c] Bastien Baldacci, Paul Jusselin, and Mathieu Rosenbaum.  How to design a derivatives market?  _arXiv preprint arXiv:1909.09257_ , 2019c.
  * Baldacci et al. [2020] Bastien Baldacci, Philippe Bergault, Joffrey Derchu, and Mathieu Rosenbaum.  On bid and ask side-specific tick sizes.  _arXiv preprint arXiv:2005.14126_ , 2020.
  * Bellia [2017] M. Bellia.  High-frequency market making: liquidity provision, adverse selection, and competition.  _SSRN preprint 3074313_ , 2017.
  * Bergault and Guéant [2021] Philippe Bergault and Olivier Guéant.  Size matters for otc market makers: general results and dimensionality reduction techniques.  _Mathematical Finance_ , 31(1):279–322, 2021\.
  * Brolley and Malinova [2013] M. Brolley and K. Malinova.  Informed trading and maker-taker fees in a low-latency limit order market.  _SSRN preprint 2178102_ , 2013.
  * Colliard and Foucault [2012] J.-E. Colliard and T. Foucault.  Trading fees and efficiency in limit order markets.  _The Review of Financial Studies_ , 25(11):3389–3421, 2012.
  * Cont and Müller [2021] Rama Cont and Marvin S Müller.  A stochastic partial differential equation model for limit order book dynamics.  _SIAM Journal on Financial Mathematics_ , 12(2):744–787, 2021.
  * Da Prato and Zabczyk [2014] G. Da Prato and J. Zabczyk.  _Stochastic equations in infinite dimensions_ , volume 45 of _Encyclopedia of mathematics and its applications_.  Cambridge University Press, 2014.
  * Dayri and Rosenbaum [2015] K. Dayri and M. Rosenbaum.  Large tick assets: implicit spread and optimal tick size.  _Market Microstructure and Liquidity_ , 1(01):1550003, 2015.
  * El Euch et al. [2021] O. El Euch, T. Mastrolia, M. Rosenbaum, and N. Touzi.  Optimal make-take fees for market making regulation.  _Mathematical Finance_ , 31:109–148, 2021.
  * Élie and Possamaï [2019] R. Élie and D. Possamaï.  Contracting theory with competitive interacting agents.  _SIAM Journal on Control and Optimization_ , 57(2):1157–1188, 2019.
  * Élie et al. [2019] R. Élie, T. Mastrolia, and D. Possamaï.  A tale of a principal and many many agents.  _Mathematics of Operations Research_ , 44(2):440–467, 2019.
  * Guéant [2016] O. Guéant.  _The financial mathematics of market liquidity: from optimal execution to market making_.  Financial Mathematics. Chapman and Hall/CRC, 2016.
  * Harris et al. [1982] M. Harris, C.H. Kriebel, and A. Raviv.  Asymmetric information, incentives and intrafirm resource allocation.  _Management Science_ , 28(6):604–620, 1982.
  * Laruelle and Lehalle [2018] S. Laruelle and C.-A. Lehalle.  _Market microstructure in practice_.  World Scientific, 2nd edition, 2018.
  * Megarbane et al. [2017] N. Megarbane, P. Saliba, C.-A. Lehalle, and M. Rosenbaum.  The behavior of high-frequency traders under different market stress scenarios.  _Market Microstructure and Liquidity_ , 3(03n04):1850005, 2017.
  * Menkveld [2013] A.J. Menkveld.  High frequency trading and the new market makers.  _Journal of Financial Markets_ , 16(4):712–740, 2013.
  * Mirrlees and Raimondo [2013] J.A. Mirrlees and R.C. Raimondo.  Strategies in the principal-agent model.  _Economic Theory_ , 53(3):605–656, 2013.
  * Øksendal and Sulem [2007] B. Øksendal and A. Sulem.  _Applied stochastic control of jump diffusions_.  Universitext. Springer-Verlag Berlin Heidelberg, 2nd edition, 2007.
  * Sannikov [2008] Y. Sannikov.  A continuous-time version of the principal-agent problem.  _The Review of Economic Studies_ , 75(3):957–984, 2008.