# Chiral EFT calculation of neutrino reactions in warm neutron-rich matter

Eunkyoung Shin<EMAIL_ADDRESS>Cyclotron Institute, Texas A&M University,
College Station, TX 77843, USA
Department of Physics and Astronomy, Texas A&M University, College Station, TX
77843, USA  Ermal Rrapaj<EMAIL_ADDRESS>NERSC, Lawrence Berkeley
National Laboratory, Berkeley, California 94720, USA
University of California, Berkeley, CA 94720, USA
RIKEN iTHEMS, Wako, Saitama 351-0198, Japan  Jeremy W. Holt
<EMAIL_ADDRESS>Cyclotron Institute, Texas A&M University, College Station,
TX 77843, USA
Department of Physics and Astronomy, Texas A&M University, College Station, TX
77843, USA  Sanjay K. Reddy<EMAIL_ADDRESS>Institute for Nuclear Theory,
University of Washington, Seattle, WA 98195, USA

###### Abstract

Neutrino scattering and absorption rates of relevance to supernovae and
neutron star mergers are obtained from nuclear matter dynamical structure
functions that encode many-body effects from nuclear mean fields and
correlations. We employ nuclear interactions from chiral effective field
theory to calculate the density, spin, isospin, and spin-isospin response
functions of warm beta-equilibrium nuclear matter. We include corrections to
the single-particle energies in the mean field approximation as well as vertex
corrections resummed in the random phase approximation (RPA), including, for
the first time, both direct and exchange diagrams. We find that correlations
included through the RPA redistribute the strength of the response to higher
energy for neutrino absorption and lower energy for antineutrino absorption.
This tends to suppress the absorption rate of electron neutrinos across all
relevant energy scales. In contrast, the inclusion of RPA correlations
enhances the electron antineutrino absorption rate at low energy and supresses
the rate at high energy. These effects are especially important at high-
density and in the vicinity of the neutrino decoupling region. Implications
for heavy element nucleosynthesis, electromagnetic signatures of compact
object mergers, supernova dynamics, and neutrino detection from galactic
supernovae are discussed briefly.

##  I Introduction

Neutrinos dominate energy, momentum, and lepton number transport in extreme
astrophysical phenomena. Their scattering and absorption rates in hot and
dense nuclear matter play a critical role in core-collapse supernovae Burrows
and Sawyer [1999], O’Connor [2015], Melson  _et al._ [2015], Roberts  _et al._
[2016], proto-neutron star cooling Pons  _et al._ [1999], Roberts and Reddy
[2017], and neutron star mergers Sekiguchi [2010], Wanajo  _et al._ [2014],
Endrizzi  _et al._ [2020], Sumiyoshi  _et al._ [2021], Cusinato  _et al._
[2022]. The effects of nuclear mean fields and correlations on neutrino
reaction rates are encoded in dynamical structure functions related to the
imaginary part of nuclear response functions. In the past, nuclear matter
response functions have been studied using a variety of nuclear interactions
and many-body approximations, including nonrelativistic and relativistic mean
field models Sawyer [1989], Reddy  _et al._ [1998, 1999], Martínez-Pinedo  _et
al._ [2012], Roberts  _et al._ [2012], Pastore  _et al._ [2012], Fermi liquid
theory Iwamoto and Pethick [1982], Burrows and Sawyer [1998, 1999], the virial
expansion Horowitz and Schwenk [2006], Horowitz  _et al._ [2017], Bedaque  _et
al._ [2018], and pseudopotentials Rrapaj  _et al._ [2015], Bedaque  _et al._
[2018]. Recent studies Martínez-Pinedo  _et al._ [2012], Roberts  _et al._
[2012], Rrapaj  _et al._ [2015] have highlighted the important role of nuclear
mean fields for calculating charged-current reactions, such as neutrino and
anti-neutrino absorption, in the supernova neutrinosphere. Here the large
asymmetry between proton and neutron densities leads to a strong splitting of
the proton and neutron mean fields that enhance neutrino absorption and
suppresses anti-neutrino absorption. This, in turn, affects the composition of
matter ejected from supernovae and neutron star mergers as well as neutrino
flavor and energy distributions that terrestrial neutrino detectors may
observe.

In addition, the use of high-precision nucleon-nucleon interactions Bacca  _et
al._ [2012], Rrapaj  _et al._ [2015], Bartl  _et al._ [2016] to calculate
neutrino scattering and reaction cross sections in hot and dense matter has
illustrated the important role of contributions beyond one-pion-exchange
together with nonperturbative resummations of the nucleon-nucleon interaction
in the particle-particle channel (see also Ref. Hanhart  _et al._ [2001] for
the role of these effects on neutrino production). At zero energy transfer,
the momentum dependence of the static density response function of neutron-
rich matter is related to the poorly known isovector gradient contribution in
nuclear energy density functionals used to model neutron star crusts Lim and
Holt [2017], Carreau  _et al._ [2019]. Recent quantum Monte Carlo simulations
of neutron matter using high-precision nuclear forces have provided
constraints on the isovector gradient contribution through studies of neutron
drops Gandolfi  _et al._ [2011] and the static density response function
Buraczynski and Gezerlis [2016].

Figure 1: Zeroth-order and first-order diagrammatic contributions to the
nuclear response function. The zeroth-order diagram is labeled (a), the first-
order direct and exchange mean field contributions are labeled (b) and (c),
and the first-order direct and exchange vertex corrections are labeled (d) and
(e). Solid lines represent nucleons, wavy lines represent the coupling to an
external probe, and dashed lines represent the nucleon-nucleon interaction.

In this study, we employ nuclear forces based on chiral effective field theory
to investigate beyond-mean-field corrections to spin and density response
functions of nuclear matter under ambient conditions typical of supernova and
neutron star merger neutrinospheres, where the nucleon number density varies
in the range $n=10^{11}-10^{13}$ g/cm3 and the temperature varies in the range
$T=5-10$ MeV. Specifically, we calculate Hartree-Fock mean field corrections
as well as resummed particle-hole vertex corrections in the random phase
approximation (RPA). In contrast to naive order-by-order perturbation theory,
the RPA with self-consistent Hartree-Fock mean fields provides a
thermodynamically consistent “conserving approximation” Baym and Kadanoff
[1961], Harris and Evans [1978] for which the computed dynamic structure
function is guaranteed to respect sum rules, detailed balance, and positivity
for all values of the energy and momentum transfer. We find that at the
densities of relevance, RPA vertex corrections can be as important as the mean
field effects studied in earlier work and should be included in a consistent
description of nuclear matter response functions for astrophysical
applications. The paper is organized as follows. In Section II, we present
expressions for density and spin response function in isospin-asymmetric
nuclear matter at finite temperature with mean field and RPA vertex
corrections. In particular, we outline an exact matrix inversion method Harris
and Evans [1978] for calculating RPA response functions, including both direct
and exchange terms for an arbitrary nucleon-nucleon potential. In Section III
we calculate neutral-current and charged-current density and spin response
functions for a range of thermodynamic conditions in beta-equilibrium matter.
Finally, we conclude with a summary and outlook in Section IV.

##  II Response functions in isospin-asymmetric nuclear matter

In the present section we derive expressions for the first-order mean field
and vertex corrections to the response functions of homogeneous nuclear matter
at nonzero temperature. We also outline the calculation of the RPA response
function employing a matrix eigenvalue method. In the region of the supernova
or neutron star merger neutrinospheres, where neutrinos decouple from nuclear
matter and their free-streaming energy spectrum is set, the proton fraction is
small $Y_{p}\approx 0.05-0.10$, the temperature is warm $T=5-10$ MeV, and the
matter is dilute $\rho\approx 10^{11}-10^{13}$ g/cm3. We assume beta
equilibrium, which provides the restriction

$\mu_{n}-\mu_{p}-\mu_{e}=0$ (1)

on the proton and neutron chemical potentials and hence their number densities

$n_{i}=\frac{2}{(2\pi)^{3}}\int d^{3}k\frac{1}{1+e^{(e_{i}(k)-\mu_{i})/T}},$
(2)

for $i=\\{n,p,e\\}$. Together with charge neutrality $\rho_{e}=\rho_{p}$, the
above equations must be solved self consistently to determine the proton,
neutron, and electron chemical potentials for a fixed baryon number density
and temperature.

In the mean field approximation, one also has to calculate the nucleon energy-
momentum dispersion relation for protons
$e_{p}(k)=k^{2}/(2M_{p})+\Sigma_{p}(k)$, and likewise for neutrons, where
$\Sigma_{p}(k)$ is the self-energy. The Hartree-Fock conserving approximation
Harris and Evans [1978] consists of computing the irreducible self-energy to
first order in perturbation theory and the response function to all orders in
the random phase approximation. The resulting perturbative approximation to
the response function is then guaranteed to satisfy properties of the exact
response function, such as sum rules and strictly positive dynamical structure
functions.

In the present work we employ nucleon-nucleon (NN) potentials derived from
chiral effective field theory (ChEFT) Entem and Machleidt [2003], Coraggio
_et al._ [2007], Marji  _et al._ [2013], which provides a systematic expansion
of nuclear two and many-body forces. In ChEFT the long-range part of the
nuclear force comes from pion-exchange processes constrained by chiral
symmetry, while the short-range part is encoded in a set of contact terms
fitted to nucleon-nucleon scattering and deuteron properties. The high-
momentum components of the ChEFT nuclear potentials employed in this study are
regulated by exponential functions with a characteristic momentum scale
$\Lambda$ and smoothness parameter $n$:

$f(p,p^{\prime})=\exp[-(p^{\prime}/\Lambda)^{2n}-(p/\Lambda)^{2n}],$ (3)

where $\vec{p}$ and $\vec{p}^{\,\prime}$ are the incoming and outgoing
relative momenta of the two nucleons. We employ ChEFT potentials with cutoff
scales $\Lambda=414,450,500$ MeV and associated smoothness parameters
$n=10,3,2$ respectively. When supplemented by the ChEFT three-body force that
appears at next-to-next-to-leading order (N2LO) in the chiral expansion, this
set of nuclear potentials has been shown to predict well the properties of
nuclear matter, such as the equation of state Holt and Kaiser [2017],
Wellenhofer  _et al._ [2014], optical potential Whitehead  _et al._ [2021],
and quasiparticle interaction in Fermi liquid theory Holt  _et al._ [2012]. In
dilute neutron-rich matter, three-body forces give negligible contribution to
single-particle and response properties of the medium. We therefore neglect
three-body forces in the present work.

###  II.1 Nuclear matter response functions at 0th order in perturbation
theory

In Figure 1 we show the zeroth-order and first-order perturbation theory
contributions to nuclear response functions. The wavy lines denote the
coupling to a $W$ or $Z$ boson defined in the non-relativistic limit by its
vector/axial vector and isoscalar/isovector nature. The dashed lines in Figure
1 represent the nucleon-nucleon interaction.

The zeroth-order contribution to the neutral-current density response function
$\chi_{\rho}$, the neutral-current spin response function $\chi_{\sigma}$, the
charged-current density response function $\chi_{\tau\rho}$ for electron-
neutrino absorption, and the charged-current spin response function for
electron-neutrino absorption $\chi_{\tau\sigma}$ in isospin-asymmetric nuclear
matter at nonzero temperature are given by

$\displaystyle\chi^{(0)}_{\rho}(\vec{q},\omega)$
$\displaystyle=\sum_{s_{1}s_{2}t_{1}t_{2}}\\!\int\\!\frac{d\vec{k}}{(2\pi)^{3}}\frac{f_{\vec{k},t_{1}}-f_{\vec{k}+\vec{q},t_{2}}}{\omega+e_{\vec{k},t_{1}}-e_{\vec{k}+\vec{q},t_{2}}+i\eta}\delta_{s_{1},s_{2}}\delta_{t_{1},t_{2}},$
(4) $\displaystyle\chi^{(0)}_{\sigma}(\vec{q},\omega)$
$\displaystyle=\sum_{s_{1}s_{2}t_{1}t_{2}}\\!\int\\!\frac{d\vec{k}}{(2\pi)^{3}}\frac{f_{\vec{k},t_{1}}-f_{\vec{k}+\vec{q},t_{2}}}{\omega+e_{\vec{k},t_{1}}-e_{\vec{k}+\vec{q},t_{2}}+i\eta}|\langle
s_{1}|\sigma_{z}|s_{2}\rangle|^{2}\delta_{t_{1},t_{2}},$
$\displaystyle\chi^{(0)}_{\tau\rho}(\vec{q},\omega)$
$\displaystyle=\sum_{s_{1}s_{2}}\\!\int\\!\frac{d\vec{k}}{(2\pi)^{3}}\frac{f_{\vec{k},n}-f_{\vec{k}+\vec{q},p}}{\omega+e_{\vec{k},n}-e_{\vec{k}+\vec{q},p}+i\eta}\delta_{s_{1},s_{2}},$
$\displaystyle\chi^{(0)}_{\tau\sigma}(\vec{q},\omega)$
$\displaystyle=\sum_{s_{1}s_{2}}\\!\int\\!\frac{d\vec{k}}{(2\pi)^{3}}\frac{f_{\vec{k},n}-f_{\vec{k}+\vec{q},p}}{\omega+e_{\vec{k},n}-e_{\vec{k}+\vec{q},p}+i\eta}|\langle
s_{1}|\sigma_{z}|s_{2}\rangle|^{2},$

where the sums are over the single-particle spin projections $s_{1},s_{2}$ and
isospin projections $t_{1},t_{2}$ of the particle-hole pair, $f_{\vec{k},t}$
is the Fermi-Dirac distribution function for a nucleon with momentum $\vec{k}$
and isospin projection $t$, and $e_{\vec{k},t}$ is the single-particle energy
for a nucleon with momentum $\vec{k}$ and isospin projection $t$.

In the case of noninteracting protons and neutrons, the density and spin
response functions in the neutral-current or charged-current channels are
identical since
$\sum_{s_{1}s_{2}}\delta_{s_{1},s_{2}}=\sum_{s_{1}s_{2}}|\langle
s_{1}|\sigma_{z}|s_{2}\rangle|^{2}=2$. In addition, the single-particle
energies in Eq. (4) are simply the free-space kinetic energies
$e_{\vec{k}}=k^{2}/(2M)$. Including effects from momentum-dependent mean
fields for protons $\Sigma_{p}(k)$ and neutrons $\Sigma_{n}(k)$, the single-
particle energies that enter in the Fermi-Dirac distribution functions and the
energy denominators of Eq. (4) become

$e_{\vec{k},n}=\frac{k^{2}}{2M}+\Sigma_{n}(k),\hskip
14.45377pte_{\vec{k},p}=\frac{k^{2}}{2M}+\Sigma_{p}(k).$ (5)

At the Hartree-Fock level (first-order perturbation theory), the proton and
neutron single-particle potentials can be well described by the effective mass
approximation

$e_{\vec{k},n}=\frac{k^{2}}{2M_{n}^{*}}+U_{n},\hskip
14.45377pte_{\vec{k},p}=\frac{k^{2}}{2M^{*}_{p}}+U_{p}.$ (6)

In the case that $M^{*}_{n,p}\simeq M$, the mean field strengths $U_{n}$ and
$U_{p}$ modify only the energy denominators in Eq. (4), which results in a
shift of the imaginary part of the response functions. Although the mean field
energies $e_{\vec{k},n}$ and $e_{\vec{k},p}$ also enter in the definition of
the Fermi-Dirac distribution functions, a simple mean field shift will be
absorbed into a redefinition of the chemical potential. The inclusion of
nuclear mean fields at the Hartree-Fock level in the calculation of the
$0^{\rm th}$-order response functions $\chi^{(0,{\rm MF})}(q,\omega)$
corresponds to iterating diagrams of type (b) and (c) in Figure 1 to all
orders in perturbation theory. For neutrino scattering and absorption, the
main effect will be a shift of the $\omega$-dependent imaginary response for
fixed momentum transfer $\vec{q}$. For additional technical details regarding
the calculation of nucleon single-particle potentials starting from realistic
nucleon-nucleon interactions, the reader is referred to Refs. Holt  _et al._
[2013], Rrapaj  _et al._ [2015].

###  II.2 Random Phase Approximation (RPA)

Figure 2: Diagrammatic contributions to the RPA response function. Solid lines
represent nucleons, wavy lines represent the coupling to an external probe,
and the dashed lines represent the nucleon-nucleon interaction.

In many-body perturbation theory, an order-by-order calculation of response
functions can lead to unphysical dynamic structure functions that do not
satisfy constraints, such as positivity and sum rules that relate the long-
wavelength response to thermodynamics. The formal solution to this problem
obtained by Baym and Kadanoff Baym and Kadanoff [1961] is to construct so-
called “conserving approximations” that relate the one-body and two-body
propagators in such a way as to maintain conservation laws of energy,
momentum, angular momentum, and particle number. In the Hartree-Fock
conserving approximation, the one-body propagator is constructed at the self-
consistent mean field level, while the two-body propagator resums to all
orders the direct and exchange particle-hole diagrams in the RPA. The RPA
leads to a linear inhomogeneous integral equation for the particle-hole vertex
function. The discretized version of this equation can be re-expressed Harris
and Evans [1978] as a linear algebraic equation that can be solved through
matrix inversion. In the rest of this section, we will outline the solution of
the RPA response function and present several benchmarks to test the method.

Figure 3: Vertex function in the random phase approximation. The wavy line
represents the weak current, the dashed line is the nucleon-nucleon
interaction carrying energy and momentum, and the dark square is the dressed
vertex. The solid lines represent propagating nucleons.

In Figure 2, we show the class of response function diagrams to be resummed in
the RPA, which includes all particle-hole bubble diagrams, ladder diagrams,
and combinations thereof with dressed intermediate-state propagators in the
mean field approximation. This infinite set of diagrammatic contributions to
nuclear matter response functions cannot be resummed directly but instead must
proceed through the resummed particle-hole density vertex function
$L(\vec{k}_{1}s_{1};\vec{q}\,\omega)$. In the case of the charged-current
density response function, the density vertex function satisfies

$\chi_{\tau\rho}(\vec{q},\omega)=\sum_{s}\int\\!\frac{d\vec{k}}{(2\pi)^{3}}L(\vec{k}s;\vec{q}\,\omega).$
(7)

The 0th-order and 1st-order vertex functions are therefore given by

$\displaystyle L^{(0)}(\vec{k}s;\vec{q}\,\omega)$
$\displaystyle=\frac{f_{\vec{k},n}-f_{\vec{k}+\vec{q},p}}{\omega+e_{\vec{k},n}-e_{\vec{k}+\vec{q},p}+i\eta},$
(8) $\displaystyle L^{(1)}(\vec{k}s;\vec{q}\,\omega)$
$\displaystyle=\frac{f_{\vec{k},n}-f_{\vec{k}+\vec{q},p}}{\omega+e_{\vec{k},n}-e_{\vec{k}+\vec{q},p}+i\eta}\sum_{s^{\prime}}\int\\!\frac{d\vec{k}^{\prime}}{(2\pi)^{3}}\frac{(f_{\vec{k}^{\prime},n}\\!-f_{\vec{k}^{\prime}+\vec{q},p})\langle\vec{k}\vec{k}^{\prime}\\!+\\!\vec{q},ss^{\prime},np|\bar{V}|\vec{k}\\!+\\!\vec{q}\,\vec{k}^{\prime},ss^{\prime},pn\rangle}{\omega+e_{\vec{k}^{\prime},n}-e_{\vec{k}^{\prime}+\vec{q},p}+i\eta}.$

Summing the higher-order particle-hole bubble and ladder diagrams, one obtains
the integral equation

$\displaystyle
L(\vec{k}s;\vec{q}\,\omega)=L_{0}(\vec{k}s;\vec{q}\,\omega)+L_{0}(\vec{k}s;\vec{q}\,\omega)$
$\displaystyle\times\sum_{s^{\prime}}\int\\!\frac{d\vec{k}^{\prime}}{(2\pi)^{3}}\langle\vec{k}\vec{k}^{\prime}\\!+\\!\vec{q},ss^{\prime},np|\bar{V}|\vec{k}\\!+\\!\vec{q}\,\vec{k}^{\prime},ss^{\prime},pn\rangle$
$\displaystyle\times L(\vec{k}^{\prime}s^{\prime};\vec{q}\,\omega),$ (9)

shown diagrammatically in Figure 3. For a spin-saturated system, the vertex
function for the density response function is independent of spin. One can
then average Eq. (II.2) over the spin to obtain

$\displaystyle
L(\vec{k};\vec{q}\,\omega)=L_{0}(\vec{k};\vec{q}\,\omega)+L_{0}(\vec{k};\vec{q}\,\omega)\int\\!\frac{d\vec{k}^{\prime}}{(2\pi)^{3}}L(\vec{k}^{\prime};\vec{q}\,\omega)$
$\displaystyle\times\left[\frac{1}{2}\sum_{ss^{\prime}}\langle\vec{k}\vec{k}^{\prime}\\!+\\!\vec{q},ss^{\prime},np|\bar{V}|\vec{k}\\!+\\!\vec{q}\,\vec{k}^{\prime},ss^{\prime},pn\rangle\right],$
(10)

where for convenience we will denote the quantity in squared brackets as
${\cal V}(\vec{k},\vec{k}^{\prime})$, suppressing the explicit dependence on
$\vec{q}$. Writing the integral in Eq. (II.2) as a summation over a discrete
set $\\{\vec{k}_{1},\vec{k}_{2},\dots\\}$ of momentum-space mesh points with
associated mesh weights $\\{w_{1},w_{2},\dots\\}$, one can rewrite Eq. (II.2)
as a matrix equation whose formal solution is

$L=\left[N^{-1}\left(E+(\omega+i\eta)\mathbb{1}\right)-{\cal V}\right]^{-1}B,$
(11)

where $L$ is a vector with elements

$L=\begin{bmatrix}L(\vec{k}_{1};\vec{q}\,\omega)\\\
L(\vec{k}_{2};\vec{q}\,\omega)\\\ \vdots\end{bmatrix},$ (12)

$N$ is a diagonal matrix with elements

$N=\begin{pmatrix}f_{\vec{k}_{1},n}-f_{\vec{k}_{1}+\vec{q},p}&0&\cdots\\\
0&f_{\vec{k}_{2},n}-f_{\vec{k}_{2}+\vec{q},p}&\cdots\\\
\vdots&\vdots&\ddots\end{pmatrix},$ (13)

$E$ is a diagonal matrix with elements

$E=\begin{pmatrix}e_{\vec{k}_{1},n}-e_{\vec{k}_{1}+\vec{q},p}&0&\cdots\\\
0&e_{\vec{k}_{2},n}-e_{\vec{k}_{2}+\vec{q},p}&\cdots\\\
\vdots&\vdots&\ddots\end{pmatrix},$ (14)

${\cal V}$ is the matrix

${\cal V}=\begin{pmatrix}w_{1}{\cal V}(\vec{k}_{1},\vec{k}_{1})&w_{2}{\cal
V}(\vec{k}_{1},\vec{k}_{2})&\cdots\\\ w_{1}{\cal
V}(\vec{k}_{2},\vec{k}_{1})&w_{2}{\cal V}(\vec{k}_{2},\vec{k}_{2})&\cdots\\\
\vdots&\vdots&\ddots\end{pmatrix},$ (15)

and $B$ is a vector whose elements are all 1.

Figure 4: Imaginary part of the charged-current density response function at
$0^{\rm th}$ order (blue dots) and computed from the RPA resummation of the
direct term analytically (black dots) and within the eigenvalue formalism
described in Section II.2 (red line). A smearing length $\epsilon=70\times
10^{-4}\,{\rm fm}^{-1}$ in Eq. (20) is employed.

Our goal will be to extract the imaginary part of the vertex function, which
is related to nuclear matter dynamical structure functions and neutrino
scattering cross sections. In order for the imaginary part of Eq. (11) to be
nonzero, the matrix $N^{-1}(E+\omega\mathbb{1})-{\cal V}$ must be singular,
which occurs when $\omega$ takes on the values defined by

$(NV-E)|\ell\rangle=\omega_{\ell}|\ell\rangle.$ (16)

In the vicinity of $\omega^{\ell}$, one can write

$\displaystyle[N^{-1}(E+(\omega+i\eta)\mathbb{1})-{\cal V}]^{-1}$
$\displaystyle=\frac{1}{\langle\ell|N^{-1}|\ell\rangle}\left[\frac{{\rm
Pr}}{\omega-\omega_{\ell}}-i\pi\delta(\omega-\omega_{\ell})\right]|\ell\rangle\langle\ell|,$
(17)

where ${\rm Pr}$ denotes the principal value. The imaginary part of the
response function is then given by

${\rm Im}\,\chi_{\tau\rho}^{\rm
RPA}(\vec{q},\omega)=-i\pi\sum_{\ell}\frac{\langle
B|\ell\rangle^{2}}{\langle\ell|N^{-1}|\ell\rangle}\delta(\omega-\omega_{\ell}).$
(18)

In terms of the discrete momentum-space mesh points and weights, we have

${\rm Im}\,\chi_{\tau\rho}^{\rm
RPA}(\vec{q},\omega)=-i\pi\sum_{\ell}\frac{\left(\sum_{i}w_{i}|\ell\rangle_{i}\right)^{2}}{\sum_{i}(w_{i}|\ell\rangle_{i})^{2}(w_{i}N_{i})^{-1}}\delta(\omega-\omega_{\ell}),$
(19)

where $|\ell\rangle_{i}$ denotes the $i^{\rm th}$ element of the eigenvector
$|\ell\rangle$. In practice, the finite number of $\delta$ functions in Eq.
(19) obtained by discretizing the integral in Eq. (II.2) must be appropriately
smeared to obtain a continuous response function. We employ the approximation

$\delta_{\epsilon}(\omega)=\frac{1}{\sqrt{2\pi\epsilon}}e^{-\omega^{2}/2\epsilon}.$
(20)

and find that it is always possible to choose a smearing length $\epsilon$
that leads to a converged result.

Figure 5: Imaginary part of the neutral-current spin response function at
$0^{\rm th}$ order (blue dots) and computed from the RPA resummation of the
direct term analytically (black dots) and within the eigenvalue formalism
described in Section II.2 (red line). A smearing length $\epsilon=70\times
10^{-4}\,{\rm fm}^{-1}$ in Eq. (20) is employed.

The above eigenvalue method is suitable to resum both the direct and exchange
RPA bubble diagrams to all orders. To benchmark the method, we consider the
simplified case of iterating just the direct part of the nuclear potential to
all orders, which for a local potential $V(q)$ can be computed analytically
through

$\chi^{\rm
RPA}(\vec{q},\omega)=\frac{\chi^{(0)}(\vec{q},\omega)}{1-V(q)\chi^{(0)}(\vec{q},\omega)}.$
(21)

In Figure 4 we show the energy dependence of the charged-current density
response function in beta-equilibrium nuclear matter at density $n=50\times
10^{11}$ g/cm3 and temperature $T=7$ MeV for a momentum transfer $q=21$ MeV
assuming a scalar-isovector interaction

$V(q)=\frac{g^{2}}{m^{2}+q^{2}}\vec{\tau}_{1}\cdot\vec{\tau}_{2},$ (22)

where we take $g=5$ and $m=700$ MeV. In Figure 4 we plot the response function
in the approximation of noninteracting particles (blue dots), the exact
analytical RPA result (black dots), and the numerical eigenvalue RPA result
(red line) with a delta function smearing length of $\epsilon=70\times
10^{-4}\,{\rm fm}^{-1}$. One sees excellent agreement between the exact and
numerical RPA resummations.

The treatment of spin response functions proceeds similarly, except that we
obtain a pair of coupled equations for the spin-up and spin-down response
functions, which effectively doubles the dimensionality of the vectors and
matrices in Eqs. (12) – (15). For the spin response functions, we also do not
use the spin-averaging approximation employed to obtain Eq. (II.2). In Figure
5 we show the energy dependence of the neutral-current spin response function
in beta-equilibrium nuclear matter at density $n=50\times 10^{11}$ g/cm3 and
temperature $T=7$ MeV for a momentum transfer $q=21$ MeV assuming an isoscalar
spin-spin interaction

$V(q)=\frac{g^{2}}{m^{2}+q^{2}}\vec{\sigma}_{1}\cdot\vec{\sigma}_{2},$ (23)

where we take $g=10$ and $m=700$ MeV. In Figure 4 we plot the response
function in the approximation of noninteracting particles (blue dots), the
exact analytical RPA result (black dots), and the numerical eigenvalue RPA
result (red line) with a delta function smearing length of $\epsilon=70\times
10^{-4}\,{\rm fm}^{-1}$. Again we find excellent agreement between the exact
and numerical RPA resummations.

###  II.3 Dynamic structure functions

Neutrino opacities are a key input to numerical simulations of core-collapse
supernovae, proto-neutron star evolution, and neutron star mergers. Both
neutral-current and charged-current weak reactions are important sources of
neutrino opacity across a wide range of densities and temperatures. Matter
effects on neutrino scattering and absorption cross sections on baryons are
encoded in dynamical structure functions related to the imaginary part of
nuclear response functions.

####  II.3.1 Neutral-current neutrino scattering

The double differential cross section for low-energy neutrinos to scatter in a
non-relativistic gas of nucleons is given by

$\displaystyle\frac{1}{V}\frac{d^{2}\sigma}{d\cos\theta\,d\omega}=\frac{G_{F}^{2}}{4\pi^{2}}(E_{\nu}-\omega)^{2}$
(24)
$\displaystyle\times\left[c_{V}^{2}(1+\cos\theta)S_{\rho}(q,\omega)+c_{A}^{2}(3-\cos\theta)S_{\sigma}(q,\omega)\right],$

where $S_{\rho}$ is the neutral-current density structure function and
$S_{\sigma}$ is the neutral-current spin structure function. The energy
transfer is given by $\omega=E_{\nu}-E_{\nu}^{\prime}$ and the momentum
transfer is given by $\vec{q}=\vec{p}_{\nu}-\vec{p}_{\nu}^{\,\prime}$ with
magnitude
$q=\sqrt{E_{\nu}^{2}+{E_{\nu}^{\prime}}^{2}-2E_{\nu}E_{\nu}^{\prime}\cos\theta}$.
The structure functions in Eq. (24) are related to the imaginary parts of the
associated response functions by

$S(q,\omega)=-\frac{2\,{\rm Im}\chi(q,\omega)}{1-e^{-\omega/T}}.$ (25)

Figure 6: Imaginary part of the charged-current density response function for
a noninteracting gas (dotted blue line) and including the effects of nuclear
mean fields in the Hartree-Fock approximation for different chiral nuclear
potentials labeled by their momentum-space cutoff value: 414 MeV, 450 MeV, and
500 MeV. Figure 7: Imaginary part of the charged-current spin response
function of beta-equilibrium nuclear matter at a density $n=0.002$ fm-3 and
temperature $T=5$ MeV. (Top panel) spin response function including mean field
“MF” corrections from the N3LO-414 chiral nuclear potential. (Bottom panel)
spin response function including random phase approximation vertex corrections
plus nuclear mean fields “RPA+MF”.

####  II.3.2 Charge-current neutrino absorption

The double differential cross section for electron neutrino absorption is
given by Roberts  _et al._ [2012], Roberts and Reddy [2017], Iwamoto and
Pethick [1982]

$\displaystyle\frac{1}{V}\frac{d^{2}\sigma}{d\cos\theta
dE_{e}}=\frac{G_{F}^{2}\cos^{2}\theta_{c}}{4\pi^{2}}p_{e}E_{e}(1-f_{e}(E_{e}))$
(26)
$\displaystyle\times[(1+\cos\theta)S_{\tau\rho}(\omega,q)+g_{A}^{2}(3-\cos\theta)S_{\tau\sigma}(\omega,q)],$

where $S_{\tau\rho}$ is the charged-current density structure function and
$S_{\tau\sigma}$ is the charged-current spin structure function. The energy
transfer is given by $\omega=E_{\nu}-E_{e}$ and the momentum transfer is given
by $\vec{q}=\vec{p}_{\nu}-\vec{p}_{e}$ with magnitude
$q=\sqrt{E_{\nu}^{2}+E_{e}^{2}-2E_{\nu}E_{e}\cos\theta}$. The structure
functions in Eq. (26) are related to the imaginary part of the associated
response functions by

$S_{\tau}(q,\omega)=-\frac{2\,{\rm
Im}\chi(q,\omega)}{1-e^{-(\omega+\mu_{n}-\mu_{p})/T}},$ (27)

where the detailed balance factor depends explicitly on the proton and neutron
chemical potentials $\mu_{p}$ and $\mu_{n}$.

Figure 8: Imaginary part of the charged-current spin response function of
beta-equilibrium nuclear matter at a density $n=0.02$ fm-3 and temperature
$T=8$ MeV. (Top panel) spin response function including mean field “MF”
corrections from the N3LO-414 chiral nuclear potential. (Bottom panel) spin
response function including random phase approximation vertex corrections plus
nuclear mean fields “RPA+MF”.

##  III Results

###  III.1 Response functions in the mean field approximation

Figure 9: Energy-dependent neutral-current (left) and charged-current (right)
dynamic structure functions of beta-equilibrium nuclear matter at density
$n=0.02$ fm-3 and temperature $T=8$ MeV for a momentum transfer of $q=24$ MeV.
We plot both the density (top) and spin (bottom) structure functions,
including the effects of nucleon mean fields (dashed lines) and RPA
correlations (black curves).

In Figure 6 we show mean field effects on the charged-current density response
function of nuclear matter for different choices of the nuclear potential:
N3LO-414, N3LO-450, and N3LO-500. In the top panel, we consider beta-
equilibrium nuclear matter at density $n=0.002\,\text{fm}^{-3}$, temperature
$T=5$ MeV, and assuming a momentum transfer of $q=15$ MeV. In the bottom
panel, we consider beta-equilibrium nuclear matter at density
$n=0.02\,\text{fm}^{-3}$, temperature $T=8$ MeV, and assuming a momentum
transfer of $q=24$ MeV. We find that mean field effects are larger for smaller
values of the momentum-space cutoff. A similar effect was observed in the
context of the nuclear equation of state Holt and Kaiser [2017] and the
nuclear single-particle potential Holt  _et al._ [2013, 2016]. Namely, low-
cutoff potentials are more perturbative and therefore generate more attraction
at first order in perturbation theory, even though the sum of first- and
second-order perturbation theory contributions to the nuclear equation of
state or single-particle potential are similar. As argued in Section II, the
largest effect of mean field corrections in the imaginary part of the neutrino
absorption charged-current response is to shift the strength by an amount
$U_{p}-U_{n}$ in the energy transfer $\omega$. At the lower value of the
density, $n=0.002\,\text{fm}^{-3}$, the mean field splitting is on the order
of $1\,\text{MeV}$, while at the larger density $n=0.02\,\text{fm}^{-3}$, the
mean field spliting is nearly $10\,\text{MeV}$ for the N3LO-414 chiral
nucleon-nucleon potential. However, such mean field splittings are still about
a factor of 2 less than those generated from typical phenomenological mean
field models Roberts  _et al._ [2012], Rrapaj  _et al._ [2016].

###  III.2 Response functions in the random phase approximation

The random phase approximation for the vertex function together with the
Hartree-Fock approximation for the single-particle energies represents a
conserving approximation that is guaranteed to preserve sum rules and the
positivity of dynamical structure functions. In addition, it is able to
capture the presence of collective oscillations such as the giant-dipole and
Gamow-Teller resonances that are known to play a role in the response of
nuclei Reddy  _et al._ [1999]. In Figure 7 we show a contour plot of the
imaginary part of the neutrino-absorption charged-current spin response as a
function of the energy $\omega$ and momentum $q$ transfer for beta-equilibrium
matter with density $n=0.002\,\text{fm}^{-3}$ and temperature $T=5$ MeV. The
top panel shows the imaginary part of the response including mean field (MF)
corrections alone, while the bottom panel shows the combined effect of RPA
correlations and mean fields (RPA + MF). We see that the MF response exhibits
a relatively broad distribution already at low momentum transfers that peaks
at a nearly constant energy $\omega\simeq-2$ MeV. In contrast, the RPA + MF
response remains sharply peaked for longer, up to a momentum transfer $q\simeq
0.04\,\text{fm}^{-1}\simeq 8$ MeV, and peaks at a value $\omega\simeq-1$ MeV
that increases slowly with $q$. The sharper structure of the RPA + MF response
is indicative of a collective mode. This feature is even more evident in
Figure 7, where we show the contour plot of the imaginary part of the
neutrino-absorption charged-current spin response for beta-equilibrium matter
with density $n=0.02\,\text{fm}^{-3}$ and temperature $T=8$ MeV. Here the
collective mode (now at positive energy) remains sharp up to a momentum
transfer $q\simeq 0.2\,\text{fm}^{-1}\simeq 40$ MeV. The shift in peak energy
of the imaginary response from $-7.5$ MeV in the mean field approximation to
$2$ MeV in the RPA + MF approximation will have an important effect on
electron neutrino absorption in dilute beta-equilibrium nuclear matter. The
shift will push the outgoing electron energy to smaller values where Pauli
blocking acts to reduce the available phase space, thereby suppressing the
absorption cross section and increasing the mean free path.

Figure 10: Differential neutrino absorption cross section per unit volume as a
function of the outgoing electron energy $E_{e}$ computed from the spin
response function of beta-equilibrium nuclear matter at density $n=0.002$ fm-3
and temperature $T=5$ MeV. The incoming neutrino energy is set to be
$E_{\nu}=3T=15$ MeV. Curves are shown for the case of noninteracting nucleons
(blue dotted line), nuclear mean fields computed from the N3LO-414 chiral
nuclear potential (blue dashed lines), and in the RPA + mean field
approximation (black solid line) also employing the N3LO-414 potential.

In the left panels of Figure 9 we plot the neutral-current density (top) and
spin (bottom) dynamic structure factors for beta equilibrium nuclear matter
under the ambient conditions $T=8$ MeV, $n=0.02$ fm-3, and momentum transfer
of $q=24$ MeV. For both the neutral-current density and spin response
functions, the mean field corrections play only a very minor role. Mean fields
change slightly the proton and neutron densities for beta equilibrium matter,
and the effective mass has only a small effect on the Fermi distribution
functions and energy denominators. The energy shifts in the single-particle
potentials are absorbed into redefinitions of the proton and neutron chemical
potentials for the Fermi distribution functions, and the mean field shifts
cancel in the response function energy denominators. The inclusion of RPA
correlations, however, is very important, enhancing the density structure
function and suppressing the spin structure function, a feature already
observed Bedaque  _et al._ [2018] including first-order vertex corrections
starting from a pseudopotential defined in terms of nucleon-nucleon scattering
phase shifts. This behavior can be traced to the large neutron-neutron
attraction in the ${}^{1}S_{0}$ partial wave.

Figure 11: Same as Figure 10 except for a density $n=0.02$ fm-3, temperature
$T=8$ MeV, and incoming neutrino energy $E_{\nu}=24$ MeV.

In the right panels of Figure 9 we plot the charged-current density (top) and
spin (bottom) dynamic structure functions in beta-equilibrium nuclear matter
under the ambient conditions $T=8$ MeV, $n=0.02$ fm-3, and momentum transfer
of $q=24$ MeV. Whereas mean fields drive absorption strength to lower energy
transfers, RPA correlations significantly shift strength to larger energy
transfers. This is due to the existence of giant dipole and Gamow-Teller
collective modes in the isovector-density and isovector-spin channels. The
combined effect of nuclear mean fields and RPA correlations is a
redistribution of strength to energies above that of the noninteracting
charged-current response functions. We conclude that the inclusion of vertex
corrections is crucial for an accurate description of nuclear matter response
functions at and above neutrinosphere densities.

###  III.3 Energy-dependent neutrino absorption cross section

Integrating Eq. (26) over the scattering angle $\theta$, we obtain the
differential energy-dependent neutrino absorption cross section

$\displaystyle\frac{1}{V}\frac{d\sigma}{dE_{e}}=\frac{G_{F}^{2}\cos^{2}\theta_{c}}{4\pi^{2}}p_{e}E_{e}(1-f_{e}(E_{e}))$
(28) $\displaystyle\times\int
d\cos\theta[(1+\cos\theta)S_{\tau\rho}(\omega,q)+g_{A}^{2}(3-\cos\theta)S_{\tau\sigma}(\omega,q)],$

where $S_{\tau\rho}$ and $S_{\tau\sigma}$ are the density and spin dynamic
structure functions. In Figures 10 and 11 we plot the differential cross
section for electron neutrino absorption, keeping only the spin dynamic
structure function in Eq. (28) for beta-equilibrium matter at two densities
$n=0.002$ fm-3 and $n=0.02$ fm-3. The incoming neutrino energy is
$E_{\nu}=3T$. We show the cross section assuming noninteracting nucleons (blue
dotted line), the cross section keeping only nuclear mean fields at the
Hartree-Fock level (blue dashed line), and finally the cross section including
RPA correlations and Hartree-Fock mean fields (black solid line). The nuclear
force is taken to be the N3LO-414 chiral two-body interaction. One finds that
the inclusion of nuclear mean fields enhances neutrino absorption since the
response is shifted to lower energy transfers $\omega$ and, therefore, higher
electron energies for which the Pauli suppression factor is reduced. However,
RPA correlations provide a stronger shift in the response function toward
higher energy transfers, leading to a reduced absorption cross section whose
peak in electron energy can be shifted below that of the non-interacting Fermi
gas.

###  III.4 Mean free path

Figure 12: Inverse neutrino and antineutrino absorption mean free paths as a
function of energy $E_{\nu}$ in beta-equilibrium nuclear matter at density
$n=0.002$ fm-3 and temperature $T=5$ MeV. Shown are results for (i)
noninteracting nucleons (dotted lines), (ii) mean field corrections (dashed
lines), and (ii) RPA + mean field corrections (solid lines) all calculated
using the N3LO-414 chiral nucleon-nucleon potential. Figure 13: Inverse
neutrino and antineutrino absorption mean free paths as a function of energy
$E_{\nu}$ in beta-equilibrium nuclear matter at density $n=0.02$ fm-3 and
temperature $T=8$ MeV. Shown are results for (i) noninteracting nucleons
(dotted lines), (ii) mean field corrections (dashed lines), and (ii) RPA +
mean field corrections (solid lines) all calculated using the N3LO-414 chiral
nucleon-nucleon potential.

We now present and discuss results for the mean free path of electron and
anti-electron neutrinos in the vicinity of the neutrinosphere due to charged
current interactions. Differences between these mean free paths directly
impact several key observable aspects of supernovae and neutron star mergers
including dynamics, nucleosynthesis, and neutrino oscillations.

The inverse of the electron and anti-electron neutrino mean free path due to
their charged current interactions is obtained by integrating the differential
absorption cross section per unit volume over the final-state lepton energy,
e.g.,

$\frac{1}{\lambda}=\int\frac{1}{V}\frac{d\sigma}{dE_{e}}dE_{e}\,.$ (29)

In Figures 12 and 13 we plot the inverse mean free paths of electron neutrino
and antineutrino absorption as a function of the incident energy for two sets
of ambient conditions $(n,T)=(0.002\,\text{fm}^{-3},5\,\text{MeV})$ and
$(0.02\,\text{fm}^{-3},8\,\text{MeV})$. The dynamic structure functions are
computed from the associated charged-current spin response functions in three
approximations. First, the inverse neutrino mean free paths neglecting
interactions between nucleons is shown by the dotted curves. Second, we show
as the dashed lines the effect of introducing proton and neutron mean fields
in the Hartree-Fock approximation employing the N3LO-414 chiral nucleon-
nucleon interaction. Finally, the solid curves show the combined effects of
nucleon mean fields and vertex corrections obtained in the random phase
approximation. We find that for the N3LO-414 potential, the mean-field effects
significantly enhance the electron neutrino absorption cross-section across
all energies considered in agreement with earlier studies. However, in
contrast, RPA correlations redistribute strength to the vicinity of the
positive-energy collective mode. This significantly reduces the outgoing
electron energy into a region where Pauli blocking suppresses the reaction.
This redistribution of strength due to a broad collective mode shifts the
response to higher energy and undoes the enhancement of the inverse mean free
path due to mean-field effects. Remarkably, correlations suppress the electron
neutrino absorption cross-sections over the entire energy range and are
especially large for low-energy neutrinos.

In Figure 14 we show the imaginary part of the spin response function for
antineutrino absorption in beta-equilibrium nuclear matter at density
$n=0.02\,\text{fm}^{-3}$ and temperature $T=8\,\text{MeV}$. For electron
antineutrinos, nuclear mean fields shift the response to higher energies. The
absorption cross-section is therefore greatly reduced because the threshold
energy to convert protons into neutrons is increased to such an extent that
there is little phase space available for the reaction. As observed in Figure
13, the corresponding mean free path increases dramatically at low energies in
agreement with earlier work. However, the presence of the negative-energy
collective mode due to RPA correlations, shown in the lower panel of Figure
14, lowers the energy required for the process. It provides a reaction pathway
even at low antineutrino energies, thereby increasing the cross-section. The
reaction at low energies can be viewed as a process involving the absorption
of a positively charged collective mode by the anti-electron neutrino to
produce a positron in the final state. At higher antineutrino energies, the
strong coupling to the negative-energy collective mode weakens the absorption
cross section through the detailed balance factor
$\left(1-e^{-(\omega+\mu_{p}-\mu_{n})/T}\right)^{-1}$. The response function
is also narrowly peaked in this region, which limits the available phase space
for final-state positron energies. Overall, at high energies the antineutrino
absorption total cross section is reduced relative to both the noninteracting
and mean field approximations as seen in Figure 13.

##  IV Conclusion

We have developed the framework to calculate the neutrino scattering and
absorption rates in a warm neutron-rich matter that consistently includes
mean-field effects and correlations through the Random Phase Approximation
(RPA). We employ nuclear interactions derived from chiral effective field
theory and, for the first time, include direct and exchange contributions to
the mean field energies and RPA vertex functions. The combination of Hartree-
Fock self energies and RPA vertex corrections to the response function
constitutes a “conserving approximation” that guarantees thermodynamic
consistency and the positivity of dynamic structure functions. The integral
equation for the RPA vertex function was discretized, leading to a matrix
eigenvalue problem that could be solved through standard diagonalization.

We find that including RPA correlations produces a broad collective mode that
shifts the strength of the charged-current neutrino response to higher energy
transfer. For a fixed incident neutrino energy, this lowers the outgoing
electron energy into a region of Pauli-blocked suppression, thereby reducing
the total cross-section. In contrast, electron-antineutrino absorption is
enhanced at low energies because the same broad collective mode lowers the
threshold energy needed to convert protons into neutrons. For higher
antineutrino energy, this enhancement is absent because of kinematic
constraints.

Figure 14: Imaginary part of the antineutrino absorption charged-current spin
response function $\bar{\chi}_{\tau\sigma}$ of beta-equilibrium nuclear matter
at a density $n=0.02$ fm-3 and temperature $T=8$ MeV. (Top panel) spin
response function including mean field “MF” corrections from the N3LO-414
chiral nuclear potential. (Bottom panel) spin response function including
random phase approximation vertex corrections plus nuclear mean fields
“RPA+MF”.

Our consistent treatment of nuclear mean fields and RPA correlations has
important implications for charged current reactions near the supernova and
neutron star merger neutrinospheres. When only mean-field effects are
included, we confirm the results of previous studies that found a large
enhancement of electron-neutrino absorption cross-section and a reduction in
antineutrino absorption cross-section. However, correlations included through
RPA qualitatively change the picture. As discussed in the previous section,
the absorption cross-sections for electron neutrinos are reduced over the
entire range of relevant energies. This would imply an increased luminosity
and average energy of electron neutrinos emitted in supernovae and mergers.
For anti-electron neutrinos, the absorption cross-section is enhanced at low
energy and suppressed at high energy. This implies that the lepton number flux
carried by neutrinos can be strongly energy-dependent when RPA correlations
are included. Since neutrino oscillations and nucleosynthesis are especially
sensitive to the neutrino lepton number flux and its energy dependence Qian
and Woosley [1996], Arcones and Thielemann [2012], our findings will likely
impact both.

Our finding that vertex corrections suppress differences between $\nu_{e}$ and
$\bar{\nu}_{e}$ absorption rates in the vicinity of the neutrinosphere implies
that the emerging spectra of $\nu_{e}$ and $\bar{\nu}_{e}$ will be more
similar than previously expected Martínez-Pinedo  _et al._ [2012], Roberts
_et al._ [2012]. This will likely inhibit the production of neutron-rich
r-process nuclei in the neutrino driven ejecta from core-collapse supernovae
and dynamical ejecta in neutron star mergers (for recent reviews see Kajino
_et al._ [2019], Cowan  _et al._ [2021]). However, collective neutrino flavor
oscillations, which are also sensitive to the energy dependence of the lepton
number flux, impact the final spectra at the nucleosynthesis site (see reviews
Chakraborty  _et al._ [2016], Tamborra and Shalgar [2021] and references
therein). To gauge if corrections to mean free paths calculated in this work
can alter nucleosynthesis it is critical to include their effect on flavor
transformation by modifying the two-point Feynman diagrams for electron
neutrino absorption and emission processes in Quantum Kinetic Equation (QKE)
treatments Vlasenko  _et al._ [2014], Richers  _et al._ [2019], and
collisional instabilities Johns and Xiong [2022], Johns [2023]. Further,
larger mean free paths for both electron and anti-electron neutrinos at higher
energy imply an increase in the total luminosity and average energy. This will
likely increase the net neutrino energy deposition behind the shock in core-
collapse supernova and aid the explosion. The results of this work can be
tabulated and included in simulations of core-collapse supernovae, neutron
star mergers, and neutron star cooling.

Finally, we mention the limitations of our study and identify directions for
future work. First, although our analysis provides a consistent treatment of
excitations above the Hartree-Fock ground state, it neglects two-body currents
and correlations beyond one-particle-one-hole RPA. At the low density, we
expect both two-body currents and two-particle-two-hole excitations to be
small because they appear at higher order in the density expansion, but more
work is needed to assess their importance at densities of relevance to the
neutrino sphere. It is well-known that the interplay between short-range
correlations and two-body currents, especially in the axial vector channel,
plays an important role in nuclear weak interactions. In addition, error
estimates for neutrino interaction rates require systematic order-by-order
calculations. These quantitive issues warrant further work before one can draw
definite conclusions about $\nu_{e}$ and $\bar{\nu}_{e}$ charged current
reactions in the neutrino sphere and their emergent spectra.

##  V Acknowledgement

The work of E. Shin and J. W. Holt is supported by the National Science
Foundation under Grant Nos. PHY1652199 and PHY2209318. Portions of the
research were conducted with the advanced computing resources provided by
Texas A&M High Performance Research Computing. The work of S. R. was supported
by the U.S. DOE under Grant No. DE-FG02-00ER41132 and by the National Science
Foundation’s Physics Frontier Center: The Network for Neutrinos, Nuclear
Astrophysics, and Symmetries.

## References

  * Burrows and Sawyer [1999] A. Burrows and R. F. Sawyer, Phys. Rev. C 59, 510 (1999).
  * O’Connor [2015] E. O’Connor, Astrophys. J. Suppl. 219, 24 (2015).
  * Melson  _et al._ [2015] T. Melson, H.-T. Janka, R. Bollig, F. Hanke, A. Marek, and B. Müller, Astrophys. J. Lett. 808, L42 (2015).
  * Roberts  _et al._ [2016] L. F. Roberts, C. D. Ott, R. Haas, E. P. O’Connor, P. Diener, and E. Schnetter, Astrophys. J. 831, 98 (2016).
  * Pons  _et al._ [1999] J. A. Pons, S. Reddy, M. Prakash, J. M. Lattimer, and J. A. Miralles, Astrophys. J. 513, 780 (1999).
  * Roberts and Reddy [2017] L. F. Roberts and S. Reddy, Phys. Rev. C 95, 045807 (2017).
  * Sekiguchi [2010] Y. Sekiguchi, Class. Quant. Grav. 27, 114107 (2010).
  * Wanajo  _et al._ [2014] S. Wanajo, Y. Sekiguchi, N. Nishimura, K. Kiuchi, K. Kyutoku, and M. Shibata, Astrophys. J. Lett. 789, L39 (2014).
  * Endrizzi  _et al._ [2020] A. Endrizzi, A. Perego, F. M. Fabbri, L. Branca, D. Radice, S. Bernuzzi, B. Giacomazzo, F. Pederiva, and A. Lovato, Eur. Phys. J. A 56, 15 (2020).
  * Sumiyoshi  _et al._ [2021] K. Sumiyoshi, S. Fujibayashi, Y. Sekiguchi, and M. Shibata, Astrophys. J. 907, 92 (2021).
  * Cusinato  _et al._ [2022] M. Cusinato, F. M. Guercilena, A. Perego, D. Logoteta, D. Radice, S. Bernuzzi, and S. Ansoldi, Eur. Phys. J. A 58, 99 (2022).
  * Sawyer [1989] R. F. Sawyer, Phys. Rev. C 40, 865 (1989).
  * Reddy  _et al._ [1998] S. Reddy, M. Prakash, and J. M. Lattimer, Phys. Rev. C 58, 013009 (1998).
  * Reddy  _et al._ [1999] S. Reddy, M. Prakash, J. M. Lattimer, and J. A. Pons, Phys. Rev. C 59, 2888 (1999).
  * Martínez-Pinedo  _et al._ [2012] G. Martínez-Pinedo, T. Fischer, A. Lohs, and L. Huther, Phys. Rev. Lett. 109, 251104 (2012).
  * Roberts  _et al._ [2012] L. F. Roberts, S. Reddy, and G. Shen, Phys. Rev. C 86, 065803 (2012).
  * Pastore  _et al._ [2012] A. Pastore, M. Martini, V. Buridon, D. Davesne, K. Bennaceur, and J. Meyer, Phys. Rev. C 86, 044308 (2012).
  * Iwamoto and Pethick [1982] N. Iwamoto and C. J. Pethick, Phys. Rev. D 25, 313 (1982).
  * Burrows and Sawyer [1998] A. Burrows and R. F. Sawyer, Phys. Rev. C 58, 554 (1998).
  * Horowitz and Schwenk [2006] C. J. Horowitz and A. Schwenk, Phys. Lett. B 642, 326 (2006).
  * Horowitz  _et al._ [2017] C. J. Horowitz, O. L. Caballero, Z. Lin, E. O’Connor, and A. Schwenk, Phys. Rev. C 95, 025801 (2017).
  * Bedaque  _et al._ [2018] P. F. Bedaque, S. Reddy, S. Sen, and N. C. Warrington, Phys. Rev. C 98, 015802 (2018).
  * Rrapaj  _et al._ [2015] E. Rrapaj, J. W. Holt, A. Bartl, S. Reddy, and A. Schwenk, Phys. Rev. C 91, 035806 (2015).
  * Bacca  _et al._ [2012] S. Bacca, K. Hally, M. Liebendorfer, A. Perego, C. J. Pethick, and A. Schwenk, Astrophys. J. 758, 34 (2012).
  * Bartl  _et al._ [2016] A. Bartl, R. Bollig, H.-T. Janka, and A. Schwenk, Phys. Rev. D 94, 083009 (2016).
  * Hanhart  _et al._ [2001] C. Hanhart, D. R. Phillips, and S. Reddy, Phys. Lett. B 499, 9 (2001).
  * Lim and Holt [2017] Y. Lim and J. W. Holt, Phys. Rev. C 95, 065805 (2017).
  * Carreau  _et al._ [2019] T. Carreau, F. Gulminelli, and J. Margueron, Eur. Phys. J. A 55, 188 (2019).
  * Gandolfi  _et al._ [2011] S. Gandolfi, J. Carlson, and S. C. Pieper, Phys. Rev. Lett. 106, 012501 (2011).
  * Buraczynski and Gezerlis [2016] M. Buraczynski and A. Gezerlis, Phys. Rev. Lett. 116, 152501 (2016).
  * Baym and Kadanoff [1961] G. Baym and L. P. Kadanoff, Phys. Rev. 124, 287 (1961).
  * Harris and Evans [1978] C. G. Harris and W. A. B. Evans, J. Phys. C: Solid State Phys. 11, 4447 (1978).
  * Entem and Machleidt [2003] D. R. Entem and R. Machleidt, Phys. Rev. C 68, 041001 (2003).
  * Coraggio  _et al._ [2007] L. Coraggio, A. Covello, A. Gargano, N. Itaco, D. R. Entem, T. T. S. Kuo, and R. Machleidt, Phys. Rev. C 75, 024311 (2007).
  * Marji  _et al._ [2013] E. Marji, A. Canul, Q. MacPherson, R. Winzer, C. Zeoli, D. R. Entem, and R. Machleidt, Phys. Rev. C 88, 054002 (2013).
  * Holt and Kaiser [2017] J. W. Holt and N. Kaiser, Phys. Rev. C 95, 034326 (2017).
  * Wellenhofer  _et al._ [2014] C. Wellenhofer, J. W. Holt, N. Kaiser, and W. Weise, Phys. Rev. C 89, 064009 (2014).
  * Whitehead  _et al._ [2021] T. R. Whitehead, Y. Lim, and J. W. Holt, Phys. Rev. Lett. 127, 182502 (2021).
  * Holt  _et al._ [2012] J. W. Holt, N. Kaiser, and W. Weise, Nucl. Phys. A 876, 61 (2012).
  * Holt  _et al._ [2013] J. W. Holt, N. Kaiser, G. A. Miller, and W. Weise, Phys. Rev. C 88, 024614 (2013).
  * Holt  _et al._ [2016] J. W. Holt, N. Kaiser, and G. A. Miller, Phys. Rev. C 93, 064603 (2016).
  * Rrapaj  _et al._ [2016] E. Rrapaj, A. Roggero, and J. W. Holt, Phys. Rev. C 93, 065801 (2016).
  * Qian and Woosley [1996] Y.-Z. Qian and S. E. Woosley, Astrophys. J. 471, 331 (1996).
  * Arcones and Thielemann [2012] A. Arcones and F.-K. Thielemann, J. Phys. G: Nucl. and Part. Phys. 40, 013201 (2012).
  * Kajino  _et al._ [2019] T. Kajino, W. Aoki, A. B. Balantekin, R. Diehl, M. A. Famiano, and G. J. Mathews, Prog. Part. Nucl. Phys. 107, 109 (2019).
  * Cowan  _et al._ [2021] J. J. Cowan, C. Sneden, J. E. Lawler, A. Aprahamian, M. Wiescher, K. Langanke, G. Martínez-Pinedo, and F.-K. Thielemann, Rev. Mod. Phys. 93, 15002 (2021).
  * Chakraborty  _et al._ [2016] S. Chakraborty, R. Hansen, I. Izaguirre, and G. Raffelt, Nucl. Phys. B 908, 366 (2016), neutrino Oscillations: Celebrating the Nobel Prize in Physics 2015\.
  * Tamborra and Shalgar [2021] I. Tamborra and S. Shalgar, Ann. Rev. Nucl. Part. Sci. 71, 165 (2021).
  * Vlasenko  _et al._ [2014] A. Vlasenko, G. M. Fuller, and V. Cirigliano, Phys. Rev. D 89, 105004 (2014).
  * Richers  _et al._ [2019] S. A. Richers, G. C. McLaughlin, J. P. Kneller, and A. Vlasenko, Phys. Rev. D 99, 123014 (2019).
  * Johns and Xiong [2022] L. Johns and Z. Xiong, Phys. Rev. D 106, 103029 (2022).
  * Johns [2023] L. Johns, Phys. Rev. Lett. 130, 191001 (2023).


# Mechanical feedback controls the emergence of dynamical memory in growing
tissue monolayers

Sumit Sinha1, Xin Li2, Rajsekhar Das2 and D. Thirumalai2 Department of
Physics, University of Texas at Austin, Austin, TX 78712, USA.  Department of
Chemistry, University of Texas at Austin, Austin, TX 78712, USA.

###### Abstract

The growth of a tissue, which depends on cell-cell interactions and
biologically relevant process such as cell division and apoptosis, is
regulated by a mechanical feedback mechanism. We account for these effects in
a minimal two-dimensional model in order to investigate the consequences of
mechanical feedback, which is controlled by a critical pressure, $p_{c}$. A
cell can only grow and divide if the pressure it experiences, due to
interaction with its neighbors, is less than $p_{c}$. Because temperature is
an irrelevant variable in the model, the cell dynamics is driven by self-
generated active forces (SGAFs) that are created by cell division. It is shown
that even in the absence of intercellular interactions, cells undergo
diffusive behavior. The SGAF-driven diffusion is indistinguishable from the
well-known dynamics of a free Brownian particle at a fixed finite temperature.
When intercellular interactions are taken into account, we find persistent
temporal correlations in the force-force autocorrelation function ($FAF$) that
extends over timescale of several cell division times. The time-dependence of
the $FAF$ reveals memory effects, which increases as $p_{c}$ increases. The
observed non-Markovian effects emerge due to the interplay of cell division
and mechanical feedback, and is inherently a non-equilibrium phenomenon.

##  I Introduction

Life around us, spanning a bewildering array of length and time scales, is
sustained through multicellular processes that are driven by non-equilibrium
events such as cell growth and cell division [1, 2]. Although known for a long
time [3], several recent experimental studies have emphasized that growth and
division in cell collectives are governed by local stresses that the cells
experience[4, 5, 6, 7]. A manifestation of coupling of growth and division to
local stress is the deviation of growth law of the cell collective from
exponential law [6]. These experiments suggest that there must exist
mechanical feedback between the local stress and cell division.

In a series of papers, we showed that single cells in a collective exhibit
anomalous dynamics due to local stress-dependent cell growth and division [8,
9, 10, 11, 12, 13, 14]. In the present study, we explore the dependence of
mechanical feedback, mediated by stress threshold $p_{c}$, on the dynamics of
single cells. A recent interesting study [15], has shown that mechanical
feedback regulates the physical properties of jammed cell collectives, which
supports experimental findings [4]. However, how the mechanical feedback
regulates individual cell migration in a collective is unknown, and is the
problem which we address in this study.

Using a two-dimensional off-lattice agent-based simulation model, we explore
the role of mechanical feedback ($p_{c}$) on single-cell dynamics. The central
results of the present study are: (1) In the absence of cell growth and
division, the dynamics is solely governed by short-ranged two body
interactions. In this limit, the cells in the long-time behave like a glass-
like solid. (2) In the presence of cell division, with systematic interactions
absent, the cells exhibit Markov dynamics resulting in diffusive motion at
long times. This finding is surprising because there is no thermal motion
(temperature is an irrelevant variable). This is different from a free
Brownian particle where temperature randomizes the particle motion. (3) When
both cell division and systematic interactions control the collective movement
in tandem, the dynamics is regulated by the mechanical feedback that is
parameterized using $p_{c}$. To quantify the dynamics, we calculated the force
auto-correlation function (FAF) inspired by works in the theory of chemical
reactions in liquids [16, 17, 18]. We show FAF increases when $p_{c}$ is
increased. The emergence of long time correlation in the FAF shows departure
from Markov dynamics, and is suggestive of memory effects in growing cell
collectives. (4) The persistence of trajectories of individual cells increases
as $p_{c}$ is increased. The trajectories are strikingly different from simple
Brownian motion. The enhanced persistence in cell dynamics, as $p_{c}$
increases, is the origin of memory in active systems. Taken together, the
present study establishes how mechanical feedback coupled with cell growth and
division leads to non-Markovian cell dynamics whose importance has not been
appreciated before.

Figure 1: Schematic of the 2D model.(a) Total force (Eq. 3) as a function of
inter-cellular distance for two cells, $i$ and $j$, with radius
$R_{i}=R_{j}=4\leavevmode\nobreak\ \mu m$. The repulsive and attractive part
of the force are given by Eqs. 1 and 2, respectively. The inset is the zoomed
in view that highlights the region in which the force is predominantly
attractive. (b) Cartoon illustrating the role of mechanical feedback. On the
left, the ‘red’ cell is dormant (cannot grow and divide) because the pressure
exerted by the neighbors exceeds $p_{c}$. The ‘green’ cell is in the growth
phase ($G$), which grows ands divides ($p<p_{c}$). The green cell from the
left gives given birth to two daughter cells (orange and cyan) when the radius
exceeds the mitotic radius $R_{m}$.

##  II Methods

We simulated the spatial and temporal dynamics of a two dimensional (2D)
growing tissue using agent-based off-lattice model in which the cells are
represented as interacting deformable disks. This simplified assumption of
representing cells as deformable disks was also used in previous studies [19],
although the details differ. In the model, the cells grow stochastically in
time and divide upon reaching a critical size ($R_{m}$), the mitotic radius.
The cell-to-cell interaction is characterized by elastic and adhesive forces.
We also consider cell-to-substrate damping as a way of accounting for the
effects of friction experienced by a moving cell by the substrate.

Physical Interactions : Each cell is modeled as a deformable disk with a time
dependent radius. A cell is characterized by physical properties such as the
radius, elastic modulus, membrane receptor and ligand concentration. In
addition, the cells attract each other through E-Cadherin mediated adhesive
interactions. This model is inspired by previous works on 3D off-lattice
multicellular tumor growth models [20, 21, 8, 9, 11, 12, 10, 13]. The elastic
(repulsive) force between two disks with radii $R_{i}$ and $R_{j}$ is given
by,

$F_{ij}^{el}(t)=\frac{h_{ij}^{3/2}(t)}{\frac{3}{4}(\frac{1-\nu_{i}^{2}}{E_{i}}+\frac{1-\nu_{j}^{2}}{E_{j}})\sqrt{\frac{1}{R_{i}(t)}+\frac{1}{R_{j}(t)}}},$
(1)

where $E_{i}$ and $\nu_{i}$, respectively, are the elastic modulus and Poisson
ratio of cell $i$. The overlap between the disks, if they interpenetrate
without deformation, is $h_{ij}$, which is given by
$\mathrm{max}[0,R_{i}+R_{j}-|\vec{r}_{i}-\vec{r}_{j}|]$ with
$|\vec{r}_{i}-\vec{r}_{j}|$ being the center-to-center distance between the
two disks.

Cell adhesion, mediated by receptors on the cell membrane, is the process by
which cells can attach to one another. For simplicity, we assume that the
receptor and ligand molecules are evenly distributed on the cell surface.
Consequently, the magnitude of the adhesive force, $F_{ij}^{ad}$, between two
cells $i$ and $j$ is expected to scale as a function of their contact line-
segment, $L_{ij}$. Keeping the 3D model as a guide [8], we calculate
$F_{ij}^{ad}$ using,

$F_{ij}^{ad}=L_{ij}f^{ad}\frac{1}{2}(c_{i}^{rec}c_{j}^{lig}+c_{j}^{rec}c_{i}^{lig}),$
(2)

where the $c_{i}^{rec}$ ($c_{i}^{lig}$) is the receptor (ligand) concentration
(assumed to be normalized with respect to the maximum receptor or ligand
concentration so that $0\leq c_{i}^{rec},c_{i}^{lig}\leq 1$). The coupling
constant $f^{ad}$ allows us to rescale the adhesion force to account for the
variabilities in the maximum densities of the receptor and ligand
concentrations. We calculate the contact length, $L_{ij}$, using the length of
contact between two intersecting circles,
$L_{ij}=\frac{\sqrt{(|4r_{ij}^{2}R_{i}^{2}-(r_{ij}^{2}-R_{j}^{2}+R_{i}^{2})^{2}|)}}{r_{ij}}$.
Here, $r_{ij}$ is the distance between cells $i$ and $j$. As before, $R_{i}$
and $R_{j}$ denote the radius of cell $i$ and $j$.

Repulsive and adhesive forces considered in Eqs.(1) and (2) act along the unit
vector ${\bf n}_{ij}$ pointing from the center of cell $j$ to the center of
cell $i$. The total force on the $i^{th}$ cell is given by the sum over its
nearest neighbors ($NN(i)$),

${\bf F}_{i}=\Sigma_{j\epsilon NN(i)}(F_{ij}^{el}-F_{ij}^{ad}){\bf n}_{ij}.$
(3)

The nearest neighbors satisfy the condition $R_{i}+R_{j}-|{\bf r}_{i}-{\bf
r}_{j}|\leavevmode\nobreak\ >\leavevmode\nobreak\ 0$. Figure 1a shows the plot
of the total force as a function of inter-cellular distance.

Three Scenarios: In order to elucidate the dramatically different dynamical
behavior, we consider three limits. (I) The collective movement arising solely
from the systematic forces, given in Eq. 3. (II) Cell movement with ${\bf
F}_{i}=0$ (no inter-cellular interactions) but allowing for cell division and
growth. Note that since inter-cellular interactions are absent, mechanical
feedback ($p_{c}$) does not play a role. In this limit, we show that the
dynamics can only arise due to active forces generated upon cell division. The
limits (I) and (II) are not relevant in describing collective movements in
Multicellular Spheroids (MCSs) [22] or evolving cell monolayers [6]. (III) In
this limit, we not only include interactions between cells (Eq. 3) but also
allow for cell growth, division and apoptosis. Most importantly, the time-
dependent growth of the tissue colony is limited by mechanical feedback, which
prohibits the biologically important process of cell growth and division if
the local stress on a cell exceeds a critical non-zero value, $p_{c}$.

Equation of Motion: The damped dynamics of the $i^{th}$ cell is computed based
on the equation of motion,

$\dot{{\bf r}}_{i}=\frac{{\bf F}_{i}}{\gamma_{i}}.$ (4)

Here, $\gamma_{i}$ is the friction coefficient of the $i^{th}$ cell. We assume
$\gamma_{i}$ to be equal to $\gamma_{o}R_{i}(t)$, where $\gamma_{o}$ is a
constant. This form of $\gamma_{i}$ is inspired from the simulations of three
dimensional models for solid tumor where $\gamma_{i}=6\pi\eta R_{i}$ with
$\eta$ being the viscosity. Note, we do not consider the effect of temperature
(set to zero in the simulations) as we assume the friction coefficient, that
in reality arises from the extracellular matrix in 3D or substrate in 2D, to
be so high [19] that thermal motion is irrelevant. The equation of motion in
Eq. 4 is similar to the case for soft granular materials where the role of
temperature is neglected [23]. However, it is crucial to note that in the
growth of the tissue colony, scenarios II and III in our case, there is a
self-generated active force (SGAF) that arises due to the biologically
important processes of cell growth and division [10].

Cell growth, division and apoptosis: In our model, cells can be either in the
dormant ($D$) or in the growth ($G$) phase depending on the local pressure
associated with a cell (Figure 1(b)). Using Irving-Kirkwood definition, we
track the pressure ($p_{i}$) experienced by the $i^{th}$ cell due to contact
with its neighbors [24]. The expression for $p_{i}$ is given by,

$p_{i}=\frac{1}{2}\Sigma_{j\epsilon NN(i)}\frac{{\bf F}_{ij}\cdot{\bf
dr}_{ij}}{A_{i}},$ (5)

where $A_{i}=\pi R_{i}^{2}$, is the area of the cell. If the local pressure,
$p_{i}$, exceeds a critical limit ($p_{c}$) the cell stops growing and enters
the dormant phase. Note, the cell can switch back to a growing phase if
$\frac{p_{i}}{p_{c}}<1$ as the tissue evolves. The critical pressure $p_{c}$,
serves as a mechanical feedback, which is known to regulate the growth of
tissues [7].

For growing cells ($\frac{p_{i}}{p_{c}}<1$), their area increases at a
constant rate $r_{A}$. The cell radius is updated from a Gaussian distribution
with the mean rate $\dot{R}=(2\pi R)^{-1}r_{A}$. Over the cell cycle time
$\tau$,

$r_{A}=\frac{\pi(R_{m})^{2}}{2\tau},$ (6)

where $R_{m}$ is the mitotic radius. The cell cycle time ($\tau$) is related
to the growth rate ($k_{b}$) by $\tau=\frac{ln\leavevmode\nobreak\ 2}{k_{b}}$.
A cell divides once it grows to the fixed mitotic radius ($R_{m}$). To ensure
area conservation, upon cell division, we use $R_{d}=R_{m}2^{-1/2}$ as the
radius of the daughter cells. The two resulting cells are placed at a center-
to-center distance $d=2R_{m}(1-2^{-1/2})$. The direction of the new cell
location is chosen randomly from a uniform distribution on the unit circle.
One source of stochasticity in the cell movement in our model is due to random
choice for the mitotic direction. In our simulations, the cells may undergo
apoptosis at the rate $k_{a}$. Throughout this work, the apoptosis rate was
fixed to $10^{-6}s^{-1}$. Table 1 depicts the parameters used in the
simulations.

Table 1: The parameters used in the simulation.

Parameters |  Values |  References
---|---|---
Timestep ($\Delta t$) |  10$\mathrm{s}$ |  This paper
Critical Radius for Division ($R_{m}$) |  5 $\mathrm{\mu m}$ |  [21, 8]
Friction coefficient ($\gamma_{o}$) |  0.1 $\mathrm{kg/(\mu m\leavevmode\nobreak\ s)}$ |  This paper
Cell Cycle Time ($\tau$) |  54000 $\mathrm{s}$ |  [25, 26, 27, 8]
Adhesive Coefficient ($f^{ad})$ |  $10^{-4}\mathrm{\mu N/\mu m}$ |  This paper
Mean Cell Elastic Modulus ($E_{i})$ |  $10^{-3}\mathrm{MPa}$ |  [28, 8]
Mean Cell Poisson Ratio ($\nu_{i}$) |  0.5 |  [21, 8]
Apoptosis Rate ($k_{a}$) |  $10^{-6}\mathrm{s^{-1}}$ |  [8]
Mean Receptor Concentration ($c^{rec}$) |  1.0 |  [8]
Mean Ligand Concentration ($c^{lig}$) |  1.0 |  [8]

Initial Conditions: We begin the simulations by placing 100 cells on a 2D
plane whose coordinates are chosen from a normal distribution with mean zero
and standard deviation $25\leavevmode\nobreak\ \mu m$. All the parameters
apart from critical pressure ($p_{c}$) are fixed.

Figure 2: Number and MSD of cells solely based on physical interactions and
cell division: (a) Number of cells, $N(t)$, solely based on physical two body
interactions. Since, cell division rate is zero, $N(t)$ is a constant. (b)
Mean-Squared Displacement, $\Delta(t)$, as a function of time when only
physical interactions are present. $\Delta(t)$ relaxes initially, governed by
elastic forces over a time-scale $\frac{\gamma}{ER_{m}}$ (red dashed vertical
line), and settles to a plateau value. (c) $N(t)$ as a function of scaled time
($\frac{t}{\tau}$) in the presence of cell division only. $N(t)$ grows
exponentially. (d) $\Delta(t)$ as a function of time in the presence of cell
division only. $\Delta(t)$ grows linearly in time and is diffusive (in black).
The red dashed line is a linear fit given by Eq. 13.

##  III Results

###  III.1 Markov dynamics in the presence of delta-correlated random force

For comparison, we briefly summarize the well-known result for a stochastic
process (over-damped Langevin equation) in one dimension for a free Brownian
particle. The equation of motion is,

$\frac{dx}{dt}=\sqrt{2D}\eta(t),$ (7)

where the random force obeys $\langle\eta(t)\rangle=0$ and
$\langle\eta(t)\eta(t^{\prime})\rangle=\delta(t-t^{\prime})$. The position of
the particle is $x$. The solution $P(x,t)$, for probability density for
finding the particle at $x$ at time $t$, is given by

$P(x,t)=\frac{1}{\sqrt{4\pi Dt}}e^{-\frac{x^{2}}{4Dt}}.$ (8)

Here, we have assumed that the particle was at the origin at $t=0$,
$P(x,0)=\delta(x)$. The moments of $P(x,t)$, which serve as the physical
observables in cell tracking experiments [22], are readily calculated. For
instance, the first moment $\langle x\rangle$ is given as

$\langle x(t)\rangle=\int_{-\infty}^{\infty}xP(x,t)=0.$ (9)

The second moment $\langle x^{2}\rangle$, also called the mean-squared
displacement (MSD), is non-zero and is given as

$\langle x^{2}(t)\rangle=\int_{-\infty}^{\infty}x^{2}P(x,t)=2Dt.$ (10)

The dynamics of a free Brownian particle is an example of a Markov process.
There is no memory because $P(\delta x(t))=\mathcal{N}(0,\sqrt{2D\delta t})$
is independent of $x(t)$, where $\delta x(t)=x(t+\delta t)-x(t)$ is the
displacement of the particle from time $t$ to $t+\delta t$ and
$\mathcal{N}(0,\sqrt{2D\delta t})$ is the normal distribution with zero mean
and variance $2D\delta t$. This example sets the stage for exploring the
emergent force auto-correlation ($FAF$), with long temporal correlation, in an
expanding tissue.

###  III.2 Role of physical interactions and cell division

We first investigate the consequences of physical interactions and cell
division when they are not coupled to one another.

(a) Physical Interactions (limiting case I): When the interactions between
cells are based only on the systematic interactions, as given in Eqs. 1 and 2
without cell division and apoptosis, the dynamics of the interacting cells is
governed by the elastic timescale $\frac{\gamma}{ER_{m}}$. In the absence of
cell growth, division and apoptosis, the number of cells, $N(t)$, is a
constant, which is confirmed in Figure 2a. Figure 2b shows the plot of mean-
square displacement, $\Delta(t)$, defined as

$\Delta(t)=\bigg{\langle}\frac{1}{N}\sum_{i=1}^{N}\big{(}{\bf r}_{i}(t)-{\bf
r}_{i}(0)\big{)}^{2}\bigg{\rangle},$ (11)

where $\langle...\rangle$ denotes the ensemble average over $20$ simulation
runs and $N$ is the initial $100$ cells. Figure 2a shows that $\Delta(t)$
relaxes rapidly to a plateau on a time scale $\approx\frac{\gamma}{ER_{m}}$.
Usually, $\Delta(t)\sim t^{\alpha}$, with $\alpha=0$, is indicative of solid-
like behavior. Figure 2b shows that in the long time limit,
$t>>\frac{\gamma}{ER_{m}}$, $\alpha=0$, and hence the cell collective behaves
as a solid in the sense there is absence of diffusion. Note, the dynamics is
performed under athermal (temperature is not relevant) open boundary
conditions and the scale of systematic interactions are short-ranged ($\approx
R_{m}$). Hence, the cells cannot move after the initial relaxation process.

(b) Effect of cell division (limiting case II): In this limit, during each
cell division event, a cell is displaced by the distance $\approx R_{m}$,
randomly in space. Hence, when the time evolution of cell colony is governed
solely by cell division (absence of apoptosis and systematic interactions) we
expect that with successive cell divisions a cell would undergo a random walk
that is uncorrelated in time and space. In other words, it would behave as a
Brownian particle due to the SGAF induced by cell division. This is purely a
non-equilibrium dynamical process. As a result, we expect that
$\Delta(t)=D_{eff}t$ at long times. Because in this limiting case $\tau$ is
the only time scale and $R_{m}$ is the only length scale, we obtain
$D_{eff}=\frac{R_{m}^{2}}{\tau}$. In the absence of systematic interactions
the pressure (Eqn. 5) on an individual cell is zero, all the cells grow and
divide independently. Hence, the number of cells, $N(t)$ increases
exponentially (see Figure 2c). In this limit,

$\frac{dN(t)}{dt}=k_{b}N(t)$ (12)

and therefore $N(t)=N_{o}e^{k_{b}t}$. Interestingly, in accord with the
arguments given above, the dynamics of individual cells is diffusive in this
limit and the mean squared displacement $\Delta(t)$ is given by

$\Delta(t)=\frac{R_{m}^{2}}{\tau}t,$ (13)

as shown in Figure 2d. Although there is no thermal motion, the scaling
behavior of $\Delta(t)$, is similar to the standard Brownian dynamics given by
Eqn. 10. The dynamics is diffusive because during every cell division the cell
is displaced by distance $R_{m}$ randomly, and hence mimics a Brownian motion
in two dimensions. Note, when only cell apoptosis is present ($k_{b}=0$), the
cells do not move, and $N(t)$ decreases exponentially to zero at the rate
$k_{a}$.

The dynamics is non-trivial when all the components (systematic forces, cell
division and apoptosis, and mechanical feedback) are included. For this case,
the growing tissue develops a core where the cells are jammed, and exhibit
glass-like dynamics [11]. In contrast, the cells in the periphery are
predominantly in the growth ($G$) phase. As a result, the cells exhibit
anomalous spatially heterogeneous dynamics with super-diffusive (sub-
diffusive) periphery (core) [11]. Our previous work has shown that this
dynamical phase separation arises due to SGAFs that arise due to local stress-
regulated cell growth and division [10]. In the model, cell division and
growth are regulated by the mechanical feedback parameter $p_{c}$, the
consequences of which is explored in the following section.

Figure 3: Emergence of highly correlated force: Plot of force auto-correlation
function (FAF) as a function of time. From top to bottom, FAF corresponds to
$p_{c}=10^{-3},10^{-4}$ and $10^{-5}$. The dashed line corresponds to the
fits. In the inset, we zoom in on the initial time regime of FAF. The order of
the plots and dashed lines are same as in the main figure. The figure shows
the emergence of FAF with two time scales: one long
($\sim\frac{1}{k_{b}-k_{a}}$) and one small (elastic time scale
$=\frac{\gamma}{ER_{m}}$).

###  III.3 Highly correlated force correlations in an expanding tissue

For a simple Brownian motion, the force is delta-correlated and hence the
dynamics is Markovian (see Eqn. 7). Therefore, a signature of such dynamics is
the fast decay of force auto-correlation function ($FAF$), in comparison to
the smallest time-scale of the problem, which in the present case is
$\frac{\gamma}{ER_{m}}$. To explore the nature of the dynamics, when both cell
division and systematic interactions are present in tandem, we calculated
$FAF$ ($t^{*}$) given as

$FAF(t^{*})=\frac{\langle{\bf F}(t+t^{*})\cdot{\bf
F}(t)\rangle_{t}}{\langle{\bf F}(t)\cdot{\bf F}(t)\rangle_{t}}$ (14)

Here, ${\bf F}(t)$ is the force on the cell at time $t$ and
$\langle...\rangle_{t}$ is the time average, and $t^{*}$ is the delay (or
waiting) time. The time averaging is performed over $\approx 2,000$ cells.
Figure 3 shows the plot of FAF for $p_{c}=10^{-3}Nm^{-1},10^{-4}Nm^{-1}$ and
$10^{-5}Nm^{-1}$. Figure 3 shows that the FAF decays on two time scales: long
($\sim\frac{1}{k_{b}-k_{a}}$) and short ($\frac{\gamma}{ER_{m}}$). In order to
extract the two time scales, we fit FAF with $Ae^{\frac{-t^{*}}{\tau_{c}}}+B$
in both the regimes.

In the short time regime (see the inset of Figure 3), for
$p_{c}=10^{-3}Nm^{-1}$, $A=0.5,\tau_{c}=\frac{1.2\gamma}{ER_{m}}$ and
$B=0.41$. For $p_{c}=10^{-4}Nm^{-1}$,
$A=0.75,\tau_{c}=\frac{0.97\gamma}{ER_{m}}$ and $B=0.16$. Lastly, for
$p_{c}=10^{-5}Nm^{-1}$, $A=0.81,\tau_{c}=\frac{0.95\gamma}{ER_{m}}$ and
$B=0.11$. As anticipated, we find that in the short time regime, the
relaxation time is approximately close to the elastic time scale
$\frac{\gamma}{ER_{m}}$, which is negligible compared to
$\frac{1}{k_{b}-k_{a}}$. However, in the long time limit, FAF shows memory
effects, especially for $p_{c}=10^{-3}Nm^{-1}$. For $p_{c}=10^{-3}Nm^{-1}$,
$A=0.41,\tau_{c}=\frac{2.2}{k_{b}-k_{a}}$ and $C=-0.06$. For $p_{c}=10^{-4}$,
$A=0.12,\tau_{c}=\frac{2.3}{k_{b}-k_{a}}$ and $B=-0.02$. Lastly, for
$p_{c}=10^{-5}Nm^{-1}$, $A=0.04,\tau_{c}=\frac{0.2}{k_{b}-k_{a}}$ and
$B\approx 0$. For $p_{c}=10^{-5}Nm^{-1}$, A is negligible which indicates the
absence of memory. In the other two cases, A for $p_{c}=10^{-3}Nm^{-1}$ is
four times larger than for $p_{c}=10^{-4}Nm^{-1}$. Larger magnitude of FAF in
the long time regime leads to higher degree of migration for
$p_{c}=10^{-3}Nm^{-1}$. The emergence of highly correlated forces elucidates
the departure from Markovian dynamics in a system comprising of many
interacting cells whose time evolution occurs under non-equilibrium
conditions, and absence of fluctuation-dissipation theorem [14].

Figure 4: Trajectories of cells for different $p_{c}$ values: (a, b, c)
Trajectories of cells for $p_{c}=10^{-5}\leavevmode\nobreak\ N/m$,
$p_{c}=10^{-4}\leavevmode\nobreak\ N/m$ and $p_{c}=10^{-3}\leavevmode\nobreak\
N/m$. For all the three panels, the color bar represents the time in units of
$(k_{b}-k_{a})t$. From the figures, it is clear that cells migration is
enhanced when $p_{c}$ is increased.

###  III.4 Persistence of trajectories increases on increasing $p_{c}$

To better understand, the significance of memory effects that are embodied in
FAF, we probed the trajectories of individual cells in the growing tissue
colony. We investigated the three cases: $p_{c}=10^{-5}\leavevmode\nobreak\
N/m,10^{-4}\leavevmode\nobreak\ N/m$ and $10^{-3}\leavevmode\nobreak\ N/m$ in
Figure 4. Note, the trajectories are for the initial $100$ cells with the
identical initial conditions and the total simulation time. Thus, the
differences in the trajectories emerge due to interplay of systematic
interactions, cell growth and division. Figure 4 shows that the nature of the
trajectories are strikingly different when $p_{c}$ values are changed. For
$p_{c}=10^{-5}\leavevmode\nobreak\ N/m$ in Figure 4a, the cells move the
smallest compared to $p_{c}=10^{-4}\leavevmode\nobreak\ N/m$ and
$p_{c}=10^{-3}\leavevmode\nobreak\ N/m$ and are not persistent in nature. This
is because the the memory effect in the $FAF$ is negligible, as illustrated in
Figure 3. The cells exhibit persistent trajectories for
$p_{c}=10^{-4}\leavevmode\nobreak\ N/m$ and $p_{c}=10^{-3}\leavevmode\nobreak\
N/m$, with the degree of persistence being higher for the latter. The
emergence of persistence in trajectories implies the presence of memory or
(non-Markovian) effects in the dynamics.

##  IV Conclusion

We have shown that in an evolving tissue colony, the dynamics of an individual
cell may be approximated as a stochastic process where the force is correlated
over many cell division times. The emergence of correlation in force is a
manifestation of memory effects, which is a hallmark of non-Markovian
dynamics. Memory effects in a growing tissue arises due to cell division and
mechanical feedback. There is no potential or energy function whose derivative
results in these processes. This immediately implies that there is no
equilibrium in the growing tissue, and hence the system is always out of
equilibrium. The present study also provides a mechanism for persistent motion
observed in many out of equilibrium active matter systems [29, 30, 31].

It is tempting to describe the simulated time-dependent force autocorrelation
using a reduced description, something like the Generalized Langevin equation
(GLE), to describe the effective dynamics of a cell in the evolving tissue. It
is natural to use such an approach in thermally controlled barrier crossing
problems, as was done decades ago in the most insightful studies [16, 17]. For
reasons stated above, construction of a similar set of equations, if it exists
at all, in any reduced variable (the analogue of reaction coordinate in
barrier crossing problems), which must also include large spatial
heterogeneity, is likely to be difficult. Cell-division and apoptosis require
energy input and dissipation, which can only be described using a physical
picture and a framework that goes beyond the usual description based on
Hamiltonians or energy functions. The work here may provide an impetus to
develop a general theoretical framework for describing feedback controlled
dynamics in active systems.

## Acknowledgement

We would like to thank Abdul N. Malmi-Kakkada and Himadri S. Samanta for
valuable comments on the manuscript. This work was supported by grants from
National Science Foundation (PHY 17-08128, PHY-1522550). Additional support
was provided by the Collie-Welch Reagents Chair (F-0019).

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] Ricard Alert and Xavier Trepat.  Physical models of collective cell migration.  Annual Review of Condensed Matter Physics, 11:77–101, 2020.
  * [2] M Reza Shaebani, Adam Wysocki, Roland G Winkler, Gerhard Gompper, and Heiko Rieger.  Computational models for active matter.  Nature Reviews Physics, 2(4):181–199, 2020.
  * [3] G Hemlinger, PA Noetti, HC Lichtenbeld, RJ Melder, and RK Jain.  Solid stress inhibits the growth of multicellular tumor spheroids.  Nat. Biotech., 15:778–783, 1997.
  * [4] Morgan Delarue, Jörn Hartung, Carl Schreck, Pawel Gniewek, Lucy Hu, Stephan Herminghaus, and Oskar Hallatschek.  Self-driven jamming in growing microbial populations.  Nature physics, 12(8):762–766, 2016.
  * [5] M. E. Dolega, M. Delarue, F. Ingremeau, J. Prost, A. Delon, and G. Cappello.  Cell-like pressure sensors reveal increase of mechanical stress towards the core of multicellular spheroids under compression.  Nature Communications, 8:14056 EP –, 01 2017.
  * [6] Alberto Puliafito, Lars Hufnagel, Pierre Neveu, Sebastian Streichan, Alex Sigal, D Kuchnir Fygenson, and Boris I Shraiman.  Collective and single cell behavior in epithelial contact inhibition.  Proceedings of the National Academy of Sciences, 109(3):739–744, 2012.
  * [7] Boris I Shraiman.  Mechanical feedback as a possible regulator of tissue growth.  Proceedings of the National Academy of Sciences, 102(9):3318–3323, 2005.
  * [8] Abdul N Malmi-Kakkada, Xin Li, Himadri S Samanta, Sumit Sinha, and Dave Thirumalai.  Cell growth rate dictates the onset of glass to fluidlike transition and long time superdiffusion in an evolving cell colony.  Physical Review X, 8(2):021025, 2018.
  * [9] Abdul Malmi-Kakkada, Xin Li, Sumit Sinha, and D Thirumalai.  Dual role of cell-cell adhesion in tumor suppression and proliferation.  arXiv:1906.11292, 2019.
  * [10] Sumit Sinha and D Thirumalai.  Self-generated persistent random forces drive phase separation in growing tumors.  The Journal of Chemical Physics, 153(20):201101, 2020.
  * [11] Sumit Sinha, Abdul N Malmi-Kakkada, Xin Li, Himadri S Samanta, and D Thirumalai.  Spatially heterogeneous dynamics of cells in a growing tumor spheroid: Comparison between theory and experiments.  Soft Matter, 16(22):5294–5304, 2020.
  * [12] Himadri S Samanta, Sumit Sinha, and D Thirumalai.  Far from equilibrium dynamics of tracer particles embedded in a growing multicellular spheroid.  arXiv:2003.12941, 2020.
  * [13] Sumit Sinha and Abdul N Malmi-Kakkada.  Inter-particle adhesion regulates the surface roughness of growing dense three-dimensional active particle aggregates.  Journal of Physical Chemistry B, 125(37):10445–10451, 2021.
  * [14] Himadri S Samanta and D Thirumalai.  Origin of superdiffusive behavior in a class of nonequilibrium systems.  Physical Review E, 99(3):032401, 2019.
  * [15] Pawel Gniewek, Carl F Schreck, and Oskar Hallatschek.  Biomechanical feedback strengthens jammed cellular packings.  Physical review letters, 122(20):208102, 2019.
  * [16] John E Straub, Michal Borkovec, and Bruce J Berne.  Calculation of dynamic friction on intramolecular degrees of freedom.  Journal of Physical Chemistry, 91(19):4995–4998, 1987.
  * [17] JE Straub, M Borkovec, and BJ Berne.  Molecular dynamics study of an isomerizing diatomic in a lennard-jones fluid.  J. Chem. Phys., 89:4833–4847, 1988.
  * [18] BJ Berne and GD Harp.  On the calculation of time correlation functions.  Advances in chemical physics, 17:63, 1970.
  * [19] DA Matoz-Fernandez, Kirsten Martens, Rastko Sknepnek, JL Barrat, and Silke Henkes.  Cell division and death inhibit glassy behaviour of confluent tissues.  Soft matter, 13(17):3205–3212, 2017.
  * [20] Dirk Drasdo and Stefan Höhme.  A single-cell-based model of tumor growth in vitro: monolayers and spheroids.  Physical biology, 2(3):133, 2005.
  * [21] Gernot Schaller and Michael Meyer-Hermann.  Multicellular tumor spheroid in an off-lattice voronoi-delaunay cell model.  Physical Review E, 71(5):051910, 2005.
  * [22] Angela M Jimenez Valencia, Pei-Hsun Wu, Osman N Yogurtcu, Pranay Rao, Josh DiGiacomo, Inês Godet, Lijuan He, Meng-Horng Lee, Daniele Gilkes, Sean X Sun, et al.  Collective cancer cell invasion induced by coordinated contractile stresses.  Oncotarget, 6(41):43438, 2015.
  * [23] Rahul N Chacko, Peter Sollich, and Suzanne M Fielding.  Slow coarsening in jammed athermal soft particle suspensions.  Physical review letters, 123(10):108001, 2019.
  * [24] Xingbo Yang, M Lisa Manning, and M Cristina Marchetti.  Aggregation and segregation of confined active particles.  Soft matter, 10(34):6477–6484, 2014.
  * [25] James P Freyer and Robert M Sutherland.  Regulation of growth saturation and development of necrosis in emt6/ro multicellular spheroids by the glucose and oxygen supply.  Cancer research, 46(7):3504–3512, 1986.
  * [26] Joseph J Casciari, Stratis V Sotirchos, and Robert M Sutherland.  Variations in tumor cell growth rates and metabolism with oxygen concentration, glucose concentration, and extracellular ph.  Journal of cellular physiology, 151(2):386–394, 1992.
  * [27] Jacques Landry, James P Freyer, and Robert M Sutherland.  Shedding of mitotic cells from the surface of multicell spheroids during growth.  Journal of cellular physiology, 106(1):23–32, 1981.
  * [28] Jörg Galle, Markus Loeffler, and Dirk Drasdo.  Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro.  Biophysical journal, 88(1):62–75, 2005.
  * [29] Gerhard Gompper, Roland G Winkler, Thomas Speck, Alexandre Solon, Cesare Nardini, Fernando Peruani, Hartmut Löwen, Ramin Golestanian, U Benjamin Kaupp, Luis Alvarez, et al.  The 2020 motile active matter roadmap.  Journal of Physics: Condensed Matter, 32(19):193001, 2020.
  * [30] Clemens Bechinger, Roberto Di Leonardo, Hartmut Löwen, Charles Reichhardt, Giorgio Volpe, and Giovanni Volpe.  Active particles in complex and crowded environments.  Reviews of Modern Physics, 88(4):045006, 2016.
  * [31] M.C. Marchetti, J. F. Joanny, S. Ramaswamy, T.B. Liverpool, J. Prost, M. Rao, and R. Aditi Simha.  Hydrodynamics of soft active matter.  Rev. Mod. Phys., 85:1143, 2013.