††thanks: Email<EMAIL_ADDRESS>Email:
<EMAIL_ADDRESS>Email<EMAIL_ADDRESS>Email:
<EMAIL_ADDRESS>Email<EMAIL_ADDRESS>Email:
<EMAIL_ADDRESS>Email<EMAIL_ADDRESS>
# Study of neutron density fluctuation and neutron-proton correlation in Au+Au
collisions using PYTHIA8/Angantyr

Zuman Zhang  School of Physics and Mechanical Electrical & Engineering, Hubei
University of Education, Wuhan 430205, China  Institute of Theoretical
Physics, Hubei University of Education, Wuhan 430205, China  Key Laboratory of
Quark and Lepton Physics (MOE), Central China Normal University, Wuhan 430079,
China  Sha Li  School of Physics and Mechanical Electrical & Engineering,
Hubei University of Education, Wuhan 430205, China  Ning Yu  School of Physics
and Mechanical Electrical & Engineering, Hubei University of Education, Wuhan
430205, China  Institute of Theoretical Physics, Hubei University of
Education, Wuhan 430205, China  Key Laboratory of Quark and Lepton Physics
(MOE), Central China Normal University, Wuhan 430079, China  Jianping Lin
School of Physics and Mechanical Electrical & Engineering, Hubei University of
Education, Wuhan 430205, China  Shuang Li  College of Science, China Three
Gorges University, Yichang 443002, China  Center for Astronomy and Space
Sciences, China Three Gorges University, Yichang 443002, China  Siyu Tang
School of Mathematical & Physical Sciences, Wuhan Textile University, Wuhan
430200, China  Daicui Zhou  Key Laboratory of Quark and Lepton Physics (MOE),
Central China Normal University, Wuhan 430079, China

###### Abstract

Utilizing the PYTHIA8 Angantyr model, which incorporates the multiple-parton
interactions (MPI) based color reconnection (CR) mechanism, we study the
relative neutron density fluctuation and neutron-proton correlation in Au+Au
collisions at $\sqrt{s_{\text{NN}}}$ = 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4,
and 200 GeV. In this study, we have not only delved into the dependence of
these two remarkable observations on rapidity, centrality, and energy, but
also presented an analysis of their interplay with the MPI and CR. Our results
have shown that the light nuclei yield ratio of proton, deuteron, and triton,
expressed by the elegant expression $N_{t}N_{p}/N_{d}^{2}$, remains unchanged
even as the rapidity coverage and collision centrality increase.
Interestingly, we have also revealed that the effect of CR is entirely
dependent on the presence of MPI; CR has no impact on the yield ratio if MPI
is off. Our findings further demonstrate that the light nuclei yield ratio
experiences a slight increase with increasing collision energy as predicted by
the PYTHIA8 Angantyr model, but it cannot describe the non-monotonic trend
observed by the STAR experiment. Based on the Angantyr model simulation
results, it is essential not to overlook the correlation between neutron and
proton fluctuations. The Angantyr model is a good baseline for studying
collisions in the absence of a Quark-Gluon Plasma (QGP) system, given its lack
of flow and jet quenching.

PACS numbers

`relativistic heavy-ion collisions, quark deconfinement,`quark-gluon plasma
production, phase transitions

##  I Introduction

The creation of a state of matter known as the quark-gluon plasma (QGP), a
mixture of deconfined quarks and gluons, is believed to occur in heavy-ion
collisions at ultra-relativistic energies at extremely high temperatures
and/or densities. Understanding the Quantum Chromodynamics (QCD) phase diagram
of strongly interacting matter is crucial in the field of nuclear physics. The
QCD phase diagram can be represented in a two-dimensional graph of temperature
(T) versus baryon chemical potential ($\mu_{\mathrm{B}}$). According to
lattice QCD calculations, the transition from the hadronic phase to the QGP
occurs as a smooth crossover at low values of $\mu_{\mathrm{B}}$
ref:QCDphaseplot ; Aoki:2006we . However, at finite values of
$\mu_{\mathrm{B}}$, it is predicted to be a first-order phase transition based
on QCD model calculations Endrodi:2011gv ; ref:LQCDCplot ; ref:lqcd03 ;
ref:lqcd05 . If these predictions hold true, there must exist a QCD critical
point, marking the endpoint of the first-order phase boundary. Despite ongoing
theoretical discussions on the location and even the existence of the QCD
critical point, relativistic heavy-ion collisions provide a controlled means
to explore the QCD phase structure and search for the critical point
ref:SPS_result ; ref:RHIC_result ; Gupta:2011wh ; Luo:2017faz ; Bzdak:2019pkr
; Aggarwal:2010wy ; Luo:2015doi ; Adamczyk:2014fia ; Adamczyk:2013dal ;
Adamczyk:2017wsl ; Adam:2020unf .

As the system approaches the QCD critical point, the correlation length and
density fluctuations experience an increase. The fluctuations of conserved
quantities, such as net-baryon, net-charge, and net-strangeness, are acutely
influenced by the correlation length. The STAR experiment has conducted
measurements of high-order cumulants and second-order off-diagonal cumulants
of net-proton, net-charge, and net-kaon multiplicity distributions RN182 ;
RN183 ; RN184 ; RN185  in Au+Au collisions across a broad range of energies,
from as low as $\sqrt{s_{\text{NN}}}$ = 7.7 GeV up to 200 GeV in the most
central collisions, a non-monotonic behavior of the fourth-order net-proton
cumulant ratio was observed, with a minimum around 19.6 GeV. The non-monotonic
behavior observed in the fourth-order net-proton cumulant ratio in the STAR
experiment cannot be explained by existing model calculations unless the
physics of the QCD critical point is considered.

Moreover, as the correlation length increases and instability spinodal domains
form, critical fluctuations and first-order phase transitions can give rise to
significant baryon density fluctuations. It is predicted that the production
of light nuclei is sensitive to the baryon density fluctuations, and thus can
be used to probe the QCD phase transition in heavy-ion collisions Sun:2017xrx
; Deng:2020zxo ; Shuryak:2019ikv ; Yu:2018kvh ; Shao:2019xpj ; RN9 . For
instance, the light nuclei yield ratio $N_{p}N_{t}/N_{d}^{2}$ of produced
proton ($p$), deuteron ($d$), and triton ($t$) can be described by the
relative neutron density fluctuation $\langle(\delta n)^{2}\rangle/\langle
n\rangle^{2}$ (also denoted as $\Delta n$ in Ref. Sun:2017xrx ).
Interestingly, it was observed that the yield ratio and the extracted neutron
density fluctuation in central Au+Au collisions measured by STAR experiment
show clear non-monotonic energy dependence, with a peak at
$\sqrt{s_{\text{NN}}}$ = $20\sim 30$ GeV Zhang:2019wun ; Zhang:2020ewj .

In this paper, we studied the light nuclei yield ratio in Au+Au collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4, and 200 GeV
from the framework of PYTHIA8 Angantyr model. Angantyr adds a framework
wherein AA-collisions can be constructed as a superposition of binary nucleon-
nucleon collisions as implemented in the PYTHIA8 event generator. Our paper is
organized as following: In Sec. II.1, we give a brief introduction to the
PYTHIA8 Angantyr model. In Sec. II.2, the relationship between the relative
neutron density fluctuation and the light nuclei yield ratio in heavy ion
collisions is given. In Sec. III, we present the behavior of neutron density
fluctuation and neutron-proton correlation ($\alpha$) . Finally, the summary
will be given in Sec. IV.

##  II Event generation and definition of light nuclei yield ratio

###  II.1 PYTHIA8 (Angantyr) model

PYTHIA Sjostrand:2006za  is an event generator that is extensively and
successfully used for the study of proton-proton and proton-lepton collisions.
In pp collisions, Multi-Parton Interaction (MPI) is generated under the
assumption that every partonic interaction is almost independent. PYTHIA8
natively does not support heavy-ion systems. Recently, PYTHIA8 Angantyr model
Bierlich:2018xfw  extrapolates the pp dynamics into heavy-ion collisions using
the PYTHIA8 event generator, enables the study of heavy nuclei collisions,
namely proton-nuclei (pA) and nuclei-nuclei (AA). The Angantyr model combines
several nucleon-nucleon collisions into one heavy-ion collision. It is a
combination of many-body physics (theoretical) models suitable for producing
hard and soft interactions, initial and final-state parton showers, particle
fragmentations, multi-partonic interactions, color reconnection (CR)
mechanisms, and decay topologies. But it doesn’t include any mechanism of QGP
medium believed to be produced in AA collisions.

In the current version of the PYTHIA8 Angantyr model Sjostrand:2008za , in a
heavy-ion collision, each projectile nucleon can interact with several target
nucleons, and the number of participant nucleons determined by Glauber model.
This model added several algorithms to distinguish different types of nucleon-
nucleon interactions, such as elastic, diffractive, and absorptive. It’s
supposed to well describe the final-state properties, such as multiplicity and
transverse momentum distributions in AA collisions Abdel:2022 ; Samsul:2022 .

We use 8.308 version of PYTHIA8 in this work. The simulation is considered
with different PYTHIA tunes using MPI and CR configurations. About one million
events for each collision energy were generated for the Au + Au collisions. In
the PYTHIA8 Angantyr model, due to the both density and density fluctuation
are dependent on system volume. To remove the system volume effect, we use two
dimensionless statistical quantities $\langle(\delta p)\rangle/\langle
p\rangle$ and $\langle(\delta n)\rangle/\langle n\rangle$ (Sec. II.2). The
nucleons are extracted for different rapidity ranges and centralities. We
define the centrality intervals which are based on the summed transverse
energy ($\sum E_{T}$) in the pseudorapidity interval [-0.5, 0.5].

###  II.2 Definition of light nuclei yield ratio

Based on the references cited as RN9 ; RN172 ; RN2022 , it has been proposed
that the creation of light nuclei occurs through the process of nucleon
coalescence. By this mechanism, protons and neutrons come together to form
deuterons and tritons. If we neglect the binding energy, the number of these
particles can be expressed as

$N_{d}\approx\frac{3}{2^{1/2}}\left(\frac{2\pi}{mT}\right)^{3/2}N_{p}\langle
n\rangle(1+\alpha)$ (1)

$N_{t}\approx\frac{3^{3/2}}{4}\left(\frac{2\pi}{mT}\right)^{3}N_{p}\langle
n\rangle^{2}(1+\Delta n+2\alpha),$ (2)

then, to eliminate the energy dependence on the local effective temperature
($T$) at coalescence, we have the light nuclei ratio

$\frac{N_{t}N_{p}}{N_{d}^{2}}=g\frac{1+\Delta n+2\alpha}{(1+\alpha)^{2}},$ (3)

where $N_{p}$, $N_{d}$, and $N_{t}$ represent the number of protons, neutrons,
and tritons respectively. $m$ stands for the mass of a nucleon. $\Delta
n=\langle(\delta n)^{2}\rangle/\langle n\rangle^{2}$ is the dimensionless
relative neutron density fluctuation. The standard deviation and mean value of
neutron production are represented by $\langle(\delta n)\rangle$ and $\langle
n\rangle$ in our work. $\alpha=\langle\delta n\delta p\rangle/(\langle
n\rangle\langle p\rangle)$ is the neutron and proton number density
correlation coefficient. $g$ = $\frac{3^{3/2}}{4}/(\frac{3}{2^{1/2}})^{2}$ =
$\frac{1}{2\sqrt{3}}$.

If we assume correlation coefficient between neutron and proton density is
small, then $\alpha$ can be approximated as zero ($\alpha\approx 0$), Eq. (3)
can be rewritten as

$\frac{N_{t}N_{p}}{N_{d}^{2}}=\frac{1+\Delta n}{2\sqrt{3}}.$ (4)

##  III Results and Discussions

Figure 1:  In the top panels of the figure, various dimensionless statistics
are displayed, including $\langle(\delta p)\rangle/\langle p\rangle$,
$\langle(\delta n)\rangle/\langle n\rangle$, and $\alpha$, which are
associated with 0-10% Au+Au collisions at center-of-mass energies of
$\sqrt{s_{\text{NN}}}=$ 14.5 and 39 GeV . Meanwhile, the bottom panels depict
the ratio of yields for light nuclei, $N_{t}N_{p}/N_{d}^{2}$, which is
calculated using the data from the top panels, represented as solid circles in
accordance with Eq. (3), and as a dash-dot line in accordance with Eq. (4).

Except for the physical parameters utilized in the model calculations
mentioned in last paragraph of Sec. II.1, we have the option to incorporate
the multiple-parton interactions (MPI) based color reconnection (CR) mechanism
by switching on/off the ColourReconnection:reconnect and PartonLevel:MPI.

The results in Figure 1, 2, 3, 4, and 6 have been investigated using the
PYTHIA8 Angantyr model with the multiple-parton interactions based color
reconnection mechanism.

The top panels depicted in Figure 1 illustrate the rapidity dependence of the
dimensionless statistics $\langle(\delta p)\rangle/\langle p\rangle$ and
$\langle(\delta n)\rangle/\langle n\rangle$ for 0-10$\%$ Au+Au collisions at
$\sqrt{s_{\text{NN}}}=$ 14.5 and 39 GeV. As the rapidity coverage is
increased, the relative fluctuations in nucleon density ($\langle(\delta
p)\rangle/\langle p\rangle$ and $\langle(\delta n)\rangle/\langle n\rangle$)
displayed in the top panels of Figure 1 appear to decrease. As the rapidity
coverage is expanded, the nucleon density fluctuations appear to converge
towards a constant value. The $\langle(\delta p)\rangle/\langle p\rangle$ and
$\langle(\delta n)\rangle/\langle n\rangle$ are consistent in Figure 1. The
correlation $\alpha$ have similar trend as relative fluctuations in nucleon
density with rapidity coverage for 0-10$\%$ Au+Au collisions at both
$\sqrt{s_{\text{NN}}}=$ 14.5 and 39 GeV. We plot the line of $1/2\sqrt{3}$,
which means both density fluctuation and correlation disappeared. From Eq.
(3), the light nuclei ratio $N_{t}N_{p}/N_{d}^{2}$ can be calculated. The
rapidity dependence of this ratio, the solid circles, is shown in the bottom
panel of Figure 1. Both the ratio for 0-10$\%$ Au+Au collisions at
$\sqrt{s_{\text{NN}}}=$ 14.5 and 39 GeV are lower than the line of
$1/2\sqrt{3}$ with rapidity distribution. One the other hand, the light nuclei
ratio with the dash-dot lines ($\alpha\approx 0$) higher than the line of
$1/2\sqrt{3}$. So, the $\alpha$ can not be neglected in the central collision
to calculate the relative neutron density fluctuation from the light nuclei
yield ratio with Eq. (4). Results from other collision energies at central
collision are similar.

Figure 2:  In the top panels of the figure, various dimensionless statistics
are displayed, including $\langle(\delta p)\rangle/\langle p\rangle$,
$\langle(\delta n)\rangle/\langle n\rangle$ and $\alpha$, which are associated
with 60-80% Au+Au collisions at center-of-mass energies of
$\sqrt{s_{\text{NN}}}=$ 14.5 and 39 GeV . Meanwhile, the bottom panels depict
the ratio of yields for light nuclei, $N_{t}N_{p}/N_{d}^{2}$, which is
calculated using the data from the top panels, represented as solid circles in
accordance with Eq. (3), and as a dash-dot line in accordance with Eq. (4).

The rapidity dependence of $\langle(\delta p)\rangle/\langle p\rangle$ and
$\langle(\delta n)\rangle/\langle n\rangle$ for 60-80$\%$ Au+Au collisions at
$\sqrt{s_{\text{NN}}}=$ 14.5 and 39 GeV in the top panels, and related light
nuclei ratio in the bottom panels are shown in Figure 2. Similar to the
results of the central collisions, $\langle(\delta p)\rangle/\langle p\rangle$
and $\langle(\delta n)\rangle/\langle n\rangle$ decrease with the rapidity
coverage is increased. The correlation $\alpha$ is dependent on rapidity
coverage and has similar trend as relative fluctuations in nucleon density. We
present the light nuclei yield ratios, calculated using Eq. (3) and Eq. (4),
in the bottom lower section of Figure 2. Our results reveal the exclusion of
the $\alpha$ parameter leads to a significant increase in the light nuclei
yield ratio. In other words, the omission of $\alpha$ results in a smaller
relative neutron density fluctuation with the light nuclei yield ratio.

Figure 3:  The dimensionless statistics $\langle(\delta p)\rangle/\langle
p\rangle$, $\langle(\delta n)\rangle/\langle n\rangle$ and $\alpha$ for Au+Au
collisions at center-of-mass energies of $\sqrt{s_{\text{NN}}}=$ 14.5 and 39
GeV, with rapidities confined to the range $|y|\leq 0.5$. In addition, the
ratio of yields for light nuclei, $N_{t}N_{p}/N_{d}^{2}$, is shown in the same
figure, represented as solid circles in accordance with Eq. (3), and as a
dash-dot line in accordance with Eq. (4).

The top panels of Figure 3 show the centrality dependent $\langle(\delta
p)\rangle/\langle p\rangle$ and $\langle(\delta n)\rangle/\langle n\rangle$
for Au+Au collisions at $\sqrt{s_{\text{NN}}}=$ 14.5 and 39 GeV with rapidity
coverage of $|y|\leq 0.5$. All two quantities are flat from central to
peripheral collision. At the bottom of Figure 3, it can be found that the two
ratios given by Eq. (3) and Eq. (4) are flat at central, mid-central, and
peripheral collisions. It is noted that the relative nucleon density
fluctuation cannot be extracted directly from the light nuclei ratio, and the
effects of neutron-proton correlation $\alpha$ must be considered in different
centralities of collisions.

Figure 4:  The PYTHIA8 Angantyr model was used to examine the collision energy
dependence of the light nuclei yield ratio $N_{t}N_{p}/N_{d}^{2}$ in Au+Au
collisions with $|y|\leq 0.5$. The results revealed solid circles and squares
for 0-10% central and 60-80% peripheral collisions, respectively.
Additionally, the dash-dot lines displayed the corresponding results with
vanished $\alpha$.

The collision energy dependence of the $N_{t}N_{p}/N_{d}^{2}$ from 0-10%
central to 60-80% peripheral Au+Au collisions with $|y|\leq 0.5$ are shown in
Figure 4. It is clear from the Figure 4 that the light nuclei yield ratio
increases slightly with increasing collision energy from the PYTHIA8 Angantyr
model. The light nuclei yield ratio of peripheral collisions is similar to
that of central collisions, and both results are lower than $1/2\sqrt{3}$. The
results of vanished $\alpha$ decreases with increasing collision energy are
also shown as dash-dot lines. In Au+Au collisions, both in central (0-10%) and
peripheral (60-80%) collisions, the yield ratio of light nuclei is not highly
consistent with each other after the disappearance of $\alpha$.

Figure 5:  The PYTHIA8 Angantyr model was used to examine the collision energy
dependence of the light nuclei yield ratio $N_{t}N_{p}/N_{d}^{2}$ in Au+Au
collisions with $|y|\leq 0.5$ in different PYTHIA8 Angantyr model tunes. The
results shown in upper and lower panel related to 0-10% central and 60-80%
peripheral collisions, respectively.

To see the effect of different PYTHIA tunes, we consider the following
configurations: MPI with CR, No CR, No MPI, and both MPI and CR off. Figure 5
shows the collision energy dependence of the light nuclei yield ratio
$N_{t}N_{p}/N_{d}^{2}$ from 0-10% central and 60-80% peripheral Au+Au
collisions with $|y|\leq 0.5$ in different PYTHIA8 Angantyr model tunes. It is
also clear from this figure that the light nuclei yield ratio increases
slightly with increasing collision energy from the PYTHIA8 Angantyr model. But
the results of vanished $\alpha$ decreases with increasing collision energy
are also shown as dash-dot lines. At Au+Au collisions, in different PYTHIA8
Angantyr model tunes, there is no effect of CR if MPI is off.

Figure 6:  The yield ratio of $N_{t}N_{p}/N_{d}^{2}$ from the PYTHIA8 Angantyr
model, which takes into account the collision energy and centrality, was
analyzed in the context of $|y|\leq 0.5$. Represented as dash lines from the
PYTHIA8 Angantyr model, and as the solid lines from AMPT model RN2022 .
PYTHIA8 Angantyr model is under the MPI and RC mode. Solid circles are the
results from STAR detector at 0-10% central Au+Au collision Zhang:2020ewj ;
2209.08058 . Additionally, open squares were used to present the findings from
NA49 in central Pb+Pb collision RN9 ; RN19 .

Figure 6 presents the experimental results of the STAR detector at 0-10%
central Au+Au collision Zhang:2020ewj ; 2209.08058  , NA49 at central Pb+Pb
collision RN9 ; RN19 , then compared with results from the AMPT model and
PYTHIA8 Angantyr model. In experimental results, It is evident that there is a
non-monotonic energy dependence. The yield ratio of light nuclei peaks at
$\sqrt{s_{\text{NN}}}=20\sim 30$ GeV, which implies the most significant
relative neutron density fluctuations in this energy range. Both AMPT model
and PYTHIA8 Angantyr model results are lower than experimental results;
moreover, the AMPT model results which includes QGP medium mechanisms and does
not consider critical physics are higher than PYTHIA8 Angantyr model results.
Due to the absence of critical physics and QGP medium mechanisms, the PYTHIA8
Angantyr model is unable to describe this non-monotonic energy dependence.

##  IV Summary

In summary, by using the PYTHIA8 Angantyr model, we explore the dependence of
the relative neutron density fluctuation, neutron-proton correlation $\alpha$,
and corresponding light nuclei yield ratio $N_{t}N_{p}/N_{d}^{2}$ on rapidity,
centrality, and collision energy. From PYTHIA8 Angantyr model, the relative
nucleon density fluctuation cannot be extracted directly from the light nuclei
ratio, and the effects of neutron-proton correlation $\alpha$ must be
considered in different centralities of collisions. The $N_{t}N_{p}/N_{d}^{2}$
does not change with increasing rapidity coverage and collision centrality,
and it increases slightly with increasing collision energy from the PYTHIA8
Angantyr model. We show the interplay between multi-parton interactions and
color reconnection on the light nuclei yield ratio; there is no effect of CR
if MPI is off. The experimental results reveal that the light nuclei yield
ratio exhibits a peak at $\sqrt{s_{\text{NN}}}=20\sim 30$ GeV, indicating a
significant fluctuation in relative neutron density. The non-monotonic energy
dependence measurements are underestimated by the PYTHIA8 Angantyr model due
to the lack of critical physics and QGP medium mechanisms. Nevertheless, it
can be utilized as a baseline in scenarios where critical physics and QGP
medium mechanisms are absent.

This work is supported in part by the Natural Science Foundation of Henan
Province (No.212300410386), Key Research Projects of Henan Higher Education
Institutions (No.20A140024), the Scientific Research Foundation of Hubei
University of Education for Talent Introduction (No. ESRC20220028 and No.
ESRC20230002), the NSFC 12005114 and NSFC key Grant 12061141008, and key
Laboratory of Quark and Lepton Physics (MOE) in Central China Normal
University (NO. QLPL2022P01).

## References

  * (1) B. J. Schaefer and M. Wagner, “On the QCD phase structure from effective models,” Progress in Particle and Nuclear Physics 62 (2009), 381-385 .
  * (2) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, “The order of the quantum chromodynamics transition predicted by the standard model of particle physics,” Nature 443, (2006) 675 .
  * (3) G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, “The QCD phase diagram at nonzero quark density,” JHEP 1104, (2011) 001.
  * (4) Alexei Bazavov $etal.$, “The equation of state in (2+1)-flavor QCD,” $Phys.~{}Rev.~{}D$ 90 (2014) 094503.
  * (5) E. V. Shuryak, “What RHIC experiments and theory tell us about properties of quark-gluon plasma ?” $Nucl.~{}Phys.~{}A$ 750 (2005), 64-83.
  * (6) F. Karsch, “Deconfinement and quarkonium suppression,” $Eur.~{}Phys.~{}J.~{}C$ 43 (2005), 35-43.
  * (7) Hans J. Specht, “Heavy Ion Physics at the CERN SPS: Roots 1974-1984 and Key Results,” CERN, September 8, 2014.
  * (8) A. Franz, $et.~{}al.$, “RHIC Experimental Evaluations,” Nuclear Physics A Volume 757, (2005) Issues 1-2.
  * (9) S. Gupta, X. Luo, B. Mohanty, H. G. Ritter and N. Xu, “Scale for the Phase Diagram of Quantum Chromodynamics,” Science 332, (2011) 1525 .
  * (10) X. Luo and N. Xu, “Search for the QCD Critical Point with Fluctuations of Conserved Quantities in Relativistic Heavy-Ion Collisions at RHIC : An Overview,” Nucl. Sci. Tech. 28, no. 8, (2017) 112.
  * (11) A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov and N. Xu, “Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan,” Phys. Rep. 853, (2020) pp. 1-87.
  * (12) M. M. Aggarwal et al. [STAR Collaboration], “Higher Moments of Net-proton Multiplicity Distributions at RHIC,” Phys. Rev. Lett. 105, (2010) 022302 .
  * (13) X. Luo, “Exploring the QCD Phase Structure with Beam Energy Scan in Heavy-ion Collisions,” Nucl. Phys. A 956, (2016) 75-82.
  * (14) L. Adamczyk et al. [STAR Collaboration], “Beam energy dependence of moments of the net-charge multiplicity distributions in Au+Au collisions at RHIC,” Phys. Rev. Lett. 113, (2014) 092301.
  * (15) L. Adamczyk et al. [STAR Collaboration], “Energy Dependence of Moments of Net-proton Multiplicity Distributions at RHIC,” Phys. Rev. Lett. 112, (2014) 032302.
  * (16) L. Adamczyk et al. [STAR Collaboration], “Collision Energy Dependence of Moments of Net-Kaon Multiplicity Distributions at RHIC,” Phys. Lett. B 785, (2018) 551.
  * (17) J. Adam et al. [STAR Collaboration], “Net-proton number fluctuations and the Quantum Chromodynamics critical point,” Phys. Rev. Lett. 126, (2021) 092301.
  * (18) M. Stephanov, “QCD phase diagram: an overview,” PoS LAT2006, (2006) 024.
  * (19) M. M. Aggarwal et al. (STAR Collaboration), “Higher Moments of Net Proton Multiplicity Distributions at RHIC,” Phys Rev Lett 105, (2010) 022302.
  * (20) L. Adamczyk et al. (STAR Collaboration), “Energy Dependence of Moments of Net-Proton Multiplicity Distributions at RHIC,” Phys Rev Lett 112, (2014) 032302.
  * (21) M. Stephanov, “Beam Energy Dependence of Moments of the Net-Charge Multiplicity Distributions in Au+Au Collisions at RHIC,” Phys Rev Lett 113, (2014) 092301.
  * (22) K. J. Sun, L. W. Chen, C. M. Ko and Z. Xu, “Probing QCD critical fluctuations from light nuclei production in relativistic heavy-ion collisions,” Phys. Lett. B 774, (2017) 103.
  * (23) X. G. Deng, and Y. G. Ma, “Light nuclei production in Au + Au collisions at $\sqrt{s_{NN}}$ = 7.7-80 GeV from UrQMD model,” Phys. Lett. B 808, (2020) 135668.
  * (24) E. Shuryak and J. M. Torres-Rincon, “Baryon preclustering at the freeze-out of heavy-ion collisions and light-nuclei production,” Phys. Rev. C 101, no.3, (2020) 034914.
  * (25) N. Yu, D. Zhang and X. Luo, “Search for the QCD Critical Point by Transverse Velocity Dependence of Anti-deuteron to Deuteron Ratio,” Chin. Phys. C 44, no. 1, (2020) 014002.
  * (26) T. Shao, J. Chen, C. M. Ko and K. Sun, “Probing QCD critical fluctuations from the yield ratio of strange hadrons in relativistic heavy-ion collisions,” Phys. Lett. B 801, (2020) 135177.
  * (27) K. J. Sun, L. W. Chen, C. M. Ko, J. Pu and Z. B. Xu, “Light nuclei production as a probe of the QCD phase diagram,” Phys. Lett. B 781, (2018) 499.
  * (28) D. Zhang, “Energy Dependence of Light Nuclei (d, t) Production at STAR,” JPS Conf. Proc. 32, (2020) 010069.
  * (29) D. Zhang, “Light Nuclei (d, t) Production in Au + Au Collisions at $\sqrt{s_{NN}}$=7.7$-$200GeV,” Nucl. Phys. A 1005, (2021) 121825.
  * (30) T. Sjostrand, S. Mrenna and P. Z. Skands, “PYTHIA 6.4 Physics and Manual,” JHEP 05 (2006), 026.
  * (31) C. Bierlich, G. Gustafson, L. Lönnblad and H. Shah, “The Angantyr model for Heavy-Ion Collisions in PYTHIA8,” JHEP 10 (2018), 134.
  * (32) C. Bierlich et al, “A comprehensive guide to the physics and usage of PYTHIA 8.3,”  arXiv:2203.11601[hep-ph].
  * (33) K. Abdel-Waged, “Investigating $\pi^{-}$ production in central ${}^{7}Be+^{9}Be$ and ${}^{40}Ar+^{45}Sc$ collisions at CERN SPS energies with PYTHIA8/Angantyr,” Eur. Phys. J. C 82 (2022), 65.
  * (34) I. M. Samsul, S. Tinku, R. Pradip, B. P. Pratim, “Charged-particles distribution in proton-proton and heavy-ion collisions using PYTHIA8 Angantyr model at LHC energies,” Eur. Phys. J. Plus, 137 12 (2022), 1327.
  * (35) H. Liu, D. Zhang, S. He, K. j. Sun, N. Yu, 326 and X. Luo, “Light nuclei production in Au+Au collisions at $\sqrt{s_{NN}}$ = 5$-$200 GeV from JAM model,” Phys. Lett. B 805, (2020) 135452.
  * (36) The STAR Collaboration, “Beam Energy Dependence of Triton Production and Yield Ratio ($\mathrm{N}_{t}\times\mathrm{N}_{p}/\mathrm{N}_{d}^{2}$) in Au+Au Collisions at RHIC,” Phys. Rev. Lett. 130, (2023) 202301.
  * (37) T. Anticic et al. (NA49 Collaboration), “Production of deuterium, tritium, and ${}^{3}{He}$ in central Pb + Pb collisions at $20A,\phantom{\rule{1.60004pt}{0.0pt}}30A,\phantom{\rule{1.60004pt}{0.0pt}}40A,\phantom{\rule{1.60004pt}{0.0pt}}80A$, and $158A$ GeV at the CERN Super Proton Synchrotron,”  Phys. Rev. C 94, (2016) 044906.
  * (38) Z. Zhang, N. Yu, H. G. Xu, “Neutron density fluctuation and neutron-proton correlation from AMPT model,” Eur. Phys. J. A 58, (2022) 240.


# Effect of Dust Evaporation and Thermal Instability on Temperature
Distribution in a Protoplanetary Disk

Ya. N. Pavlyuchenkov, V. V. Akimkin, A. P. Topchieva, E. I. Vorobyov
Institute of Astronomy of the Russian Academy of Sciences, Moscow, 119017
Russia

(Received 16.01.2023; revised 22.03.2023; accepted 27.03.2023
<EMAIL_ADDRESS>
###### Abstract

The thermal instability of accretion disks is widely used to explain the
activity of cataclysmic variables, but its development in protoplanetary disks
has been studied in less detail. We present a semi-analytical stationary model
for calculating the midplane temperature of a gas and dust disk around a young
star. The model takes into account gas and dust opacities, as well as the
evaporation of dust at temperatures above 1000 K. Using this model, we
calculate the midplane temperature distributions of the disk under various
assumptions about the source of opacity and the presence of dust. We show that
when all considered processes are taken into account, the heat balance
equation in the region $r<1$ au has multiple temperature solutions. Thus, the
conditions for thermal instability are met in this region. To illustrate the
possible influence of instability on the accretion state in a protoplanetary
disk, we consider a viscous disk model with $\alpha$-parameterization of
turbulent viscosity. We show that in such a model the disk evolution is non-
stationary, with alternating phases of accumulation of matter in the inner
disk and its rapid accretion onto the star, leading to an episodic accretion
pattern. These results indicate that this instability needs to be taken into
account in evolutionary models of protoplanetary disks.

DOI:10.1134/S1063772923050086

##  1 Introduction

The thermal structure of a protoplanetary disk of gas and dust around a young
star is inextricably linked with many key processes that affect the evolution
of the disk itself. Temperature determines the physical, chemical and
ionization structure of the protoplanetary disk, the location of ice
evaporation fronts and dead zones. The observed properties in the infrared
range depend on the temperature distribution in the disk. The details of
heating and cooling processes determine the development of a number of
instabilities, in particular, gravitational, convective, and thermal
instabilities, see e.g. Armitage (2022); Bae et al. (2022). In this regard,
the problem of describing the self-consistent thermal structure of the disk in
numerical simulations can hardly be overestimated.

The main mechanisms of disk heating are the absorption of stellar and
interstellar radiation, dissipative processes associated with turbulence and
the presence of a magnetic field, and the work of gas pressure forces. The
cooling of the disk is mainly caused by diffusion and release of infrared
radiation. The combination of these processes leads to the formation of a
complex temperature disk structure with vertical stratification and radial
gradients. The main source of absorption in the disk is dust, so it is
important to know the parameters of the dust itself and how they change during
evolution. At relatively small distances from the star, the temperature in the
disk can be high enough for the evaporation of dust. Under such conditions,
gas becomes the main source of absorption of stellar radiation and, at the
same time, the main coolant. The opacity of the gas, in turn, is provided by a
wide variety of absorption and emission processes (vibrational, rotational and
electronic transitions of various atoms and molecules).

The strong dependence of the gas absorption coefficient on temperature, in
particular, due to the ionization of hydrogen, leads to conditions for the
development of thermal instability. Under the conditions of circumstellar
disks, thermal instability appears when several values of the equilibrium
temperature are possible at a fixed surface density, i.e. the disk can be in
thermal and hydrostatic equilibrium in one of two phases, which can be called
a cold and a hot phase. Assuming that temperature affects the turbulent
viscosity of the disk, thermal instability leads to dynamic (viscous)
instability of accretion disks. The thermal instability of accretion disks is
widely used to explain the activity of dwarf novae, X-ray novae, and other
cataclysmic variables, see reviews by Hameury (2020); Lasota (2001).

Study of thermal structure and thermal instability in protoplanetary disks is
important for understanding the nature of FU Ori and EX Lupi type young
stellar objects. These objects are outbursting low-mass protostars
characterized by a sharp increase in luminosity by tens and hundreds of times
(Audard et al., 2014; Connelley & Reipurth, 2018). The nature of these
outbursts is still not clear, but the explanation of this phenomenon is
necessary for any self-consistent theory of star formation. Despite the fact
that the number of known FU Ori and EX Lupi type objects is estimated at
several dozen (Vorobyov et al., 2021), recent observations of young stellar
objects in the optical and near infrared ranges have shown that most young
protostars exhibit radiation variability on time scales from several months to
several years (Flaherty et al., 2016; Contreras Peña et al., 2017; Rigon et
al., 2017). While short luminosity outbursts in EX Lupi type objects can be
explained by the interaction of the inner boundary of the protoplanetary disk
with the protostellar magnetosphere (D’Angelo & Spruit, 2012), longer FU Ori
type outbursts are most likely caused by processes occurring directly in the
protoplanetary disk (Kravtsova et al., 2007). The latter may include both
magneto-rotational instability (Armitage, Livio & Pringle, 2001),
gravitational fragmentation of the disk (Vorobyov & Basu, 2015), and thermal
instability in the inner disk (Bell & Lin, 1994). These instabilities are
essentially sensitive to the thermal structure of the protoplanetary disk,
which makes the study of thermal processes an important astrophysical problem.

Thermal instability in a protoplanetary disk in the context to explain the
phenomenon of FU Ori type objects has been modeled in a number of studies. In
particular, the work by Kawazoe & Mineshige (1993) presented periodic
luminosity curves obtained in a one-dimensional disk model. The key elements
of their model are the simple functions of disk heating and cooling, the
postulated switching of the viscosity coefficient between ionized and neutral
gas, and convective heat transfer in the vertical direction. Bell & Lin (1994)
studied the influence of the model parameters on the characteristics of the
emerging periodic accretion regime in more detail, adopting similar one-
dimensional disk model. Kley & Lin (1999) presented the results of two-
dimensional hydrodynamic simulations of thermal instability and investigated
the transition from the quiet (cold) to active (hot) phase of the evolution of
the inner part of the disk. Two-dimensional hydrodynamic calculations of the
disk evolution with magneto-rotational and gravitational instabilities by
Kawazoe & Mineshige (1993) also indicated the development of thermal
instability in the inner regions of the disk, but the their authors concluded
that thermal instability alone is insufficient to provide the outbursts.

The purpose of this work is to calculate and analyze the distribution of the
midplane temperature of the protoplanetary disk based on a model with a more
detailed description of a number of physical processes. In particular, we use
realistic absorption coefficients by gas and dust, and explicitly include dust
evaporation. In addition, along with viscous heating of the disk, we take into
account heating by stellar radiation, while the heating and cooling functions
are appropriate for arbitrary optical thicknesses. Based on the presented
semi-analytical model, we study the possibility of developing thermal
instability in protoplanetary disks.

##  2 Model description

We consider a stationary axially symmetric Keplerian circumstellar disk with
the following radial distribution of gas surface density

$\Sigma^{\text{gas}}=\Sigma^{\text{gas}}_{\text{au}}\left(\dfrac{R}{R_{\text{au}}}\right)^{-1},$
(1)

where surface density is measured from the midplane, i.e.,
$\Sigma^{\text{gas}}=\int\limits_{0}^{\infty}\rho(R,z)dz$,
$\Sigma^{\text{gas}}_{\text{au}}$ is surface density at $R_{\text{au}}$=1 au,
$R$ is the distance between the disk element and the star. The inner and outer
disk edges are chosen to be 0.1 and 100 au, respectively. The mass and
luminosity of the central star are equal to those of the Sun. In the general
case, two sources of disk heating are taken into account: stellar radiation
and viscous heating. Stellar radiation is assumed to be blackbody with a
temperature of 6000 K. The viscous heating of the disk is calculated in the
stationary accretion approximation with a given accretion rate $\dot{M}$. The
disk is cooled via release of thermal radiation.

###  2.1 Method for calculating midplane temperature

The midplane temperature $T_{\text{mid}}$ at each distance $R$ from the star
is found from the balance between heating and cooling:

$\Lambda_{\text{IR}}=\Gamma_{\text{star}}+\Gamma_{\text{vis}},$ (2)

where $\Lambda_{\text{IR}}$ is the cooling rate of the disk midplane layers
due to escape of IR radiation, $\Gamma_{\text{star}}$ is the heating rate of
the midplane layers by stellar radiation, $\Gamma_{\text{vis}}$ is the heating
rate due to viscous dissipation:

$\displaystyle\Lambda_{\text{IR}}=\frac{4\tau_{\text{P}}\sigma
T_{\text{mid}}^{4}}{1+2\tau_{\text{P}}\left(1+\dfrac{3}{4}\tau_{\text{R}}\right)}$
(3) $\displaystyle\Gamma_{\text{star}}=$ $\displaystyle=\frac{\mu
F_{0}\tau_{\text{P}}\left[2\left(1-e^{-\tau_{\text{uv}}}\right)+3\mu\dfrac{\tau_{\text{R}}}{\tau_{\text{uv}}}\left(1-e^{-\tau_{\text{uv}}}-\tau_{\text{uv}}e^{-\tau_{\text{uv}}}\right)+\frac{\tau_{\text{uv}}}{\tau_{\text{P}}}e^{-\tau_{\text{uv}}}\right]}{1+2\tau_{\text{P}}\left(1+\dfrac{3}{4}\tau_{\text{R}}\right)}$
(4) $\displaystyle\Gamma_{\text{vis}}=\frac{3}{8\pi}\frac{GM\dot{M}}{R^{3}}.$
(5)

The detailed derivation of the above functions is presented in the Appendix.
In these equations, $F_{0}$ is the radiation flux from the star reaching the
disk surface, $\mu$ is the cosine of the angle between the direction to the
star and the normal to the disk surface (chosen to be 0.05),
$\tau_{\text{uv}}$ is optical depth to stellar radiation, $\tau_{\text{P}}$
and $\tau_{\text{R}}$ are Planck and Rosseland mean optical depths to disk
thermal radiation, $\sigma$ is the Stefan-Boltzmann constant, $G$ is the
gravitational constant, $M$ is the mass of the central star, $\dot{M}$ is the
accretion rate over the disk. Optical depths are defined as follows:

$\displaystyle\tau_{\text{uv}}=\dfrac{1}{\mu}(\kappa_{\text{F}}^{\text{gas}}\Sigma^{\text{gas}}+\kappa_{\text{F}}^{\text{dust}}\Sigma^{\text{dust}})$
(6)
$\displaystyle\tau_{\text{P}}=\kappa_{\text{P}}^{\text{gas}}\Sigma^{\text{gas}}+\kappa_{\text{P}}^{\text{dust}}\Sigma^{\text{dust}}$
(7)
$\displaystyle\tau_{\text{R}}=\kappa_{\text{R}}^{\text{gas}}\Sigma^{\text{gas}}+\kappa_{\text{R}}^{\text{dust}}\Sigma^{\text{dust}},$
(8)

where $\Sigma^{\text{gas}}$ and $\Sigma^{\text{dust}}$ are gas and dust
surface densities, $\kappa_{\text{F}}^{\text{gas}}$,
$\kappa_{\text{P}}^{\text{gas}}$, and $\kappa_{\text{R}}^{\text{gas}}$ are
two-temperature mean, Planck mean, and Rosseland mean opacity coefficients for
gas, respectively, $\kappa_{\text{F}}^{\text{dust}}$,
$\kappa_{\text{P}}^{\text{dust}}$, and $\kappa_{\text{R}}^{\text{dust}}$ are
dust opacity coefficients. Opacity coefficients are generally functions of
density and temperature. Note that the expression (8) is valid in the gray
body approximation, but is incorrect in the general case, i.e. Rosseland mean
optical depths are not additive. This is due to the non-linear nature of the
averaging procedure. However, we use this formula because it provides
acceptable accuracy and significantly simplifies the subsequent analysis.

Equations (3)–(4) are derived under the assumption that the opacity
coefficients are constant along the vertical direction. This is obviously a
rather rough approximation, since both temperature and density change
significantly with height (distance from the midplane). Nevertheless, we use
these expressions because they allow to describe the thermal structure of the
disk in the first-order approximation. Temperature and density in the midplane
are used as arguments for calculating the opacity coefficients. The gas
opacity coefficients, as will be shown below, depend on gas density. At a
fixed midplane temperature, the midplane density is found from the condition
of hydrostatic equilibrium of a vertically isothermal disk:

$\rho_{\text{mid}}=\frac{1}{\sqrt{2\pi}}\frac{2\Sigma^{\text{gas}}}{H}$ (9)

where $H=c_{s}/\Omega$ is the disc scale height,
$c_{s}=\sqrt{k_{\text{B}}T_{\text{mid}}/m}$ is isothermal sound speed, $m$ is
molecular mass, $k_{\text{B}}$ is Boltzmann constant, $\Omega=\sqrt{GM/R^{3}}$
is Keplerian angular velocity at radius $R$.

Since the opacity coefficients are functions of temperature, the equation (2)
is a non-linear equation with respect to $T_{\text{mid}}$. We find its
solution graphically by calculating the sign of the difference between the
left and right sides of equation (2) for a sequence of trial values of
$T_{\text{mid}}$. The areas of sign change of this difference correspond to
the roots of the equation.

###  2.2 Gas absorption coefficients

Gas opacity coefficients $\kappa_{\text{P}}^{\text{gas}}$,
$\kappa_{\text{R}}^{\text{gas}}$, $\kappa_{\text{F}}^{\text{gas}}$ as
functions of temperature and gas density (see Fig. 1) are adopted from Malygin
et al. (2014)111https://vizier.cds.unistra.fr/viz-
bin/VizieR?-source=J/A+A/568/A91. To calculate these coefficients, Malygin et
al. (2014) used the DFSYNTHE code (Castelli, 2005; Kurucz, 1970), originally
oriented for calculation of stellar atmospheres. The coefficients we use
correspond to the solar metallicity, $\kappa_{\text{F}}$ is taken for a
stellar temperature of 6000 K. We note that at relatively low temperatures
($T<1000$ K) the main contribution to $\kappa_{\text{P}}^{\text{gas}}$ and
$\kappa_{\text{R}}^{\text{gas}}$ comes from the absorption lines of various
molecules. At higher temperatures ($T>3000$ K), absorption by atomic hydrogen,
the H- ion, various metals, and other processes become important. Note also
that the distribution of $\kappa_{\text{F}}^{\text{gas}}(T,\rho)$ differs
significantly from the distributions of
$\kappa_{\text{P}}^{\text{gas}}(T,\rho)$ and
$\kappa_{\text{R}}^{\text{gas}}(T,\rho)$.

When calculating the thermal structure of the disk for arbitrary
$T_{\text{mid}}$ and $\rho_{\text{mid}}$, we use linear interpolation of the
coefficients between neighboring nodes of the input grid and extrapolation of
the opacity coefficients by edge values in case of going beyond the original
grid in terms of temperature and density.

Figure 1: Planck mean (top left panel), Rosseland mean (top right panel), and
two-temperature mean (bottom panel) gas opacity coefficients from Malygin et
al. (2014). The white curve in the top right panel corresponds to the
equilibrium parameters of the disk in the instability region from the M5
model, see Fig. 3.

###  2.3 Dust absorption coefficients and evaporation model

The adopted Planck and Rosseland mean opacitis as functions of temrerature are
shown in the left panel of Fig. 2. We calculate these coefficients from the
frequency-dependent absorption and scattering coefficients for spherical
silicate dust grains. We chose forsterite Mg2SiO4 as a model of silicate dust.
The spectral absorption and scattering coefficients themselves were calculated
using the Mie theory, with the size distribution of dust grains taken as a
power-law $n(a)\propto a^{-3.5}$ with the minimum and maximum dust grain radii
$a_{\rm min}=5\times 10^{-7}$ cm and $a_{\rm max}=10^{-4}$ cm. Note that the
opacity is shown in Fig. 2 up to $10^{4}$ K, although the dust should
evaporate at much lower temperatures (see below). However, the calculation of
the dust opacity coefficients (per unit mass of dust) and the fraction of
evaporated dust are independent procedures, the combination of which will give
the optical depth necessary for modeling the thermal structure based on the
relations (6)–(8).

Figure 2: Left panel: Planck and Rosseland mean dust opacity versus
temperature. Right panel: fraction of unevaporated dust as a function of
midplane temperature for three values of the disk surface density at 1 au .

At high temperatures the dust evaporates. We account for this process by
calculating the surface density of the dust as follows:

$\Sigma^{\text{dust}}=\xi(T_{\text{mid}})\,\mu_{\text{dg}}\,\Sigma^{\text{gas}},$
(10)

where $\mu_{\text{dg}}=0.01$ is the dust-to-gas mass ratio we assume in the
absence of evaporation, $\xi(T)$ is the fraction of unevaporated dust at
temperature $T$. Thermal evaporation of forsterite is a complex process, it is
accompanied by the formation of a number of components in the gas phase, such
as Si, O, Mg, MgO, O2, SiO, SiO2. We calculate $\xi(T)$ assuming that dust is
in thermodynamic equilibrium between the solid and the gas phases. In this
approximation, the coexistence of phases is possible provided that the partial
pressures of gases from the dust components are equal to the pressures of
their saturated vapors at a given temperature. The calculation of equilibrium
partial pressures is carried out by the methods of chemical thermodynamics and
is presented, for example, in the paper by Duschl, Gail & Tscharnuter (1996).

We consider an acceptable simplification that all silicon in the gaseous phase
is in the form of a SiO molecule. To calculate the function $\xi(T)$ we use
the following expression:

$\xi(T)=1-\dfrac{n_{\text{vap}}}{n_{\text{tot}}},$ (11)

where $n_{\text{vap}}$ is number density of silicon (in the form of SiO) in
saturated vapor at temperature $T$, $n_{\text{tot}}$ is total number density
of silicon (in solid and gaseous phases). Let us express the number density
$n_{\text{vap}}$ in terms of saturated vapor pressure $P_{\text{vap}}$:

$n_{\text{vap}}=\frac{P_{\text{vap}}}{k_{\text{b}}T}.$ (12)

We take the approximation of $P_{\text{vap}}$ temperature dependence for
silicon oxide vapor from Duschl, Gail & Tscharnuter (1996) as:

$P_{\text{vap}}=f_{\text{a}}\,\exp\left(x_{1}/T+x_{2}+x_{3}T+x_{4}T^{2}+x_{5}T^{3}\right),$
(13)

with the parameters $x_{1}=-6.28\times 10^{4}$ K, $x_{2}=1.80\times 10^{1}$,
$x_{3}=3.59\times 10^{-4}$ K-1, $x_{4}=-3.72\times 10^{-7}$ K-2,
$x_{5}=6.53\times 10^{-11}$ K-3, and the coefficient $f_{\text{a}}=10^{6}$ dyn
cm-2 which takes into account the conversion of pressure to CGS units.

The total number density of silicon $n_{\text{tot}}$ is considered
proportional to the number density of the matter in the midplane:

$n_{\text{tot}}=\frac{X_{\text{Si}}\,\rho_{\text{mid}}}{\mu_{\text{Si}}\,m_{\text{a}}},$
(14)

where $X_{\text{Si}}=3.5\times 10^{-5}$ is mass fraction of silicon in the
interstellar medium, $\mu_{\text{Si}}=28$ is atomic weight of silicon,
$m_{\text{a}}$ is atomic unit of mass. Since $n_{\text{tot}}$ is determined by
the midplane density $\rho_{\text{mid}}$, which depends on a number of
parameters (see equation (9)), then the dependence $\xi(T)$ will be in general
unique for each $R$ in the disk.

At a sufficiently high temperature, $n_{\text{vap}}$ may turn out to be
greater than $n_{\text{tot}}$, and $\xi$ becomes formally negative according
to the formula (11). This means that at a given temperature, the silicon
available in a given volume is insufficient to saturate the vapors, and phase
equilibrium is impossible. In this case, all silicon should pass into the
gaseous state, i.e. the dust is completely evaporated. To avoid physical
uncertainties in the numerical model, we additionally restrict the smallest
value of $\xi$ by the parameter $\xi_{\text{min}}=10^{-5}$, i.e. we assume
that under the condition $n_{\text{vap}}>n_{\text{tot}}$ a small fraction of
dust grains does not evaporate.

The right panel of Fig. 2 shows the function $\xi(T_{\text{mid}})$ for three
values of the surface density of the disk at 1 au ($\Sigma=$1, $10^{2}$,
$10^{4}$ g cm-2) around a solar mass star. It can be seen that at low
temperatures, the function $\xi(T)$ is close to unity, and when the region of
intense evaporation is reached, the function $\xi(T)$ decreases very rapidly.

The adopted dust model is quite primitive and is used only to study the
physical effect of evaporation. In fact, protoplanetary dust can include a
carbon component, refractory organic matter, polycyclic aromatic hydrocarbons,
ice mantles. The dust particles themselves can have a complex and fractal
shape.

###  2.4 Considered models

The “M1” model considers a sparse low-mass ($M_{\text{disk}}$ $\approx$
$10^{-5}M_{\odot}$) disk heated only by stellar irradiation, and only takes
into account dust opacity, while the the dust is assumed to be non-evaporable.
In the “M2” model, the surface density of the disk is increased by four orders
of magnitude relative to the “M1” model, such a massive disk
($M_{\text{disk}}$ $\approx$ 0.1 $M_{\odot}$) corresponds to the initial
phases of the evolution of a protostellar system. In the “M3” model, viscous
heating corresponding to the accretion rate of $10^{-6}M_{\odot}$/yr is added
to “M2” as an additional source of heating. In the “M4” model, the restriction
on the impossibility of dust evaporation is removed. Finally, the “M5” model
takes into account, in addition to everything else, the contribution of the
gas to absorption and emission. Designations and parameters of considered
models are given in Table 1.

Model | Dust | Dust | Gas |  $\Sigma_{0}$, g/cm2 |  $\dot{M}$, $M_{\odot}$/yr
---|---|---|---|---|---
| Opacity | Evaporation | Opacity |  |
M1 | $+$ | $-$ | $-$ | $0.1$ | 0
M2 | $+$ | $-$ | $-$ | $1000$ | 0
M3 | $+$ | $-$ | $-$ | $1000$ | $10^{-6}$
M4 | $+$ | $+$ | $-$ | $1000$ | $10^{-6}$
M5 | $+$ | $+$ | $+$ | $1000$ | $10^{-6}$
Table 1: Considered disk models and their parameters

##  3 Simulation results

Fig. 3 shows the results of calculating the temperature distributions for the
considered models. The equilibrium temperature lays at the border between the
red and blue regions, the color is determined by the sign of the difference
between the right and left parts of the equation (2).

Figure 3: Calculated disk thermal structure for models M1–M5. The equilibrium
temperature lays at the border between the red and blue regions, the color is
determined by the sign of the difference between the right and left parts of
the equation (2). The dashed line shows the distribution according to equation
(15). The white curve highlights the multi-valued solution that is responsible
for the instability.

In the M1 model, the temperature gradually decreases with distance up to a
radius of 7 AU, after which a “hump” forms on the distribution. The monotonic
inner portion of the distribution corresponds to the region of the disk opaque
to stellar radiation and is well described by the relation (45):

$aT^{4}_{\text{mid}}=\dfrac{2\mu F_{0}}{c},$ (15)

which depends only on the parameters of the stellar radiation flux. Dependence
(15) is shown in the figure by a dashed line. The disk in the M1 model is so
rarefied that the outer parts of this disk ($R>7$ AU) are transparent to
stellar radiation. The formation of a temperature hump is associated with a
smooth transition to the limit (44) for a disk that is optically thin to
stellar radiation:

$aT^{4}_{\text{mid}}=\frac{\kappa_{\text{F}}}{\kappa_{\text{P}}}\frac{F_{0}}{c},$
(16)

where the thermal structure depends not only on the radiation flux, but also
on the opacity ratio $\kappa_{\text{F}}/\kappa_{\text{P}}$.

The optical depth to stellar radiation in the M2 massive disk model is large
for the entire disk, so the distribution of the midplane temperature is
monotonous and is determined only by the dilution of stellar radiation in
accordance with the formula (15).

In the M3 model of a massive disk with viscous heating, the temperature
distribution depends not only on the stellar flux parameters, but also on the
dissipation rate and the opacity of the medium. The thermal structure of such
a disk is described by the expression (46):

$aT^{4}_{\text{mid}}=\dfrac{2\mu
F_{0}}{c}+\frac{\Gamma_{\text{vis}}}{c}\left[\frac{1+2\tau_{\text{P}}(1+\frac{3}{4}\tau_{\text{R}})}{\tau_{\text{P}}}\right].$
(17)

With the adopted model parameters, the temperature changes from $4\times
10^{4}$ K to 20 K when moving from the inner to the outer boundary of the
disk. The uneven distribution of the equilibrium temperature is related to the
features of the $\kappa_{\text{P}}^{\text{dust}}(T)$ and
$\kappa_{\text{R}}^{\text{dust}}(T)$ dependences. In particular, these
features are associated with the nonmonotonic behavior of the absorption
coefficient of silicate dust in the vicinity of 10 microns. Note that the
temperature inside 2 au turns out to be much higher than the dust evaporation
temperature (which does not exceed 2000 K, see Fig. 2). Therefore, the M3
model is certainly physically inconsistent, but we present it for
methodological purposes to show what the considered approximations lead to
when calculating the thermal structure of a disk.

In the M4 model, which takes into account both viscous heating and dust
evaporation, the temperature in the inner region of the disk ($R<2$ au) is
significantly lower than in the “M3” model. The temperature in this region is
$T_{\text{mid}}\approx$2000 K. When such temperatures are reached in the
specified region, the abundance of unevaporated dust becomes equal to the
minimum value postulated in the model $\xi_{\text{min}}=10^{-5}$. The value of
$10^{-5}$ was chosen in a way that the corresponding optical depths
$\tau_{\text{P}}$ and $\tau_{\text{R}}$ in this region are close to unity,
thereby ensuring the maximum cooling rate $\Lambda_{\text{IR}}$. In this case,
the cooling rate $\Lambda_{\text{IR}}$ fully compensates for the heating rate
$\Gamma_{\text{star}}$ and $\Gamma_{\text{vis}}$, providing the equilibrium
value of $T_{\text{mid}}$. A plateau with a slightly lower temperature
$T_{\text{mid}}\approx 1600$ K is also described in the paper by D’Alessio et
al. (1999) (see their Fig. 1 and description on page . 895, third paragraph of
section 2.2). D’Alessio et al. (1999) also take dust sublimation into account,
while their gas opacity is small for these conditions, which makes their model
conceptually close to our M4 model. Differences in plateau temperatures are
probably due to different models of dust evaporation.

In the M5 model, which takes into account the contribution of the gas to the
opacity of the medium, the optical depths in the inner part of the disk
increase significantly, which generally leads to an increase in temperature
compared to the M4 model. The key feature of this model is that the
equilibrium temperature within 0.15–0.3 au has multiple solutions, which is
caused by the strong dependence of gas opacity on temperature. The red color
on the distributions corresponds to heating dominating over cooling, while the
blue color corresponds to the predominance of cooling. With this in mind, the
arrival at equilibrium can be considered in this diagram as an upward movement
(heating leads to an increase in temperature) in the red region and a downward
movement (cooling leads to a decrease in temperature) in the blue region to
the boundary of the regions. Thus within 0.15–0.3 AU stable equilibrium is
possible at temperatures of 2–3 thousand K or 10–20 thousand K. Taking into
account the ambiguity of the solution for the equilibrium temperature, it can
be expected that in the inner zone of the considered disk, the conditions for
thermal instability are fulfilled. The actual temperature will be determined
by the side of the distribution from which the disk comes to the equilibrium.

It is also useful to represent the results obtained for the M5 model as the
function of surface density $T_{\text{mid}}(\Sigma)$ shown in the lower right
panel of Fig. 3. The resulting S-shaped distribution of the equilibrium
temperature $T(\Sigma)$ inevitably leads to an analogy with the shape of the
$\Sigma-T_{\text{eff}}$ dependence for classical thermal instability in
accretion disks of cataclysmic variables (Lasota, 2001), where the formation
of an equilibrium temperature bend is associated with hydrogen ionization,
which leads to a strong dependence of opacity on temperature. In Fig. 1 the
white curve shows the values of midplane density and temperature corresponding
to the region of thermal instability. One can see that the white curve crosses
the area of strong gradients in the opacity distribution. This result allows
us to say that when the opacity of gas in the inner regions of protoplanetary
disks is taken into account, the development of instability is possible. This,
among other things, can lead to the formation of morphological features in the
inner regions of the disks and/or a nonstationary (episodic) character of
accretion. The episodic accretion will naturally appear in the model of
viscous disk evolution if the $\alpha$-parameterization of turbulent viscosity
is used. With such a parameterization, there is a positive feedback between
the disk temperature and the accretion rate, which, in the presence of a jump
in the equilibrium temperature, leads to an accumulative mode of disk
evolution, followed by a rapid accretion of disk matter onto the star. An
example of such a model is discussed in the next section.

##  4 Evolutionary disk model

Figure 4: Results of the calculation of the viscous disk evolution. Left
panel: surface density distributions for three moments. Middle panel: midplane
temperature distributions for the same time moments. Right panel: disk-to-
stellar accretion rate change over time. The zero moment of time is chosen
arbitrarily and corresponds to 15.4 thousand years of evolution from the
initial state.

To illustrate possible influence of the discussed instability on the nature of
accretion in a protoplanetary disk, we consider a model of a viscous disk, the
evolution of which is described using the following system of equations:

$\displaystyle\frac{\partial\Sigma}{\partial
t}=\frac{3}{R}\frac{\partial}{\partial R}\left(R^{1/2}\frac{\partial}{\partial
R}\left(R^{1/2}\nu\Sigma\right)\right)$ (18)
$\displaystyle\frac{\partial\varepsilon}{\partial
t}=\Gamma_{\text{star}}+\Gamma_{\text{vis}}-\Lambda_{\text{IR}}+\frac{3}{R}\frac{\partial}{\partial
R}\left(R^{1/2}\frac{\partial}{\partial
R}\left(R^{1/2}\nu\varepsilon\right)\right),$ (19)

where $\varepsilon=c_{\text{v}}T_{\text{mid}}\Sigma$ is thermal energy per
unit disk surface, $c_{\text{v}}=\dfrac{k_{\text{B}}}{m(\gamma-1)}$,
$\gamma=7/5$ is adiabatic exponent, $\Gamma_{\text{star}}$,
$\Gamma_{\text{vis}}$ are the rates of disk heating by stellar radiation and
viscous dissipation, calculated by the equations (4) and (5), respectively,
$\Lambda_{\text{IR}}$ is cooling rate due to IR emission (3), $\nu$ is
turbulent viscosity coefficient. The equation (18) for surface density is the
classical Pringle equation.

Equation (19) for the evolution of thermal energy along with the rates of
heating and cooling includes radial transfer of thermal energy (the last term
on the right side). The transfer of thermal energy is treated here similarly
to the transfer of individual components in the diffusion accretion-decretion
disk approximation. In this approximation, it is assumed that turbulence leads
to efficient mixing of matter, i.e. to turbulent diffusion. As a result, the
dynamics of all individual components of matter (for example, impurities) is
described by the same diffusion equation, similar to the Pringle equation, see
equation (13) in Pavlyuchenkov & Dullemond (2007). Assuming that thermal
energy is inextricably linked with the matter itself, the equation for its
transfer in this approximation is similar to the Pringle equation. At zero
rates of heating and cooling, such a system of equations leads to the fact
that the initially isothermal disk will remain so in the course of evolution.
We note that for the case of a purely viscous disk, the spatial dynamics of
individual components (and thermal energy) must be described using the
advection equation. The relationship between the viscosity and diffusion
coefficients (Schmidt number) in protoplanetary disks is a debatable issue,
see e.g. Pavlyuchenkov & Dullemond (2007), and therefore the choice in favor
of one or another approximation is not obvious. Note, however, that the use of
the diffusion approximation to calculate the thermal energy transfer greatly
simplifies the numerical model and, at the same time, reflects the main
features of the spatial redistribution of energy.

The viscosity coefficient is given in the framework of the classical $\alpha$
parameterization of Shakura & Sunyaev (1973): $\nu=\alpha c_{\text{s}}H$,
where $c_{\text{s}}$ is speed of sound in the midplane, $H$ is typical disk
scale height. The scale height is calculated from the vertical hydrostatic
equilibrium condition: $H/R=c_{\text{s}}/v_{\text{k}}$. We use the constant
value of $\alpha=10^{-3}$. The surface density distribution from the M5 model
is taken as the initial one.

The integration of equations (18)– (19) at each time step is carried out in
two stages. At the first stage, we solve the problem of diffusion of surface
density and thermal energy. At the second stage, the change in thermal energy
due to heating and cooling sources is calculated. For both stages, implicit
methods are used (the tridiagonal matrix algorithm for solving a system of
linear algebraic equations when calculating diffusion and the bisection method
for updating the temperature due to the heating and cooling, respectively),
which allows one to significantly weaken the constraint on the time step.

Fig. 4 shows the calculation results of the considered disk evolution model
after 15.4 thousand years from the initial distribution. In the inner region
of the disk ($R<0.5$ au), the disk structure is periodically rearranged: when
a sufficient amount of matter is accumulated, the IR radiation is “locked”,
which leads to an increase in temperature and a transition to the upper branch
of the quasi-equilibrium solution. High temperature, in turn, leads to
increased viscous heating, which maintains a high-energy regime until the
moment when a significant part of the matter from the inner region accretes
onto the star as a result of increased viscosity. This regime leads to an
episodic accretion pattern with a period of about 1000 years, which, however,
will change as the disk is depleted. We note that the characteristics of flare
activity in this model strongly depends on the viscosity coefficient. The
maximum accretion rate decreases with $\alpha$, while the flare period
increases, and the flare shape also changes.

##  5 Conclusions

The model of the thermal structure of a protoplanetary disk presented in this
article is mainly illustrative due to a number of rather rough approximations.
Its main goal was to demonstrate that conditions for thermal instability can
indeed be fulfilled in a protoplanetary disk. Undoubtedly, a full study of the
evolution of the protoplanetary disk, taking into account dust evaporation and
using realistic gas extinction coefficients, should be carried out on the
basis of a more consistent dynamical model. In such a model it is particularly
necessary to abandon the approximation of thermodynamic equilibrium between
the solid and gaseous phases when calculating the fraction of evaporated dust.
Instead, it should be taken into account that the dust evaporates and
condenses in a finite time, and the characteristic times of these processes
can differ significantly. The dynamic model must take into account not only
the motion of the gas, but also the drift, settling, and growth of dust
grains, since these processes have a strong influence on the thermal and
spatial structure of the disk. Nevertheless, the expressions for the rates of
heating and cooling presented in this paper, together with the formalism of
calculating optical depths taking into account gas opacity and dust
evaporation, can be used for more detailed modeling of the long-term evolution
of the disk, for example, within the model of a viscous self-gravitating disk
from the studies by Vorobyov & Basu (2015); Vorobyov et al. (2021).

One of the key issues related to thermal instability in gas and dust disks, in
our opinion, is the question of its actual effect on disk morphology and the
nature of accretion onto the star. The periodic nature of accretion,
illustrated in the last section, arises when using the
$\alpha$-parameterization of turbulent viscosity, which provides a positive
feedback between the accretion rate and temperature. Meanwhile, the source of
viscosity in protoplanetary disks has not yet been reliably established. Note
that when using the $\beta$-parameterization of turbulent viscosity, where
there is no dependence on temperature, no periodic accretion will occur within
the viscous disk model,see Huré, Richard & Zahn (2001). Thus, the effect of
thermal instability should be studied together with the question of the source
of turbulent viscosity in the disk.

## Funding

The study was supported by the Russian Science Foundation grant No.
22-72-10029, https://rscf.ru/project/22-72-10029/.

## ACKNOWLEDGMENTS

The authors are grateful to the referee for valuable comments and suggestions
for improving the article. We also express our gratitude to Yu.A. Fadeev, D.A.
Semenov, L.I. Mashonkina, T.M. Sitnova for discussing the problem of
calculating gas opacity.

##  Appendix A Derivation of the formula for the midplane temperature of the
circumstellar disk

We consider a circumstellar disk in a state of thermal equilibrium. In the
plane-parallel approximation, the thermal structure of such a disk in the
vertical direction can be described by a system of moment transfer equations
for thermal radiation:

$\displaystyle\frac{dF}{dz}=c\rho\kappa_{\text{P}}(B-E)$ (20)
$\displaystyle\frac{c}{3}\frac{dE}{dz}=-\rho\kappa_{\text{R}}F,$ (21)

where $F$ is radiation flux, $E$ is density of radiation energy, $B=aT^{4}$ is
radiation energy density at thermodynamic equilibrium, $a$ is radiation
density constant, $T$ is medium temperature, $c$ is speed of light, $z$ is
vertical coordinate measured from the midplane, $\rho$ is medium density,
$\kappa_{\text{P}}$ and $\kappa_{\text{R}}$ are Planck and Rosseland mean
opacity coefficients. The equation (20) describes the change in flux due to
the difference between the emission and absorption of radiation energy. The
equation (21) relates the radiation flux to the energy density in the
Eddington approximation. The system of equations (20)–(21) is closed by the
equation:

$\frac{dF}{dz}=\rho S,$ (22)

according to which the thermal radiation flux is generated by some heating
source $\rho S$, where $S$ [erg s-1 r${}^{-1}]$ is defined as the heating
power per unit mass. Let us rewrite these equations using the surface density
$\Sigma=\int\limits_{0}^{z}\rho(z^{\prime})dz^{\prime}$ as a variable:

$\displaystyle c\kappa_{\text{P}}(B-E)=S$ (23)
$\displaystyle\frac{c}{3}\frac{dE}{d\Sigma}=-\kappa_{\text{R}}F$ (24)
$\displaystyle\frac{dF}{d\Sigma}=S.$ (25)

We will assume that the heating power $S$ originates from two processes:
absorption of stellar radiation and viscous dissipation:

$S=S_{\text{star}}+S_{\text{vis}}.$ (26)

Heating power due to viscous dissipation of gas in the stationary
approximation can be found as:

$S_{\text{vis}}=\frac{\Gamma_{\text{vis}}}{\Sigma_{0}}=\frac{3}{8\pi}\frac{GM\dot{M}}{R^{3}}\bigg{/}\Sigma_{0},$
(27)

where $M$ is the stellar mass, $\dot{M}$ is the accretion rate through the
disk, $R$ is the distance from the star to the disk element in question,
$\Sigma_{0}$ is surface density from the midplane to the upper boundary of the
disk, $G$ is gravitational constant. The use of (27) is also based on the
assumption that the rate of viscous dissipation per unit volume is
proportional to the density of the medium. The heating of the disk by stellar
radiation is found using the formula:

$S_{\text{star}}=\kappa_{\text{F}}F_{0}\,e^{-\dfrac{\kappa_{\text{F}}(\Sigma_{0}-\Sigma)}{\mu}},$
(28)

where $\kappa_{\text{F}}$ is the extinction coefficient averaged over the
stellar spectrum, $F_{0}=\dfrac{L}{4\pi R^{2}}$ is the radiation flux from the
star reaching the surface of the disk, $L$ is stellar luminosity, $\mu$ is
cosine of the angle between the direction to the star and the normal to the
disk surface. The formula (28) is derived from a formal solution of the
radiative transfer equation under the assumption that the absorption
coefficient $\kappa_{\text{F}}$ is constant along the vertical direction. In
this case, we also neglect stellar radiation from the opposite surface of the
disk. Accounting for disk heating by stellar radiation by introducing the
$S_{\text{star}}$ function into the (23)–(25) system is based on the
assumption that the disk radiates weakly in the visible range, i.e. this range
weakly intersects with the range of thermal radiation of the disk itself. We
introduce the notation

$\tau_{\text{uv}}=\dfrac{\kappa_{\text{F}}\Sigma_{0}}{\mu},$ (29)

which is the optical depth of the medium to stellar radiation up to the
current position in the disk. Integration of the equation (25), taking into
account the expressions (26)–(28) and the condition that thermal radiation
flux in the midplane is zero (due to the symmetry of the problem), gives:

$F=\mu
F_{0}\,e^{-\tau_{\text{uv}}}\left(e^{\dfrac{\tau_{\text{uv}}\Sigma}{\Sigma_{0}}}-1\right)+S_{\text{vis}}\Sigma.$
(30)

In particular, on the disk surface, the thermal radiation flux is equal to:

$F(\Sigma_{0})=\mu
F_{0}\,\left(1-e^{\tau_{\text{uv}}}\right)+S_{\text{vis}}\Sigma_{0}.$ (31)

By substituting the equation (30) into the equation (24) and integrating the
resulting equation from the midplane to the upper boundary of the disk, we can
obtain a relation between the radiation energy density on the surface
$E(\Sigma_{0})$ and in the disk midplane $E(0)$:

$E(\Sigma_{0})-E(0)=-\dfrac{3\kappa_{\text{R}}\mu^{2}F_{0}}{c\kappa_{\text{F}}}\left(1-\tau_{\text{uv}}e^{-\tau_{\text{uv}}}-\tau_{\text{uv}}\right)-\dfrac{3\kappa_{\text{R}}S_{\text{vis}}}{2c}\Sigma_{0}^{2}.$
(32)

When obtaining the relation (32), it was assumed that $\kappa_{\text{R}}$ is
constant along the vertical direction. As a boundary condition on the disk
surface, we can use the relation:

$F(\Sigma_{0})=\eta\,cE(\Sigma_{0}),$ (33)

where the coefficient $\eta$ depends on the assumed anisotropy of the outgoing
thermal radiation. The value $\eta=0.5$ corresponds to isotropy over the upper
hemisphere, while $\eta=1$ describes the case of strictly vertical radiation
output. In what follows, we will set $\eta=0.5$. Combining the equations (31),
(32) and (33), one can obtain an expression for the radiant energy of thermal
radiation in the midplane:

$E(0)=\frac{F_{0}}{c}\left[2\mu(1-e^{-\tau_{\text{uv}}})+3\mu^{2}\frac{\kappa_{\text{R}}}{\kappa_{\text{F}}}(1-e^{-\tau_{\text{uv}}}-\tau_{\text{uv}}e^{-\tau_{\text{uv}}})\right]+\\\
+\frac{2S_{\text{vis}}\Sigma_{0}}{c}\left(1+\frac{3}{4}\tau_{\text{R}}\right),$
(34)

where the Rosseland optical depth to thermal radiation is introduced:

$\tau_{\text{R}}=\kappa_{\text{R}}\Sigma_{0}.$ (35)

The desired midplane temperature $T_{\text{mid}}$ is found from:

$B(0)=aT_{\text{mid}}^{4}.$ (36)

The value of $B(0)$, in turn, is expressed through the equation (23), which,
taking into account the value of the source function in the midplane
(26)–(28), takes the form:

$B(0)=E(0)+\frac{\kappa_{\text{F}}F_{0}\,e^{-\tau_{\text{uv}}}}{c\kappa_{\text{P}}}+\frac{S_{\text{vis}}}{c\kappa_{\text{P}}}.$
(37)

Combining the equations (34), (36) and (37), we get:

$aT_{\text{mid}}^{4}=\dfrac{\mu
F_{0}}{c}\left[2\left(1-e^{-\tau_{\text{uv}}}\right)+3\mu\dfrac{\tau_{\text{R}}}{\tau_{\text{uv}}}\left(1-e^{-\tau_{\text{uv}}}-\tau_{\text{uv}}e^{-\tau_{\text{uv}}}\right)\right]+\\\
+\dfrac{\mu
F_{0}\tau_{\text{uv}}}{\tau_{\text{P}}c}e^{-\tau_{\text{uv}}}+\dfrac{\Gamma_{\text{vis}}}{c}\left[\dfrac{1+2\tau_{\text{P}}\left(1+\dfrac{3}{4}\tau_{\text{R}}\right)}{\tau_{\text{P}}}\right],$
(38)

where the Planck mean optical depth with respect to thermal radiation is
introduced:

$\tau_{\text{P}}=\kappa_{\text{P}}\Sigma_{0}.$ (39)

The (38) equation can also be usefully expressed as follows:

$\Lambda_{\text{IR}}=\Gamma_{\text{star}}+\Gamma_{\text{vis}},$ (40)

where $\Lambda_{\text{IR}}$ [erg c-1 cm-2] is cooling rate of midplane disk
layers, $\Gamma_{\text{star}}$ [erg c-1 cm-2] is heating rate by stellar
radiation for the midplane layers, $\Gamma_{\text{vis}}$ [erg c-1 cm-2] is
heating rate due to viscous dissipation:

$\displaystyle\Lambda_{\text{IR}}=\frac{4\tau_{\text{P}}\sigma
T_{\text{mid}}^{4}}{1+2\tau_{\text{P}}\left(1+\dfrac{3}{4}\tau_{\text{R}}\right)}$
(41) $\displaystyle\Gamma_{\text{star}}=\hskip 207.7052pt$
$\displaystyle=\frac{\mu
F_{0}\tau_{\text{P}}\left[2\left(1-e^{-\tau_{\text{uv}}}\right)+3\mu\dfrac{\tau_{\text{R}}}{\tau_{\text{uv}}}\left(1-e^{-\tau_{\text{uv}}}-\tau_{\text{uv}}e^{-\tau_{\text{uv}}}\right)+\frac{\tau_{\text{uv}}}{\tau_{\text{P}}}e^{-\tau_{\text{uv}}}\right]}{1+2\tau_{\text{P}}\left(1+\dfrac{3}{4}\tau_{\text{R}}\right)}$
(42) $\displaystyle\Gamma_{\text{vis}}=\frac{3}{8\pi}\frac{GM\dot{M}}{R^{3}}.$
(43)

Let us analyze the behavior of (38) in the absence of viscous heating when
$S_{\text{vis}}=0$. For small optical depths with respect to stellar radiation
$\tau_{\text{uv}}\ll 1$ we get:

$aT_{\text{mid}}^{4}=\frac{\kappa_{\text{F}}}{\kappa_{\text{P}}}\frac{F_{0}}{c},$
(44)

where the temperature of the medium is determined by the ratio of the opacity
of the medium to the stellar and thermal radiation. If the disk is optically
thick to stellar radiation $\tau_{\text{uv}}\gg 1$ and the ratio
$\kappa_{\text{R}}/\kappa_{\text{F}}$ is small (which is usually the case),
then:

$aT_{\text{mid}}^{4}=\dfrac{2\mu F_{0}}{c},$ (45)

where midplane temperature depends only on the total flux of stellar radiation
entering the disk. If the disk is optically thick to stellar radiation and
there is viscous heating $S_{\text{vis}}\neq 0$, then:

$aT_{\text{mid}}^{4}=\dfrac{2\mu
F_{0}}{c}+\frac{\Gamma_{\text{vis}}}{c}\left[\frac{1+2\tau_{\text{P}}(1+\frac{3}{4}\tau_{\text{R}})}{\tau_{\text{P}}}\right].$
(46)

Figure 5: Dependences of the midplane temperature on the disk surface density
for different accretion rates $\dot{M}$ = 0, $10^{-8}$ and $10^{-6}$
$M_{\odot}$/year.

As an example, Fig. 5 shows the midplane temperature distributions obtained
using the equation (38), depending on the surface density of the disk. The
following parameters were used: $M=1M_{\odot}$, $L=1L_{\odot}$, $R=1$ au,
$\mu=0.05$, $\kappa_{\text{P}}=\kappa_{\text{R}}=1$ cm${}^{2}/g$,
$\kappa_{\text{F}}=100$ cm${}^{2}/g$. These dependences illustrate the limits
obtained in (44) and (45), and also show the effect of viscous heating.

## References

  * Armitage (2022) Armitage P. J., 2022, arXiv e-prints, arXiv:2201.07262
  * Armitage, Livio & Pringle (2001) Armitage P. J., Livio M., Pringle J. E., 2001, Mon. Not. R. Astron., 324, 705
  * Audard et al. (2014) Audard M. et al., 2014, in Protostars and Planets VI, Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds., pp. 387–410
  * Bae et al. (2022) Bae J., Isella A., Zhu Z., Martin R., Okuzumi S., Suriano S., 2022, arXiv e-prints, arXiv:2210.13314
  * Bell & Lin (1994) Bell K. R., Lin D. N. C., 1994, ApJ, 427, 987
  * Castelli (2005) Castelli F., 2005, Memorie della Societa Astronomica Italiana Supplementi, 8, 34
  * Connelley & Reipurth (2018) Connelley M. S., Reipurth B., 2018, ApJ, 861, 145
  * Contreras Peña et al. (2017) Contreras Peña C. et al., 2017, Mon. Not. R. Astron., 465, 3011
  * D’Alessio et al. (1999) D’Alessio P., Calvet N., Hartmann L., Lizano S., Cantó J., 1999, ApJ, 527, 893
  * D’Angelo & Spruit (2012) D’Angelo C. R., Spruit H. C., 2012, Mon. Not. R. Astron., 420, 416
  * Duschl, Gail & Tscharnuter (1996) Duschl W. J., Gail H. P., Tscharnuter W. M., 1996, Astron. and Astrophys., 312, 624
  * Flaherty et al. (2016) Flaherty K. M., DeMarchi L., Muzerolle J., Balog Z., Herbst W., Megeath S. T., Furlan E., Gutermuth R., 2016, ApJ, 833, 104
  * Hameury (2020) Hameury J. M., 2020, Advances in Space Research, 66, 1004
  * Huré, Richard & Zahn (2001) Huré J. M., Richard D., Zahn J. P., 2001, Astron. and Astrophys., 367, 1087
  * Kawazoe & Mineshige (1993) Kawazoe E., Mineshige S., 1993, Publ. Astron. Soc. Japan, 45, 715
  * Kley & Lin (1999) Kley W., Lin D. N. C., 1999, ApJ, 518, 833
  * Kravtsova et al. (2007) Kravtsova A. S., Lamzin S. A., Errico L., Vittone A., 2007, Astronomy Letters, 33, 755
  * Kurucz (1970) Kurucz R. L., 1970, SAO Special Report, 309
  * Lasota (2001) Lasota J.-P., 2001, New Astronomy Reviews, 45, 449
  * Malygin et al. (2014) Malygin M. G., Kuiper R., Klahr H., Dullemond C. P., Henning T., 2014, Astron. and Astrophys., 568, A91
  * Pavlyuchenkov & Dullemond (2007) Pavlyuchenkov Y., Dullemond C. P., 2007, Astron. and Astrophys., 471, 833
  * Rigon et al. (2017) Rigon L., Scholz A., Anderson D., West R., 2017, Mon. Not. R. Astron., 465, 3889
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, Astron. and Astrophys., 24, 337
  * Vorobyov & Basu (2015) Vorobyov E. I., Basu S., 2015, ApJ, 805, 115
  * Vorobyov et al. (2021) Vorobyov E. I., Elbakyan V. G., Liu H. B., Takami M., 2021, Astron. and Astrophys., 647, A44