# Expected geoneutrino signal at JUNO using local integrated 3-D refined
crustal model

Ran HAN corresponding author<EMAIL_ADDRESS>Beijing Institute of
Spacecraft Environment Engineering, Beijing 100094, China.  ZhiWei LI
corresponding author<EMAIL_ADDRESS>Beijing Institute of Spacecraft
Environment Engineering, Beijing 100094, China.  Ruohan GAO corresponding
author<EMAIL_ADDRESS>China University of Geosciences (Beijing),
Beijing, 100083, China.  Yao SUN  China University of Geosciences (Beijing),
Beijing, 100083, China.  Ya XU  Key Laboratory of Petroleum Resource Research,
Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing,
100029, China.  Yaping CHENG  Beijing Institute of Spacecraft Environment
Engineering, Beijing 100094, China.  Guangzheng JIANG  State Key Laboratory
for Oil and Gas Reservoir Geology and Exploitation, Chengdu University of
Technology, Chengdu, 610059, China.  Jie PANG  School of Nuclear Science &
Engineering, North China Electric Power University, Beijing 102206, China.
Fengcheng LIU  School of Nuclear Science & Engineering, North China Electric
Power University, Beijing 102206, China.  Andong WANG  East China University
of Technology, Nanchang, 330013, China.  Yufei XI  Institute of Hydrogeology
and Environmental Geology, Chinese Academy of Geological Sciences,
Shijiazhuang, 050061, China.  Liangjian WEN  Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing, 100049, China.  Jun CAO  Institute of
High Energy Physics, Chinese Academy of Sciences, Beijing, 100049, China.
YuFeng LI  Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing, 100049, China.

###### Abstract

Geoneutrinos are a unique tool that brings to the surface information about
our planet, in particular, its radiogenic power, which, however, can also
provide insight into its formation and chemical composition. To date, only the
KamLAND and Borexino experiments observed geoneutrinos, with the former
characterized by low concentration of heat-producing elements in the Earth in
contrast to the latter that sets tight upper limits on the power of a
georeactor hypothesized to be present inside the Earth. With respect to the
results yielded therefrom, a small discrepancy has been identified. On this
account, next generation experiments like JUNO are needed if it is to provide
definitive results with respect to the Earth’s radiogenic power, and to fully
exploit geoneutrinos to better understand the natural radioactivity of our
planet.

An accurate a priori prediction of the crustal contribution plays an important
role in enabling the translation of a particle physics measurement into above
geo-scientific questions. The existing GIGJ model of JUNO only focused on
constructing a geophysical model of the local crust, without local geochemical
data. Another existing JULOC includes both data, but only able to be achieved
for the top layer of the upper crust, not in deep vertical.

This paper reports on the development of JUNO’s first 3-D integrated crustal
model, JULOC-I, which combines seismic, gravity, rock sample data and thermal
flow data with new building method, solved the problem of geochemical
distribution in vertical depth. JULOC-I estimates the total lithosphere
geoneutrino signal to be $44.4^{+5.1}_{-4.4}$ TNU. Results shows higher than
expected geoneutrino signals are mainly attributable to higher uranium (U) and
thorium (Th) in southern China than that found elsewhere on Earth. Moreover,
the high level of accuracy of the JULOC-I model, complemented by 10 years of
experimental data, indicates that JUNO has an opportunity to test different
mantle models. Predictions by JULOC-I can be tested after JUNO goes online and
higher accuracy local crustal model continue to play an important role to
improve mantle measurements precision.

##  1 Introduction

Collaboration between geologists and physicists over the past 150 years or so
has yielded many new scientific disciplines. However, the 150-year-old
conundrum posed by Lord Kelvin about the origin of the Earth and its power
remains unresolved. An understanding of the Earth’s energy budget requires an
assessment of the relative contributions of mantle cooling and radiogenic heat
production, and calculation of the rate of heat dissipation. Geothermal, as
the heat generated within the Earth, promotes the development and evolution of
Earth in entirety and also controls the activity of the crust and mantle,
formation of volcanoes, magmatic activity, earthquakes, and plate tectonics
[1, 2]. Among the Earth’s energy sources, mantle energy is extremely
important. Understanding the process by which Earth generates energy, as well
as the present condition of Earth’s energy sources, is essential information
for all fields of the Earth sciences.

Seismology has revealed the physical properties of the Earth, showing the
crust-mantle-core layer structure, Meanwhile, cosmochemistry has successfully
modeled the bulk chemical composition of the Earth. In combination with the
physical properties of the Earth, cosmochemistry has facilitated the
determination of the crystal structure and chemical constitution of rocks in
each layer of the Earth. However, since the deepest part of the Earth like
mantle is essentially unreachable, no direct observation of the chemical
composition of the whole Earth can be made.

Fortunately, geoneutrinos, the antineutrinos emitted from nuclei undergoing
$\beta^{-}$ decay of 238U, 232Th and 40K inside the Earth, can provide
information on the concentrations of heat-producing elements of the whole
planet. The flux of geoneutrinos at any point on the Earth is a function of
the abundance and distribution of radioactive elements within our planet. This
flux has been successfully detected by the 1-kt Kamioka Liquid Scintillator
Antineutrino Detector(KamLAND) and 0.3-kt Borexino detectors[4, 6, 7, 8, 9,
10], though their measurements have been limited by their low statistics.
Jiangmen Underground Neutrino Observatory(JUNO) will join the family of
geoneutrino experiments and will provide an opportunity to measure geoneutrino
more accurately [11, 12, 13, 14, 15]. Within the first year of its operation,
JUNO will record more geoneutrino events than all other detectors will have
accumulated to that point in time. Given sufficient statistics, JUNO will be
able to separate the relative contributions of U and Th.

Geoneutrino signal prediction by JUNO is of great significance to the
translation of particle physics measurements into meaningful geoscientific
ones. Comparisons of predictions and real data will not only verify different
crustal models, but also disentangle the mantle signal and U/Th ratio. The
first published JUNO sensitivity studies[12, 40] in relation to the
composition of the Earth[22] globally showed an expected geoneutrino signal of
$39.7^{+6.5}_{-5.2}$ TNU (Terrestrial Neutrino Units), whereby a TNU
corresponds to one event per $10^{32}$ free protons per year at 100%
efficiency. In this model, a crustal model based on reflection and refraction
data “CRUST 2.0” [17, 18], a global shear-velocity model of the crust and
upper mantle “CUB 2.0” [19] and a high-resolution map of the depth of the Moho
(crust-mantle boundary) based on gravity field data-’GEMMA’ [20, 21].

To test different compositional models of the mantle, an accurate a priori
estimate of the crustal geoneutrino flux is required for the purpose of
reducing the uncertainty of prediction. However, the regional crust, i.e. the
crust within 500 km of JUNO, contributes more than 50% to the total
signal[23]. So, local refined geological models are needed for more accurate
estimation of the crustal signal. Local refined geological models of the South
China Block have recently been published with Reguzzoni et al. [14] outlining
a GOCE Inversion for Geoneutrinos at JUNO(GIGJ) model, a 3-D crustal model
built by inverting GOCE gravimetric data for the $6^{\circ}\times 4^{\circ}$
area centered on the location of JUNO experiment. This model only focused on
constructing a geophysical and gravimetric model of the local crust, the heat-
producing element(HPE) abundances were still based on the homogenous global
model. Meanwhile, Gao et al. [15] built the first refined 3-D JUno LOCal
model(JULOC). JULOC is the local crust model of the closest $10^{\circ}\times
10^{\circ}$ grid surrounding JUNO, using local geological, geochemical, and
geophysical data. The JULOC model has provided a way to build a crustal
composition model cell-by-cell to lessen the uncertainty in estimating
geoneutrino flux in the crust. However, due to limitations with respect to
geophysical and geochemical data, this was only able to be achieved for the
top layer of the upper crust, with homogenous local data being used for the
remaining basement layer of the upper crust. Using such a cell-by-cell
composition model has enabled the uncertainty of the top layer to be reduced
to 5%. In contrast, the U and Th abundances of the whole basement layer have
been estimated based on the composition of low grade metamorphic rocks and
considered to be uniform. Thus, the errors of different cells are correlated,
resulting in an uncertainty as high as 45%. In the JULOC model, not only is
there a large difference in uncertainty between the top layer and base layers
of the upper crust, but also the boundary of the top layer and the basement
layer have been arbitrarily assigned to the uppermost 5-kilometers. To solve
these issues with JULOC, an integrated updated JULOC crustal model –
JULOC-I(JUno LOCal - Integrate model) has been developed following the
collection of more sample data and the design of a new building method.

This paper introduces JULOC-I model, including its data sources, building
method and expected geoneutrino signal and associated uncertainty. The
expected geoneutrino signals of the four crustal models are also summarized by
combining the Global, JULOC, JULOC-I and GIGJ models, with the potential of
JUNO measurements being used to separate these models.

##  2 JULOC-I: An integrated local 3-D refined crustal model

This paper introduces an updated version of the JULOC model, JULOC-I, which
offers improvements in data and methodology compared to its predecessor(see
Fig.1).

  * 1)

With respect to the geophysics model, the JULOC model used 450 permanent
broadband seismic stations in South China, with the JULOC-I model using an
additional 150 seismic stations. In terms of model construction, the JULOC
model only used seismic ambient noise tomography methodology, whereas the
JULOC-I model uses joint tomographic inversion methodology with seismic
surface wave dispersion and gravity data for obtaining more reliable crustal
3-D shear wave structures(S wave).

  * 2)

With respect to the geochemical model, JULOC-I used data sourced from
literature reviews with a sampling survey subsequently undertaken to complete
information for cells that had missing data or missing rock type. A different
model building method was also used in JULOC-I. To ascertain the proportion of
different rocks in each cell, JULOC used a geological map to determine rock
distribution; however,this method only effective for the surface, and was very
arbitrary with respect to depth. In contrast, JULOC-I uses the crustal 3-D S
wave structures determined in step 1, and then combines a priori probability
of lithology to produce a posteriori probability lithology map for each cell
for every 1 km of depth. The U and Th concentrations of each rock type plus
the probability lithology map, provide a 3-D composition probabilistic map.

  * 3)

The probability map and U/Th abundance derived from the geochemical model are
then constrained by 182 thermal flow data points from around JUNO to produce
the final probability map for geoneutrino flux calculation.The U, Th and K
data for each layer, along with its heat production and one-dimensional heat
transfer function, are used to test the reliability of the geochemical model.

Figure 1: Flow chart of how the JUNO-I model was built. The blue one
components represent the data using, while those in red one pertain to the
model we built.

###  2.1 Research area and sampling survey

The JUNO detector is located on the coastline of South China
($N22.12^{\circ},E112.52^{\circ}$). The JULOC-I model is based on the same
research area as JULOC, focuses on the $10^{\circ}\times 10^{\circ}$ grid
immediately surrounding the detector, i.e. $N17^{\circ}\sim
27^{\circ},E108^{\circ}\sim 118^{\circ}$. The research area is mainly
comprised of the South China Block to the northwest and the South China Sea to
the southeast and has two parts – the northern continental margin and a small
portion of the oceanic basin. The oceanic basin is widely accepted to have
been formed by Paleogene continental rifting and subsequent oceanic floor
spreading during the Oligo-Miocene age (e.g [24, 25]). The oceanic data were
mainly sourced from literature due to limitations posed by challenging survey
conditions. With respect to the northern continental margin, data from the
literature were supplemented by six sampling surveys during which more 1,002
rock samples were collected.

A goal of the sampling surveys was to collect enough data samples in each cell
of the research area. In the JULOC model, a cell-by-cell calculation method
was used to reduce the error in the Upper Crust (UC) layer, which is the main
contributor of geoneutrino signals. However,the accuracy of this methodology
depends upon having enough samples in each cell to reduce the associated error
between different cells and thereby deliver a more precise estimate of
geoneutrino flux [26, 15]. So,to achieve this goal, the sampling range within
a 200-kilometer radius of the JUNO center was divided into $0.4^{\circ}\times
0.4^{\circ}$,whereas in the 200-500-kilometer range, the sampling cell size
was $1^{\circ}\times 1^{\circ}$. In each cell, 10 to 20 representative
lithology samples were collected, though this step was skipped if sufficient
data could be sourced from the literature.

Another sampling survey goal was to collect the lithology samples that the
JULOC model lacked.The continental crust around JUNO is mainly composed of
crystalline basement including meta-sedimentary and meta-igneous rocks,
Precambrian to Quaternary sedimentary rocks with minor volcanic rocks, and
granitic intrusions. The survey included the following upper crust samples :
Yanshanian granite (especially the rock mass in Guangdong province),
Caledonian granite (Guangning County and Deqing County, Zhaoqing City,
Guangdong province), Variscan period granite (Pubei County, Qinzhou City,
Guangxi Province), Carboniferous-Devonian limestone (in Guangxi and Guangdong
provinces), Quaternary sediments, Permian-Triassic limestone and clastic rocks
(in Guangxi province). GPS location and geological data were recorded for each
sample. Most of the samples were collected from fresh outcrops, representative
of the geological formation, and placed in sealed bags. Fig.2 outlines the
locations from where the 1,002 rock samples were collected.

Figure 2: Locations within a 500-kilometer radius of JUNO from where 1002 rock
samples were collected. Every sample was collected from the upper crust layer
of fresh outcrops, and projected onto a geologic map.

###  2.2 Joint tomographic inversion with seismic and gravity data

JULOC was the first JUNO 3-D crustal with seismic ambient noise tomography,
using one-year continuous seismic ambient noise waveform data from 450
permanent broadband seismic stations in South China, with 150 additional
seismic stations being used by JULOC-I.Ambient noise tomography is a popular
method for imaging the 3-D shear wave structure of the Earth’s crust. However,
a single type of geophysical data is sensitive to a few physical parameters in
a specific regions, leading to an incomplete data set. This problem can be
overcome however by combining different geophysical data sets in the inverse
problem, with gravity data a significant part of such a multidata set joint
inversion method, especially in the area of joint inversion seismic and
gravity data. The JULOC-I model uses a newly developed joint tomographic
inversion method of surface wave dispersion and Bouguer gravity anomaly data
to obtain more reliable crustal 3-D shear wave structures[27]. The gravity
data used in this study came from the global satellite Bouguer gravity anomaly
model WGM2012 [28], which is based on the global satellite gravity model
EGM2008, using the DTU10 digital elevation model to determine $1^{\circ}\times
1^{\circ}$ terrain correction results[29](see Fig.3).

Figure 3: Distribution of seismic stations in South China from which data were
collected(red and green triangles represent permanent seismic stations of the
China Earthquake Networks Center,while the yellow triangles represent
temporary seismic stations), and distribution of gravity anomalies in the
study area (WGM2012 sourced data).

With respect to 3-D density, the same methodology as used in JULOC has been
kept, using empirical formulae transfer S-wave velocity to density [31]. Fig.4
shows the 3-D density distributions at 3 km, 9 km,30 km and 50 km depths in
the South China Block.

Figure 4: 3-D density structures in the South China Block, derived from the
3-D velocity model using a gravity-seismology joint inversion method and based
on the empirical relationship between the density and velocity of crustal
materials.

The crustal layering is the same as in the JULOC model, the UC and MC of
CRUST1.0 model were adopted. As for the Moho depth, the model derived from
Bouguer gravity inversion under the constraints of $\leavevmode\nobreak\ 120$
deep seismic sounding profiles is preferred [32], as it provides better
resolution and coverage of the research area than other available models.

###  2.3 3-D U/Th concentrations map

Following the geophysical 3-D crustal model above, here we attempt to
construct a 3-D model for the chemical (U and Th) abundances by combining U/Th
abundance of each rock type with the 3-D lithology proportion map.

More than 500 rocks samples were collected to fill in the data gaps in the
literature(e.g. for sedimentary rocks and carbonate rocks) or to rectify
uneven data point distribution in cells. The samples were tested for trace
elements by means of the Thermo ScientificTM Element XRTM High Resolution ICP-
MS system. The measurement process was undertaken as Sun et.al(2022). Table.1
shows the minimum, median, and maximum values of U and Th concentrations from
those laboratory analyses.

The median values of U (7.48 ppm) and Th (24.3 ppm) for 213 GRA (Granite-
Granodiorite) samples were higher than for other types of rock, with a U
abundance of up to 296 ppm. These results were similar to those found in the
literature at 7.49 ppm and 26.0 ppm, respectively. MGW(siltstone, sandstone,
conglomerate and shale) had low median abundances of 2.75 ppm and 12.9 ppm for
U and Th, respectively.However, there was wide variation for both U and Th, at
0.07 to 39.1 ppm and 0.10 to 120 ppm, respectively. The distribution of U and
Th abundances in the sample is similar to that reported in the
literature,though the reported abundances of U and Th in MGW in the literature
are 2.34 ppm and 9.21 ppm, respectively, both of which are slightly lower than
the observed data. PHY(Phyllite) and QSC(Schist) have similar U-Th median and
variation ranges. Compared with other rock types, SLT(Slate) had the lowest
median U and Th values at 3.79 ppm and 8.65 ppm, respectively.

Table 1: Minimum (min), median (median) and maximum (max) U-Th concentrations of rock types. Rock type | Concentration | Mini | Median | Max | Sample size
---|---|---|---|---|---
GRA | U(ppm) | 0.56 | 7.48 | 296 | 213
| Th(ppm) | 0.99 | 24.3 | 85.5
MGW | U(ppm) | 0.07 | 2.75 | 39.1 | 246
| Th(ppm) | 0.10 | 12.9 | 120
PHY | U(ppm) | 1.70 | 3.30 | 6.08 | 13
| Th(ppm) | 6.26 | 17.2 | 26.3
QSC | U(ppm) | 1.07 | 4.13 | 8.05 | 8
| Th(ppm) | 1.33 | 17.5 | 22.0
SLT | U(ppm) | 1.40 | 3.79 | 33.7 | 19
| Th(ppm) | 0.68 | 8.65 | 41.6

The method of building a 3-D lithology proportion map by means of 3-D S-wave
velocity has previously been used by Takeuchi et al.(2019) [26]. This method
firstly assumes a priori probability map of lithology as $P^{x}_{(i)}$, where
’i’ is the type of rock and ’x’ is the location of that rock. This priori
probability allows local lithological features to be considered. An a
posteriori probability lithology map is then developed by establishing a
mathematical connection between the crustal P-wave velocity $V_{P}^{obs}$ and
rock type proportion, as described in Appendix A. Following the method
described in Appendix B, 3-D PDFs of U and Th concentrations were constructed.
Fig.5 shows the average abundance of U and Th in each cell of the research
area, the depth from top to bottom is 1-km,5-km,10-km,15-km and 20-km
respectively.

Figure 5: Abundance of U and Th in each $0.4^{\circ}\times 0.4^{\circ}$ cell
of the research area,the depth from top to bottom is 1-km,5-km,10-km,15-km and
20-km respectively.

##  3 Geoneutrino signal at JUNO as predicted by JULOC-I

The joint geophysical and geochemical 3-D local crustal models conducted
around JUNO can be used to calculate the geoneutrino flux produced by the
Earth’s crust. However, predicting the geoneutrino signal that could be
detected at the JUNO detector also requires detection information of the
particle by means of inverse beta decay (IBD) reaction. The detailed
mathematics of geoneutrino signal calculation is outlined in R.Gao et.al(2020)
[15], the parameters setting and error estimation method of this paper are
described in Appendix C. Summary of expected geoneutrino signals and
uncertainty (1 standard deviation) in TNU for JUNO from local
crust($10^{\circ}\times 10^{\circ}$ as JULOC-I research area) as reported in
Table 2

Table 2: Predicted geoneutrino signal (S) from the local crust(LOC) in the research area around JUNO and uncertainty ($\sigma$) in TNU, including the contributions from the Continental Crust(includes UC, Middle Crust (MC), Lower Crust (LC) and Oceanic Crust (OC). |  $S_{\text{U}}\pm\sigma$ |  $S_{\text{Th}}\pm\sigma$ |  $S_{\text{U+Th}}\pm\sigma$
---|---|---|---
UC |  $18.6\pm 1.8$ |  $4.7\pm 0.3$ |  $23.3\pm 1.9$
MC |  $5.0^{+4.0}_{-2.2}$ |  $1.5^{+1.8}_{-0.8}$ |  $7.0^{+4.0}_{-2.5}$
LC |  $1.2^{+0.3}_{-0.2}$ |  $0.2\pm 0.1$ |  $1.5\pm 0.3$
Continental Crust |  $25.3^{+3.9}_{-3.4}$ |  $6.5^{+1.2}_{-1.0}$ |  $32.2\pm 3.8$
OC |  $0.2\pm 0.05$ |  $0.1\pm 0.01$ |  $0.3\pm 0.05$
Total Crust |  $25.5^{+3.9}_{-3.4}$ |  $6.6^{+1.2}_{-1.0}$ |  $32.5\pm 3.8$

The expected total geoneutrino signal at JUNO can be expressed as the sum of
four independent contributions, signals from the bulk crust, the continental
lithospheric mantle(CLM), the depleted mantle(DM) and the enriched mantle(EM).
The bulk crust signal comprises the signal from two parts of the local crust
(LOC) and the signal from the far filed crust (FFC). The total geoneutrino
signal prediction at JUNO is shown in Table.3.

Table 3: Summary of expected geoneutrino signals and uncertainty (1 standard deviation) in TNU for JUNO, from uranium ($S_{\text{U}}$), thorium ($S_{\text{Th}}$), and the sum thereof ($S_{TOT}$). The bulk crust contribution includes that of the LOC and FFC. The lithosphere contribution is the sum of the bulk crust and CLM contributions. The total signal is the sum of that of the lithosphere, the DM and EM. All data is taken from the JULOC model other than the LOC contribution. |  $S_{\text{U}}\pm\sigma$ |  $S_{\text{Th}}\pm\sigma$ |  $S_{\text{U+Th}}\pm\sigma$
---|---|---|---
LOC |  $25.5^{+3.9}_{-3.4}$ |  $6.6^{+1.2}_{-1.0}$ |  $32.5\pm 3.8$
FFC |  $7.6\pm 1.6$ |  $2.3\pm 0.2$ |  $9.8\pm 1.7$
Bulk crust |  $33.1^{+4.2}_{-3.7}$ |  $8.9^{+1.2}_{-1.0}$ |  $42.3\pm 4.2$
CLM |  $1.3^{+2.4}_{-0.9}$ |  $0.4^{+1.0}_{-0.3}$ |  $2.1^{+3.0}_{-1.3}$
Lithosphere |  $34.4^{+4.6}_{-4.3}$ |  $9.3\pm 1.0$ |  $44.4^{+5.1}_{-4.4}$
DM |  $4.2$ |  $0.8$ |  $4.9$
EM |  $2.9$ |  $0.9$ |  $3.8$
Total |  $41.5^{+4.6}_{-4.3}$ |  $11.0\pm 1.0$ |  $53.1^{+5.1}_{-4.4}$

Table.3 shows the LOC contributes around $76.8\%$ of the total crustal
geoneutrino signal, more than 50% comes from the UC layer and MC layer
contributed 17%. Thus, to improve model precision, except UC layer, reduce the
uncertainty of MC layer are also important. The uncertainty of the estimated
geoneutrino signal is mainly attributable to crustal model inputs (i.e.crustal
layer thickness, density and U-Th abundances), the cross-section of the
geoneutrino detection and the electron antineutrino survival probability.The
latter two inputs are expected to introduce uncertainty of about $1\%$ [34],
which is excluded from the calculation. The uncertainty associated with the
crustal layer thickness and density are also negligible, as any uncertainty
from different cells could be cancelled out during the calculation, for
example when the uncertainty of different cells are uncorrelated, which is the
case when crustal layer thickness and density are obtained by means of seismic
ambient noise tomography [15]. In the end, only the uncertainty attributable
to U and Th abundances was considered, with the method of determining how
geoneutrino uncertainty is propagated from abundances being the same as used
for JULOC. If each cell of UC layer had enough rock samples to effectively
describe the abundance distribution, then any uncertainty was uncorrelated,
i.e. rock samples from different cells are not correlated. With respect to the
MC and LC layers, there were not enough samples, so the uncertainty is
considered to be correlated, resulting in the UC having a smaller uncertainty
than the MC and LC.

##  4 Comparison of models and conclusions

The existing local refined models are GIGJ (GOCE Inversion for Geoneutrinos at
JUNO) [14] and JULOC, with the GIGJ model being a refined geophysical model
and the JULOC being the first comprehensive 3-D high resolution geophysical
and geochemical model. Using the JULOC model as foundation, JULOC-I was
developed as an integrated model, incorporating seismic ambient noise,
gravity, rock sample and thermal flow data, and with these data being
constrained to each other.

The total lithosphere geoneutrino signal estimated by JULOC-I is
$44.4^{+5.1}_{-4.4}$ TNU, while the JULOC model generated an expected signal
of $40.4^{+5.6}_{-5.0}$ TNU.It can be seen that the updated model that has
incorporated additional data points has been able to reduce the total expected
geoneutrino signal error by 3%, from 14% to 11%. Compares with
$(30.9^{+6.5}_{-5.2})$ TNU[13] based on the global crustal model and 29.9 TNU
based on the GIGJ model, the expected geoneturino signal are significantly
high. The left hand plot in Fig.6 shows expected lithosphere geoneutrino
signals using by Global(Crust2.0), GIGJ, JULOC and JULOC-I models.The error
bars represent error mainly attributable to U and Th abundance. It can be seen
from this plot that the JULOC-I and JULOC models have higher expected
geoneutrino signals than the GIGJ and Global models on the basis of the
different local composition data. The significant increase in signals is
mainly attributable to the upper crust of the South China Block which is
richer in U/Th than the global average. This finding is consistent with the
widespread distribution of high U-Th granite intrusions in the area around
JUNO[41], and is also able to be confirmed by sampling survey results and
thermal flow data.

Figure 6: The plot on the left shows the expected lithosphere geoneutrino
signals $S_{\text{U}}$, $S_{\text{Th}}$ and $S_{\text{U+Th}}$ according to
different models, from left to right, JULOC-I, JULOC, GIGJ and Global models
(Crust2.0). The error bars represent error mainly attributable to U and Th
abundance. The plot on the right shows total expected geoneutrino signal, i.e.
with the addition of the contribution of the mantle to the plot on the left.
The shaded areas represent the precision of geoneutrino measurements at JUNO
over a one-year period. The mean value of the red-shaded area is derived from
the JULOC-I model plus mantle signal, while the blue-shaded area is derived
from the Global model plus mantle signal. The experiment precision is 13%[40].

The right-hand plot in Fig.6 compares the experimental precision level and the
total expected signal. The total expected signal is the sum of the lithosphere
crustal signal and the mantle signal. The mantle contribution has been taken
from the published paper on expected geoneutrino signals at JUNO, with it
considered to be a geochemical mantle model[38]. The precision level of
geoneutrino experiment measurements at JUNO are the red- and blue-shaded
areas. The mean value of the red-shaded area is based on the JULOC-I model
plus mantle contribution, while the blue-shaded area is based on the Global
model plus mantle contribution. Both are taken as 13% experiment precision for
one year data from R.Han et.al(2016)[40], though it should be noted that
different expected geoneutrino numbers will cause this value to differ
slightly. The expected signal error is the same as the crustal error in the
left hand plot. It can be seen from the right hand plot, based on one year of
data, and using the same Mantle model, the measurement precision at JUNO is
sufficient to distinguish different crustal models.

Figure 7: The total expected geoneutrino signal at JUNO in TNU based on
different mantle assumptions, i.e. cosmochemical, geochemical and
geodynamical. For each mantle model, the lithosphere contribution to total
expected signal has been estimated by means of the JULOC-I crustal model. The
10-year precision level of the geoneutrino measurement at JUNO are plotted as
shaded ranges, with the pink-shaded area representing the cosmochemical model,
the red-shaded area on the geochemical model and the blue on the geodynamical
model. The 10-year experiment precision is 5% [40].

As previously noted, the total measurement of geoneutrino signals at JUNO
comprises two parts, a lithosphere contribution and a mantle contribution,
with the mantle contribution being predicted on the basis of various
assumptions in relation to the distribution of HPE. All published expected
geoneutrino measurements at JUNO use values that are based on geochemical
assumptions. However, there are three classes of compositional
estimates–cosmochemical[35, 36], geochemical[37, 38], and geodynamical[39],
all of which give different abundances of HPEs, and which are also denoted as
Low-Q (low radiogenic power), Middle-Q (middle radiogenic power) and High-Q
(High radiogenic power) [16]. In Table.3, the mantle contribution is mainly
based on the geochemical model. The mantle results of the cosmochemical and
geodynamical models(from Table 4 in Šrámek et.al.(2013)[16]) were also added
to the lithosphere contribution, while the total expected geoneutrino signals
at JUNO using different mantle models are shown in Fig.7 with the crustal
model being JULOC-I. In this paper,uncertainty was not included in the mantle
assumptions, with the error bar being derived only from the lithosphere parts
as outlined in Table.3. The 10-year precision level of geoneutrino measurement
at JUNO is also plotted as shaded areas in Fig.7. The mean value of the pink-
shaded area is based on the cosmochemical mantle model, the red-shaded area on
the geochemical mantle model and the blue on the geodynamical mantle model.
The 10-year experiment precision is 5% [40].

As Fig.7 shows, the different mantles also produce different expected
geoneutrinos at JUNO. So, to distinguish between the different mantle models
sufficient experiment data is required, as is enhanced accuracy of crustal
models. With 10 years of data, JUNO will have collected sufficient data to
improve the precision level of measurement to 5%, and the JULOC-I model gives
the lowest error 11% currently. Both of these provide JUNO with an opportunity
to distinguish different mantle models. Looking into the future, a higher
accuracy crustal model will still be required, not only to test different
mantle models, but also to constrain the U/Th ratio.

In conclusion, the JUNO detector will provide an exciting opportunity to
obtain a statistically significant and precise measurement of geoneutrinos. An
accurate a priori estimation of the crustal contribution is an important and
necessary step to translate particle physics measurements into geo-scientific
goals. The first accurate, integrated local crustal model based on seismic,
gravity, abundance and thermal flow data has been developed and which has been
proven to predict geoneturino flux with the highest accuracy of the four
existing crustal models. So, within its first year of operation, JUNO should
be able to distinguish the four different local crustal models. Within 10
years of data, accuracy local crustal models should offer the possibility of
extracting the mantle component, which is much sought-after goal in the geo-
sciences.

## List of abbreviations used

KamLAND: Kamioka Liquid Scintillator Antineutrino Detector
JUNO: Jiangmen Underground Neutrino Observatory
GIGJ: GOCE Inversion for Geoneutrinos at JUNO
JULOC: JUno LOCal model
JULOC-I: JUno LOCal - Integrate model
LOC: Local Crust
FFC: Far Field Crust
CLM: Continental Lithospheric Mantle
DM: Depleted Mantle
EM: Enriched Mantle
UC: Upper Crust
MC: Middle Crust
LC: Lower Crust
OC: Oceanic Crust
IBD: Inverse Beta Decay
TNU: Terrestrial Neutrino Units
SED: Sediments
BAS: Basalt
GRA: Granite-Granodiorite
BGN-GGN: Granite Gneiss, Biotite Gneiss
MGW: Siltstone, Sandstone, Conglomerate and Shale
SLT-PHY-QSC: Slate, Phyllite, Schist
BGN-GGN: Granite Gneiss, Biotite Gneiss
AMP: Amphibolite
IGR: Intermediate Granulite
MGR: Mafic Granulite
GAB: Gabbro
MPE: Metapelite
MTL: Mantle Materials

## Acknowledgement

This study is supported by National Natural Science Foundation of China
(U1865206 and 42074092), National Key Research and Development Program of
China (2018YFA0404100).

## Appendices

##  Appendix A Estimation of lithology proportion by means of 3-D shear-wave
velocity

The method of building a lithology proportion map by means of 3-D S-wave
velocity has previously been used by Takeuchi et al.(2019) [26] and Sun
et.al(2022). This method firstly assumes a priori probability map of lithology
as $P^{x}_{(i)}$, where ’$i$’ is the type of rock and ’$x$’ is the location of
that rock. This priori probability allows local lithological features to be
considered.

In the prior probability map, the lithology of the UC is estimated from the
geological map, while the lithology of the MC and LC are estimated based on
the average crustal composition of East China[3]. Table.4 shows the a priori
probability map, there are SED(Sediments), BAS(Basalt), GRA, BGN-GGN(Granite
Gneiss, Biotite Gneiss), MGW, and SLT-PHY-QSC(Slate, Phyllite, Schist) in UC
layer. In MC layer there are BGN-GGN(Granite Gneiss, Biotite Gneiss) and
AMP(Amphibolite). IGR(Intermediate granulite),MGR(Mafic Granulite),
GAB(Gabbro), MPE (Metapelite) and MTL(Mantle Materials) are in LC layer.

Table 4: The input a priori probability of continental crust rock composition in the research area. UC |  MC |  LC
---|---|---
SED 0.08 |  BGN-GGN 0.2 |  IGR 0.6
BAS 0.02 |  AMP 0.8 |  MGR 0.15
BGN-GGN 0.1 |  |  MPE 0.15
GRA 0.4 |  |  GAB 0.05
MGW 0.25 |  |  MTL 0.05
SLT-PHY-QSC 0.15 |  |

An a posteriori probability lithology map is then developed by establishing a
mathematical connection between the crustal P-wave velocity $V_{P}^{obs}$ and
rock type proportion. The $V_{P}^{obs}$ velocity is calculated from S-wave
velocity $V_{S}^{obs}$ by means of the empirical equation[31],as
eq.LABEL:Eq:vstovp, with the S-wave velocity having been derived from the 3-D
velocity model constructed by means of the gravity-seismology joint inversion
method as described earlier, and

$\displaystyle V_{P}^{obs}(km/s)$
$\displaystyle=0.9409+2.0947V_{S}^{obs}-0.8206(V_{S}^{obs})^{2}$ (1)
$\displaystyle+0.2683(V_{S}^{obs})^{3}-0.0251(V_{S}^{obs})^{4}$

where $V_{P}^{obs}$ and $V_{S}^{obs}$ are wave velocities at location ’x’ as
derived from the geophysical model. The a posteriori probability of a rock
type can be calculated:

$P(i|V_{P}^{obs})=\frac{P(V_{P}^{obs}|i)P(x|i)}{(\Sigma_{i=1}^{n}P(V_{P}^{obs}(x|i)P(x|i)}$
(2)

where $P(x|i)$ refers to the a priori probability and $P(V_{P}^{obs}|i)$
represents the probability of $V_{P}^{obs}$ when a certain rock occurs,
$P(V_{P}^{obs}(x|i))$is sampled from a gaussian distribution with mean and
standard deviation summarized in in Tab.4 [33, 22].

The 3-D geological proportion map shows that SED and MGW are dominant in the
top 1 km of the research area, accounting for 50% and 80%, respectively. GRA
is mainly distributed in the southwest and northeast of the research area,
accounting for less than 40%. At a depth of 2 km, SED has almost disappeared
with GRA becoming more prevalent, accounting for up to 70%, in some places,
though MGW still dominates overall. At a depth of 3 to 5 km, GRA begins to
surpass MGW as the main component of crustal rocks, but trends downward from
60% to about 50%. At a depth of 3 km, in most areas, the proportion of MGW
decreases sharply to below 30%. SLT-PHY-QSC and BGN-GGN begin to appear at a
depth of 4 km, occupying about 20% and 10% of the crust, respectively.

From a depth of 6 km, GRA and BGN-GGN are grouped together due to their
similar P wave velocities. The proportion of GRA-BGN-GGN remains at about 50%,
but it gradually decreases after 15 km and disappears completely at a depth of
30 km. SLT-PHY-QSC are the major constituents of the upper crustal basement,
increasing from 40% at 6 km to 50% at 15 km. MGW decreases with increasing
depth and completely disappears at a depth of 15 km.

The results are consistent with previous studies on the composition of the MC,
which suggests that AMP is the main components. Starting at 20 km, the main
components are IGR, MGR, MPE, and GAB. Within the range of 20–30 km, IGR and
MGR accounted for more than 60% and 20%, respectively, with MPE and GAB
represnting only a small proportion. As the depth increases from 31 km to 45
km, the lithology of the crust in the research area becomes dominated by MPE
and GAB, with mantle material becoming the main component thereafter. Fig.8
shows the proportions of different rock types at a depth of 5 km.

Figure 8: Proportions of volume of different rock types at 5 km depth.

##  Appendix B Building of U/Th concentrations cell-by-cell

Within a 200-kilometer radius of the JUNO center was divided into
$0.4^{\circ}\times 0.4^{\circ}\times$1 km cell size, whereas in the
200-500-kilometer range, the cell size was $1^{\circ}\times 1^{\circ}\times$1
km. Each cell was further subdivided into ranges according to the 2.1
mentioned rock types. U/Th abundance of each rock type combined with the 3-D
lithology proportion map derived from seismic wave velocity, the average U/Th
abundance of each cell was subsequently calculated.

The processor to calculate U/Th abundance of each rock type as below: if there
were more than five data points available for a rock type in a cell, the data
were subsequently used to calculate abundance and related uncertainty. Given
this criterion, when the available sample size was less than 10, expected
abundance and uncertainty were calculated by using the mean and standard
deviation, respectively. For sample size exceeding 10, expected abundance and
uncertainty were determined by means of normal,logarithmic or gamma
distribution fitting.

The fitting function selected depended on the shape of the curve for the
abundance distribution in a particular rock type. When the shape was
symmetric, i.e. the arithmetic mean and median are equal, normal distribution
was selected. When the geometric mean was greater than the median, a log
normal distribution was chosen. Takeuchi et.al(2019)[26] demonstrated that
when a distribution is a highly skewed log normal distribution, the median may
be smaller than the arithmetic mean to an order of magnitude, increasing the
error in the fitted estimate. However, a gamma distribution does not cause a
deviation from the mean, which would seem to render it a better choice than
normal and lognormal distributions(see Fig.9).

Figure 9: The distribution of GRA and MGW U/Th with different fitting
function.

When there were less than five samples of a certain rocks in the cell, the
expected abundance and error were estimated using all the data for that kind
of rock in the study area.

The 3-D PDFs of U/Th concentrations in each cell were constructed by
$A(i)\times P(i)$, where $A(i|x)$ and $P(i|x)$ are the abundance and
posteriori probability of i rock type,respectively. These 3-D PDFs were
subsequently constrained by thermal data. A total of 182 thermal flow data
sets around JUNO were used for comparison. The U, Th and K data from the 3-D
model were used to determine the heat production of each layer, with one-
dimensional heat transfer then being used to test the reliability of the 3-D U
and Th concentrations model.

##  Appendix C Predicting the geoneutrino signal based on crustal model

The joint geophysical and geochemical 3-D local crustal models conducted
around JUNO can be used to calculate the geoneutrino flux produced by the
Earth’s crust. However, predicting the geoneutrino signal that could be
detected at the JUNO detector also requires detection information of the
particle by means of inverse beta decay (IBD) reaction. The detailed
mathematics of geoneutrino signal calculation is outlined in JULOC paper R.Gao
et.al [15]. All parameters were largely kept the same, for example, the
detector efficiency $\epsilon(E)$ was assumed to be 1, $\sigma(E)$ was taken
from Strumia and Vissani[42],the geoneutrino energy spectrum $f_{i}(E)$ was
taken from Sanshiro[44]. With respect to survival probability, for simplicity
sake, the baseline averaged probability, which approximates as $P_{ee}=0.55$
was used in JULOC model. Along with the improvement for the experimental
accuracy of the oscillation parameters and current and future geoneutrino
measurements, this simple approximation cannot be taken for granted. In
JULOC-I, the precision oscillation effects was taking into account during the
crustal geoneutrino calculations. The three-flavor survival probability for
electron antineutrino $\bar{{\nu}_{e}}$ propagating in the vacuum can be
described as

$P_{ee}=1-P_{0}-P_{\ast};$ (3)

with

$P_{0}=sin^{2}2\theta_{12}cos^{4}\theta_{13}sin^{2}\Delta_{21},\\\ $ (4)
$P_{\ast}=\frac{1}{2}sin^{2}2\theta_{13}(1-cos\Delta_{\ast}cos\Delta_{21})\\\
+cos2\theta_{12}sin\Delta_{\ast}sin\Delta_{21},$ (5)

$\Delta_{ij}=1.27(\Delta
m^{2}_{ij}L)/E_{\bar{\nu}},\Delta_{\ast}=\Delta_{31}+\Delta_{32}$ (6)

where $\theta_{12}$ and $\theta_{13}$ are the mixing angles, $\Delta
m^{2}_{ij}$ is the square mass difference of antineutrinos in $eV^{2}$,
$E_{\bar{\nu}}$ is the antineutrino energy in MeV, and L is the propagation
distance in meters. In the calculation, we assume the mass ordering to be
normal ($m1<m2<m3$) and take the oscillation parameters provided by the
Particle Data Group (PDG). At last, Density and the abundance for each cell
were taken from the new JULOC-I model.

When calculating the geoneutrino flux produced by the local crust(LOC)
surrounding JUNO, the regional crust was divided into three layers, with each
layer being further divided into $0.4^{\circ}\times 0.4^{\circ}\times 1$ km
cells and spatial, geophysical, and geochemical attributes being assigned to
each cell. The total LOC contribution to the geoneutrino signal at JUNO can be
obtained by summing the signal of each cell, layer by layer(see table2). The
JULOC-I integrated model only updated the JULOC model in terms of the
contribution made by the continetal crust, with very small contribution made
by the Oceanic crust being kept the same, as JULOC-I did not involve any new
Oceanic crust sampling. Fig.10 shows the 1,000 times geoneutrino signal
distribution as calculated from 1,000 random samples chosen in accordance with
the statistical distribution of input parameters, namely the density, layer
thickness and U-Th abundances of each cell. Uncertainties associated with the
input random numbers is propagated through the geoneutrino signal calculation
to the distribution of the geoneutrino signal.

Figure 10:  The 1,000-time geoneutrino signal distribution from 238U(left),
232Th(middle) and sum of them(right) of the UC (upper),MC (upper middel), LC
(lower middel) and the sum(bottom ). The corresponding geoneutrino signals are
the mean values and $1\sigma$ errors in the plot.

## References

  * [1] C. P. Conrad and C. Lithgow-Bertelloni, “How mantle slabs drive plate tectonics,” _Science_ , vol. 298, no. 5591, pp. 207–209, 2002.
  * [2] R. White and D. McKenzie, “Magmatism at rift zones: the generation of volcanic continental margins and flood basalts,” _Journal of Geophysical Research: Solid Earth_ , vol. 94, no. B6, pp. 7685–7729, 1989.
  * [3] S. Gao, T.-C. Luo, B.-R. Zhang, H.-F. Zhang, Y.-w. Han, Z.-D. Zhao, and Y.-K. Hu, “Chemical composition of the continental crust as revealed by studies in east china,” _Geochimica et cosmochimica acta_ , vol. 62, no. 11, pp. 1959–1975, 1998.
  * [4] T. Araki, S. Enomoto, K. Furuno, Y. Gando, K. Ichimura, H. Ikeda, K. Inoue, Y. Kishimoto, M. Koga, Y. Koseki _et al._ , “Experimental investigation of geologically produced antineutrinos with kamland,” _Nature_ , vol. 436, no. 7050, pp. 499–503, 2005.
  * [5] S. Abe, T. Ebihara, S. Enomoto, K. Furuno, Y. Gando, K. Ichimura, H. Ikeda, K. Inoue, Y. Kibe, Y. Kishimoto _et al._ , “Precision measurement of neutrino oscillation parameters with kamland,” _Physical Review Letters_ , vol. 100, no. 22, p. 221803, 2008.
  * [6] A. Gando, Y. Gando, H. Hanakago, H. Ikeda, K. Inoue, K. Ishidoshiro, H. Ishikawa, M. Koga, R. Matsuda, S. Matsuda _et al._ , “Reactor on-off antineutrino measurement with kamland,” _Physical Review D_ , vol. 88, no. 3, p. 033001, 2013.
  * [7] G. Bellini, J. Benziger, S. Bonetti, M. B. Avanzini, B. Caccianiga, L. Cadonati, F. Calaprice, C. Carraro, A. Chavarria, F. Dalnoki-Veress _et al._ , “Observation of geo-neutrinos,” _Physics Letters B_ , vol. 687, no. 4-5, pp. 299–304, 2010.
  * [8] G. Bellini, J. Benziger, D. Bick, G. Bonfini, D. Bravo, M. B. Avanzini, B. Caccianiga, L. Cadonati, F. Calaprice, P. Cavalcante _et al._ , “Measurement of geo-neutrinos from 1353 days of borexino,” _Physics Letters B_ , vol. 722, no. 4-5, pp. 295–300, 2013.
  * [9] M. Agostini, S. Appel, G. Bellini, J. Benziger, D. Bick, G. Bonfini, D. Bravo, B. Caccianiga, F. Calaprice, A. Caminata _et al._ , “Spectroscopy of geoneutrinos from 2056 days of borexino data,” _Physical Review D_ , vol. 92, no. 3, p. 031101, 2015.
  * [10] M. Agostini, K. Altenmüller, S. Appel, V. Atroshchenko, Z. Bagdasarian, D. Basilico, G. Bellini, J. Benziger, D. Bick, G. Bonfini _et al._ , “Comprehensive geoneutrino analysis with borexino,” _Physical Review D_ , vol. 101, no. 1, p. 012009, 2020.
  * [11] Y.-F. Li, J. Cao, Y. Wang, and L. Zhan, “Unambiguous determination of the neutrino mass hierarchy using reactor neutrinos,” _Physical Review D_ , vol. 88, no. 1, p. 013008, 2013.
  * [12] F. An, G. An, Q. An, V. Antonelli, E. Baussan, J. Beacom, L. Bezrukov, S. Blyth, R. Brugnera, M. B. Avanzini _et al._ , “Neutrino physics with juno,” _Journal of Physics G: Nuclear and Particle Physics_ , vol. 43, no. 3, p. 030401, 2016.
  * [13] V. Strati, M. Baldoncini, I. Callegari, F. Mantovani, W. F. McDonough, B. Ricci, and G. Xhixha, “Expected geoneutrino signal at juno,” _Progress in Earth and Planetary Science_ , vol. 2, no. 1, pp. 1–7, 2015\.
  * [14] M. Reguzzoni, L. Rossi, M. Baldoncini, I. Callegari, P. Poli, D. Sampietro, V. Strati, F. Mantovani, G. Andronico, V. Antonelli _et al._ , “Gigj: A crustal gravity model of the guangdong province for predicting the geoneutrino signal at the juno experiment,” _Journal of Geophysical Research: Solid Earth_ , vol. 124, no. 4, pp. 4231–4249, 2019.
  * [15] R. Gao, Z. Li, R. Han, A. Wang, Y. Li, Y. Xi, J. Liu, X. Mao, Y. Sun, and Y. Xu, “Juloc: A local 3-d high-resolution crustal model in south china for forecasting geoneutrino measurements at juno,” _Physics of the Earth and Planetary Interiors_ , vol. 299, p. 106409, 2020.
  * [16] O. Šrámek, W. F. McDonough, E. S. Kite, V. Lekić, S. T. Dye, and S. Zhong, “Geophysical and geochemical constraints on geoneutrino fluxes from earth’s mantle,” _Earth and Planetary Science Letters_ , vol. 361, pp. 356–366, 2013.
  * [17] C. Bassin, “The current limits of resolution for surface wave tomography in north america,” _EOS Trans. AGU. 81: Fall Meet. Suppl., Abstract_ , 2000\.
  * [18] G. Laske, G. Masters, and C. Reif, “Crust 2.0: A new global crustal model at 2$\times$ 2 degrees,” _Institute of Geophysics and Planetary Physics, The University of California, San Diego, website http://mahi. ucsd. edu/Gabi/rem. dir/crust/crust2. html_ , 2001.
  * [19] N. Shapiro and M. Ritzwoller, “Monte-carlo inversion for a global shear-velocity model of the crust and upper mantle,” _Geophysical Journal International_ , vol. 151, no. 1, pp. 88–105, 2002.
  * [20] M. Reguzzoni and N. Tselfes, “Optimal multi-step collocation: application to the space-wise approach for goce data analysis,” _Journal of Geodesy_ , vol. 83, no. 1, pp. 13–29, 2009.
  * [21] M. Negretti, M. Reguzzoni, and D. Sampietro, “A web processing service for goce data exploitation,” in _First International GOCE Solid Earth workshop, Enschede, The Netherlands_ , 2012, pp. 16–17.
  * [22] Y. Huang, V. Chubakov, F. Mantovani, R. L. Rudnick, and W. F. McDonough, “A reference earth model for the heat-producing elements and associated geoneutrino flux,” _Geochemistry, Geophysics, Geosystems_ , vol. 14, no. 6, pp. 2003–2029, 2013.
  * [23] F. Mantovani, L. Carmignani, G. Fiorentini, and M. Lissia, “Antineutrinos from earth: A reference model and its uncertainties,” _Physical Review D_ , vol. 69, no. 1, p. 013001, 2004.
  * [24] N. Holloway, “North palawan block, philippines—its relation to asian mainland and role in evolution of south china sea,” _AAPG Bulletin_ , vol. 66, no. 9, pp. 1355–1383, 1982.
  * [25] P. Clift and J. Lin, “Preferential mantle lithospheric extension under the south china margin,” _Marine and Petroleum Geology_ , vol. 18, no. 8, pp. 929–945, 2001.
  * [26] N. Takeuchi, K. Ueki, T. Iizuka, J. Nagao, A. Tanaka, S. Enomoto, Y. Shirahata, H. Watanabe, M. Yamano, and H. K. Tanaka, “Stochastic modeling of 3-d compositional distribution in the crust with bayesian inference and application to geoneutrino observation in japan,” _Physics of the Earth and Planetary Interiors_ , vol. 288, pp. 37–57, 2019.
  * [27] N. Du, Z. Li, T. Hao, X. Xia, Y. Shi, and Y. Xu, “Joint tomographic inversion of crustal structure beneath the eastern tibetan plateau with ambient noise and gravity data,” _Geophysical Journal International_ , vol. 227, no. 3, pp. 1961–1979, 2021.
  * [28] G. Balmino, N. Vales, S. Bonvalot, and A. Briais, “Spherical harmonic modelling to ultra-high degree of bouguer and isostatic anomalies,” _Journal of Geodesy_ , vol. 86, no. 7, pp. 499–520, 2012.
  * [29] N. K. Pavlis, S. A. Holmes, S. C. Kenyon, and J. K. Factor, “The development and evaluation of the earth gravitational model 2008 (egm2008),” _Journal of geophysical research: solid earth_ , vol. 117, no. B4, 2012.
  * [30] H. CHEN, Z. LI, R. CHU, J. XIE, J. CHONG, and Q. YE, “3d s-wave velocity and radial anisotropy structure of the crust and uppermost mantle in yunnan, se tibet,” _Chinese Journal of Geophysics_ , vol. 61, no. 8, pp. 3219–3236, 2018.
  * [31] T. M. Brocher, “Empirical relations between elastic wavespeeds and density in the earth’s crust,” _Bulletin of the seismological Society of America_ , vol. 95, no. 6, pp. 2081–2092, 2005.
  * [32] T.-Y. HAO, W.-J. HU, J. XING, L.-T. HU, Y. XU, J.-X. QIN, S.-H. LIU, and S.-M. LEI, “The moho depth map (1: 5000000) in the land and seas of china and adjacent areas and its geological implications,” _Chinese Journal of Geophysics_ , vol. 57, no. 12, pp. 3869–3883, 2014.
  * [33] N. I. Christensen and W. D. Mooney, “Seismic velocity structure and composition of the continental crust: A global view,” _Journal of Geophysical Research: Solid Earth_ , vol. 100, no. B6, pp. 9761–9788, 1995.
  * [34] X. Mao, R. Han, and Y.-F. Li, “Non-negligible oscillation effects in the crustal geoneutrino calculations,” _Physical Review D_ , vol. 100, no. 11, p. 113009, 2019.
  * [35] H. S. C. O’Neill and H. Palme, “Collisional erosion and the non-chondritic composition of the terrestrial planets,” _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 366, no. 1883, pp. 4205–4238, 2008.
  * [36] H. Palme and H. S. C. O’Neill, “Cosmochemical estimates of mantle composition,” _Treatise on geochemistry_ , vol. 2, p. 568, 2003.
  * [37] R. Arevalo Jr, W. F. McDonough, and M. Luong, “The k/u ratio of the silicate earth: Insights into mantle composition, structure and thermal evolution,” _Earth and Planetary Science Letters_ , vol. 278, no. 3-4, pp. 361–369, 2009\.
  * [38] W. F. McDonough and S.-S. Sun, “The composition of the earth,” _Chemical geology_ , vol. 120, no. 3-4, pp. 223–253, 1995.
  * [39] D. L. Turcotte and G. Schubert, _Geodynamics_. Cambridge university press, 2002.
  * [40] R. Han, Y.-F. Li, L. Zhan, W. McDonough, J. Cao, and L. Ludhova, “Potential of geo-neutrino measurements at juno,” _Chinese Physics C_ , vol. 40, p. 33003, 11 2015.
  * [41] Y. Wang, “Heat flow pattern and lateral variations of lithosphere strength in china mainland: constraints on active deformation,” _Physics of the Earth and Planetary Interiors_ , vol. 126, no. 3, pp. 121–146, 2001, thermal Studies of the Earth’s Structure and Geodynamics .IUGG 99 Assembly at Birmingham , UK . [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0031920101002515
  * [42] A. Strumia and F. Vissani, “Precise quasielastic neutrino/nucleon cross-section,” _Physics Letters B_ , vol. 564, no. 1, pp. 42–54, 2003. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0370269303006166
  * [43] F. Capozzi, G. L. Fogli, E. Lisi, A. Marrone, D. Montanino, and A. Palazzo, “Status of three-neutrino oscillation parameters, circa 2013,” _Phys. Rev. D_ , vol. 89, p. 093018, May 2014. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevD.89.093018
  * [44] E. Sanshiro, “geoneutrino spectrum,” https://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/index.html.


# Predictor-Corrector(PC) Temporal Difference(TD) Learning (PCTD)

Caleb Bowyer
Department of Electrical and Computer Engineering
University of Florida
Gainesville, FL
<EMAIL_ADDRESS>

###### Abstract

Using insight from numerical approximation of ODEs and the problem formulation
and solution methodology of TD learning through a Galerkin relaxation, I
propose a new class of TD learning algorithms. After applying the improved
numerical methods, the parameter being approximated has a guaranteed order of
magnitude reduction in the Taylor Series error of the solution to the ODE for
the parameter $\theta(t)$ that is used in constructing the linearly
parameterized value function. Predictor-Corrector Temporal Difference (PCTD)
is what I call the translated discrete time Reinforcement Learning(RL)
algorithm from the continuous time ODE using the theory of Stochastic
Approximation(SA). Both causal and non-causal implementations of the algorithm
are provided, and simulation results are listed for an infinite horizon task
to compare the original TD(0) algorithm against both versions of PCTD(0).

###### Contents

  1. 1 Introduction
  2. 2 Numerical Method Motivation
  3. 3 PCTD(0) Derivation
    1. 3.1 Causal PCTD(0)
    2. 3.2 Non-causal PCTD(0)
    3. 3.3 Extension to PCTD($\lambda$)
  4. 4 Prediction Task
  5. 5 Simulation Results
  6. 6 Conclusion and Future Work
  7. 7 Broader Impact

##  1 Introduction

This novel family of reinforcement learning algorithms will stay within the
realm of Eulerian approximations to the parameter ODE which is used to
construct a TD-like learning algorithm, but will show to be a strict
improvement over the explicit Eulerian based algorithms found in [1] by an
order of magnitude in the step size of the algorithm.

Firstly, I will rederive "TD(0) learning" for the discounted cost objective.
Discounted cost is the only objective considered in this paper although the
ideas carry over to other cost settings. Additionally, a finite state space
$X=\\{1,2,...,n\\}$ will be used for all developments. An uncontrolled
irreducible aperiodic Markov chain samples points from this state space. To
understand the TD(0) algorithm, one must first understand the discounted
Poisson equation:

$h(x)=c(x)+\beta Ph(x)=c(x)+\beta E[h(X_{k+1})|X_{k}=x]$ (1)

where $h(x)$ is our relative cost value function and $c(x)$ is the cost
incurred by being in state $x$. Lastly, $P$ is the transition matrix for the
Markov chain that also functions as the conditional expectation on the
function $h(x)$.

Rewriting, we have

$0=E[-h(X_{k})+c(X_{k})+\beta h(X_{k+1})|\mathcal{F}_{k}],$

where $\mathcal{F}_{k}=\\{\sigma{(X_{i}):i\leq k}\\}$ is our sigma algebra or
history up to time $k$. The above step follows from the Markov property.

TD learning is based on a linearly parameterized family of value functions
$h^{\theta}:X\mapsto\mathbb{R}$ with a Galerkin relaxation applied to
approximate the conditional expectation. For the linear parameterization, let

$h^{\theta}(x)=\theta^{T}\phi(x)=\sum_{i=1}^{d}\theta_{i}\phi_{i}(x),$ (2)

where $d$ is some integer, usually much smaller than the state space size and
$\phi_{i}$ for $i=1,2,...,d$ are the basis functions. As in TD(0), I use
$\zeta_{k}=\phi(X_{k})\in\mathbb{R}^{d}$ for our eligibility vectors. Notation
is now introduced for the Bellman error as
$\epsilon(X_{k},X_{k+1})=-h(X_{k})+c(X_{k})+\beta h(X_{k+1})$. It is a known
fact that a Galerkin relaxation can be applied to the above conditional
expectation:

$E[\epsilon(X_{k},X_{k+1})\zeta_{k}]=0,$ (3)

for any bounded $\mathcal{F}_{k}$ measurable R.V. $\zeta_{k}$.

Now, using 2 and 3 in 1, one gets:

$\displaystyle E[-h^{\theta}(X_{k})+c(X_{k})+\beta
h^{\theta}(X_{k+1})|\mathcal{F}_{k}]$ $\displaystyle=$ $\displaystyle
E[-\phi^{T}(X_{k})\theta+c(X_{k})+\beta\phi^{T}(X_{k+1})\theta|\mathcal{F}_{k}]$
$\displaystyle=$ $\displaystyle
E[[-\phi^{T}(X_{k})\theta+\beta\phi^{T}(X_{k+1})\theta+c(X_{k})]\zeta_{k}]|_{\theta=\theta^{*}}$
$\displaystyle=$ $\displaystyle 0,$

where the last equality is achieved at the optimal parameter $\theta^{*}$
which is the root of the equation.

From this, I define the parameterized Bellman error as

$\epsilon^{\theta}(X_{k},X_{k+1})=[(\beta\phi^{T}(X_{k+1})-\phi^{T}(X_{k}))\theta+c(X_{k})].$

This expression is also called the temporal difference sequence:

$d_{k+1}=(\beta\phi^{T}(X_{k+1})-\phi^{T}(X_{k}))\theta_{k}+c(X_{k}).$

The TD(0) algorithm is often presented as:

$\theta_{k+1}=\theta_{k}+\alpha_{k}\zeta_{k}d_{k+1}$ (4)

without any of the steps to the derivation presented above.

Now, consider how the TD(0) algorithm can be derived from the analysis of a
particular ODE and leads to a root finding problem with expectation being
taken in steady-state. It is a two part process:

Step 1: Consider an ODE we wish to solve

$\frac{d\theta(t)}{dt}=\alpha\bar{f}(\theta(t)),$ (5)

$\alpha>0$ where

$\bar{f}(\theta(t))=E_{\infty}[\zeta_{k}(-h^{\theta}(X_{k})+c(X_{k})+\beta
h^{\theta}(X_{k+1}))]$ (6)

What has been done for the last 30 years or so has been to take this fixed
point equation and apply explicit or direct Euler method to the ODE.

Step 2a) Translate to discrete time

$\theta_{k+1}=\theta_{k}+\alpha_{k}\bar{f}(\theta_{k})$ (7)

Step 2b) The expectation is dropped when applying SA

$\theta_{k+1}=\theta_{k}+\alpha_{k}[\zeta_{k}(-h^{\theta}(X_{k})+c(X_{k})+\beta
h^{\theta}(X_{k+1}))]$ (8)

This has led to the well known family of TD learning algorithms and later
development of the control algorithms: SARSA and Q-learning. For more details
on the ODE design and connection to the classical discrete time algorithms see
[2] and the connection with Q-learning as well as more function approximation
considerations in the asymptotic case [3]. In this paper, I am proposing to
make all of the above algorithms more efficient in learning the optimal value
function faster by improving the updates to the parameter $\theta_{k}$ with a
better translation step from the algorithm design.

Furthermore, I will provide conditions upon which the algorithm is stable in
the convergence of the parameter $\theta_{k}$ to $\theta^{*}$, that results in
learning the optimal value function. In a separate research direction, some
work on nonlinear parameterized value functions has been done, but it is often
dismissed because of divergent counterexamples found in the literature. E.g.,
[1] and more recently a study in [4] investigate this type of function
approximation. For linear value function learning in a non-asymptotic setting
see [5] and for the control learning algorithms of this flavor see [6]. Now, I
will reformulate the basic problem of TD(0) from a modern viewpoint before
deriving the PCTD(0) algorithms.

The goal of RL can be stated as trying to solve i.e., find the zeros of this
equation:

$\bar{f}(\theta)=E_{\infty}[\zeta_{k}\epsilon^{\theta}(X_{k},X_{k+1})]$ (9)

The goal is to find a $\theta^{*}$ s.t. $\bar{f}(\theta^{*})=0$. In other
words, at the fixed point $\theta^{*}$ of the ODE:

$\frac{d\theta(t)}{dt}=A\theta(t)+b|_{\theta(t)=\theta^{*}}=0$ (10)

and so the deterministic recursion after applying SA produces a stable
algorithm provided A is negative definite.

From TD(0) learning in the discounted case, one can obtain:

$A=-E[\phi(X_{k})\phi(X_{k})^{T}]+\beta E[\phi(X_{k})\phi(X_{k+1})]$

and

$b=-E[\phi(X_{k})c(X_{k})]$

which forms our linear equation for $\theta^{*}$ of the form
$\bar{f}(\theta^{*})=0\Longrightarrow A\theta^{*}+b=0$ implies the solution is
easily obtained from solving the system of equations. Provided A is full rank
implies a sufficient condition for learning the optimal parameter vector
$\theta^{*}$. Alternatively, the first matrix $E[\phi(X_{k})\phi(X_{k})^{T}]$
must form a complete linear independent basis, and because of the contraction
term with $\beta$, that ensures stability of TD(0).

A note on step-size of the algorithm: The step-size sequence in the final
version of the translated algorithms will take the form of
$\alpha_{k}=\frac{g}{k+n_{0}}$ with $g,n_{0}>0$. This is sometimes done to
tune the step size to ensure optimal asymptotic covariance of the parameter
estimate. The gain sequence, $\alpha_{k}$ whether it is decayed or a fixed
learning rate $\alpha$ must be tuned for the problem at hand. See [7] for a
complete exposition. The form of the above varying step-size has been shown to
handle transients and steady-state performance through appropriate choice of
$n_{0}$ and $g$.

##  2 Numerical Method Motivation

In these pages, I propose and derive a new family of RL algorithms that come
from a better study of the numerical methods that lead from the translation of
the continuous time ODE to the final discrete time RL algorithm. Step 1:
Create your ODE as in 5 which will be used for the rest of this paper,
although the methods and algorithms clarified and elaborated in this paper can
be adapted to more complicated ODEs or cost objectives for the RL solution of
choice. Step 2 is what will be improved from previous research. Denote

$f(\theta)=[\zeta_{n}(-h^{\theta}(X_{n})+c(X_{n})+\beta h^{\theta}(X_{n+1}))]$
(11)

From the numerical ODE literature it is known that the predictor-corrector
method performs better than explicit Euler because it cancels out the leading
order Taylor Series error (T.E.) First let us review explicit Euler and
Modified Euler Methods. The Taylor Series(T.S.) expansion of

$\dot{\theta}=f(\theta,t)$

is

$\theta(t)=\theta_{0}+\theta_{0}^{\prime}(t-t_{0})+\frac{1}{2!}\theta_{0}^{\prime\prime}(t-t_{0})^{2}+...$

For small $\alpha=(t_{1}-t_{0})$ if we just keep the $\mathcal{O}(\alpha)$
term then we obtain

$\displaystyle\theta(t_{0}+\alpha)$ $\displaystyle=$
$\displaystyle\theta(t_{0})+\alpha\theta_{0}^{\prime}$

This implies that

$\displaystyle\theta_{n+1}$ $\displaystyle=$
$\displaystyle\theta_{n}+\alpha\theta_{n}^{\prime}+\mathcal{O}(\alpha^{2})=\theta_{n}+\alpha
f(t_{n},\theta_{n})+\mathcal{O}(\alpha^{2})\forall n=0,1,2,...$

Euler’s method is explicit. It is easy to implement, but inaccurate.

We see that the direct Euler’s method has only a first order accuracy. The
next development before I can set up the new algorithm is to review the Euler
Implicit method: Evaluate the slope $f$ at $(t_{n+1},\theta_{n+1})$

$\frac{\theta_{n+1}-\theta_{n}}{\alpha}+\alpha\theta_{n+1}^{\prime\prime}=f(t_{n+1},\theta_{n+1})$
$\Longrightarrow\theta_{n+1}=\theta_{n}+\alpha
f_{n+1}-\alpha^{2}\theta_{n+1}^{\prime\prime}=\theta_{n}+\alpha
f_{n+1}+\mathcal{O}(\alpha^{2})$

The local T.E. in the explicit method is
$0.5\theta_{n}^{\prime\prime}\alpha^{2}$, and the local T.E. in the implicit
method is $-0.5\theta_{n+1}^{\prime\prime}\alpha^{2}$. Note the two terms are
opposite to each other. Major problems with both Implicit and Explicit Euler
methods: 1) In general $f(\theta_{n+1},t_{n+1})$ depends on $\theta_{n+1}$
which is an unknown at time $n$. 2) If $f(\theta,t)$ is a non-linear function
of $\theta$, the procedure for obtaining $\theta_{n+1}$ can be complicated
using the implicit method. In both methods, the accuracy is poor. The T.E. is
$\mathcal{O}(\alpha)$. This is where the modified Euler predictor-corrector
method improves the accuracy of $\theta_{n}$.

##  3 PCTD(0) Derivation

###  3.1 Causal PCTD(0)

PCTD(0) uses the improved numerical methods to improve the approximation to
the optimal value function. In my notation for the derivation of PCTD(0) I
will reserve the symbol h for the relative cost value function. For step-size,
the symbol $\alpha$ will be used. One can combine the Euler Explicit and
Implicit methods to cancel out the leading order T.E.

Consider slope:

$\displaystyle\frac{\theta_{n+1}-\theta_{n}}{\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{2}[f(t_{n},\theta_{n})+f(t_{n+1},\theta_{n+1})]$
$\displaystyle\equiv$ $\displaystyle\textbf{avg. slope over }[t_{n},t_{n+1}]$

$\Longrightarrow\theta_{n+1}=\theta_{n}+\frac{1}{2}\alpha(f_{n}+f_{n+1})+\mathcal{O}(\alpha^{3})$

In this way, the local T.E. is reduced to $\mathcal{O}(\alpha^{3})$. To obtain
$f_{n+1}$ with $\theta_{n+1}$ yet to be determined, estimate $\theta_{n+1}$
for the next time period using the Explicit Euler at every step. Predictor:

$\bar{\theta}_{n+1}=\theta_{n}+\alpha f_{n}+\mathcal{O}(\alpha^{2})$ (12)

(As a preliminary predicted value of $\theta$ at $t_{n+1}$)
Corrector:

$\theta_{n+1}=\theta_{n}+\frac{\alpha}{2}[f(t_{n},\theta_{n})+f(t_{n+1},\bar{\theta}_{n+1})]$
(13)

Qualitative Error Analysis: When calculating $f_{n+1}$ in the corrector step,
$\theta_{n+1}$ is replaced by $\bar{\theta}_{n+1}$ The modified Euler
(Predictor-Corrector) method has $\mathcal{O}(\alpha^{2})$ error.

To denote the linear value function’s parameterization by $\theta$, it will be
written as $h^{\theta}$. Note:

$h^{\theta}(x)=\phi^{T}(x)\theta$ (14)

After substituting,

$f(\theta_{n})=[(\zeta_{n}(-h^{\theta_{n}}(X_{n})+c(X_{n})+\beta
h^{\theta_{n}}(X_{n+1})))]|_{\theta=\theta_{n}}$ (15)

and

$f(\bar{\theta}_{n+1})=[\zeta_{n}(-h^{\bar{\theta}_{n+1}}(X_{n})+c(X_{n})+\beta
h^{\bar{\theta}_{n+1}}(X_{n+1}))].$ (16)

The new improved TD learning algorithm called Causal PCTD(0) causal is:

$\theta_{n+1}=\theta_{n}+\alpha_{n}\frac{1}{2}[f(t_{n},\theta_{n})+f(t_{n+1},\bar{\theta}_{n+1})].$
(17)

Futhermore the predictor is the same as in 12.

From the definitions we can arrive at after substitution and reduction to the
form

$\theta_{n+1}=\theta_{n}+A\theta_{n}+b,$ (18)

where

$A=\alpha\zeta_{n}[\beta\phi^{T}(X_{n+1})-\phi^{T}(X_{n})]+\frac{\alpha^{2}}{2}[\zeta_{n}[\beta\phi^{T}(X_{n+1})-\phi^{T}(X_{n})]]^{2}$

and

$b=\alpha
c(X_{n})\zeta_{n}+\frac{\alpha^{2}}{2}\zeta_{n}(\beta\phi^{T}(X_{n+1})-\phi^{T}(X_{n}))\zeta_{n}c(X_{n}).$

After expanding, we have six terms remaining to analyze in the matrix $A$.
Splitting $A$ into two matrices, it is easy to see what terms are involved:
$A$ = $A_{1}$ \+ $A_{2}.$

$A_{1}=\alpha\beta\phi(X_{n})\phi^{T}(X_{n+1})-\alpha\phi(X_{n})\phi^{T}(X_{n})$

and

$A_{2}=\frac{\alpha^{2}\beta^{2}}{2}H^{2}-\frac{\alpha^{2}\beta}{2}GH-\frac{\alpha^{2}\beta}{2}HG+\frac{\alpha^{2}}{2}G^{2}$

where

$G=\phi(X_{n})\phi^{T}(X_{n+1})$ and $H=\phi(X_{n})\phi^{T}(X_{n}).$

$A_{1}$’s contribution to the dynamics are exactly what TD($0$) learning
contributes to the learning algorithm. The remaining four terms come from the
application of the Predictor-Corrector method to the ODE and are the
additional terms contributed by Causal PCTD(0).

###  3.2 Non-causal PCTD(0)

Note that the time support for $\theta$ and $X$ are not exactly the same nor
do they have to share the same time support entirely, only be relatively
close. For example, at time $n+1$, $\theta_{n+1}$ will depend on both
$(X_{n},X_{n+1})$ the previous state and the present state. However, a minimal
non-causal filter may improve the updates to $\theta_{n+1}$ drastically by
using information from the past, present and future in the update of the
parameter. I say the update is minimally non-causal in the sense that it only
depends on the state information one time-step into the future.

Now, I present the Non-causal PCTD(0) derivation taking this extra information
into account so we can see the final update form:

Starting from 13 with definitions 12 and 15 the same, we redefine 16 to be:

$f(\bar{\theta}_{n+1})=[\zeta_{n+1}(-h^{\bar{\theta}_{n+1}}(X_{n+1})+c(X_{n+1})+\beta
h^{\bar{\theta}_{n+1}}(X_{n+2}))].$ (19)

From these definitions, it is natural to see that a better update may result
from using the triple: $(X_{n},X_{n+1},X_{n+2})$ in this form of the update,
rather than using only past and present information. Practically, to apply
this algorithm in real-time, you would only have to wait one time-step to use
the updated improved parameter.

After substituting and reducing, we can again arrive at the same standard
update form as in 18,

where $A=A_{1}+A_{2}.$ For Non-causal PCTD(0), I split $A$ into two parts.
$A_{1}$ will hold the causal terms, and $A_{2}$ will hold the non-causal
terms.

$A_{1}=\frac{\alpha}{2}[-\zeta_{n}\phi^{T}_{n}-\zeta_{n+1}\phi^{T}_{n+1}+\beta\zeta_{n}\phi^{T}_{n+1}+\alpha\zeta_{n+1}\phi^{T}_{n+1}\zeta_{n}\phi^{T}_{n}-\beta\alpha\zeta_{n}\phi^{T}_{n+1}\zeta_{n}\phi^{T}_{n+1}]$

and

$A_{2}=\frac{\alpha}{2}[\beta\zeta_{n+1}\phi^{T}_{n+2}-\alpha\beta\zeta_{n+1}\phi^{T}_{n+2}\zeta_{n}\phi^{T}_{n}+\alpha\beta^{2}\zeta_{n+1}\phi^{T}_{n+2}\zeta_{n}\phi^{T}_{n+1}]$

and

$b=\frac{\alpha}{2}[\zeta_{n}c_{n}+\zeta_{n+1}c_{n+1}]-\frac{\alpha^{2}}{2}\zeta_{n+1}\phi^{T}_{n+1}\zeta_{n}c_{n}+\frac{\alpha^{2}\beta}{2}\zeta_{n+1}\phi^{T}_{n+2}\zeta_{n}c_{n}.$

Comments: this leads to a direct eigenvalue test for stability. Using the law
of large numbers, the exact parameter $\theta^{*}$ can be computed only
assuming a full rank condition on the basis functions.

###  3.3 Extension to PCTD($\lambda$)

Futhermore, we can generalize these results to a generalized
TD($\lambda$)-like algorithm, which I call PCTD($\lambda$). To obtain the
PCTD($\lambda$) general versions of the new RL algorithms, one needs to
substitute for the eligibility vector
$\zeta_{n}=\phi_{n}+\beta\lambda\zeta_{n-1}$. This generalizes the algorithms
further for any $\lambda\in[0,1)$, and stability is guaranteed under the same
mild set of assumptions.

##  4 Prediction Task

Figure 1: Optimal Value Function $h^{*}(x)$ from sparse costs in the
Prediction Task

For testing PCTD(0) varieties vs TD(0) learning I will be implementing an
infinite state MDP random walk on 7 states. Two states are "terminal": the
leftmost and rightmost states, i.e., once the agent reaches either x=0 or x=6
the agent directly returns to state x=3 in the next step with probability one.
All independent games start in state x=3, and with probability $\frac{1}{2}$
the agent transitions left or right. The learning environment is inspired by
the random walk example from the introductory RL book [8] Read Ch.6 on
Temporal Difference Learning. However, there it is formulated as episodic,
which means the game ends once the agent reaches either x=0, or x=6. Then, in
the next episode, the agent starts at x=3 whereas in my MDP, the play is
unending. In the simulation section, the transient performance of the
algorithms will be compared. Additionally, comparisons will be made between
using a decayed learning rate and a fixed learning rate for both Causal and
Non-causal PCTD(0) algorithms. The true optimal value function in 1 was
obtained by running value iteration until the Bellman error reached zero with
a $\beta=0.95$ discount factor.

##  5 Simulation Results

Figure 2: Causal PCTD(0) for N=100 independent trials for decayed learning
rates Figure 3: Non-Causal PCTD(0) for N=100 independent trials for decayed
learning rates

In this section we will test the performance in practice for the two versions
of PCTD(0), both Causal and Non-causal to determine which leads to lower MSE
during early learning for the value-function prediction task. For TD(0) and
both versions of PCTD(0) a band of learning rates with decayed learning rates
was found which resulted in stable behavior of the iterations. Below an
initial learning rate of $1e^{-5}$ for the decayed learning rate sequences,
stability was not guaranteed for any of the algorithms using our fixed
polynomial basis of linear, quadratic, and cubic features of the form
$[x,\frac{1}{2}x^{2},\frac{1}{3}x^{3}]$. Hence $d=3$ was the number of
features.

Within the range of $\alpha_{0}=1e^{-5}-1e^{-6}$ performance improved
monotonically with decreasing initial learning rates. Graphed in both figures
2 and 3 is the MSE averaged over all states for the linearly parameterized
value function averaged over 100 independent trials each of time extent
$T=100$ steps for both versions of PCTD(0) against TD(0). The performance of
TD(0) vs. Causal PCTD(0) and TD(0) vs. Non-causal PCTD(0) on the sparse cost
Prediction Task is illustrated. The causal version of PCTD(0) performed just
as well as TD(0) on the prediction task whereas the Non-causal version of
PCTD(0) offered a noticeable advantage at early learning of the optimal value
function. These experiments were repeated using fixed learning rates, and
essentially the same results were obtained. Lastly, I will show how the
PCTD(0) and TD(0) value functions’ estimates compare to the optimal value
function by showing confidence bounds using $N=100$ independent experiments
with $g=0.01$, and $n0=200$ for the same Horizon $T=100$ to show clearly the
early learning advantages of Non-causal PCTD(0) over TD(0) version at learning
the optimal-value function $h^{*}(x)$.

Figure 4: $h(x)$ over $N=100$ independent runs for $T=100$

##  6 Conclusion and Future Work

The Causal PCTD(0) performs slightly worse than TD(0) across the band of
stable learning rates. Meanwhile, the Non-causal PCTD(0) variant appears to
have faster transient learning capabilities for learning the optimal value
function during the early stages of learning as compared to the TD(0)
algorithm in terms of utilizing the past, present, and future information in
its updates. It should be noted that if the stability of the algorithms could
be made to improve with larger initial learning rates, the Non-causal
PCTD(0)’s MSE would most likely decrease even faster than standard TD(0)’s MSE
based on being more accurate in the parameter estimate by a factor of
$\mathcal{O}(\alpha)$. To improve Non-causal PCTD(0)’s faster learning gains
more, different types of feature vector’s may show other advantages for this
algorithm over standard TD(0) that increases its effectiveness than those
features employed in the simulation section of this paper.

Furthermore, the improvement in MSE of Non-causal PCTD(0) over Causal PCTD(0)
suggests that Non-causal PCTD(0)’s derivation leads to the correct form of the
update based on the Predictor-Corrector numerical ODE method being applied to
the original ode 5. This implies Non-causal PCTD(0) has the correct
interpretation of the time support for $\theta$ and $X$ involved in the
expression on the right hand side of 19. Non-causal PCTD(0)’s form of
algorithm seems to be able to admit further extensions that will be explored
in future work such as incorporating additional future and past feature terms
for every update of the parameter $\theta_{k+1}$. Additional numerical methods
of higher order will also lead to further RL algorithm developments. In a
future paper, I will seek to tie together different temporal difference
learning algorithm around a common thread of improving numerical methods, such
as the Runge-Kutta methods or even Richardson extrapolation. This will create
a new family of temporal difference learning algorithms. Lastly, the PCTD
algorithms developed in this paper will be studied in a deterministic setting
later where the noise comes from deterministic probing signals, such as well
structured mixtures of sine waves. In this future paper, the state space will
be continuous, and the control space will be as well. The proof of convergence
developed in this future work may shed light on the probabilistic setting and
lead to alternative proofs of stability and convergence rates.

##  7 Broader Impact

Because this research is largely theoretical in nature, no material
applications are foreseen. Hence, this section is not applicable for the
current work.

## References

  * [1] John N Tsitsiklis and Benjamin Van Roy.  An analysis of temporal-difference learning with function approximation.  IEEE transactions on automatic control, 42(5):674–690, 1997.
  * [2] Vivek S Borkar and Sean P Meyn.  The ode method for convergence of stochastic approximation and reinforcement learning.  SIAM Journal on Control and Optimization, 38(2):447–469, 2000.
  * [3] Francisco S Melo, Sean P Meyn, and M Isabel Ribeiro.  An analysis of reinforcement learning with function approximation.  In Proceedings of the 25th international conference on Machine learning, pages 664–671, 2008.
  * [4] Joshua Achiam, Ethan Knight, and Pieter Abbeel.  Towards characterizing divergence in deep q-learning.  arXiv preprint arXiv:1903.08894, 2019.
  * [5] Jalaj Bhandari, Daniel Russo, and Raghav Singal.  A finite time analysis of temporal difference learning with linear function approximation.  In Conference On Learning Theory, pages 1691–1692. PMLR, 2018.
  * [6] Shaofeng Zou, Tengyu Xu, and Yingbin Liang.  Finite-sample analysis for sarsa with linear function approximation.  arXiv preprint arXiv:1902.02234, 2019.
  * [7] AM Devraj, A Bušic, and S Meyn.  Fundamental design principles for reinforcement learning algorithms.  Handbook on Reinforcement Learning and Control. Springer, 2020.
  * [8] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.