# Turbulent Transport of Dust Particles in Protostellar Disks: The Effect of
Upstream Diffusion

Tingtao Zhou(周廷弢)  Department of Mechanical and Civil Engineering, California
Institute of Technology, Pasadena, CA 91125, USA  Hong-Ping Deng(邓洪平)
Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai
200030, China  Yi-Xian Chen(陈逸贤)  Department of Astrophysics, Princeton
University  Douglas N. C. Lin(林潮)  Department of Astronomy & Astrophysics,
University of California, Santa Cruz, CA 95064, USA  Institute for Advanced
Studies, Tsinghua University, Beijing 100086, China

###### Abstract

We study the long-term radial transport of micron to mm-size grain in
protostellar disks (PSDs) based on diffusion and viscosity coefficients
measured from 3D global stratified-disk simulations with a Lagrangian
hydrodynamic method. While gas-drag tend to transport dust species radially
inwards, stochastic diffusion can spread a considerable fraction of dust
radially outwards (upstream) depending on the nature of turbulence. In
gravitationally unstable disks, we measure a high radial diffusion coefficient
$D_{\rm r}\sim H^{2}\Omega$ with little dependence on altitude. This leads to
strong and vertically homogeneous upstream diffusion in early PSDs. In the
solar nebula, the robust upstream diffusion of $\mu$m to mm size grains not
only efficiently transports highly refractory $\mu$m-size grains (such as
those identified in the samples of comet 81P/Wild 2) from their regions of
formation inside the snow line out to the Kuiper Belt, but can also spread mm-
size CAI formed in the stellar proximity to distances where they can be
assimilated into chondritic meteorites. In disks dominated by
magnetorotational instability (MRI), the upstream diffusion effect is
generally milder, with a separating feature due to diffusion being stronger in
the surface layer than the midplane. This variation becomes much more
pronounced if we additionally consider a quiescent midplane with lower
turbulence and larger characteristic dust size due to non-ideal MHD effects.
This segregation scenario helps to account for dichotomy of two dust
populations’ spatial distribution as observed in scattered light and ALMA
images.

##  1 Introduction

The omnipresence of protostellar disks (PSDs) (Hartmann, 1998) provides birth
sites and incubators for emergence of planets (Cameron, 1978; Hayashi et al.,
1985; Lin & Papaloizou, 1985; Ida & Lin, 2004). The temperature distribution
in these disks inferred from their observed spectral energy distribution
(Adams et al., 1987; Hartmann, 1998) matches well with those from theoretical
thermal equilibrium models where the cooling rate from the disk surface is
balanced by both intrinsic viscous dissipation and stellar irradiation (Chiang
& Goldreich, 1997; Garaud & Lin, 2007). This provides a natural division for
the condensation sequence of terrestrial, asteroids, gas and ice giant planets
(Ruden & Lin, 1986).

In such a setting, refractory and volatile grains are expected to condense in
the stellar proximity (within $\sim 0.1$ au) and outside the snow line
(several aus) respectively. Yet, $\mu$m-size crystalline refractory grains are
found in abundance among the Stardust return samples from comet 81P/Wild 2
(Brownlee et al., 2006). Their presence is also inferred from the IR spectra
of other comets (Wooden et al., 2007) as well as extended regions (out to
$\sim 10$ au) of PSDs (Bouwman et al., 2003, 2008; van Boekel et al., 2005;
Watson et al., 2009).

The crystalline structure requires the melt-down and re-condensation of
amorphous silicate grains at temperatures $\gtrsim 1000$K, depending on
composition (Bockelée-Morvan et al., 2002), while their implantation into
comets occurs below the icy condensation temperature $\lesssim 160$K.
Therefore, very efficient mechanisms need to be invoked for these grains to be
transported outwards after re-condensation in the stellar proximity.
Explanations for this thermal paradox include large scale diffusion (Morfill &
Voelk, 1984; Clarke & Pringle, 1988), spreading of the solar nebula (Safronov
& Ruzmaikina, 1985; Ciesla, 2007), meridional circulation (Urpin, 1984; Kley &
Lin, 1992), radiation-pressure driven outflow in the disk surface (Takeuchi &
Lin, 2003), fluctuating X-wind (Shu et al., 2000, 2001; Shang et al., 2000),
and shock heating (Weidenschilling et al., 1998; Wooden et al., 2007; Desch et
al., 2012; Gong et al., 2019).

Optical and NIR maps of scattered light also indicate that in certain PSDs,
the micron-size particles have larger radial extent and thickness than those
of larger mm-size particles inferred from ALMA images of reprocessed radiation
(Villenave et al., 2020; Benisty et al., 2022). This suggests that dust of
different characteristic size, occupying different vertical layers of the
disk, may have distinct capabilities of radial transport.

In this paper, we revisit the diffusion hypothesis. The structure and
evolution of PSDs including the solar nebula are regulated by turbulent
viscosity. Some causes of turbulence include the magneto-rotational
instability (MRI) (Balbus & Hawley, 1991) and gravitational instability (GI)
(Safronov, 1960; Toomre, 1964) in non-magnetised (see, e.g., Durisen et al.,
2007) and magnetised disks (Riols & Latter, 2019; Deng et al., 2020). The
chaotic turbulent motion is correlated and leads to Reynold and Maxwell stress
and effective torque which transport gas angular momentum outwards, while
inducing mass spreading (Pringle, 1981; Balbus & Hawley, 1998) and diffusion
of passive contaminants including small grains (Morfill & Voelk, 1984; Clarke
& Pringle, 1988) 111An analogous turbulent-diffusion mechanism has also been
invoked to account for the shallow abundance gradient in disk galaxies
(Clarke, 1989; Yang & Krumholz, 2012; Krumholz et al., 2018). A key issue is
whether turbulent diffusion of crystalline grains or heavy-element tracers can
diffuse against the inward gas flow and spread widely over extended regions of
the disk. We present in §2 three sets of numerical simulation for MRI, GI in
the absence of magnetic fields (GI-HD), and GI of disks with magnetic fields
(GI-MHD). Based on measurements from state-of-the-art simulations in §3, we
construct simplified 1D models of layered accretion disks and examine the
global diffusion process of dust elements in turbulent disks in §4, focusing
on the evolution of their concentration distribution. In §5, we summarize our
results and discuss their implications.

##  2 Simulations and analysis

The simulations of GI-HD, MRI and GI-MHD turbulence are drawn from Deng et al.
(2020) dedicated to studying the interplay between magnetic fields and spiral
density waves. The simulations were carried out with a Godunov-type Lagrangian
method, the meshless finite mass (MFM) method in the GIZMO code 222The code is
publicly available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html.
(Hopkins, 2015) which shows excellent conservation property in GI disk
simulations (Deng et al., 2017). At sufficiently high resolution, MFM can also
capture the MRI and the ensuing turbulence showing good agreement with grid-
code simulations (Deng et al., 2019, 2020). In all the simulations, the disk
mass is about 0.07 the stellar mass and the disk is resolved by $>35$M
particles/elements, i.e., a Jupiter mass object is resolved by half a million
computational elements, among the highest resolution GI disk simulations.333We
note that thanks to the adaptive nature of Lagrangian methods, even a 1M
particle simulation produces converged results comparable to some ultra-high
resolution grid-code simulations (up to 2 billion effective cells) in GI disk
modeling. See the Wengen code comparison project test4, including one
Lagrangian code, Gasoline and four grid codes (some of which with mesh
refinement), Enzo, FLASH, CHYMERA and IU code at
http://users.camk.edu.pl/gawrysz/test4/. The initial surface density is
inversely proportional to the radius spanning 5–25 au. The code units are 1
au, 1 solar mass and $G=1$.

Figure 1: The mid-plane volume density in quasi-equilibrium turbulence of MRI,
GI-MHD and GI-HD simulations. The colors, from black to yellow, show the
density ranging from $10^{-15}$ to $10^{-9}$ g/cm-3 in a logarithmic scale.

In the following analysis, the MRI and GI-MHD turbulence corresponds to the
model MRI-hr and grvmhd1 simulation. We then turn off the magnetic fields in
the grvmhd1 model and run it for 400 additional years and this simulation is
named grvhd1 here. We show representative mid-plane density maps 444This is
similar to the mid-row of figure 3 in Deng et al. (2020) and the high density
inner edge in that figure is due to artificial effects during plotting. in the
quasi-steady state in Figure 1. The MRI turbulence is mild with small density
contrast while GI turbulence, regardless of the presence of magnetic fields,
show strong density fluctuations. The flocculent spirals in the grvmhd1 model
become more coherent after turning off the magnetic fields and saturating to
grvhd1 resembling previous hydrodynamic simulations of gravitational
instability (Deng et al., 2017; Béthune et al., 2021). This in turn confirmed
the important back reaction of magnetic fields on the spiral density wave
morphology (for details, see Deng et al., 2020).

In the well-coupled regime, dusts are tied to the Lagrangian fluid elements so
that we can directly analyse the fluid elements’ trajectories to uncover
particle motions. This is different from injecting and analysing tracers of
different sizes into either Lagrangian (Rice et al., 2004) or Eulerian
hydrodynamic simulations (Haghighipour & Boss, 2003; Boley & Durisen, 2010;
Zhu et al., 2015; Riols et al., 2020b; Baehr & Zhu, 2021; Hu & Bai, 2021). To
our knowledge, this work amounts to be the first study of fluid tracers in
global Lagrangian simulations of magnetised/self-gravitating accretion discs.

Figure 2: The trajectories of particles in the MRI simulation ($R-z$ plane)
initially at 12.5 au within five local orbits.

We trace fluid elements in the quasi-equilibrium turbulence for 250 years (two
rotation periods at 25 au), which is sufficient to identify the particles
diffusion coefficients. We choose particles in small annuli centered at 10 au
and 12.5 au with a width of 0.001 au. In each annulus, there are more than
2000 particles. For example, we follow the trajectories of fluid elements in a
narrow annulus of width 0.001 au in the MRI-hr simulation and plot the
trajectories of 100 randomly chosen particles out of the $>$2000 particles in
Figure 2.

The three simulations have distinct thermophysics and vertical structures
(see, Figure 2 of Deng et al., 2020) and thus, for ease in comparison, we need
to normalise the statistics to an appropriate length scale. Specifically, the
MRI-hr simulation is adiabatic and the accretion heats up the disk gradually;
the grvmhd1 and grvhd1 saturates to a state where the heating is balanced by a
prescribed cooling at a timescale of the orbital period (Deng et al., 2020).
In addition, the disk vertical structure is strongly fluctuating the GI
simulations so that the traditional definition of the vertical scale-height is
not adequate.

We adopt the dynamic scale-height ($H$) definition of Ogilvie (2018) where $H$
is the square root of the second vertical moment of the density field, i.e.,

$H=\sqrt{\frac{\int\rho z^{2}dz}{\int\rho dz}}.$ (1)

It naturally reduce to the ratio between the speed of the sound ($c_{s}$) and
the local rotational frequency ($\Omega$) for a vertical Gaussian profile. We
take the time (250 years) and azimuthal average the dynamic scale-heights for
the three simulations and show them in Figure 3.

Figure 3: The averaged dynamic scale height in the quasi-equilibrium
turbulence.

The adiabatic MRI-hr simulation is hot with a typical $H/R=0.1$ (due to the
initial condition). The inner disc is relatively hotter due to a faster
development (small $\Omega^{-1}$) of MRI turbulence than the outer region in
the adiabatic simulation. The grvmhd1 and grvhd1 have significantly
thinner/colder discs with $H/R$ around 0.03 and 0.02. In further analysis, we
normalise displacements to $H$ to enable comparison across different models
and with previous grid-code simulations (Zhu et al., 2015; Riols et al.,
2020b). However, for simplicity, a thermal profile with constant $H/R$ is
adopted for GI, GI-MHD and MRI in the simplified 1-D dust transport model in
§4.

Runaway types of super-diffusion processes have been reported for certain
highly supersonic turbulence, in which there is no well-defined diffusion
coefficient (Colbrook et al., 2017). In contrast, our global simulations
corroborate with previous local simulations (Zhu et al., 2015; Shi et al.,
2016) to show that GI and MRI, as turbulence typically on the sonic level
(characteristic length scale comparable to $H$), induce regular
Gaussian/Brownian diffusion of passive tracers. We can measure diffusion
coefficients through the displacement of these tracers (see Figure 2), as
elaborated in §3. Our approach is similar to the method applied in grid-code
simulations to measure the diffusion coefficients of Lagrangian particles with
finite Stokes number $T_{S}$. As we mentioned, since particles represent gas
itself in our simulations, we directly measure the diffusion coefficients of
passive particles well-coupled with gas (the limit of small $T_{S}\rightarrow
0$), however, diffusivity measured in this limit are proven to be valid for
small dust Stokes number up to $T_{S}\lesssim 0.1$ (Zhu et al., 2015; Shi et
al., 2016), consistent with analytical scaling from Youdin & Lithwick (2007).

The measured relevant diffusion coefficient discussed below are summarised in
Table 1. Although $D\sim\alpha H^{2}\Omega$ is often applied in grid-based
dust-gas hydrodynamic simulations (Zhu et al., 2011, 2012; Chen et al., 2021),
we show that this approximation can be off by up to an order of magnitude.
$D_{r}$ is the overall radial diffusion coefficient and $D_{z}$ is the overall
vertical diffusion coefficient. To explore possibility of layered radial
diffusion, we also divide disks into the disk midplane within $|z|<H$ and the
disk atmosphere $|z|>H$, and measure the radial diffusion coefficients $D_{\rm
mid}$ and $D_{\rm sur}$ for each layer. We generally expect that with a large
$D_{z}$, the vertical mixing of dust/gas at different altitudes becomes very
efficient such that $D_{\rm mid}\sim D_{\rm sur}$ whereas they may differ more
significantly when $D_{z}$ is small.

Table 1: dust diffusion coefficients  Run | Physics | $D_{r}$ | $D_{\rm mid}$ | $D_{\rm sur}$ | $D_{z}$ | $\alpha$
---|---|---|---|---|---|---
grvhd1 | GI | 2. | 1.9 | 2.1 | 0.2 | 0.1
grvmhd1 | GI-MHD | 0.5 | 0.49 | 0.5 | 0.1 | 0.23
MRI-hr | MHD | 0.045 | 0.033 | 0.07 | 0.013 | 0.02
Table 1:  The diffusion coefficients (given in a unit of $H^{2}\Omega$) by
linear regression of the dust mean squared displacement (MSD). Note, for self-
gravitating disk, only the early stage vertical MSD evolution reflects the
fast vertical motion (see text). Here $D_{\rm mid}$ and $D_{\rm sur}$ are
measured dust radial diffusion coefficients for particles with initial
position $|z|<H$ and $|z|>H$, respectively. In the unit of $H^{2}\Omega$,
$\nu$ corresponds to the value of $\alpha$ in the widely adopted $\alpha$
prescription for effective viscosity Shakura & Sunyaev (1973), and $D/\nu$ are
referred to as Schmidt numbers. We see that while radial $D_{r}/\nu$ span a
wide range, the vertical $D_{z}/\nu$ are always of order-unity. Schmidt
numbers are also sometimes used to refer to the ratio between dust diffusivity
and gas diffusivity, which we assume to be 1 in this study for small $T_{S}$
(Youdin & Lithwick, 2007). We will not be using this term in the rest of this
paper to avoid confusion.

##  3 Dust diffusion in real turbulence

###  3.1 gravitationally unstable disks

Young circumstellar disks are relatively massive to their hosts and can be
unstable due to gas self-gravity forming spiral density waves (see, e.g.
reviews by Durisen et al., 2007; Kratter & Lodato, 2016). The spiral density
waves are critical for angular momentum transport and thus disk accretion. The
spiral density wave also excites strong vertical circulation and likely powers
a larges scale dynamo amplifying any weak magnetic fields to near thermal
strength (Riols & Latter, 2019; Deng et al., 2020; Löhnert & Peeters, 2022).
The strong fields facilitates accretion in early discs significantly and may
impact the population of planet formed via disk instability (Deng et al.,
2021).

Figure 4: The vertical motion of randomly chosen particles in the annuls
centered at 12.5 au. .

In an gravitational unstable disks, the vertical circulation (Boley & Durisen,
2006; Deng et al., 2020) besides spirals efficiently mixes small particles at
different altitudes above the disk mid-plane in local shearing box simulations
(Riols et al., 2020b; Baehr & Zhu, 2021). In order to characterize the
difference in turbulent diffusion in disks regulated by MRI and GI, we plot
the vertical motion of 15 randomly chosen particles in the annulus centered at
12.5 au in Figure 4. Particles in the MRI turbulence show minor vertical
motion (model MRI-hr), impliying a slow vertical diffusion. This weak coupling
also enables particles in the disk mid-plane and surface to diffuse at
different rates. In contrast, the vertical motion over a dynamical time scale
extend to the scale of $H$ in GI turbulence (models grvmhd1 and grvhd1). This
random motion leads to efficient particle exchanges in the vertical direction
and inhibits a significant differential radial diffusion between the disk mid-
plane and surface.

Figure 5: The evolution (given in units of the local dynamical time) of the
mean squared displacement along the vertical direction for particles in annuli
centered at 12.5 au initially. Here the MSD is normalised to the dynamical
scale-height at 12.5 au. .

The vertical diffusion process in grvmhd1 and grvhd1 can be quantified by the
vertical mean-squared displacement (MSD). In Figure 5, we plot the averaged
vertical MSD for particles in the disk mid-plane ($|z|<H$) and surface
($|z|>H$) versus the all-particle average MSD around 12.5 au, as functions of
diffusing time. In both grvmhd1 and grvhd1, the vertical MSD experience a fast
increase followed by a slower but steady growth. We measure a mean vertical
diffusion rate of $D_{z}=0.2H^{2}\Omega$ and $0.1H^{2}\Omega$ for the grvhd1
and grvmhd1, comparable to their viscosity $\nu=\alpha
H^{2}\Omega=0.1H^{2}\Omega$ and $0.23H^{2}\Omega$ (Deng et al., 2020), also
summarised in Table 1. After the particles in the narrow annulus (initially
confined) spreads out along a large scale eddy, the growth of vertical MSD
become slow in grvmhd1 or negligible in grvhd1. The initial fast increase in
MSD is typical for high altitude particles ($|z|>H$), suggesting faster flow
in the disk surface layer (dashed lines). However, over a few dynamical
timescales the MSD growth quickly saturate towards a lower pace that converges
with the average diffusion rate (dotted lines), since particles’ movement is
hindered by their rapid circulation across lower altitudes ($|z|\sim H$),
where diffusion is much slower. In principle, we can measure somewhat
different vertical diffusion coefficients $D_{z,mid}\ll D_{z,sur}$ from the
initial trend of dashed and dotted lines in Figure 5. Nevertheless, since
vertical mixing is efficient suggesting particles commute across different
layers frequently (as shown in Figure 4), it would be more meaningful to
simply apply the average diffusion coefficient in following discussions for
disks with gravitational instability.

Moreover, due to efficient vertical mixing, the radial MSD shows little
altitude dependence $D_{r}\approx D_{\rm mid}\approx D_{\rm sur}$. As a
result, we only plot the average all-particle MSD (used to measure $D_{r}$)
for all particles centered at 10 au and 12.5 au in Figure 6. At the two
different radii, the MSD evolves similarly so that our results apply to
general GI and GI-MHD turbulence. Similar to the vertical MSD evolution, the
radial MSD also increase rapidly at the early stage when the particles are
moved around by the vertical circulation. However, particles can hop onto new
eddies and be radially transported over large distances. As a result, the
radial MSD continues to grow. In the grvmhd1 model, the radial diffusion rate
is $D_{r}=0.5H^{2}\Omega$, twice the viscosity. In the grvhd1 model, the
radial diffusion rate is $2H^{2}\Omega$ which is much larger than its
viscosity $0.1H^{2}\Omega$, as summarised in the first 2 rows of Table 1.

Figure 6: The evolution (given in units of the local dynamical time) of the
mean squared displacement for particles. Here MSD is normalised to the
dynamical scale-height at 12.5 au or 10 au. .

###  3.2 MRI turbulence

Magnetised accretion disks are prone to the MRI provided sufficient ionization
(Balbus & Hawley, 1991), which further leads to vigorous turbulence that
transports angular momentum outwards enabling gas accretion. In ideal MHD
simulations, MRI turbulence has a typical effective viscosity $\nu=\alpha
H^{2}\Omega$, with $\alpha\sim 0.02$ in grid-code (Simon et al., 2012;
Beckwith et al., 2011) as well as Lagrangian simulations (Deng et al., 2019,
2020; Wissing et al., 2022). However, realistic PSDs are poorly ionized so
that only the surface layer may be MRI active and turbulent or the whole
vertical extent is laminar with winds launched at high altitude (Gammie, 1996;
Okuzumi & Hirose, 2011; Bai & Stone, 2013; Riols et al., 2020a). In that case,
the constrast between $D_{\rm sur}$ and $D_{\rm mid}$ would be even more
pronounced. Our ideal MHD (active at all altitudes) model is admittedly
oversimplified and it offers the most conservative illustration for the
layered dust diffusion tendency in MRI turbulence.

Figure 7: The evolution (given in units of the local dynamical time) of the
mean squared displacement for particles initially at 12.5 au for the MRI-hr
model. Here MSD is normalised to the dynamical scale-height at 12.5 au. .

In Figure 7, we plot the MSD evolution for all particles in the annulus
centered at 12.5 au and particles initially in the disk body with $|z|<H$ and
disk surface with $|z|>H$ separately. There are some 1700 and 600 particles in
the disk body and surface respectively. By the end of the simulation about 200
particles crossed the $z=\pm H$ boundary. If we restrict our analysis to those
particles confined in the same (mid-plane or surface) portion of the disk
(see, e.g. Zhu et al., 2015), we would find almost identical results to those
in Figure 7. The particles show a mild uniform vertical diffusion rate about
$D_{z}=0.013H^{2}\Omega$ which is smaller than the viscous $\alpha$ parameter
in units of $H^{2}\Omega$. A similar analysis for the particles in an annulus
centered at 10 au reveals a vertical diffusion rate about $0.017H^{2}\Omega$
close to $\alpha$. However, in Table 1, we restrict ourselves to the reference
value measured at 12.5 au since the latter is less likely affected by the
inner boundary (see Figure 3).

In the MRI context, the diffusion rates start to show dependence on the
vertical altitude. Particles in the disk mid-plane diffuse at a rate $D_{\rm
mid}=0.033H^{2}\Omega$ whereas those at the disk surface diffuse at a faster
rate of $D_{\rm sur}=0.07H^{2}\Omega$. This introduces a factor of two
difference in their characteristic radial diffusion “velocity” $\sim D_{\rm
r}/R$. The analysis for particles centered at 10 au show similar radial
diffusion rates within a few percent deviation.

Since both $\alpha$ and $D_{z}\approx\alpha$ are small, vertical mixing in the
MRI-prone disks is much less efficient than the GI cases. Under such
circumstance, the considerable difference between $D_{\rm mid}<D_{\rm sur}$
introduces potentials for a layered diffusion process. Moreover, to go beyond
our ideal MHD simulation, in realistic disks the moderate turbulence rendering
$\alpha$ and $D_{\rm sur}$ as measured in our simulations may only be active
in the surface layer (Gammie, 1996), while residue $\alpha$ and $D_{\rm mid}$
is much smaller in the dead zone below the atmosphere.

To study the long-term evolution of grain in fore-mentioned environments, we
construct a 1D two-fluid disk model to further explore the long-term
implications of the diffusion process with or without layered structure.

##  4 Dust transport models

In the case of GI or GI-MHD turbulence, the radial diffusion coefficients are
similar and independent of the vertical position in the disk, as in Table 1.
We write down the evolution for the dust concentration $C$, assuming
cylindrical symmetry so that there will be no explicit dependence on the
azimuthal angle $\theta$.

$\displaystyle\Sigma\frac{\partial C}{\partial t}=$
$\displaystyle\frac{1}{R}\frac{\partial}{\partial
R}\left(RD\Sigma\frac{\partial C}{\partial R}\right)$
$\displaystyle-\frac{1}{R}\frac{\partial}{\partial R}\left(R\Sigma v_{\rm
d}C\right)+C\frac{1}{R}\frac{\partial}{\partial R}\left(R\Sigma v_{g}\right)$
(2) $\displaystyle\frac{\partial\Sigma}{\partial t}=$
$\displaystyle\frac{3}{R}\frac{\partial}{\partial
R}\left(R^{1/2}\frac{\partial}{\partial R}(\Sigma\nu R^{1/2})\right)$ (3)

where $\Sigma$ is the gas density, $C$ is the dust-to-gas surface density
ratio, and Equation 3 describes the evolution of the gas density (Ruden & Lin,
1986). In Equation 2 we generalize the formulation of Morfill & Voelk (1984);
Clarke & Pringle (1988), where we include the radial velocities of dust
$v_{d}$, which could differ from the radial velocity of the gas $v_{g}$. We
note that the dust back reaction is often ignored in the gas nebula evolution
(Equation 3) given the low dust concentration level of a few percent (see,
e.g., Liu et al., 2022). The gas advective velocity in the radial direction
follows (Lynden-Bell & Pringle, 1974; Armitage, 2019)

$\displaystyle v_{g}=$ $\displaystyle-\frac{3}{\Sigma
R^{1/2}}\frac{\partial}{\partial R}(\nu\Sigma R^{1/2})$ (4) $\displaystyle=$
$\displaystyle-\frac{3\nu}{R}\left(\frac{\partial\ln\Sigma}{\partial\ln
R}+\frac{\partial\ln\nu}{\partial\ln R}+\frac{1}{2}\right).$ (5)

For dust completely coupled with gas $v_{d}=v_{g}$, we can recover Equation
2.1.4 of Clarke & Pringle (1988). The motion of the dust particles are
influenced by the hydrodynamic drag by the disk gas. Considering dust
particles with finite inertia, their radial velocity (Takeuchi & Lin, 2002,
2005; Zhu et al., 2012) is

$\displaystyle v_{d}=\frac{v_{g}T_{S}^{-1}-\eta v_{K}}{T_{S}+1/T_{S}}$ (6)

where $\eta=-h^{2}d\ln P/d\ln R$ with $h=H/R$ is a measure for the sub-
Keplerian rotation of gas and $v_{K}=\sqrt{GM/R}$ is the Keplerian velocity.
The Stokes number, or dimensionless stopping time in the Epstein regime (e.g.
Birnstiel et al., 2010)

$\displaystyle T_{S}$ $\displaystyle=\dfrac{\pi}{2}\frac{\rho_{d}s}{\Sigma}$
(7)
$\displaystyle=0.0015\times\left(\frac{s}{1\mathrm{~{}mm}}\right)\left(\frac{\rho_{d}}{1\mathrm{~{}g}\mathrm{~{}cm}^{-3}}\right)\left(\frac{100\mathrm{~{}g}\mathrm{~{}cm}^{-2}}{\Sigma}\right),$

depends on individual dust grains’ radius $s$, internal physical density
$\rho_{d}$, and the gas volume density $\rho_{g}$.

In the case of MRI turbulence, the diffusion coefficients measured for the
middle plane and surface of the disk are different, and we correspondingly
treat the dust distribution as two stratified layers, with the middle plane
dust concentrations $C_{\rm
mid}(R,\theta,t)=\int_{0}^{H}\rho_{d}dz/\int_{0}^{H}\rho_{g}dz$ and surface
layer concentration $C_{\rm
sur}(R,\theta,t)=\int_{H}^{\infty}\rho_{d}dz/\int_{0}^{H}\rho_{g}dz$ where
$\rho_{g}$ is the gas density and $\rho_{d}$ is the dust density. Locally if
the two layers have different densities, there would be a linear mass transfer
rate $\lambda=D_{z}/H^{2}$ between the two layers. We can write down the
coupled two-layer diffusion problem as

$\displaystyle\phi\Sigma\frac{\partial C_{\rm mid}}{\partial t}=$
$\displaystyle\frac{1}{R}\frac{\partial}{\partial R}\left(RD_{\rm
mid}\phi\Sigma\frac{\partial C_{\rm mid}}{\partial
R}\right)-\frac{1}{R}\frac{\partial}{\partial R}\left(R\phi\Sigma v_{\rm
d,mid}C_{\rm mid}\right)+C_{\rm mid}\frac{1}{R}\frac{\partial}{\partial
R}\left(R\phi\Sigma v_{g}\right)+F_{z}$ (8)
$\displaystyle(1-\phi)\Sigma\frac{\partial C_{\rm sur}}{\partial t}=$
$\displaystyle\frac{1}{R}\frac{\partial}{\partial R}\left(RD_{\rm
sur}(1-\phi)\Sigma\frac{\partial C_{\rm sur}}{\partial
R}\right)-\frac{1}{R}\frac{\partial}{\partial R}\left(R(1-\phi)\Sigma v_{\rm
d,sur}C_{\rm sur}\right)+C_{\rm sur}\frac{1}{R}\frac{\partial}{\partial
R}\left(R(1-\phi)\Sigma v_{g}\right)-F_{z}$ (9)

where $\phi$ is the fraction of gas mass in the middle plane. $F_{z}=\Sigma
D_{z}/H^{2}((1-\phi)C_{\rm sur}-\phi C_{\rm mid})$ is the vertical diffusion
term. $v_{\rm d,mid}$ and $v_{\rm d,sur}$ are the radial velocity of dust in
the mid-plane and surface layer.

Figure 8:  Coupled evolution of gas and dust in an early disk with GI
turbulence. The gas density is initialized by a bell curve centered around
$R_{i}=10$, and it evolves according to Eq. 3, with an absorbing inner
boundary. $\alpha=0.1$, $D_{\rm r}=2~{}H^{2}\Omega$. (a) The gas density
evolution. (b) the gas drift velocity (black lines) and the dust velocities of
the middle plane (blue lines) and the surface layer (red lines). (c) The
evolution of dust concentration (blue lines). In all three plots the times are
indicated by linestyles: dotted for $t=100$, dashed for $t=1000$ and solid for
$t=10000$.

For simplicity, we take $\phi=0.7$ approximating a Gaussian vertical profile,
assuming a constant $h=H/R=0.05$, and therefore $\eta=(1-d\ln\Sigma d\ln
R)h^{2}$ and $\nu=\alpha h^{2}\sqrt{GMR}\propto R^{0.5}$.

We choose different $\alpha$ for different turbulence models according to
Table 1 unless otherwise specified. The inner boundary is set at a normalizing
length scale $R_{\rm in}=1.0$ with absorbing boundary condition for the dust
concentration $C(R=R_{\rm in})=0$. For the GI-dominated disk as shown in Fig.
8 the inner boundary is absorbing for the gas density $\Sigma(R=R_{\rm
in})=0$, while for the other cases presented below we use Dirichlet boundary
condition for gas density (to maintain a steady power law). The outer boundary
is set at $R_{\rm out}=130$ with a Dirichlet boundary condition for the gas
density and absorbing condition for the dust concentration(s). The time unit
is set to $\Omega_{\rm in}^{-1}$, the Keplerian dynamical timescale at $R_{\rm
in}$. Focusing on sub-mm size to micron size grains, we assume a
characteristic dust size corresponding to $T_{S}=0.001$, which is well-coupled
with gas (see Equation 7) in most of our models with vertically-homogeneous
$\alpha$, neglecting vertical differentiation of dust characteristic size due
to sedimentation (Takeuchi & Lin, 2003; Garaud et al., 2004). This can be
justified by considering the maximum dust Stokes number allowed by
fragmentation (Ormel & Cuzzi, 2007; Birnstiel et al., 2010; Chen et al., 2020;
Li et al., 2019):

$T_{\rm S,max}\approx\dfrac{u_{f}^{2}}{2\alpha c_{s}^{2}},$ (10)

for fragmentation velocity $u_{f}\sim 1$m/s, typical sound speed $c_{s}\sim
1$km/s and $\alpha>0.02$, the characteristic dust species’ size growth by
coagulation will be halted at $T_{S}$ smaller than a few 0.001, even in the
midplane. Because these dust species already couple fairly well with dust, the
diffusion of $\mu$m crystalline grains with even smaller $T_{s}$ will be even
better coupled to gas. We will discuss possible effects of having a dead zone
with much smaller midplane turbulence in §4.3, where characteristic dust
species may be allowed to grow to larger than mm-size.

In the following sections we apply one- and two-layered diffusion model to
PSDs at different stages of evolution with the dust diffusion parameter
extracted from the turbulence simulations as summarised in Table 1.

Figure 9:  Coupled evolution of gas and dust in a GI-MHD disk perturbed from
the steady-state power-law distribution. $\alpha=0.2$, $D_{\rm
r}=0.5~{}H^{2}\Omega$. (a) The gas density evolution and (b) the gas drift
velocity (black lines) and the dust velocity (blue lines) almost overlap
throughout the evolution, due to small grain sizes with $T_{S}=0.001$. (c) The
evolution of dust/gas ratio from solving Equations 2 and 3. In all three plots
the times are indicated by linestyles: dotted for $t=100$, dashed for $t=1000$
and solid for $t=10000$.  Figure 10: The maximum upstream fractions $P(R)$
(dimensionless) at $t=1\times 10^{4}~{}\Omega_{\rm in}^{-1}$, due to ideal MRI
turbulence, computed for the initial dust ring positioned at (a) $R_{0}=4$ and
(b) $R_{0}=10$ with $\alpha=0.02$, $D_{\rm mid}=0.033~{}H^{2}\Omega$, $D_{\rm
sur}=0.07~{}H^{2}\Omega$, $D_{z}=0.013~{}H^{2}\Omega$. $T_{\rm S,mid}=T_{\rm
S,sur}=0.001$. Solid lines represent the surface layer, and dashed lines are
for the middle plane layer. The inset shows the inferred disk concentration
images, with the left half representing the middle plane layer and the right
half representing the surface layer. (c) The steady state power-law gas
density (dotted black line, label on the left), gas radial velocity (black
solid line) and the dust radial drag velocity (red line, label on the right).
Figure 11: The maximum upstream fractions $P(R)$ at $t=1\times
10^{4}~{}\Omega_{\rm in}^{-1}$, due to non-ideal MRI turbulence, computed for
the initial dust ring positioned at (a) $R_{0}=4$ and (b) $R_{0}=10$ with
$\alpha_{\rm mid}=0.003$, $\alpha_{\rm sur}=0.02$, $D_{\rm
mid}=D_{z}=0.003~{}H^{2}\Omega$, $D_{\rm sur}=0.07~{}H^{2}\Omega$. $T_{\rm
S,mid}=T_{\rm S,sur}=0.001$. Solid lines represent the surface layer, and
dashed lines are for the middle plane layer. The inset shows the inferred disk
concentration images, with the left half representing the middle plane layer
and the right half representing the surface layer. (c) The steady state power-
law gas density (dotted black line), gas radial velocities of the middle plane
(dashed black line) and the surface layer (solid black line), and the dust
radial drag velocities of the middle plane (blue line) and the surface layer
(red line).  Figure 12: The maximum upstream fractions $P(R)$ at $t=1\times
10^{4}~{}\Omega_{\rm in}^{-1}$, due to non-ideal MRI turbulence and
differential stopping times, computed for the initial dust ring positioned at
(a) $R_{0}=4$ and (b) $R_{0}=10$ with $\alpha_{\rm sur}=0.02$, $\alpha_{\rm
mid}=0.003$, $D_{\rm sur}=0.07~{}H^{2}\Omega$, $D_{\rm
mid}=D_{z}=0.003~{}H^{2}\Omega$. The stopping times are chosen as $T_{\rm
S,mid}=0.01,~{}T_{\rm S,sur}=0.001$. Solid lines represent the surface layer,
and dashed lines are for the middle plane layer. The inset shows the inferred
disk concentration images, with the left half representing the middle plane
layer and the right half representing the surface layer. (c) The steady state
power-law gas density (dotted black line), gas radial velocities of the middle
plane (dashed black line) and the surface layer (solid black line), and the
dust radial drag velocities of the middle plane (blue line) and the surface
layer (red line)..

###  4.1 Coupled dust-gas evolution in a gravitationally unstable disk

As the PSD evolves in time, the disk mass and accretion rate gradually
decline, regions of gravitational instability will recede and MRI turbulence
will dominate in the inner disk (Nakamoto & Nakagawa, 1994; Hueso & Guillot,
2005). We first examine the coupled evolution of dust and gas densities in the
youngest GI-dominated disks (Xu, 2022), by solving Equations 2 and 3 for a
homogeneous layer of dust with $T_{s}=0.001$, using $\alpha=0.1$ and
$D_{r}=2$. Considering the early disk to be very dynamic and GI dominated, we
demonstrate an initially bell-shaped gas density using a Gamma distribution,
as shown in Figure 8(a)

$\displaystyle\Sigma(R)=\Sigma_{0}\frac{6(R-R_{\rm in})^{k-1}e^{-(R-R_{\rm
in})/\theta}}{\Gamma(k)\theta^{k}}$ (11)

where $k=15$ and $\theta=0.6$ are the shape parameters and $\Sigma_{0}$ is the
gas surface density unit. The gas density evolution induces a separatrix
around the peak of the initial distribution $R_{\rm peak}\sim 10$, where gas
drift velocity is inward for $R<R_{\rm peak}$ and outward for $R>R_{\rm
peak}$, as seen in Figure 8(b). The factor $6$ is chosen such that the peak
value of gas density is $\approx\Sigma_{0}$. For the initial dust density, we
choose a narrow ring around an arbitrary radius $R_{0}=8$ analogous to the
delta function adopted by Clarke & Pringle (1988). After 10000 dynamical
timescales at the inner boundary, significant outward/upstream diffusion of
the dust concentration is seen in Figure 8(c). By the end of the simulation
(solid lines), both the gas and dust density profiles have significantly
flattened, while the gas and dust advection velocities have coupled with each
other, conserving the location of the velocity separatrix. Due to large values
of $D_{r}=2.0$, a considerable portion of grains can diffuse out to $R\gg
R_{0}$ although its peak will still drift towards the inner boundary at
velocities close to $v_{g}$, which is negative for $R\lesssim R_{\rm peak}$.

At a later stage, significant GI and MHD turbulence co-exist to influence disk
evolution. Reading off Table 1, we apply viscosity coefficient $\alpha=0.2$
and radial diffusion coefficient $D_{r}=0.5$. We find that the Bell curve
initial condition (Equation 11) follows a similar diffusion process, only
slower in time by a factor of 3-4 due to change in $D_{r}$. To mimic a
relatively more stable disk profile, we then consider another initial gas
density by superimposing a Gaussian bump on a steady-state power-law. The gas
bump is also centered at around $R_{\rm peak}=10$.

$\displaystyle\Sigma=\Sigma_{0}R^{-0.5}+0.2\Sigma_{0}e^{-\frac{(R-R_{\rm
peak})^{2}}{2\sigma_{g}^{2}}}$ (12)

with the Gaussian spread chosen as $\sigma_{g}=2.0$, as shown in Figure. 9(a)
with dotted lines. The gas/dust radial velocities are shown in Figure 9(b),
and the typical dust concentrations at times $t=100,1000,10000$ are shown in
Figure 9(c). Throughout the simulation, the dust and gas velocities are
tightly coupled with each other. At the end of the simulation, the gas density
relaxes towards a power law, and the dust concentration always spread out
beyond the initial ring position $R_{0}$, although its peak still migrates
inward due to a generally negative advection velocity profile. For GI-MHD
turbulence we also tried to start with a smooth steady-state gas profile
without the Gaussian bump (with $v_{d}\approx v_{g}$ being negative throughout
the simulation), and found the evolution of dust concentration is very similar
to Figure 9(c), confirming that outward radial transport is dominated by
diffusion rather than advection as $D_{r}>\alpha$.

###  4.2 Upstream diffusion of dust in a steady-state MRI disk with ideal MHD
turbulence

For quasi-steady MRI turbulence applicable to a later stage of PSD evolution,
in the simulation we measure a radial accretion flow with a velocity
$\sim-10^{-4}v_{K}$. With $H/R\lesssim 0.1$ and $\alpha\sim 2\times 10^{-2}$,
this value is on the order of the radial velocity in the steady solution (e.g.
Clarke & Pringle, 1988)

$v_{g}=-\frac{3\nu}{2R}.$ (13)

Here we examine dust evolution also based on the steady-state model for gas
accretion flow. We demonstrate the possibility of upstream diffusion with the
corresponding steady-state $\Sigma=\Sigma_{0}R^{-0.5}$ profile . Since the gas
evolution equation gives $\partial\Sigma/\partial t=0$, we only solve the dust
evolution equations, but we distinguish two layers of dust with different
radial diffusion coefficients $D_{\rm mid}=0.033$ and $D_{\rm sur}=0.07$ from
Table 1. We further apply $\alpha=0.02$ and make use of $D_{z}=0.013$. The
dust evolution is then solved using Equations 8 & 9, neglecting feedback on
the gas profile, where initially a ring of dust with $T_{S}=0.001$ is located
at $R_{0}$. The $F_{z}$ term represents vertical diffusion between two layers.
The time-independent gas density (black) and the radial advection velocities
(red) are shown in Figure 10(c). To quantify the extent of outward grain
diffusion in our MHD case and better compare two layers, instead of directly
showing the radial evolution of $C$, which is much more quiescent than Figure
8(c) and Figure 9(c), we compute the maximum fraction of dust mass that
reaches the region $R>R_{0}$ beyond the initial radial location of the dust
ring, defined as

$\displaystyle P_{i}(R)=\max_{0<t<T}\frac{\int_{R}^{\infty}2\pi R\Sigma
C_{i}dR}{\int_{R_{\rm in}}^{\infty}2\pi R\Sigma C_{i}dR}$ (14)

where $R_{\rm in}=1.0$ is taken as the position of the inner absorbing
boundary, $T=1\times 10^{4}~{}\Omega^{-1}$ the total simulation time, and
$i={\rm mid,sur}$ for the midplane layer and surface layer respectively.
Figure 10(a)(b) show the maximum fractions of dust mass reaching different
radius $R>R_{0}$, for two typical values of $R_{0}=4$ and $R_{0}=10$. The
insets of Figure 10(a)(b) show axisymmetric 2D distribution of $P(R)$ from the
upstream fraction curves. In both cases, the surface layer display stronger
upstream diffusion compared to the middle plane due to an order-unity larger
$D_{\rm sur}$, although the midplane dust can also diffuse outwards
significantly in the ideal MHD setup.

###  4.3 Dust evolution in a steady-state MRI disk with non-ideal MHD
turbulence

We then consider the effect of non-ideal MHD turbulence, where the midplane of
the disk is much less ionized and less coupled to the MRI turbulence (Gammie,
1996; Kretke & Lin, 2010), with effective $\alpha_{\rm mid}\sim 10^{-3}$.
Hence, the diffusion in the middle plane and the vertical direction may be
constrained by $\alpha_{\rm mid}$ and much less effective. We test this effect
by going beyond simulation measurements and reducing coefficients to $D_{\rm
mid}=D_{z}=0.003$. The midplane layer’s $v_{g}$ and $v_{d}$ are also reduced
according to a smaller $\alpha_{\rm mid}=0.003$. In 3D non-ideal MHD
simulations, gas in certain regions of the upper layer $|z|\gtrsim H$ can have
positive radial velocity, albeit significant outflow ($v_{g}\sim v_{K}$) only
occurs at $|z|\gg H$ where gas density is negligible (Bai, 2017). In our 2D
treatment we still make the simplification to match a characteristic negative
$v_{g}$ of the MRI-active layer with the inward accretion rate to ensure
steady-state (Gammie, 1996), proportional to the characteristic $\alpha$ that
we assume in that layer. Note that if our upper layer has a significant
positive gas velocity, dust outward transport will be aided by a combination
of outward advection and diffusion, which will appear even more significant.
Nevertheless, currently we only focus on the effect of the latter. The results
are presented in Figure 11 similar to Figure 10. The midplane now shows much
smaller upstream fractions, as seen in (a)(b). The gas inflow velocity is much
larger near the surface than the mid-plane, due to $\alpha_{\rm
sur}\gg\alpha_{\rm mid}$, and the dust inward-drift speed near the surface is
also faster than that in the mid-plane (c). Nevertheless, upstream diffusion
is much more pronounced near the surface than the mid-plane, due to $D_{\rm
sur}\gg D_{\rm mid}$.

Additionally, we consider the dust sizes in the midplane being larger (super-
mm) than the surface layer (sub-mm or micron) due to relaxation of the
fragmentation barrier (Equation 10), producing a larger characteristic Stokes
number for the midplane characteristic dust species $T_{\rm s,mid}=0.01$.
Results from this setup is shown in Figure 12. The difference in dust
velocities for the two layers are significantly enhanced, as seen in Figure
12(c), while segregation effect becomes stronger and the surface layer can
diffuse upstream significantly while the midplane is dominated by slow inward
radial drift, as seen in Figure 12 (a) and (b). Changing the stopping time for
mid-plane layer to $T_{\rm s,mid}=0.1$ could enhance this segregation effect
even further.

###  4.4 Dust evolution in an ideal MRI disk with stellar radiation wind

At the final stages of gaseous disk evolution, accretion heating becomes
subdominant and the central star illuminates the outer disk (Nakamoto &
Nakagawa, 1994; Hueso & Guillot, 2005), which could cause significant
evaporation of the gas density (Takeuchi & Lin, 2003). To study the sole
effect of radiation wind, based on the ideal MHD model (§4.2), we modify
Equation 3 to include a sink term

$\displaystyle\dfrac{\partial\Sigma}{\partial t}=$
$\displaystyle\dfrac{3}{R}\dfrac{\partial}{\partial_{R}}\left(R^{1/2}\dfrac{\partial}{\partial_{R}}(\Sigma\nu
R^{1/2})\right)-\dfrac{\Sigma}{\tau_{\rm wind}}f(R)$ (15)

where $f(R)=\left(\tanh[(R-R_{w})/\sigma_{w}]+1\right)/2$ is a smoothing
function.

We solve the modified gas evolution equation together with Equations 8 and 9,
choosing fiducial parameters within a range around $\tau_{\rm wind}\sim 100$,
$R_{w}\sim 15$ and $\sigma_{w}\sim 5$. All other parameters are similar to
that of Figure 10. The reduction of gas density beyond $R_{w}$ induces an
outward gas velocity, on the same order of magnitude as the steady-state
inward gas velocity. The resulting dust evolution is similar to the ideal MHD
case presented in Figure 10, and both the surface layer and middle plane show
considerable upstream diffusion.

##  5 Summary

Based on measurements from numerical simulations, we constructed 1D accretion
disk models of gas and dust to determine outcomes of long-term diffusion
process. In the simulations, we separate the disk into the midplane
$|z|\lesssim H$ section and the surface section $|z|\gtrsim H$. Characteristic
radial diffusion coefficients for each layer $D_{mid}$ and $D_{sur}$ are
measured from Lagrangian particle diffusion, while exchange of particles
between two layers is controlled by an overall vertical diffusion coefficient
$D_{z}$. In the GI and GI-MHD runs, we found $D_{z}\sim\alpha\sim 0.1-0.2$ is
quite efficient and the radial diffusion has little vertical dependence
$D_{mid}\approx D_{\rm sur}\sim D_{\rm r}$, therefore in §4.1 we only consider
one homogeneous layer of $\mu$m size dust with characteristic Stokes number
$T_{s}=0.001$. The early stage of a PSD is GI or GI-MHD dominated and highly
dynamical, therefore we apply impulsive initial conditions with an evolving
gas profile. We found very efficient diffusion of dust with $T_{S}\lesssim
0.001$ due to $D_{\rm r}$ much larger than the gas accretion coefficient
$\alpha$. This implies that not only micron-size crystalline (highly
refractory) grains condensed interior to the main belt of asteroid region can
be transported outwards to Heliocentric distance of the Kuiper Belt region and
be incorporated into comets (such as those retrieved from comet 81P Wild 2 by
the Stardust mission), but up to mm-size CAI condensed in the sub-au inner
regions are also prone to marginal gravitational instability, and can spread
out to the super-au main-belt regions to be assimilated into the planetesimal
parent bodies of chondritic meteorites.

At a later stage, we consider a relatively quiescent MRI dominated disk with a
fixed density profile; We first apply a diffusion model with two layers of
dust of different diffusion coefficient $D_{\rm sur}>D_{\rm mid}$ according to
the MHD simulation. Although diffusion coefficients are much smaller than in
the gravito-turbulence case, we still observe upstream diffusion of dust in
both layers. That is to say, a considerable fraction of dust may reach the
region beyond a its initial radial location through the effect of
concentration diffusion. While the difference between two layers is small for
our ideal-MHD model (§4.2), we also explore the influence of the existence of
dead zone (§4.3), where the turbulence parameter $\alpha$ is significantly
smaller and characteristic dust Stokes number $T_{\rm S,mid}$ can be
significantly larger than the surface layer $T_{\rm S,sur}$. In this case, the
upstream diffusion of mm or cm grains in the dead zone is much less vigorous
than grains in the surface. This effect attributes to the difference in disk
radii inferred from NIR emission by mm-size dust in the midplane, or scattered
light by micron-size grains near the disk surface (Villenave et al., 2020;
Benisty et al., 2022). We lastly consider photo-evaporating winds in the outer
regions of the disk, and found its sole effect can mildly enhance the upstream
diffusion, albeit its efficiency is limited.

We thank William Bethune, Hans Baehr, Zhaohuan Zhu and Wenrui Xu for useful
conversations. We thank the anonymous referee for helpful suggestions. We have
benefited from many past discussions with Willy Kley and would like to
dedicate this paper to his memory. Y. X. C. thanks Man Hoi Lee and Hong Kong
University for their hospitality during the completion of this manuscript.
T.Z. is supported by the Cecil and Sally Drinkward postdoc fellowship during
this work.

## References

  * Adams et al. (1987) Adams, F. C., Lada, C. J., & Shu, F. H. 1987, ApJ, 312, 788, doi: 10.1086/164924
  * Armitage (2019) Armitage, P. J. 2019, in From Protoplanetary Disks to Planet Formation (Springer), 1–150
  * Baehr & Zhu (2021) Baehr, H., & Zhu, Z. 2021, ApJ, 909, 136, doi: 10.3847/1538-4357/abddb4
  * Bai (2017) Bai, X.-N. 2017, ApJ, 845, 75, doi: 10.3847/1538-4357/aa7dda
  * Bai & Stone (2013) Bai, X.-N., & Stone, J. M. 2013, The Astrophysical Journal, 769, 76
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Balbus & Hawley (1998) —. 1998, Reviews of Modern Physics, 70, 1, doi: 10.1103/RevModPhys.70.1
  * Beckwith et al. (2011) Beckwith, K., Armitage, P. J., & Simon, J. B. 2011, Monthly Notices of the Royal Astronomical Society, 416, 361
  * Benisty et al. (2022) Benisty, M., Dominik, C., Follette, K., et al. 2022, arXiv preprint arXiv:2203.09991
  * Béthune et al. (2021) Béthune, W., Latter, H., & Kley, W. 2021, Astronomy & Astrophysics, 650, A49
  * Birnstiel et al. (2010) Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79, doi: 10.1051/0004-6361/200913731
  * Bockelée-Morvan et al. (2002) Bockelée-Morvan, D., Gautier, D., Hersant, F., Huré, J. M., & Robert, F. 2002, A&A, 384, 1107, doi: 10.1051/0004-6361:20020086
  * Boley & Durisen (2006) Boley, A., & Durisen, R. 2006, The Astrophysical Journal, 641, 534
  * Boley & Durisen (2010) Boley, A. C., & Durisen, R. H. 2010, The Astrophysical Journal, 724, 618
  * Bouwman et al. (2003) Bouwman, J., de Koter, A., Dominik, C., & Waters, L. B. F. M. 2003, A&A, 401, 577, doi: 10.1051/0004-6361:20030043
  * Bouwman et al. (2008) Bouwman, J., Henning, T., Hillenbrand, L. A., et al. 2008, ApJ, 683, 479, doi: 10.1086/587793
  * Brownlee et al. (2006) Brownlee, D., Tsou, P., Aléon, J., et al. 2006, Science, 314, 1711, doi: 10.1126/science.1135840
  * Cameron (1978) Cameron, A. G. W. 1978, Moon and Planets, 18, 5, doi: 10.1007/BF00896696
  * Chen et al. (2020) Chen, Y.-X., Li, Y.-P., Li, H., & Lin, D. N. C. 2020, ApJ, 896, 135, doi: 10.3847/1538-4357/ab9604
  * Chen et al. (2021) Chen, Y.-X., Wang, Z., Li, Y.-P., Baruteau, C., & Lin, D. N. C. 2021, ApJ, 922, 184, doi: 10.3847/1538-4357/ac23d7
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368, doi: 10.1086/304869
  * Ciesla (2007) Ciesla, F. J. 2007, Science, 318, 613, doi: 10.1126/science.1147273
  * Clarke (1989) Clarke, C. J. 1989, MNRAS, 238, 283, doi: 10.1093/mnras/238.2.283
  * Clarke & Pringle (1988) Clarke, C. J., & Pringle, J. E. 1988, MNRAS, 235, 365, doi: 10.1093/mnras/235.2.365
  * Colbrook et al. (2017) Colbrook, M. J., Ma, X., Hopkins, P. F., & Squire, J. 2017, MNRAS, 467, 2421, doi: 10.1093/mnras/stx261
  * Deng et al. (2021) Deng, H., Mayer, L., & Helled, R. 2021, Nature Astronomy, 5, 440
  * Deng et al. (2020) Deng, H., Mayer, L., & Latter, H. 2020, The Astrophysical Journal, 891, 154
  * Deng et al. (2019) Deng, H., Mayer, L., Latter, H., Hopkins, P. F., & Bai, X.-N. 2019, The Astrophysical Journal Supplement Series, 241, 26
  * Deng et al. (2017) Deng, H., Mayer, L., & Meru, F. 2017, The Astrophysical Journal, 847, 43
  * Desch et al. (2012) Desch, S. J., Morris, M. A., Connolly, H. C., & Boss, A. P. 2012, , 47, 1139, doi: 10.1111/j.1945-5100.2012.01357.x
  * Durisen et al. (2007) Durisen, R. H., Boss, A. P., Mayer, L., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 607.  https://arxiv.org/abs/astro-ph/0603179
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355, doi: 10.1086/176735
  * Garaud et al. (2004) Garaud, P., Barrière-Fouchet, L., & Lin, D. N. C. 2004, ApJ, 603, 292, doi: 10.1086/381385
  * Garaud & Lin (2007) Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606, doi: 10.1086/509041
  * Gong et al. (2019) Gong, M., Zheng, X., Lin, D. N. C., et al. 2019, ApJ, 883, 164, doi: 10.3847/1538-4357/ab3e70
  * Haghighipour & Boss (2003) Haghighipour, N., & Boss, A. P. 2003, The Astrophysical Journal, 583, 996
  * Hartmann (1998) Hartmann, L. 1998, Accretion Processes in Star Formation
  * Hayashi et al. (1985) Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews, 1100–1153
  * Hopkins (2015) Hopkins, P. F. 2015, Monthly Notices of the Royal Astronomical Society, 450, 53
  * Hu & Bai (2021) Hu, Z., & Bai, X.-N. 2021, MNRAS, 503, 162, doi: 10.1093/mnras/stab542
  * Hueso & Guillot (2005) Hueso, R., & Guillot, T. 2005, A&A, 442, 703, doi: 10.1051/0004-6361:20041905
  * Ida & Lin (2004) Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388, doi: 10.1086/381724
  * Kley & Lin (1992) Kley, W., & Lin, D. N. C. 1992, ApJ, 397, 600, doi: 10.1086/171818
  * Kratter & Lodato (2016) Kratter, K., & Lodato, G. 2016, Annual Review of Astronomy and Astrophysics, 54, 271
  * Kretke & Lin (2010) Kretke, K. A., & Lin, D. 2010, The Astrophysical Journal, 721, 1585
  * Krumholz et al. (2018) Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716, doi: 10.1093/mnras/sty852
  * Li et al. (2019) Li, Y.-P., Li, H., Ricci, L., et al. 2019, ApJ, 878, 39, doi: 10.3847/1538-4357/ab1f64
  * Lin & Papaloizou (1985) Lin, D. N. C., & Papaloizou, J. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews, 981–1072
  * Liu et al. (2022) Liu, B., Johansen, A., Lambrechts, M., Bizzarro, M., & Haugbølle, T. 2022, Science Advances, 8, eabm3045, doi: 10.1126/sciadv.abm3045
  * Löhnert & Peeters (2022) Löhnert, L., & Peeters, A. 2022
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603
  * Morfill & Voelk (1984) Morfill, G. E., & Voelk, H. J. 1984, ApJ, 287, 371, doi: 10.1086/162697
  * Nakamoto & Nakagawa (1994) Nakamoto, T., & Nakagawa, Y. 1994, ApJ, 421, 640, doi: 10.1086/173678
  * Ogilvie (2018) Ogilvie, G. I. 2018, Monthly Notices of the Royal Astronomical Society, 477, 1744
  * Okuzumi & Hirose (2011) Okuzumi, S., & Hirose, S. 2011, The Astrophysical Journal, 742, 65
  * Ormel & Cuzzi (2007) Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413, doi: 10.1051/0004-6361:20066899
  * Pringle (1981) Pringle, J. E. 1981, ARA&A, 19, 137, doi: 10.1146/annurev.aa.19.090181.001033
  * Rice et al. (2004) Rice, W., Lodato, G., Pringle, J., Armitage, P., & Bonnell, I. A. 2004, Monthly Notices of the Royal Astronomical Society, 355, 543
  * Riols & Latter (2019) Riols, A., & Latter, H. 2019, Monthly Notices of the Royal Astronomical Society, 482, 3989
  * Riols et al. (2020a) Riols, A., Lesur, G., & Menard, F. 2020a, Astronomy & Astrophysics, 639, A95
  * Riols et al. (2020b) Riols, A., Roux, B., Latter, H., & Lesur, G. 2020b, Monthly Notices of the Royal Astronomical Society, 493, 4631
  * Ruden & Lin (1986) Ruden, S. P., & Lin, D. N. C. 1986, ApJ, 308, 883, doi: 10.1086/164559
  * Safronov (1960) Safronov, V. S. 1960, Annales d’Astrophysique, 23, 979
  * Safronov & Ruzmaikina (1985) Safronov, V. S., & Ruzmaikina, T. V. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews, 959–980
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shang et al. (2000) Shang, H., Shu, F. H., Lee, T., & Glassgold, A. E. 2000, Space Sci. Rev., 92, 153, doi: 10.1023/A:1005234909848
  * Shi et al. (2016) Shi, J.-M., Zhu, Z., Stone, J. M., & Chiang, E. 2016, MNRAS, 459, 982, doi: 10.1093/mnras/stw692
  * Shu et al. (2000) Shu, F. H., Najita, J. R., Shang, H., & Li, Z. Y. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell, 789–814
  * Shu et al. (2001) Shu, F. H., Shang, H., Gounelle, M., Glassgold, A. E., & Lee, T. 2001, ApJ, 548, 1029, doi: 10.1086/319018
  * Simon et al. (2012) Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, Monthly Notices of the Royal Astronomical Society, 422, 2685
  * Takeuchi & Lin (2005) Takeuchi, T., & Lin, D. 2005, The Astrophysical Journal, 623, 482
  * Takeuchi & Lin (2002) Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344, doi: 10.1086/344437
  * Takeuchi & Lin (2003) —. 2003, ApJ, 593, 524, doi: 10.1086/376496
  * Toomre (1964) Toomre, A. 1964, The Astrophysical Journal, 139, 1217
  * Urpin (1984) Urpin, V. A. 1984, Soviet Ast., 28, 50
  * van Boekel et al. (2005) van Boekel, R., Min, M., Waters, L. B. F. M., et al. 2005, A&A, 437, 189, doi: 10.1051/0004-6361:20042339
  * Villenave et al. (2020) Villenave, M., Ménard, F., Dent, W., et al. 2020, Astronomy & Astrophysics, 642, A164
  * Watson et al. (2009) Watson, D. M., Leisenring, J. M., Furlan, E., et al. 2009, ApJS, 180, 84, doi: 10.1088/0067-0049/180/1/84
  * Weidenschilling et al. (1998) Weidenschilling, S. J., Marzari, F., & Hood, L. L. 1998, Science, 279, 681, doi: 10.1126/science.279.5351.681
  * Wissing et al. (2022) Wissing, R., Shen, S., Wadsley, J., & Quinn, T. 2022, Astronomy & Astrophysics, 659, A91
  * Wooden et al. (2007) Wooden, D., Desch, S., Harker, D., Gail, H. P., & Keller, L. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 815
  * Xu (2022) Xu, W. 2022, ApJ, 934, 156, doi: 10.3847/1538-4357/ac7b94
  * Yang & Krumholz (2012) Yang, C.-C., & Krumholz, M. 2012, ApJ, 758, 48, doi: 10.1088/0004-637X/758/1/48
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588, doi: 10.1016/j.icarus.2007.07.012
  * Zhu et al. (2012) Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, The Astrophysical Journal, 755, 6
  * Zhu et al. (2011) Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47, doi: 10.1088/0004-637X/729/1/47
  * Zhu et al. (2015) Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, The Astrophysical Journal, 801, 81


# Generalized Covariant Entropy Bound in Lanczos-Lovelock Gravity

Ming Zhang1<EMAIL_ADDRESS>Jie Jiang2 Corresponding author.
<EMAIL_ADDRESS>1Department of Physics, Jiangxi Normal University,
Nanchang 330022, China  2College of Education for the Future, Beijing Normal
University, Zhuhai 519087, China

###### Abstract

In this paper, we investigate the generalized covariant entropy bound in the
theory where the Einstein gravity is perturbed by the higher-order Lovelock
terms. After replacing the Bekenstein-Hawking entropy with the Jacobson-Myers
entropy and introducing two reasonable physical assumptions, we showed that
the corresponding generalized covariant entropy bound is satisfied under a
higher-order approximation of the perturbation from the higher-order Lovelock
terms. Our result implies that the Jacobson-Myers entropy strictly obeys the
entropy bound under the perturbation level, and the generalized second law of
Lanczos-Lovelock gravity is also satisfied when the Einstein gravity is
perturbed by the higher-order Lovelock terms.

##  I Introduction

The investigation of black hole thermodynamics has led to some interesting
entropy bounds that should be observed to guarantee theoretical consistency.
Bekenstein Bekenstein:1980jp  has conjectured that the entropy $S$ and energy
$E$ of any stable gravitational thermodynamic system satisfies a universal
bound

$\displaystyle\begin{aligned} S\leq 2\pi ER\,,\end{aligned}$ (1)

in which $R$ is defined as the circumferential radius of the sphere
surrounding the thermodynamical system. This bound is called the Bekenstein
bound and it can be indicated by the generalized second law (GSL) of black
holes. The Bekenstein bound has been confirmed in many weakly gravitational
systems with finite size. In a strongly gravitational system, it is hard to
define the energy $E$ and radius $R$ locally. Counterexamples can be found in
the process of gravitational collapse Bousso:1999xy . For a spherical system
in Einstein gravity, the Bekenstein bound can be simplified as

$\displaystyle\begin{aligned} S\leq\frac{A}{4}\,,\end{aligned}$ (2)

in which $S$ and $A$ are the entropy and area of the system. It is worth
noting that this bound is not well defined in a strongly gravitational system
since the area $A$ is dependent on the choice of the spacelike region in the
system and it can always be selected to be arbitrarily small by an almost null
hypersurface. It was shown that this bound can be violated in the system for
large volume Fischler:1998st .

To find a covariant version of the entropy bound, Bousso considered the
entropy cross a light sheet and proposed a covariant entropy bound, called the
Bousso bound Bousso:1999xy , which can be well formulated in arbitrarily
curved spacetime. Consider a $(D-2)$-dimensional compact spacelike surface $B$
with area $A(B)$. Let $L$ be a null hypersurface generated by the null
geodesics which starts at $B$ and is orthogonal to $B$. Assume that the
expansion of the null congruence is nonpositive (i.e., $L$ is a light sheet)
and $L$ is not terminated until a caustic point is reached. Then, the entropy
$S_{L}$ passing through the light sheet $L$ is bounded by the quarter of
$A(B)$, i.e.,

$\displaystyle\begin{aligned} S_{L}\leq\frac{A(B)}{4}\,.\end{aligned}$ (3)

This is the covariant bound proposed by Bousso and it is conjectured to be
valid in any strongly gravitational system with arbitrary large regions. This
bound is shown to be hold in various cases Flanagan:1999jp ; Gao:2004mc ;
Gao:2005bk ; Gao:2008zy ; Bousso:2003kb ; Pesci:2007rp ; Pesci:2008yy  and it
can be regarded as a formulation of holographic principle in spacetime.

Note that the above conjecture requires that the light sheet $L$ ends at a
caustic point. Flanagan $et.\,al.$ Flanagan:1999jp  extended this bound in
which the light sheet can be terminated at another $(D-2)$-dimensional spatial
surface $B^{\prime}$ before reaching a caustic. Then the entropy bound is
modified as

$\displaystyle\begin{aligned}
S_{L}\leq\frac{1}{4}\left|A(B^{\prime})-A(B)\right|\,,\end{aligned}$ (4)

in which $A(B^{\prime})$ is the area of the spatial surface $B^{\prime}$. This
is called the generalized covariant entropy bound or generalized Bousso bound
and it has been proved in Einstein gravity under some physical assumptions
Flanagan:1999jp ; Bousso:2003kb ; Strominger:2003br .

General relativity (GR) is not a complete theory of gravity due to the lack of
a definitive quantum gravity theory and it can only be regarded as an
effective theory in a certain region of scale. After considering the quantum
effect or string modification, the higher-curvature terms are often added to
Einstein-Hilbert action to modify the effective action of the gravitational
theory A13 ; A14 ; A15 ; A16 . In these cases, the Einstein gravity is
perturbed by the higher-curvature terms. In this paper, we focus on the
Lanczos-Lovelock gravity, which is the only natural generalization of Einstein
gravity to higher-dimensional spacetime if we demand that the equations of
motion are second-order differential equations of metric Lanczos ; Lovelock .
Moreover, unlike most higher-curvature gravitational theories, Lanczos-
Lovelock gravity is a ghost-free theory and admits consistent initial value
formulation Lovelock2 ; Kovacs:2020ywu . As mentioned above, the generalized
covariant entropy bound (4) is only valid in Einstein gravity. It is natural
to ask whether the higher-curvature corrections can affect the entropy bound
of the gravitational theory. Recently, Matsuda $et.\,al.$ Matsuda:2020yvl
extended the generalized covariant entropy bound into the modified
gravitational theory by replacing the quarter of area $A(B)/4$ with some
appropriate black hole entropy $S_{\text{bh}}(B)$, such as the Wald entropy
A17 ; A18  or Jacobson-Myers (JM) entropy Jacobson:1993xs . Under two
reasonable assumptions, they proved the entropy bound for Wald entropy in
$f(R)$ gravity and canonical scalar-tensor theory. Moreover, they also showed
that the bound using JM entropy holds for the GR branch of spherically
symmetric configurations in Einstein-Gauss-Bonnet gravity. In the following,
we would like to extend their discussion into the case where the Einstein
gravity is perturbed by the higher-order Lovelock terms and show that the JM
entropy can give a reasonable entropy bound in this theory.

The outline of this paper is as follows. In Sec. II, we briefly review the
Lanczos-Lovelock gravity and discuss the features of Wald entropy and JM
entropy. In Sec. III, we introduce the generalized entropy bound in Lanczos-
Lovelock gravity and show the physical assumptions as well as the key point
for proving this bound. In Sec. IV, we prove the generalized entropy bound in
the theory where the Einstein gravity is perturbed by the higher-order
Lovelock terms and show that the JM entropy strictly obeys the entropy bound
under the perturbation level. Finally, the conclusion and discussion are
presented in Sec. V.

##  II Lanczos-Lovelock gravity

In this paper, we consider the Lanczos-Lovelock gravitational theory with some
minimally coupled matter fields. The action of this theory in $D$-dimensional
spacetime is given by

$\displaystyle\begin{aligned} I=\frac{1}{16\pi}\int
d^{D}x\sqrt{g}\left(\sum_{k=0}^{k_{\text{max}}}\frac{a_{k}}{2^{k}}\mathcal{L}^{(k)}+\mathcal{L}_{\text{mat}}\right),\end{aligned}$
(5)

in which $\mathcal{L}_{\text{mat}}$ is the Lagrangian density of the matter
fields, $g_{ab}$ is the Minkowski metric of the spacetime, and

$\displaystyle\begin{aligned} \mathcal{L}^{(k)}=\delta_{c_{1}d_{1}\cdots
c_{k}d_{k}}^{a_{1}b_{1}\cdots a_{k}b_{k}}R^{c_{1}d_{1}}_{a_{1}b_{1}}\cdots
R^{c_{k}d_{k}}_{a_{k}b_{k}}\end{aligned}$ (6)

is the $k$-order Lovelock term. Here
$k_{\text{max}}=\left[(D-1)/2\right]$111The Square brackets [x] denotes the
integer part of x. and

$\displaystyle\begin{aligned} \delta_{c_{1}d_{1}\cdots
c_{k}d_{k}}^{a_{1}b_{1}\cdots
a_{k}b_{k}}=(2k)!\delta^{[a_{1}}_{c_{1}}\delta^{b_{1}}_{d_{1}}\cdots\delta^{a_{k}}_{c_{k}}\delta^{b_{k}]}_{d_{k}}\end{aligned}$
(7)

is the generalized Kronecker tensor. The equation of motion is given by

$\displaystyle\begin{aligned} &E_{ab}=8\pi T_{ab}\,,\end{aligned}$ (8)

in which $T_{ab}$ is the stress-energy tensor of the matter fields, and

$\displaystyle\begin{aligned}
E_{ab}=-\sum^{k_{\text{max}}}_{k=0}\frac{a_{k}}{2^{k+1}}\delta_{ac_{1}d_{1}\cdots
c_{k}d_{k}}^{ba_{1}b_{1}\cdots a_{k}b_{k}}R^{c_{1}d_{1}}_{a_{1}b_{1}}\cdots
R^{c_{k}d_{k}}_{a_{k}b_{k}}\end{aligned}$ (9)

is the generalized Einstein tensor of Lanczos-Lovelock gravity. Employing the
Noether charge method of Iyer and Wald A17 ; A18 , the Wald entropy of
Lanczos-Lovelock gravity can be obtained and it is given by

$\displaystyle\begin{aligned}
S_{W}=-2\pi\int_{s}d^{D-2}x\sqrt{\gamma}P^{abcd}\hat{\bm{\epsilon}}_{ab}\hat{\bm{\epsilon}}_{cd}\,,\end{aligned}$
(10)

where we have denoted

$\displaystyle\begin{aligned}
P_{ab}^{cd}=\frac{1}{16\pi}\sum^{k_{\text{max}}}_{k=0}\frac{ka_{k}}{2^{k}}\delta_{cdc_{2}d_{2}\cdots
c_{k}d_{k}}^{aba_{2}b_{2}\cdots a_{k}b_{k}}R^{c_{2}d_{2}}_{a_{2}b_{2}}\cdots
R^{c_{k}d_{k}}_{a_{k}b_{k}}\,.\end{aligned}$ (11)

Here $s$ is a cross-section of event horizon, $\gamma_{ab}$ is the induced
metric on $s$, and $\hat{\bm{\epsilon}}_{ab}$ is the binormal to $s$. The Wald
entropy gives the correct first law in the stationary black holes. However, as
discussed in Refs. A25 ; A26 ; A27 ; A28 , the Wald entropy of the Lanczos-
Lovelock gravity does not obey the linearized second law and we need to focus
on the JM entropy, i.e.,

$\displaystyle
S_{\text{JM}}=\frac{1}{4}\int_{s}d^{D-2}x\sqrt{\gamma}\rho_{\text{JM}}$ (12)

with

$\displaystyle\begin{aligned}
\rho_{\text{JM}}=\sum^{k_{\text{max}}}_{k=1}\frac{ka_{k}}{2^{k-1}}\delta_{c_{2}d_{2}\cdots
c_{k}d_{k}}^{a_{2}b_{2}\cdots
a_{k}b_{k}}\hat{R}^{c_{2}d_{2}}_{a_{2}b_{2}}\cdots\hat{R}^{c_{k}d_{k}}_{a_{k}b_{k}}\,.\end{aligned}$
(13)

in which $\hat{R}^{cd}_{ab}$ is the Riemann tensor of the induced metric
$\gamma_{ab}$ on the cross-section $s$. In the stationary black hole, JM
entropy and Wald entropy give the same result, and therefore the JM entropy
also obeys the first law. Considering the relationship between the generalized
covariant entropy bound and the generalized second law of the black holes, it
is natural to apply the JM entropy to discuss the entropy bound in the
Lanczos-Lovelock gravity.

##  III Generalized covariant entropy bound

In this section, we first introduce the basic setups of the generalized
covariant entropy bound in Lanczos-Lovelock gravity. Let $L$ be a null
hypersurface generalized by null geodesics, which starts at a compact
$(D-2)$-dimensional spatial surface $B_{0}$ and ends at another compact
$(D-2)$-dimensional spatial surface $B_{1}$. Let $k^{a}=(\partial/\partial
u)^{a}$ be the tangent vector field of the null geodesics, in which $u$ is an
affine parameter of the null geodesics such that the spatial surfaces $B_{0}$
and $B_{1}$ are given by $u=0$ and $u=1$, separately. Any spatial surface $B$
determined by the same $u$ is called the cross-section of the null
hypersurface. Then, we can choose $(u,x)$ to a coordinate system on the null
hypersurface $L$, in which $x=\\{x^{1},\cdots x^{D-2}\\}$ denotes the
coordinate of the cross-section and every geodesic is determined by a constant
$x$. Then, the covariant entropy bound in Lanczos-Lovelock gravity demands
that the entropy $S_{L}$ passing through the null hypersurface $L$ should
satisfy

$\displaystyle\begin{aligned}
S_{L}\leq|S_{\text{JM}}(B_{0})-S_{\text{JM}}(B_{1})|\,,\end{aligned}$ (14)

in which $S_{\text{JM}}(B)$ is evaluated by the JM entropy formula (12) on the
cross-section $B$.

To prove the entropy bound, we first define the generalized expansion $\Theta$
of the JM entropy as the change of entropy per unit area, i.e.,

$\displaystyle\begin{aligned}
\frac{dS_{\text{JM}}}{du}=\frac{1}{4}\int_{B}d^{D-2}x\sqrt{\gamma}\Theta\,.\end{aligned}$
(15)

Noting that the JM entropy is a purely spatial quantity in the
$(D-2)$-dimensional slice $B$, i.e., it is determined by the induced metric
$\gamma_{ab}$, we can regard $S_{\text{JM}}$ as an action on the cross-section
$B$. Then, the Lie-derivative $\mathcal{L}_{k}=\partial_{u}$ can be seen as a
variation on $S_{\text{JM}}$. After assuming that $B$ is compact and dropping
the surface terms in $\partial_{u}S_{\text{JM}}$, we can get

$\displaystyle\begin{aligned}
\partial_{u}S_{\text{JM}}&=-\frac{1}{4}\sum^{k_{\text{max}}}_{k=1}ka_{k}\int_{B}d^{D-2}x\sqrt{\gamma}[\hat{E}^{(k-1)}]^{ab}\partial_{u}\gamma_{ab}\\\
&=-\frac{1}{2}\sum^{k_{\text{max}}}_{k=1}ka_{k}\int_{B}d^{D-2}x\sqrt{\gamma}[\hat{E}^{(k-1)}]^{ab}K_{ab}\end{aligned}$
(16)

in which

$\displaystyle\begin{aligned}
{}^{b}_{a}=-\frac{1}{2^{k+1}}\delta_{ac_{1}d_{1}\cdots
c_{k}d_{k}}^{ba_{1}b_{1}\cdots
a_{k}b_{k}}\hat{R}^{c_{1}d_{1}}_{a_{1}b_{1}}\cdots\hat{R}^{c_{k}d_{k}}_{a_{k}b_{k}}\,,\end{aligned}$
(17)

and

$\displaystyle K_{ab}=\frac{1}{2}\partial_{u}\gamma_{ab}$ (18)

is the extrinsic curvature associated with $k^{a}$. These results imply that

$\displaystyle\begin{aligned}
\Theta=K_{b}^{a}\sum^{k_{\text{max}}}_{k=1}\frac{ka_{k}}{2^{k-1}}\delta_{ac_{2}d_{2}\cdots
c_{k}d_{k}}^{ba_{2}b_{2}\cdots
a_{k}b_{k}}\hat{R}^{c_{2}d_{2}}_{a_{2}b_{2}}\cdots\hat{R}^{c_{k}d_{k}}_{a_{k}b_{k}}\end{aligned}$
(19)

after neglecting the total-derivative terms.

Choose $u$ to be an affine parameter of the null geodesics. Using the equation
of motion, we can write the change of $\Theta$ as

$\displaystyle\begin{aligned}
\frac{d\Theta}{du}=-8\pi\mathcal{T}+\mathcal{F}\,,\end{aligned}$ (20)

in which

$\displaystyle\begin{aligned} \mathcal{T}&=T_{ab}k^{a}k^{b}\,,\\\
\mathcal{F}&=H_{ab}k^{a}k^{b}+k^{a}\nabla_{a}\Theta\,.\end{aligned}$ (21)

This can be regarded as the Raychaudhuri equation in Lanczos-Lovelock gravity.

In the thermodynamic limit, there exists an entropy flux vector field $s^{a}$
such that the entropy passing through the null hypersurface $L$ can be written
as

$\displaystyle\begin{aligned}
S_{L}=\int_{L}d^{D-2}xdu\sqrt{\gamma}s\end{aligned}$ (22)

with the entropy density

$\displaystyle s=-k_{a}s^{a}\,.$ (23)

Analogies to the assumptions in Einstein gravity Bousso:2003kb ;
Strominger:2003br , Ref. Matsuda:2020yvl  made two following assumptions in
the modified gravitational theories,

$\displaystyle\begin{aligned} &\text{(i)}\quad\partial_{u}s(x,u)\leq
2\pi\mathcal{T}(x,u)\,,\\\ &\text{(ii)}\quad
s(x,0)\leq-\frac{1}{4}\Theta(x,0)\,\end{aligned}$ (24)

on the null hypersurface $L$. The first assumption is from the requirement
that the change rate of the entropy flux is not large than the energy flux and
it can also be regarded as the consequence of the version of Bekenstein bound
Matsuda:2020yvl . The second assumption is just an initial choice of the
hypersurface such that the entropy bound is valid at the beginning of $L$.

With the above setups and assumptions, it is not hard to get

$\displaystyle\begin{aligned}
s(x,u)&=s(x,0)+\int_{0}^{u}du\partial_{u}s(x,u)\\\ &\leq
s(x,0)+2\pi\int_{0}^{u}du\mathcal{T}(x,u)\,,\end{aligned}$ (25)

in which we have used the first assumption at the last step. Then, using Eq.
(20) and together with the second assumption in Eq. (LABEL:assumption), we
have

$\displaystyle\begin{aligned}
s(x,u)\leq&s(x,0)-\frac{1}{4}\Theta(x,\lambda)+\frac{1}{4}\Theta(x,0)+\frac{1}{4}\int_{0}^{u}d\tilde{u}\mathcal{F}(x,\tilde{u})\\\
\leq&-\frac{1}{4}\Theta(x,\lambda)+\frac{1}{4}\int_{0}^{u}d\tilde{u}\mathcal{F}(x,\tilde{u})\,.\\\
\end{aligned}$ (26)

Finally, after integrating the above identity over $L$, we have

$\displaystyle\begin{aligned} S_{L}\leq
S_{\text{JM}}(B_{0})-S_{\text{JM}}(B_{1})+\frac{1}{4}\int_{0}^{1}du\int
d^{D-2}x\sqrt{\gamma}\Delta(x,u)\,,\quad\end{aligned}$ (27)

in which we have denoted

$\displaystyle\begin{aligned}
\Delta(x,u)=\int_{0}^{u}d\tilde{u}\mathcal{F}(x,\tilde{u})\,.\end{aligned}$
(28)

From the above results, we can see that the key point to examining the
generalized covariant entropy bound is to judge the sign of
$\mathcal{F}(x,u)$. If we have $\mathcal{F}(x,u)\leq 0$ on $L$, the inequality
(27) reduces to

$\displaystyle\begin{aligned} S_{L}\leq
S_{\text{JM}}(B_{0})-S_{\text{JM}}(B_{1})\,,\end{aligned}$ (29)

which is the entropy bound given by Eq. (14). For the Einstein gravity, Eq.
(20) is just the Raychaudhuri equation and we have $\mathcal{F}\leq 0$, which
gives the proof of the generalized covariant entropy bound in Einstein
gravity. In the following, we would like to judge the sign of $\mathcal{F}$ in
the Lanczos-Lovelock gravity.

##  IV Proof of the entropy bound with higher-curvature corrections

From the perspective of quantum corrections and string theory, it is natural
to consider the models of gravity where the Einstein gravity is perturbed by
higher curvature terms. Therefore, in the following, we consider the Lanczos-
Lovelock gravity where the higher-order Lovelock terms are treated as small
corrections to the Einstein gravity, i.e., we consider the Lovelock theory
with $a_{0}=-2\Lambda$, $a_{1}=1$ and $a_{k}=\lambda\alpha_{k}$ for $k\geq 2$,
in which $\lambda$ is a small quantity which describes the perturbation from
the higher-curvature terms. Then, we have

$\displaystyle\begin{aligned}
H_{a}^{b}=G_{a}^{b}-\Lambda\delta_{a}^{b}-\lambda\sum^{k_{\text{max}}}_{k=2}\frac{\alpha_{k}}{2^{k+1}}\delta_{ac_{1}d_{1}\cdots
c_{k}d_{k}}^{ba_{1}b_{1}\cdots a_{k}b_{k}}R^{c_{1}d_{1}}_{a_{1}b_{1}}\cdots
R^{c_{k}d_{k}}_{a_{k}b_{k}}\,,\quad\end{aligned}$ (30)

and

$\displaystyle\begin{aligned} \rho_{\text{JM}}=1+\lambda\rho\end{aligned}$
(31)

with

$\displaystyle\rho=\sum^{k_{\text{max}}}_{k=2}\frac{k\alpha_{k}}{2^{k-1}}\delta_{c_{2}d_{2}\cdots
c_{k}d_{k}}^{a_{2}b_{2}\cdots
a_{k}b_{k}}\hat{R}^{c_{2}d_{2}}_{a_{2}b_{2}}\cdots\hat{R}^{c_{k}d_{k}}_{a_{k}b_{k}}.$
(32)

After considering the higher-curvature corrections, the solution in the theory
will depend on the small parameter $\lambda$, i.e., $g_{ab}(\lambda)$, in
which $\lambda=0$ describes the solution of Einstein gravity.

To evaluate $\mathcal{F}(x,u,\lambda)$ on the null hypersurface $L$, we
introduce the Gaussian null coordinate system $\\{z,u,x\\}$, in which the line
element can be expressed as

$\displaystyle\begin{aligned} ds^{2}(\lambda)=2(dz+z^{2}\alpha
du+z\beta_{i}dx^{i})du+\gamma_{ij}dx^{i}dx^{j}\,,\quad\quad\end{aligned}$ (33)

in which the null hypersurface $L$ is given by $z=0$, and $\alpha$,
$\beta_{i}$ and $\gamma_{ij}$ are the function of $u,z,x,\lambda$. Here the
index $i,j,k,l$ denotes the coordinate of the cross-section $B$. The null
generator of $L$ is given by $k^{a}=(\partial/\partial u)^{a}$. Using this
line element, the nonvanishing component of the Christoffel symbol on $L$ can
be further obtained,

$\displaystyle\begin{aligned}
&\Gamma^{k}{}_{ij}=\hat{\Gamma}^{k}{}_{ij}\,,\quad\Gamma^{j}{}_{ui}=K^{j}_{i}\,,\quad\Gamma^{j}{}_{zi}=\bar{K}^{j}_{i}\,,\quad\Gamma^{1}{}_{uz}=\frac{1}{2}\beta^{i}\,,\\\
&\Gamma^{u}{}_{ij}=-\bar{K}_{ij}\,,\,\,\Gamma^{u}{}_{ui}=-\frac{1}{2}\beta_{i}\,,\,\,\Gamma^{z}{}_{ij}=-\bar{K}_{ij}\,,\,\,\Gamma^{z}{}_{zi}=\frac{1}{2}\beta_{i}\,,\end{aligned}$
(34)

in which $\hat{\Gamma}^{k}{}_{ij}$ is the Christoffel symbol of the induced
metric $\gamma_{ij}$, and

$\displaystyle\begin{aligned}
K_{ij}=\frac{1}{2}\partial_{u}\gamma_{ij}\,,\quad\quad\bar{K}_{ij}=\frac{1}{2}\partial_{z}\gamma_{ij}\end{aligned}$
(35)

are the extrinsic curvature associated with the null vectors
$(\partial/\partial u)^{a}$ and $(\partial/\partial z)^{a}$ separately.
Further calculation gives

$\displaystyle\begin{aligned}
&R_{ij}^{kl}=\hat{R}_{ij}^{kl}-4K_{[i}^{[k}\bar{K}_{j]}^{l]}\,,\quad
R^{zj}_{ui}=-\partial_{u}K_{i}^{j}-K_{i}^{k}K^{j}_{k}\,,\\\
&R_{ui}^{jk}=-2D^{[j}K_{i}^{k]}+K_{i}^{[j}\beta^{k]}\,,\quad
R^{zi}_{jk}=-2D_{[j}K^{i}_{k]}+K^{i}_{[j}\beta_{k]}\end{aligned}$ (36)

on the hypersurface $L$. Using the above results and considering the symmetry
of the generalized Kronecker tensor, it is not difficult to get

$\displaystyle\begin{aligned}
E_{ab}k^{a}k^{b}=R_{u}^{z}+\lambda\sum^{k_{\text{max}}}_{k=2}\alpha_{k}[E^{(k)}]^{z}_{u}\end{aligned}$
(37)

on the null hypersurface $L$. Considering the antisymmetry of the generalized
Kronecker tensor and using Eq. (LABEL:expR), it is not hard to get

$\displaystyle\begin{aligned}
{}^{z}_{u}=&\frac{k}{2^{k-1}}R^{zj}_{ui}\delta^{ji_{2}j_{2}\cdots
i_{k}j_{k}}_{il_{2}m_{2}\cdots l_{k}m_{k}}R_{i_{2}j_{2}}^{l_{2}m_{2}}\cdots
R_{i_{k}j_{k}}^{l_{k}m_{k}}\\\
&+\frac{k(k-1)}{2^{k-1}}R^{l_{1}m_{1}}_{uj_{1}}R^{zm_{2}}_{i_{2}j_{2}}\delta^{j_{1}\,i_{2}\,j_{2}\,i_{3}\,j_{3}\,\cdots\,i_{k}\,j_{k}}_{l_{1}m_{1}m_{2}l_{3}m_{3}\cdots
l_{k}m_{k}}R_{i_{3}j_{3}}^{l_{3}m_{3}}\cdots
R_{i_{k}j_{k}}^{l_{k}m_{k}}\end{aligned}$ (38)

for $k\geq 2$.

Then, using the result

$\displaystyle\begin{aligned}
\partial_{u}\hat{R}_{ab}^{cd}=K^{[c}_{e}R_{ab}^{d]e}-2D_{[a}D^{[c}K_{b]}^{d]}\,,\end{aligned}$
(39)

and together with Eq. (19), we can further obtain

$\displaystyle\begin{aligned}
&\partial_{u}\Theta=\partial_{u}\theta+\lambda\partial_{u}K_{b}^{a}\sum^{k_{\text{max}}}_{k=2}\frac{k\alpha_{k}}{2^{k-1}}\delta_{ac_{2}d_{2}\cdots
c_{k}d_{k}}^{ba_{2}b_{2}\cdots
a_{k}b_{k}}\hat{R}^{c_{2}d_{2}}_{a_{2}b_{2}}\cdots\hat{R}^{c_{k}d_{k}}_{a_{k}b_{k}}\\\
&+\lambda
K_{b}^{a}\partial_{u}\hat{R}^{c_{2}d_{2}}_{a_{2}b_{2}}\sum^{k_{\text{max}}}_{k=2}\frac{k(k-1)\alpha_{k}}{2^{k-1}}\delta_{ac_{2}d_{2}c_{3}d_{3}\cdots
c_{k}d_{k}}^{ba_{2}b_{2}a_{3}b_{3}\cdots
a_{k}b_{k}}\hat{R}^{c_{3}d_{3}}_{a_{3}b_{3}}\cdots\hat{R}^{c_{k}d_{k}}_{a_{k}b_{k}}.\end{aligned}$
(40)

Combing the above results, we have

$\displaystyle\begin{aligned} \mathcal{F}&=E_{u}^{z}+\partial_{u}\Theta\\\
&=-K_{ab}K^{ab}+\lambda(\hat{H}^{b}_{a}-H^{b}_{a})\partial_{u}K_{b}^{a}-\lambda
K_{b}^{c}K^{a}_{c}H_{a}^{b}\\\
&+\lambda(2D^{d}K_{a}^{e}-K_{a}^{d}\beta^{e})(2D_{b}K_{c}^{f}-K^{f}_{b}\beta_{c})P_{abc}^{def}\\\
&+\lambda
K_{a}^{d}(K^{e}_{\tilde{e}}R_{bc}^{f\tilde{e}}-2D_{b}D^{e}K_{c}^{f})\hat{P}_{abc}^{def}\,,\end{aligned}$
(41)

in which we have denoted

$\displaystyle\begin{aligned}
H^{b}_{a}&=\sum^{k_{\text{max}}}_{k=2}\frac{k\alpha_{k}}{2^{k-1}}\hat{\delta}_{ac_{2}d_{2}\cdots
c_{k}d_{k}}^{ba_{2}b_{2}\cdots a_{k}b_{k}}R^{c_{2}d_{2}}_{a_{2}b_{2}}\cdots
R^{c_{k}d_{k}}_{a_{k}b_{k}}\,,\\\
\hat{H}^{b}_{a}&=\sum^{k_{\text{max}}}_{k=2}\frac{k\alpha_{k}}{2^{k-1}}\hat{\delta}_{ac_{2}d_{2}\cdots
c_{k}d_{k}}^{ba_{2}b_{2}\cdots
a_{k}b_{k}}\hat{R}^{c_{2}d_{2}}_{a_{2}b_{2}}\cdots\hat{R}^{c_{k}d_{k}}_{a_{k}b_{k}}\,,\\\
P_{abc}^{def}&=\sum^{k_{\text{max}}}_{k=2}\frac{k(k-1)\alpha_{k}}{2^{k-1}}\hat{\delta}_{abcc_{3}d_{3}\cdots
c_{k}d_{k}}^{defa_{c}b_{d}\cdots a_{k}b_{k}}R^{c_{3}d_{3}}_{a_{3}b_{3}}\cdots
R^{c_{k}d_{k}}_{a_{k}b_{k}}\,,\\\
\hat{P}_{abc}^{def}&=\sum^{k_{\text{max}}}_{k=2}\frac{k(k-1)\alpha_{k}}{2^{k-1}}\hat{\delta}_{abcc_{3}d_{3}\cdots
c_{k}d_{k}}^{defa_{c}b_{d}\cdots
a_{k}b_{k}}\hat{R}^{c_{3}d_{3}}_{a_{3}b_{3}}\cdots\hat{R}^{c_{k}d_{k}}_{a_{k}b_{k}}\,,\end{aligned}$
(42)

in which

$\displaystyle\begin{aligned} \hat{\delta}_{a_{1}\cdots a_{i}}^{b_{1}\cdots
b_{i}}=i!\gamma_{[a_{1}}^{b_{1}}\cdots\gamma_{a_{i}]}^{b_{i}}\end{aligned}$
(43)

is the $i$th-order generalized Kronecker tensor on the cross-section $B$.

In the following, we would like to judge the sign of $\mathcal{F}(\lambda)$
when the coupling constant $\lambda$ is regarded as a small parameter. If we
consider the solution $g_{ab}(\lambda)$ which is an analytic function of
$\lambda$, then we can expand $\mathcal{F}(\lambda)$ by $\lambda$,

$\displaystyle\begin{aligned}
\mathcal{F}(\lambda)=\mathcal{F}+\lambda\delta\mathcal{F}+\frac{\lambda^{2}}{2}\delta^{2}\mathcal{F}+\cdots\,,\end{aligned}$
(44)

in which we have introduced the notation

$\displaystyle\begin{aligned}
\delta^{i}\eta(x)=\left.\frac{\partial^{i}\eta(x,\lambda)}{\partial\lambda^{i}}\right|_{\lambda=0}\end{aligned}$
(45)

to denote the $i$th-order variation of the quantity $\eta(x,\lambda)$, and the
symbol without $\lambda$ denotes its counterpart of $\lambda=0$.

Zeroth-order approximation:

First, we consider the zeroth-order approximation of $\lambda$. From Eq. (41),
we can obtain

$\displaystyle\mathcal{F}(\lambda)=-K^{ab}K_{ab}+\mathcal{O}(\lambda)\,.$ (46)

Considering the fact that $K_{ab}$ is a spatial tensor on $B$, we have
$K_{ab}K^{ab}\geq 0$ and therefore $\mathcal{F}(\lambda)\leq 0$ under the
zeroth-order approximation of $\lambda$, which implies that the covariant
entropy bound is satisfied under the zeroth-order approximation. This result
is straightforward because the theory with $\lambda=0$ is just the Einstein
gravity. However, it is worthy noting that there exists a zeroth-order optimal
condition with $K_{ab}=0$ (here we denote $K_{ab}=K_{ab}|_{\lambda=0}$) on the
null hypersurface $L$ such that the zeroth-order term of
$\mathcal{F}(\lambda)$ vanishes and the sign of $\mathcal{F}(\lambda)$ cannot
be determined by the zeroth-order approximation of $\lambda$. In this case, we
need to consider the higher-order approximation of $\lambda$. Thus we next
consider the first-order approximation of $\lambda$ under the zeroth-order
optimal condition $K_{ab}=0$.

First-order approximation:

From Eq. (41), the first-order variation of $\mathcal{F}(\lambda)$ gives

$\displaystyle\begin{aligned} \delta\mathcal{F}&=-2K_{b}^{a}\delta
K_{a}^{b}+(\hat{H}^{b}_{a}-H^{b}_{a})\partial_{u}K_{b}^{a}-K_{b}^{c}K^{a}_{c}H_{a}^{b}\\\
&+(2D^{d}K_{a}^{e}-K_{a}^{d}\beta^{e})(2D_{b}K_{c}^{f}-K^{f}_{b}\beta_{c})P_{abc}^{def}\\\
&+K_{a}^{d}(K^{e}_{\tilde{e}}R_{bc}^{f\tilde{e}}-2D_{b}D^{e}K_{c}^{f})\hat{P}_{abc}^{def}\,.\end{aligned}$
(47)

Noting that the zeroth-order optimal condition $K_{ab}=0$ on $L$ also implies

$\displaystyle\partial_{u}K^{b}_{a}=D_{c}K^{b}_{a}=K_{a}^{b}=0$ (48)

on $L$, we have $\delta\mathcal{F}=0$ under the zeroth-order approximation of
$\lambda$. Then, we have

$\displaystyle\begin{aligned}
\mathcal{F}(\lambda)=\mathcal{F}+\lambda\delta\mathcal{F}+\mathcal{O}(\lambda^{2})=\mathcal{O}(\lambda^{2})\,,\end{aligned}$
(49)

which indicates that the generalized covariant entropy bound is satisfied
under the first-order approximation of $\lambda$. Moreover, this also implies
that the sign of $\mathcal{F}(\lambda)$ should be determined by the higher-
order terms of $\lambda$.

Moreover, from Eq. (LABEL:expR), we can see that

$\displaystyle\begin{aligned} R_{ij}^{kl}=\hat{R}_{ij}^{kl}\end{aligned}$ (50)

on the null hypersurface under the zeroth-order optimal condition
$K_{i}^{j}=0$, which implies that $H^{a}_{b}=\hat{H}^{a}_{b}$ and therefore we
have $\delta\mathcal{F}=0$. Then, we have

$\displaystyle\begin{aligned}
\mathcal{F}(\lambda)=\mathcal{F}+\lambda\delta\mathcal{F}+\mathcal{O}(\lambda^{2})=\mathcal{O}(\lambda^{2})\,,\end{aligned}$
(51)

which indicates that the generalized covariant entropy bound is satisfied
under the first-order approximation of $\lambda$. Moreover, this also implies
that the sign of $\mathcal{F}(\lambda)$ should be determined by the higher-
order terms of $\lambda$.

Second-order approximation:

Considering the zeroth-order optimal condition
$K^{b}_{a}=\partial_{u}K^{b}_{a}=D_{c}K^{b}_{a}=0$ on $L$, the second-order
variation of $\mathcal{F}(\lambda)$ can be further obtained

$\displaystyle\begin{aligned} \delta^{2}\mathcal{F}&=-2\delta K_{ab}\delta
K^{ab}+(\hat{H}^{b}_{a}-H^{b}_{a})\partial_{u}\delta
K_{b}^{a}\,.\end{aligned}$ (52)

Moreover, from Eq. (LABEL:expR), we can see that

$\displaystyle\begin{aligned} R_{ij}^{kl}=\hat{R}_{ij}^{kl}\end{aligned}$ (53)

on the null hypersurface under the zeroth-order optimal condition, which
implies that $H^{a}_{b}=\hat{H}^{a}_{b}$ on $L$ and the second-term of Eq.
(52) vanishes. Then, we have

$\displaystyle\begin{aligned} \mathcal{F}(\lambda)=-\lambda^{2}\delta
K_{ab}\delta K^{ab}+\mathcal{O}(\lambda^{3})\,.\end{aligned}$ (54)

This result shows that $\mathcal{F}(\lambda)\leq 0$ under the second-order
approximation of $\lambda$. Similarly, there also exists a second-order
optimal condition with $\delta K_{ab}=0$ on the null hypersurface $L$ such
that the higher-order approximation should be further considered to judge the
sign of $\mathcal{F}(\lambda)$.

Third-order approximation:

Under the zeroth-order and second-order optimal conditions, it is not hard to
see

$\displaystyle\begin{aligned}
K_{a}^{b}=\partial_{u}K_{a}^{b}=&D_{c}K_{a}^{b}=0\,,\quad\delta
K_{a}^{b}=\partial_{u}\delta K_{a}^{b}=D_{c}\delta K_{a}^{b}=0\,,\\\
&H_{a}^{b}=\hat{H}_{a}^{b}\,,\quad\quad\delta
H_{a}^{b}=\delta\hat{H}_{a}^{b}\,\end{aligned}$ (55)

on the null hypersurface $L$. Then, we can further find that
$\delta^{3}\mathcal{F}=0$, and therefore we have
$\mathcal{F}(\lambda)=\mathcal{O}(\lambda^{4})$, which implies that
$\mathcal{F}(\lambda)$ vanishes under the third-order approximation and we
need further consideration for the higher-order approximation.

$n$th-order approximation

In the following, we would like to prove $\mathcal{F}(\lambda)\geq 0$ under
the $n$th-order approximation when the $(n-1)$th-order approximation of
$\mathcal{F}(\lambda)$ vanishes based on the mathematical induction. By
concluding the first three-order results, it is equivalent to proving the
following proposition: under the $n$th-order approximation of $\lambda$, we
have

$\displaystyle\mathcal{F}(\lambda)=\frac{\lambda^{n}}{n!}\delta^{n}\mathcal{F}\geq
0\,,$ (56)

when the $(n-1)$th-order approximation of $\mathcal{F}(\lambda)$ vanishes,
i.e., $\delta^{i}\mathcal{F}=0$ for $i\leq n-1$. Then, $n$th-order
approximation of $\mathcal{F}(\lambda)$ vanishes demands

$\displaystyle\begin{aligned} \delta^{i}K_{a}^{b}=0\,,\quad\quad\text{for
$i\leq[n/2]$}\,,\end{aligned}$ (57)

on the null hypersurface $L$, where $[n/2]$ denotes the integer part of $n/2$.

Proof. Obviously, the proposition is satisfied for $k=0,1,2$.

(Case of $n=2m$.) We assume that the proposition is satisfied when $n=2m$ for
$m\geq 0$, i.e., we have

$\displaystyle\begin{aligned} \delta^{2m}\mathcal{F}\geq 0\,,\end{aligned}$
(58)

when $\delta^{i}\mathcal{F}=0$ for $i\leq 2m-1$, and $\mathcal{F}(\lambda)$
vanishes under the $2m$th-order approximation (i.e., $\delta^{i}\mathcal{F}=0$
for $i\leq 2m$) demands

$\displaystyle\begin{aligned} \delta^{i}K_{a}^{b}=0\,,\quad\quad\text{for
$i\leq m$}\end{aligned}$ (59)

on $L$, which implies that $K_{a}^{b}(\lambda)=0$ vanishes under the $m$th-
order approximation of $\lambda$ and thus
$H_{a}^{b}(\lambda)=\hat{H}_{a}^{b}(\lambda)$ under the $m$th-order
approximation. Therefore, we have

$\displaystyle\begin{aligned}
&K_{a}^{b}(\lambda)=\lambda^{m+1}\tilde{K}_{a}^{b}\,,\quad\quad\hat{H}^{b}_{a}(\lambda)-H^{b}_{a}(\lambda)=\lambda^{m+1}\mathcal{H}^{b}_{a}\,,\end{aligned}$
(60)

in which

$\displaystyle\begin{aligned}
\tilde{K}_{a}^{b}&=\sum^{\infty}_{i=m+1}\frac{\lambda^{i-m-1}}{i!}\delta^{i}K_{a}^{b}\,,\\\
\mathcal{H}_{a}^{b}&=\sum^{\infty}_{i=m+1}\frac{\lambda^{i-m-1}}{i!}(\delta^{i}\hat{H}_{a}^{b}-\delta^{i}H_{a}^{b})\,.\end{aligned}$
(61)

With a straightforward calculation, we have

$\displaystyle\begin{aligned}
\mathcal{F}(\lambda)&=-\lambda^{2m+2}\tilde{K}_{ab}\tilde{K}^{ab}+\lambda^{2m+3}\mathcal{H}^{b}_{a}\partial_{u}\tilde{K}_{b}^{a}-\lambda^{2m+3}\tilde{K}_{b}^{c}\tilde{K}^{a}_{c}H_{a}^{b}\\\
&+\lambda^{2m+3}(2D^{d}\tilde{K}_{a}^{e}-\tilde{K}_{a}^{d}\beta^{e})(2D_{b}\tilde{K}_{c}^{f}-\tilde{K}^{f}_{b}\beta_{c})P_{abc}^{def}\\\
&+\lambda^{2m+3}\tilde{K}_{a}^{d}(\tilde{K}^{e}_{\tilde{e}}R_{bc}^{f\tilde{e}}-2D_{b}D^{e}\tilde{K}_{c}^{f})\hat{P}_{abc}^{def}\,.\end{aligned}$
(62)

This result shows that $\mathcal{F}(\lambda)=0$ under the $(2m+1)$th-order
approximation of $\lambda$ under the optimal condition (59), this also means
that $\delta^{2m+1}\mathcal{F}=0$ would not lead to additional conditions than
the $2m$th-order case, which actually gives the saturation condition (57) with
$n=2m+1$. These show that the proposition with $n=2m+1$ is satisfied.

(Case of $n=2m+1$.) Assume that the proposition is satisfied when $n=2m+1$ for
$m\geq 0$, and the $(2m+1)$th-order approximation of $\mathcal{F}(\lambda)$
vanishes demanding

$\displaystyle\begin{aligned} \delta^{i}K_{a}^{b}=0\,,\quad\quad\text{for
$i\leq m$}\,,\end{aligned}$ (63)

which also leads to the results in Eq. (LABEL:KH) and therefore
$\mathcal{F}(\lambda)$ is given by Eq. (62). That is to say, under the
$(2m+2)$th-order approximation, we have

$\displaystyle\begin{aligned}
\mathcal{F}(\lambda)=-\frac{\lambda^{2m+2}}{[(m+1)!]^{2}}\delta^{m+1}K_{a}^{b}\delta^{m+1}K^{a}_{b}+\mathcal{O}(\lambda^{2m+3})\,,\end{aligned}$
(64)

which implies that $\mathcal{F}(\lambda)\leq 0$ $(\delta^{2m+2}\mathcal{F}\leq
0)$ under the $(2m+2)$th-order approximation. Together with the condition
(63), the saturation of $\delta^{2m+2}\mathcal{F}$ demands

$\displaystyle\begin{aligned} \delta^{i}K_{a}^{b}=0\,,\quad\quad\text{for
$i\leq m+1$}\,.\end{aligned}$ (65)

This is actually the proposition with $n=2m+2$, i.e., we have completed the
proof. $\Box$

The above result shows $\mathcal{F}(\lambda)$ is always nonpositive and thus
the generalized covariant entropy bound associated with the JM entropy is
valid under any higher-order approximation of $\lambda$.

##  V Conclusion and discussion

In this paper, we consider the generalized covariant entropy bound for the
theory in which the Einstein gravity is perturbed by the higher-order Lovelock
terms and introduce a small parameter $\lambda$ to characterize these
perturbations. After considering the linearized second law of black holes in
Lanczos-Lovelock gravity, the entropy bound in this theory is naturally
proposed by replacing the Bekenstein-Hawking entropy with the JM entropy.
Then, we showed that the key point to examine the validity of covariant
entropy bound is to judge the sign of the quantity
$\mathcal{F}=H_{ab}k^{a}k^{b}+k^{a}\nabla_{a}\Theta$, and the entropy bound is
satisfied if $\mathcal{F}\leq 0$. After assuming two physical assumptions and
that the metric $g_{ab}$ is an analytic function of $\lambda$, we illustrated
that the dominant term of $\mathcal{F}(\lambda)$ is always nonpositive based
on the mathematical induction, i.e., the generalized covariant entropy bound
is valid under any higher-order approximation of $\lambda$. This indicates
that the entropy bound using the JM entropy is strictly satisfied under the
perturbation level of the higher-order Lovelock terms. From the discussion in
Sec. V D of Ref. Matsuda:2020yvl , we can see that the above result also
indicates the validity of the generalized second law under the higher-order
approximation of $\lambda$ for the theory where the Einstein gravity is
perturbed by the higher-order Lovelock terms, this is a different result from
the linearized second law of Lanczos-Lovelock gravity.

From the calculations presented in this paper, it is not hard to check that if
we replace the Bekenstein-Hawking with Wald entropy formula instead of JM
entropy formula in the entropy bound, we cannot show the nonnegativity of
$\mathcal{F}$ only using the assumptions given by the paper. This implies that
the covariant entropy bound might be used to select the black hole entropy of
the gravitational theory. Moreover, it is worth noting that our result is only
suitable for the case where the higher-order Lovelock terms are regarded as
some small corrections to Einstein gravity, what about the non-perturbation
cases? From the discussion in Sec. III, the key point to examine the entropy
bound is also to check the sign of $\mathcal{F}$ given by Eq. (41). However,
due to the complexity of the expression, it is difficult for us to judge its
sign directly only based on the setups and assumptions in our paper. One of
our future work is going to consider these cases. Furthermore, it is also
interesting to extend the discussion into the theory in which the Einstein
gravity is perturbed by other higher-curvature terms.

## Acknowledgement

Jie Jiang is supported by the GuangDong Basic and Applied Basic Research
Foundation with Grant no. 217200003 and the Talents Introduction Foundation of
Beijing Normal University with Grant no. 310432102. Ming Zhang is supported by
the National Natural Science Foundation of China with Grant No. 12005080.

## References

  * (1) J. D. Bekenstein, “A Universal Upper Bound on the Entropy to Energy Ratio for Bounded Systems,” Phys. Rev. D 23, 287 (1981).
  * (2) R. Bousso, “A Covariant entropy conjecture,” JHEP 07, 004 (1999).
  * (3) W. Fischler and L. Susskind, “Holography and cosmology,” arXiv:hep-th/9806039 [hep-th].
  * (4) S. Gao and J. P. S. Lemos, “The Covariant entropy bound in gravitational collapse,” JHEP 04, 017 (2004).
  * (5) S. Gao and J. P. S. Lemos, “Local conditions for the generalized covariant entropy bound,” Phys. Rev. D 71, 084010 (2005).
  * (6) S. Gao and X. Wu, “Proof of the entropy bound on dynamical horizons,” JHEP 08, 005 (2008).
  * (7) A. Pesci, “From Unruh temperature to generalized Bousso bound,” Class. Quant. Grav. 24, 6219-6226 (2007).
  * (8) A. Pesci, “On the statistical-mechanical meaning of Bousso bound,” Class. Quant. Grav. 25, 125005 (2008).
  * (9) E. E. Flanagan, D. Marolf and R. M. Wald, “Proof of classical versions of the Bousso entropy bound and of the generalized second law,” Phys. Rev. D 62, 084035 (2000).
  * (10) R. Bousso, E. E. Flanagan and D. Marolf, “Simple sufficient conditions for the generalized covariant entropy bound,” Phys. Rev. D 68, 064001 (2003).
  * (11) A. Strominger and D. M. Thompson, “A Quantum Bousso bound,” Phys. Rev. D 70, 044007 (2004).
  * (12) B. Zwiebach, “Curvature Squared Terms and String Theories,” Phys. Lett. 156B, 315 (1985).
  * (13) D. J. Gross and E. Witten, “Superstring Modifications of Einstein’s Equations,” Nucl. Phys. B 277, 1 (1986).
  * (14) A. Sen, “Black Hole Entropy Function, Attractors and Precision Counting of Microstates,” Gen. Rel. Grav. 40, 2249 (2008).
  * (15) A. Dabholkar and S. Nampuri, “Quantum black holes,” Lect. Notes Phys. 851, 165 (2012).
  * (16) C. Lanczos, “A Remarkable property of the Riemann Christoffel tensor in four dimensions,” Annals Math. 39, 842 (1938).
  * (17) D. Lovelock, “The Einstein tensor and its generalizations,” J. Math. Phys. 12, 498 (1971).
  * (18) B. Zwiebach, “ Curvature squared terms and string theories,” Phys. Lett. B156, 315 (1985).
  * (19) A. D. Kovacs and H. S. Reall, “Well-posed formulation of Lovelock and Horndeski theories,” Phys. Rev. D 101, 124003 (2020).
  * (20) T. Matsuda and S. Mukohyama, “Covariant entropy bound beyond general relativity,” Phys. Rev. D 103, no.2, 024002 (2021).
  * (21) R. M. Wald, “Black hole entropy is the Noether charge,” Phys. Rev. D 48, R3427 (1993).
  * (22) V. Iyer and R. M. Wald, “Some properties of Noether charge and a proposal for dynamical black hole,” Phys. Rev. D 50, 846 (1994).
  * (23) T. Jacobson and R. C. Myers, “Black hole entropy and higher curvature interactions,” Phys. Rev. Lett. 70, 3684-3687 (1993).
  * (24) A. Chatterjee and S. Sarkar, “Physical Process First Law and Increase of Horizon Entropy for Black Holes in Einstein Gauss-Bonnet Gravity,” Phys. Rev. Lett. 108, 091301 (2012).
  * (25) S. Kolekar, T. Padmanabhan, and S. Sarkar, “Entropy increase during physical processes for black holes in Lanczos-Lovelock gravity,” Phys. Rev. D 86, 021501 (2012).
  * (26) S. Sarkar and A. C. Wall, “Generalized second law at linear order for actions that are functions of Lovelock densities,” Phys. Rev. D 88, 044017 (2013).
  * (27) J. Jiang and M. Zhang, “Entropy increases at linear order in scalar-hairy Lovelock gravity,” JHEP 04, 148 (2020)