text:
\section{Introduction}
\label{sec:introduction}

What boundary conditions may be consistently imposed on gauge and gravitational
fields, both asymptotically or at finite distances, is an important problem that
has recently attracted a lot of attention because it is directly related not
only to asymptotic and residual gauge symmetries but also to the number and the
nature of the degrees of freedom of the system.

In his report on quantum fields in curved space (\cite{DeWitt:1975ys}, section
2.4), DeWitt explains that we know ``from years of experiment and years of model
building'' what a conductor is and what consistent boundary conditions to impose
for electric and magnetic fields. He then goes to infer meaningful boundary
conditions and boundary terms for the case of a massless scalar field.

In this work, we follow the spirit of DeWitt's approach and impose the analog of
``perfectly conducting'' boundary conditions on free massless spin 2 fields. We
will show that the partition function of gravitons in such a Casimir box can be
computed in a straightforward way and is identical to that of photons in the
original Casimir set-up, and also to that of a massless scalar field with
periodic boundary conditions in a box twice as big. This is not entirely
surprising since these are results on the level of free field theories.

At this stage, we will limit ourselves to deriving the exact analytic and
modular covariant result that allows one to access the qualitatively very
different high and low temperature expansions, in agreement with recent results
on a Cardy formula and modular properties in higher dimensions
\cite{Cappelli:1988vw,Cardy:1991kr,Shaghoulian:2015kta}. We will not address
whether one may in principle confine gravitons to a box so as to achieve thermal
equilibrium
\cite{Smolin1984,garfinkle1985possibility,Dell1987,Padmanabhan_2003}. Neither
will we speculate on the physical meaning of perfect conductors for gravitons
nor comment on the relation to a recent study of the gravitational Casimir
effect at zero temperature with non-idealized boundary conditions
\cite{Quach:2015qwa}. As compared to the closely related computations in
\cite{Viaggiu:2017zax}, where a spherical box as a model for a Schwarzschild
black hole is considered, the use of a slab geometry simplifies a detailed
discussion of the contribution of various degrees of freedom to the partition
function and thus to the entropy, in much the same way than such a geometry
makes computations and conceptual issues relatively straightforward in the
original Casimir analysis \cite{casimir1948attraction}, whereas the spherically
symmetric case \cite{Boyer:1968uf,doi:10.1063/1.1665021} is much more involved.
Finally, we will not discuss here of how such computations in linearized gravity
could be relevant for black hole physics in full non-linear general relativity
or compare to results in the context of cosmological models
\cite{Brandenberger:1992jh}.

The paper is organized as follows. As a useful preparation, we re-discuss in
appendix \ref{sec:reduced-phase-space-2}, our analysis
\cite{Barnich:2018zdg,Barnich2019} of photons with Casimir boundary conditions
by using polarization vectors adapted to $E$ and $H$ modes. We also briefly
review in appendix \ref{sec:reduced-phase-space} the well-known reduced phase
space of linearized gravity in empty space. In the main text, we start in
section \ref{sec:casim-bound-cond}, by defining Casimir-type boundary conditions
with the associated mode expansions for metric perturbations and their canonical
momenta. This allows us in section \ref{sec:reduced-phase-space-1} to work out
the physical degrees of freedom in this precise set-up and to isolate the sector
with vanishing transverse momentum, or in other words, the sector of
gravitational waves propagating parallel to the plates. The next step in section
\ref{sec:scal-field-form} consists in organizing all physical degrees of
freedom, first in terms of two massless scalar fields with Dirichlet and Neumann
boundary conditions, respectively and then in terms of a single scalar field. In
order to do so, we follow recent work on the finite temperature electromagnetic
Casimir effect and modular invariance in this context \cite{Alessio:2020okv} to
show that this scalar field should satisfy periodic boundary on an interval of
double the length of the separation of the plates. This allows us in section
\ref{sec:extend-part-funct} to simply infer the exact result for the partition
function for gravitons, with a chemical potential for suitably modified spin
angular momentum turned on, from the well-studied case of a massless scalar
field with linear momentum in the compact direction turned on
\cite{Cappelli:1988vw,Cardy:1991kr,Shaghoulian:2015kta,Alessio2020}, together
with modular properties and consistent high and low temperature expansions.

Arguably the most interesting aspect of the exact analytic result is the fact
that the low temperature/small box expansion is qualitatively very different
from the more standard high temperature/large volume expansion. The latter is
dominated by the standard black body result, up to finite size corrections. At
low temperature however, the leading contribution to the partition function
scales like inverse temperature and is responsible for the zero temperature
gravitational Casimir force. It does not contribute to the entropy. The
sub-leading contribution to the partition function is then the leading
contribution to the entropy. It is proportional to the area of the plates and
our analysis shows that it originates from gravitons that propagate parallel to
those plates. The dynamics of those gravitons is equivalent to that of photons
in two spatial dimensions, which in turn is described by a massless scalar field
in 2 spatial dimensions. Let us stress here that these degrees of freedom are
due to the non-trivial boundary conditions imposed on the fields but, except for
the particle present in the spectrum, there is ultimately no interpretation in
terms of residual gauge degrees of freedom nor a direct relation to large gauge
transformations, contrary to what is implied in
\cite{Barnich:2018zdg,Barnich2019}.

\section{Casimir boundary conditions and mode expansions}
\label{sec:casim-bound-cond}

We want to promote the empty space analysis of massless spin 2 fields in the
Hamiltonian approach, as reviewed in appendix \ref{sec:reduced-phase-space}, to
a slab geometry with two infinite parallel planes at $x^3=-\frac{d}{2}$ and
$x^3=\frac{d}{2}$. Let $a=1,2$, $i=a,3$, $V=L_1L_2d$ the volume of the system
with $L_a$ taken large in the end, and $d$ the separation of the plates. Let
$k_3=\frac{\pi n_3}{d}$, $k_a=\frac{2\pi n_a}{L_a}$ (with no summation over $a$).

In analogy with the electromagnetic case reviewed in appendix
\ref{sec:reduced-phase-space-2}, we define Casimir-type boundary conditions by
requiring that $h_{ab},\pi^{ab}$ and $h_{33},\pi^{33}$ satisfy Dirichlet
conditions, while $h_{a3},\pi^{a3}$ satisfy Neumann conditions. More explicitly,
at the boundary $x^3=-\frac{d}{2}$ and $x^3=\frac{d}{2}$, we require
\begin{equation}
  \label{eq:78}
  h_{ab}=0=\pi^{ab},\quad h_{33}=0=\pi^{33},
  \quad (\partial_3 h_{a3})=0=(\partial_3 \pi^{a3}).
\end{equation}
Let us use $T^{ij}$ for either $h_{ij}$ or $\pi^{ij}$ and define, as in
\eqref{eq:14}, a complete set of scalar functions $\psi^H_k$ with $n_a,n_3>0$
adapted to Dirichlet boundary conditions, and $\psi^E_k$ with $n_a,n_3\geq 0$
adapted to Neumann boundary conditions.

The boundary conditions on the canonically conjugate variables are thus
implemented through the mode expansions \footnote{The factor $(-i)$ in front of
  the expansion of $(h_{a3},\pi^{a3})$ is chosen for later convenience.}
\begin{equation}
  \label{eq:2a}
  T^{ab}(x)=\sum_{n_a,n_3>0}T^{ab}_{k}\psi^H_k,\quad
  T^{33}(x)=\sum_{n_a,n_3>0}T^{33}_{k}\psi^H_{k},\quad
  T^{a3}(x)=(-i)\sum_{n_a,n_3\geq 0}T^{a3}_{k}\psi^E_{k},
\end{equation}
where we take $T^{ab}_{k_a,0}=0=T^{33}_{k_a,0}$. Reality conditions
are
\begin{equation}
  \label{eq:3a}
  T^{ab}_{k_a,k_3}=T^{*ab}_{-k_a,k_3},\quad T^{33}_{k_a,k_3}=T^{*33}_{-k_a,k_3},\quad
  T^{a3}_{k_a,k_3}=-T^{*a3}_{-k_a,k_3},
\end{equation}
while parity conditions are
\begin{equation}
  \label{eq:79}
  T^{ab}_{k_a,k_3}=-T^{ab}_{k_a,-k_3},\quad T^{33}_{k_a,k_3}=-T^{33}_{k_a,-k_3},\quad
  T^{a3}_{k_a,k_3}=T^{a3}_{k_a,-k_3}.
\end{equation}
Odd variables do not have modes at $n_i=0$, $T^{ab}_{0,0}=0=T^{33}_{0,0}$, but
there are two such modes in the even sector, $T^{a3}_{0,0}\neq 0$.

Note that these boundary conditions have been constructed in such a way that the
left hand sides of the constraints have well defined expansions, either in terms
of sines or in terms of cosines. If one requires the lapse and the shift to
impose all constraints in a first order variational principle and not to
introduce spurious variables, it follows from taking into account the mode
expansion of the constraints \eqref{eq:57}, that the lapse $n$ and $n^a$ should
satisfy Dirichlet while $n^3$ should satisfy Neumann conditions,
\begin{equation}
  \label{eq:103}
  n=\sum_{n_a,n_3>0}n_k\psi_k^H,\quad n^a=\sum_{n_a,n_3>0}n^a_k\psi_k^H,\quad
  n^3=\sideset{}{'}\sum_{n_a,n_3\geq 0}n^3_k\psi_k^E,
\end{equation}
so that $n,n^a$ are odd, while $n^3$ is even under $k_3\to -k_3$.

\section{Reduced phase space sectors}
\label{sec:reduced-phase-space-1}

In terms of these decomposition,
the kinetic term is given by
\begin{equation}
  \label{eq:81}
  \int \text{d}^3x\ \partial_0 h_{ij}\pi^{ij}=\sum_{n_a,n_3\geq 0}\partial_0
  h_{ij}\pi^{*ij},
\end{equation}
where we sometimes omit the $k$ dependence of the Fourier coefficients for
notational simplicity. The Pauli-Fierz Hamiltonian \eqref{eq:56} is given by
\begin{equation}
  \label{eq:80}
  H_{\rm PF}=\sum_{n_a,n_3\geq 0}\Big[\pi^{ij}\pi^*_{ij}-\frac{1}{2}\pi\pi^*
  +\frac{1}{4}k^2 h_{ij}h^{*ij}-\frac{1}{2}k_i h^{ij}k^lh_{lj}^*+\frac{1}{2}k^ik^jh h_{ij}^*
  -\frac{1}{4}k^2 hh^*\Big].
\end{equation}

\subsection{Particle sector}
\label{sec:particle-sector}

We start with the particle sector, which has to be treated separately since the
change of variables that we will use below does not apply to this sector. For
$n_i=0$, we define $h_{a3,0,0}=\sqrt{2}q^a$ and $\pi^{a3}_{0,0}=\frac{1}{\sqrt
  2}p_a$. From the expansion of the kinetic term, these variables are
canonically conjugate. None of these variables is involved in the constraints nor
affected by a gauge transformation generated by these constraints. Their
contribution to the Pauli-Fierz Hamiltonian is that of two particles of unit
mass,
\begin{equation}
  \label{eq:82}
  H_{n_i=0}=\frac{1}{2}p^ap_a.
\end{equation}

\subsection{Adapted polarization tensors}
\label{sec:adapt-polar-tens}

For $n_i\neq 0$, one may use the polarization tensors as defined in
\eqref{eq:63}-\eqref{eq:66} to decompose the components. Furthermore, it will
turn out to be convenient to use the polarization tensors built out of the
adapted polarization vectors ${e_\lambda}^i$, $\lambda=(H,E)$ given in
\eqref{eq:9} rather than to use generic transverse vectors ${e_\alpha}^i$.

As a consequence,
\begin{equation}
  \label{eq:96} \left\{ \begin{array}{ll} {e_{TT+}}^{ab}=\frac{1}{\sqrt
2k_\perp^2k^2}{(k^2\epsilon^{ac}k_c\epsilon^{bd}k_d-k^2_3k^ak^b)}, &
{e_{TT\times}}^{ab}=\frac{k_3}{\sqrt
2k_\perp^2k}{(\epsilon^{ac}k_ck^b+k^a\epsilon^{bc}k_c)},\\
{e_{TT+}}^{a3}=\frac{1}{\sqrt 2k^2}k^ak^3,& {e_{TT\times}}^{a3}=-\frac{1}{\sqrt
2k}\epsilon^{ac}k_c,\\ {e_{TT+}}^{33}=-\frac{1}{\sqrt 2k^2}k_\perp^2, &
{e_{TT\times}}^{33}=0.
\end{array}\right.
\end{equation}
All variables now have the same reality conditions
\begin{equation}
  \label{eq:3a}
  T^{TTs}_{k_a,k_3}=T^{*TTs}_{-k_a,k_3},\quad T^{T}_{k_a,k_3}=T^{*T}_{-k_a,k_3},\quad
  T^{LT\lambda}_{k_a,k_3}=T^{*LT\lambda}_{-k_a,k_3},\quad T^{LL}_{k_a,k_3}=T^{*LL}_{-k_a,k_3},
\end{equation}
while parity conditions become
\begin{equation}
  \label{eq:79}
  \begin{split}
    & T^{TT\times}_{k_a,k_3}=T^{TT\times}_{k_a,-k_3},\quad T^{LTE}_{k_a,k_3}=T^{LTE}_{k_a,-k_3},\\
    & T^{TT+}_{k_a,k_3}=-T^{TT+}_{k_a,-k_3},\quad
    T^{T}_{k_a,k_3}=-T^{T}_{k_a,-k_3},\quad
    T^{LTH}_{k_a,k_3}=-T^{LTH}_{k_a,-k_3},\quad
    T^{LL}_{k_a,k_3}=-T^{LL}_{k_a,-k_3}.
\end{split}
\end{equation}

If $\Xi=(TTs,T,LT\lambda,LL)$, the kinetic term is
\begin{equation}
  \label{eq:84}
  \int \text{d}^3x\ \partial_0 h_{ij}\pi^{ij=}\sideset{}{'}\sum_{n_a,n_3\geq 0}\ \partial_0h_{\Xi}\pi^{*\Xi},
\end{equation}
so that the Poisson brackets are
\begin{equation}
  \label{eq:86}
  \{h_{\Xi k},\pi^{*\Gamma}_{k'}\}=\delta^\Gamma_\Xi\prod_{i=1}^3\delta_{n_i,n'_i}.
\end{equation}
The Hamiltonian constraints \eqref{eq:57} become
\begin{equation}
  \label{eq:83}
  \begin{split}
    \mathcal H_a&=-2i\sum_{n_a,n_3>0}[k^b\pi_{ab}+k^3\pi_{a3}]{\psi^H_k}
    =-2i\sum_{n_a,n_3>0}k[{e^\lambda}_a\frac{\pi_{LT
      \lambda}}{\sqrt 2}
    +{e^\parallel}_a\pi_{LL}]{\psi^H_k},\\
    \mathcal H_3&=(-2)\sideset{}{'}\sum_{n_a,n_3\geq 0}
    [k^b\pi_{b3}+k^3\pi_{33}]\psi^E_k
    =-2\sideset{}{'}\sum_{n_a,n_3\geq 0}
    k[{e^\lambda}_3\frac{\pi_{LT \lambda}}{\sqrt 2}+{e^\parallel}_3\pi_{LL}]\psi^E_k, \\
    \mathcal H_{\perp}&=-\sum_{n_a,n_3>0}k^2 \sqrt 2  h_T\psi^H_k.
\end{split}
\end{equation}
The Pauli-Fierz Hamiltonian \eqref{eq:56} is given by
\begin{equation}
  \label{eq:91}
  H_{\rm PF}=\sideset{}{'}\sum_{n_a,n_3\geq 0}\big[\pi^{TTs}\pi^*_{TTs}+
  \pi^{LT\lambda}\pi^*_{LT\lambda}+\frac{1}{2}\pi^{LL}\pi^*_{LL}
  +\frac{1}{4}k^2(h_{TTs}h^{*TTs}-h_Th^{*T})\big].
\end{equation}

\subsection{Modes propagating parallel to the plates}
\label{sec:massl-scal-field}

For $n_i\neq 0$, let us start the analysis with the sector at $n_3=0$. It is
useful to consider the dyad
\begin{equation}
  {e_\perp}^b(k_a,0)=\frac{\epsilon^{bc}k_c}{k_\perp},\quad {e_\parallel}^c(k_a,0)
  =\frac{k^c}{k_\perp}, \label{eq:88}
\end{equation}
in terms of which a two-dimensional vector $T^b_{k_a,0}$ may be decomposed as
\begin{equation}
T^b_{k_a,0}=T^\perp_{k_a,0}{e_\perp}^b(k_a,0)+T^\parallel_{k_a,0} {e_\parallel}^b(k_a,0).\label{eq:92}
\end{equation}

When using $T^{ab}_{k_a,0}=0=T^{33}_{k_a,0}$ with only $T^{a3}_{k_a,0}\neq 0$,
it is straightforward to check using the explicit form of the polarization
vectors that
\begin{equation}
  \label{eq:87}
  \begin{split}
    T^{LL}_{k_a,0}=0,\quad T^T_{k_a,0}=0,\quad T^{LTH}_{k_a,0}=0, \quad T^{TT+}_{k_a,0}=0,\\
    T^{LTE}_{k_a,0}= -\frac{2}{\sqrt 2} T^{\parallel 3}_{k_a,0},\quad
    T^{TT\times}_{k_a,0}=-\frac{2}{\sqrt 2} T^{\perp 3}_{k_a,0}.
\end{split}
\end{equation}

At $n_3=0$, $n_a\neq 0$, the constraints reduce to
\begin{equation}
\pi^{LTE}_{k_a,0}=0,\label{eq:104}
\end{equation}
or equivalently to $\pi^{\parallel 3}_{k_a,0}=0$. On the constraint surface, the
contribution from the $n_3=0$ modes to the first order action Pauli-Fierz action
then reduces to the contribution of the canonical pair
\begin{equation}
h^{TT\times}_{k_a,0}=-\frac{2}{\sqrt
  2}h^{\perp 3}_{k_a,0},\quad \pi^{TT\times}_{k_a,0}
=-\frac{2}{\sqrt 2}\pi^{\perp 3}_{k_a,0},\label{eq:93}
\end{equation}
with Hamiltonian
\begin{equation}
  \label{eq:90}
  H_{{\rm PF},n_3=0}\approx \sideset{}{'}\sum_{n_a}(\pi^{TT\times}_{k_a,0}\pi^{*TT\times}_{k_a,0}
  +\frac{1}{4}k_\perp^2 h^{TT\times}_{k_a,0}h^{*TT\times}_{k_a,0}).
\end{equation}

After the canonical transformation $\pi^{TT\times}_{k_a,0}\to \frac{1}{\sqrt
  2}\pi^{TT\times}_{k_a,0}$, $h^{TT\times}_{k_a,0}\to \frac{1}{\sqrt
  2}h^{TT\times}_{k_a,0}$, it is obvious that this collection of harmonic
oscillators is the same than one would get from the expansion of a massless
scalar field in $2$ spatial dimensions with periodic boundary conditions
(omitting its particle mode). That is the reason why we sometimes refer to this
sector as the lower dimensional scalar field. This sector is somewhat unexpected
when comparing to the empty space analysis in appendix
\ref{sec:reduced-phase-space} and is due to the perfectly conducting boundary
conditions discussed in the previous section. From the explicit discussion of
the structure of the constraints and the Hamiltonian, we have thus shown the
following:

{\em With Casimir-type boundary conditions \eqref{eq:78}, gravitational waves
  propagating parallel to the plates are equivalently described by
  electromagnetic waves in $2$ spatial dimensions, whose physical degrees of
  freedom are in turn described by a massless scalar field in $2$ spatial
  dimensions.}

\subsection{Modes with a normal component of propagation}
\label{sec:modes-prop-norm}

When $n_3\neq 0$, the constraints $\mathcal H_i=0$ in \eqref{eq:83} are equivalent
to
\begin{equation}
  \label{eq:95}
  {e^\lambda}_i\frac{\pi_{LT\lambda}}{\sqrt 2}+{e^\parallel}_i \pi_{LL}=0.
\end{equation}
After contracting with ${e_A}^i$, the constraints \eqref{eq:83} are then
equivalent to the vanishing of specific components,
\begin{equation}
  \label{eq:94}
\pi_{LT\lambda}=0=\pi_{LL},\quad h_{T}=0.
\end{equation}
On the constraint surface, the only variables that remain in the first order
Hamiltonian action are $h_{TTs},\pi^{TTs}$ with
\begin{equation}
  \label{eq:89}
  H_{{\rm PF},n_3>0}\approx \sum_{n_a,n_3> 0}(\pi_{TTs}\pi^{*TTs}
  +\frac{1}{4}k^2 h_{TTs}h^{*TTs}).
\end{equation}

Note that, if desired, one may also work with the gauge fixing conditions
\begin{equation}
\triangle \pi-\partial_i\partial_j\pi^{ij}=0,\quad -2\partial_jh^{ij}=0.\label{eq:97}
\end{equation}
and consider Dirac brackets.

\section{Scalar field formulations}
\label{sec:scal-field-form}

\subsection{Bromwich-Borgnis fields}
\label{sec:bromw-borgn-fields}

Disregarding particle zero modes, the reduced phase space consists of two scalar
fields and their momenta, one pair satisfying Dirichlet and the other Neumann
boundary conditions,
\begin{equation}
  \label{eq:99}
  \begin{split}
  h^H=\sum_{n_a,n_3> 0}h^{TT+}_{k}\psi^H_k,\quad \pi^H=\sum_{n_a,n_3> 0}\pi^{TT+}_{k}\psi^H_k,\\
  h^E=\sideset{}{'}\sum_{n_a,n_3\geq 0}h^{TT\times}_k\psi^E_k,\quad
  \pi^E=\sideset{}{'}\sum_{n_a,n_3\geq 0}\pi^{TT\times}_k\psi^E_k.
\end{split}
\end{equation}
On account of the form of the Pauli-Fierz Hamiltonian, \eqref{eq:89} and
\eqref{eq:90}, there is thus a direct gravitational analogue of the electromagnetic results
reviewed in appendix \ref{sec:borgnis-fields}:

{\em With Casimir-type boundary conditions \eqref{eq:78} the dynamics of
  physical gravitons is equivalent to that of two massless scalar fields, one
  with Dirichlet and one with Neumann boundary conditions, with action
\begin{equation}
  \label{eq:100}
  S[h^\lambda]=-\frac{1}{2}\int\text{d}^4 x\ \partial_\mu h^\lambda\partial^\mu h_\lambda,\quad \lambda=(H,E).
\end{equation}
}

Note that the massless scalar field in $2$ spatial dimensions discussed in
section \ref{sec:massl-scal-field} corresponds to the modes at $k_3=0$ contained
in $h_E$.

The set-up follows closely the electromagnetic case discussed in appendix
\ref{sec:reduced-phase-space-2}, with the following substitutions:
\begin{equation}
  \label{eq:107}
  A^E_k\to \frac{1}{\sqrt{2}} h^{TT\times}_k,\quad \pi^E_k\to \sqrt 2  \pi ^{TT\times}_k,\quad
  A^H_k\to \frac{1}{\sqrt{2}} h^{TT+}_k,\quad \pi^H_k\to \sqrt 2  \pi ^{TT+}_k.
\end{equation}
This defines the oscillator variables $a^\lambda_k,a^{*\lambda}_k$ in
\eqref{eq:13} in terms of the Fourier components of the physical gravitational
fluctuations with standard Poisson brackets \eqref{eq:23}. The Pauli-Fierz
Hamiltonian without the particle contributions, given by the sum of the pieces
in \eqref{eq:90} and of \eqref{eq:89}, is equivalent to
\begin{equation}
  \label{eq:105}
  H'_{\rm PF}=\sideset{}{'}\sum_{\lambda,n_a,n_3\geq 0}
  \frac{1}{2} (\pi^\lambda_{k}\pi^{*\lambda}_k+
  k^2h_k^\lambda h^{*\lambda}_k)
  =\sideset{}{'}\sum_{\lambda,n_a,n_3\geq
    0}\frac{k}{2} (a^\lambda_ka^{*\lambda}_k+ a^{*\lambda}_ka^\lambda_k).
\end{equation}

{\bf Remark:}

In momentum space, the shift can be further expanded as
$n^i_k=n^\lambda_k{e_\lambda}^i+n^\parallel_k {e_\parallel}^i$, with
$n^H_k,n^\parallel_k$ even and $n^E_k$ odd. In a variational principle with
constraints included, variation with respect to these variables then directly
impose the constraints \eqref{eq:104} and \eqref{eq:94}. These variables may be
lumped into fields
\begin{equation}
  \label{eq:101}
  n^{H,\parallel}=\sum_{n_a,n_3> 0}n^{H,\parallel}_{k}\psi^H_k,\quad
  n^{E}=\sideset{}{'}\sum_{n_a,n_3\geq 0}n^{E}_{k}\psi^E_k,
\end{equation}
together with the expansion of the shift in \eqref{eq:103}.

In a Hamiltonian BFV-BRST approach, spatial diffeomorphism ghosts and the ghost
associated to deformations normal to the surface that goes with the Hamiltonian
constraint follow the same pattern, with the same expansion holding for their
momenta. In an extended phase space, there are further canonically conjugate
pairs of scalar fields,
\begin{equation}
  \label{eq:102}
  \begin{split}
    & \phi^{T,LTH,LL} = \sum_{n_a,n_3> 0}h^{T,LTH,LL}_{k}\psi^H_k,\quad
    \pi^{T,LTH,LL} = \sum_{n_a,n_3> 0}\pi^{T,LTH,LL}_{k}\psi^H_k,\\
    & \phi^{LTE} = \sideset{}{'}\sum_{n_a,n_3\geq 0}h^{LTE}_{k}\psi^E_k,\quad
    \pi^{LTE} = \sideset{}{'}\sum_{n_a,n_3\geq 0}\pi^{LTE}_{k}\psi^E_k.
  \end{split}
\end{equation}
Quantization of all these variables should be done in terms of quartets so as to
have achieve equivalence with the reduced phase space quantization described
above. After integration over all momenta in the BRST gauge fixed path integral,
this then gives rise to the BRST gauge fixed path integral in the Lagrangian
approach, including the correct boundary conditions on all fields. Such an
extended setting is useful in order to set up computations in covariant gauges
(see e.g.~\cite{Henneaux:1992ig}) for a review), and also when including
interactions. We will not pursue this further here, but focus on reduced phase
space quantization.

\subsection{Additional observable}
\label{sec:addit-observ}

In order to achieve full modular covariance of the partition function below, an
additional observable is of interest. Its expression is most transparent in
terms of generalized vector calculus operations that feature prominently in the
context of the Hamiltonian approach to duality invariance
\cite{Deser:1976iy,Henneaux:2004jw,Deser:2004xt,Deser:2005sz,Barnich:2008ts}.
First, one considers the generalized vector product between a vector $v^k$ and a
symmetric tensor $T_{ij}$
\begin{equation}
  \label{eq:109}
  (\vec v\times T)_{ij}=\frac{1}{2}(\epsilon_{ilm}v^lT^m_j+\epsilon_{jlm}v^lT^m_i).
\end{equation}
This is then implies that the generalized curl for a symmetric spacetime tensor is,
\begin{equation}
  \label{eq:108}
  (\mathcal O T)_{ij}=(\vec \nabla\times T)_{ij}=\frac{1}{2}(\epsilon_{ilm}\partial^lT^m_j+\epsilon_{jlm}
  \partial^lT^m_i).
\end{equation}
It rotates electric and magnetic components and projects out both the trace part
$(T)$ and the longitudinal part of the longitudinal part $(LL)$. Then there is
the operator that projects out the transverse components of the longitudinal
part $(LT)$,
\begin{equation}
  (\mathcal P T)_{ij}=-\triangle T_{ij}+\partial_i \partial_m T^m_j
  +\partial_j \partial_m T^m_i.
\label{eq:117}
\end{equation}
Combining these two operations gives a projector onto the transverse-traceless
part,
\begin{equation}
  \label{eq:118}
  (\mathcal P^{TT} T)_{ij}=-(\mathcal O \mathcal P T)_{ij}=-(\mathcal P \mathcal O T)_{ij}.
\end{equation}
Details on how these operator act are provided in appendix
\ref{sec:gener-curl-proj}.

For simplicity, we start by working in the reduced phase space with
$T^T_{k}=T^{LT\lambda}_k=T^{LL}_k=0$. In terms of space-time fields, we thus
define $T_{TT}^{ij}$ through
\begin{equation}
  \label{eq:116}
  \begin{split}
    &  T^{ab}_{TT}=\sum_{n_a,n_3>0}{e_{TTs}}^{ab}\,T^{TTs}_k\psi_k^H,\quad
    T^{33}_{TT}=\sum_{n_a,n_3>0}{e_{TTs}}^{ab}\,T^{TTs}_k\psi_k^H,\\
    & T^{a3}_{TT}=-i\sideset{}{'}\sum_{n_a,n_3\geq 0}{e_{TTs}}^{a3}\,T^{TTs}_k\psi_k^E.
  \end{split}
\end{equation}

In analogy with \eqref{eq:25}, we then consider spin angular momentum in linearized
gravity, defined on the reduced phase space by
\begin{equation}
  \label{eq:98}
  J_i=\int_V\text{d}^3 x\ \pi^{mn}_{TT}(\vec e_i\times h^{TT})_{mn},
\end{equation}
where $\vec e_i$ is the unit vector along $x^i$ in position space. When using
that
\begin{equation}
  \label{eq:110}
  {e_{TT+}}^{ij}{e_{TT\times}}^{mn}-{e_{TT\times}}^{ij}{e_{TT+}}^{mn}=
  \frac{1}{2}\big[\epsilon^{inl}(\delta^{jm}-\frac{k^jk^m}{k^2})
  +\epsilon^{jml}\left(\delta^{in}-\frac{k^ik^n}{k^2}\right)\big]\frac{k_{l}}{k},
\end{equation}
which follows from \eqref{eq:106}, one finds in terms of modes that
\begin{equation}
  \label{eq:111}
 J_3=\sum_{n_a,n_3>0}\frac{k_3}{k}(h^{TT+}_k\pi^{*TT\times}_k-h^{TT\times}_k\pi^{*TT+}_k)
    = \sum_{n_a,n_3>0}\frac{ik_3}{k}(a^{*H}_ka^{E}_k-a^{*E}_ka^{H}_k).
\end{equation}

The observable we are interested in is
\begin{equation}
  \label{eq:112}
  P_3^S=-\int_V \text{d}^3x\ \partial_3(\mathcal O h)_{mn}
  \frac{1}{\sqrt{-\triangle}^3}(\mathcal P\pi)^{mn},
\end{equation}
which is gauge invariant without the need to put unphysical components to zero
by hand. When using the actions of $\mathcal O$ and $\mathcal P$ established in
appendix \ref{sec:gener-curl-proj}, together with the orthonormality properties
\eqref{eq:67}, it is straightforward to verify that its mode expansion differs
from that of $J_3$ by multiplying each term in momentum space by $k$,
\begin{equation}
  \label{eq:113}
  P^S_3=\sum_{n_a,n_3>0}{k_3}(h^{TT+}_k\pi^{*TT\times}_k-h^{TT\times}_k\pi^{*TT+}_k)
    = \sum_{n_a,n_3>0}{ik_3}(a^{*H}_ka^{E}_k-a^{*E}_ka^{H}_k).
\end{equation}

\subsection{Single scalar field formulation}
\label{sec:single-scalar-field-1}

In the slab geometry considered here, one can go one step further and combine
the Fourier components of the two scalar fields with Neumann and Dirichlet
boundary into those of a single massless free scalar field with periodic
boundary conditions on an interval of length $2d$.

More precisely, from the expression of the Pauli-Fierz Hamiltonian
\eqref{eq:105} and the observable $P_3^S$ \eqref{eq:113} in terms of modes and
oscillators, one may show directly that these coincide with the standard
Hamiltonian, omitting its particle mode, and linear momentum in the $x^3$
direction, of a massless scalar field with periodic boundary conditions in a
rectangular box with sides $L_1,L_2,2d$. The proof is identical to that in the
electromagnetic case reviewed in appendix \ref{sec:single-scalar-field} when
taking into account the identifications \eqref{eq:107}. In summary,

{\em With Casimir-type boundary conditions \eqref{eq:78}, the dynamics of
  physical gravitons is equivalent to that of a massless scalar field with two
  large dimensions and with periodic boundary conditions on an interval of
  length $2d$.}

\section{Partition function}
\label{sec:extend-part-funct}

From the above equivalence, the partition function for gravitons with
Casimir-type boundary conditions may be obtained from that of a massless scalar
field on $\mathbb S^1_{2d}\times \mathbb R^2$. In turn, as shown in
\cite{Alessio:2020okv} and reviewed in appendix \ref{sec:reduced-phase-space-2},
the same equivalence holds in the electromagnetic case. In the absence of a
chemical potential, $\alpha=0$, up to a subtlety related to whether the black
body result is subtracted or not, the result can thus also be obtained from the
literature on the finite temperature electromagnetic Casimir effect
\cite{Mehra:1967wf,Brown:1969na} (see also
\cite{Dowker_1976,BALIAN1978165,10.1143/PTP.75.262,Ambjorn:1981xw,Ambjorn:1981xv,%
  Plunien:1987fr,Lutken:1988ge,Ford:1988gt,PhysRevD.40.4191} and
\cite{Plunien:1986ca,Nesterenko:2005xv,Bordag:2009zzd}) for reviews). It is
however more economical to derive the result, whether the chemical potential is
turned on or not, directly from the scalar field formulation. With $\alpha\neq
0$, this has been discussed instance in \cite{Cappelli:1988vw}. The detailed
derivation by a variety of different methods of the expressions provided below
can be found in \cite{Alessio2020}.

\subsection{Real analytic Eisenstein series from functional approach}
\label{sec:expr-terms-analyt}

The partition function
\begin{equation}
  \label{eq:119}
  Z(\beta,\alpha)={\rm Tr}\ e^{-\beta \hat H+i\alpha \hat P^S_3},
\end{equation}
may be conveniently expressed in terms of the modular
parameter defined by
\begin{equation}
  \label{eq:37}
  \tau=\frac{\alpha+i\beta}{L_3},\quad L_3=2d.
\end{equation}
A further ingredient is the real analytic Eisenstein series,
\begin{equation}
  \label{eq:38}
   f_2(\tau,\bar\tau)=\sum_{(n_{3},n_{4})\in \mathbb Z^2/(0,0)}
  \frac{\mathfrak{Im}(\tau)^2}{\abs{n_{4}+n_3\tau}^{4}},
\end{equation}
which is invariant under modular transformations
\begin{equation}
  \label{eq:39}
  \tau'= \frac{a\tau +b}{c\tau +d}, \quad a,b,c,d\in \mathbb Z,\quad
  ad-bc=1.
\end{equation}

Starting from the Lagrangian path integral formulation and its relation to the
appropriate zeta function, the partition function $\mathcal
Z(\tau,\bar\tau)=Z(\beta,\alpha)$ can be compactly written as
\begin{equation}
  \label{eq:40}
\ln{\mathcal Z}(\tau,\bar\tau)=\frac{1}{2\pi^{2}}
  \frac{L_1L_2}{L_3^{2}}
  \frac{f_{2}(\tau,\bar\tau)}{\mathfrak{Im}(\tau)}.
\end{equation}
Since $f_2(\tau,\bar\tau)$ is modular
invariant and $\mathfrak{Im}(\tau)$ transforms as
\begin{equation}
  \label{eq:41}
  \mathfrak{Im}(\tau')=\frac{\mathfrak{Im}(\tau)}{\abs{c\tau+d}^2},
\end{equation}
the partition function \eqref{eq:40} transforms as
\begin{equation}
  \label{eq:42}
  \ln{\mathcal Z}(\tau',\bar\tau')=\abs{c\tau+d}^2
     \ln{\mathcal Z}(\tau,\bar \tau).
\end{equation}

If the chemical potential vanishes, ${\alpha}=0$, $\tau=it$, with
$t=\frac{\beta}{L_3}$ and ${\rm Z}(t)=\mathcal Z(it,-it)=Z(\beta,0)$,
the result can be written in terms of an Epstein zeta function,
\begin{equation}
  \label{eq:59}
  \mathfrak Z(2;t^2,1)=\sideset{}{'}\sum_{n_{4},n_3}\frac{1}{(n_3^2t^2+n_{4}^2)^2},
\end{equation}
as
\begin{equation}
  \label{eq:58}
  \ln{\rm Z}(t)=\frac{1}{2\pi^{2}}
  \frac{L_1L_2t}{L_3^{2}}\mathfrak Z(2;t^2,1).
\end{equation}
If $a=0=d$, $b=1=-c$, the modular transformation reduces to temperature
inversion $t\to \frac{1}{t}$ with
\begin{equation}
  \label{eq:60}
  \ln {\rm Z}(\frac{1}{t})=t^{2}\ln {\rm Z}(t).
\end{equation}

\subsection{Low and high temperature expansions from canonical approach}
\label{sec:low-high-temperature}

Let
\begin{equation}
  \label{eq:44}
  \xi(d)=\frac{\Gamma(\frac{d}{2})\zeta(d)}{\pi^{\frac{d}{2}}},\quad \xi(4)=\frac{\pi^2}{90},
  \quad \xi(3)=\frac{\zeta(3)}{2\pi}.
\end{equation}
The computation of the Casimir energy of the system yields the same result than
for electromagnetism and can be performed in several ways. The result is
\begin{equation}
  \label{eq:45}
  E_0(0)=-
      \xi(4)\frac{L_1L_2}{L_3^{3}}.
\end{equation}
The partition function can then be directly computed using a quantum statistical
approach in a Hilbert space basis that diagonalizes the Hamiltonian. The result
is
  \begin{multline}
    \label{eq:48}
     \ln{Z}(\beta,\alpha)=-\beta
  E_0+\xi(3)\frac{L_1L_2}{\beta^{2}}\\
 +2\frac{L_1L_2}{L_3^{\frac{3}{2}}\beta^{\frac{1}{2}}}
  \sideset{}{'}\sum_{n_3}\sum_{l\in\mathbb{N}^*}(\frac{\abs{n_3}}{l})^{\frac{3}{2}}
  K_{\frac{3}{2}}(2\pi
      l\abs{n_3}\frac{\beta}{L_3})e^{2\pi i l n_3 \frac{ \alpha }{L_3}},
  \end{multline}
where $K_{\frac{3}{2}}(2\pi m \frac{\beta}{L_3})$ is a modified Bessel function
of the second kind. This result can also be written as
\begin{multline}
  \label{eq:52}
  \ln{Z}(\beta,\alpha)=\frac{L_1L_2}{L_3^2}\Big[\xi(4)\frac{\beta}{L_3}
  +\xi(3)(\frac{L_3}{\beta})^2\\
 +2(\frac{L_3}{\beta})^{\frac{1}{2}}
  \sideset{}{'}\sum_{n_3}\sum_{l\in\mathbb{N}^*}(\frac{\abs{n_3}}{l})^{\frac{3}{2}}
  K_{\frac{3}{2}}(2\pi
      l\abs{n_3}\frac{\beta}{L_3})e^{2\pi i l n_3 \frac{ \alpha }{L_3}}\Big].
\end{multline}

At low temperature/small distance, $\frac{\beta}{L_3}\gg 1$, the leading term in
the expansion to the partition function is directly related to the Casimir
energy. The leading correction is the contribution of the modes with spatial
frequencies $n_3=0$ and with normal ordered Hamiltonian. In this context, these
are the gravitons that propagate parallel to the plates, which is the sector of
the theory explicitly discussed in section \ref{sec:massl-scal-field}. It is
independent of the small dimension $L_3$ and coincides with the black body
result of electromagnetism in $2$ spatial dimensions or with a single massless
scalar field in $2$ spatial dimensions. On account of the asymptotic expansion
\begin{equation}
  \label{eq:47}
  K_\nu(z)=\sqrt{\frac{\pi}{2z}}e^{-z}(1+O(z^{-1})),
  \quad (\abs{{\rm arg}\ z}<\frac{3}{2}\pi),
\end{equation}
for large $\abs{z}$, all other terms are exponentially suppressed.
It follows that in the low-temperature expansion of the entropy,
\begin{equation}
  \label{eq:49}
  S(\beta,\alpha)=(1-\beta\partial_\beta)\ln Z(\beta,\alpha)
  =3\xi(3)\frac{L_1L_2}{\beta^{2}}+\dots,
\end{equation}
the first term in \eqref{eq:48} proportional to the Casimir energy drops out
since it is linear in $\beta$. The dots denote exponentially suppressed terms.
The leading contribution now comes from the lower dimensional
electromagnetic/scalar field, i.e., the modes with $n_3=0$, and scales with the
area.

When expressed in terms of inverse temperature and chemical potential, the
generating set of transformations of the modular group become
\begin{equation}
  \label{eq:50}
  \tau'=\tau+1\iff \left\{\begin{array}{l} \beta'=\beta\\
                              \alpha'=\alpha+L_3
                            \end{array}
                          \right.\quad
                          \tau'=-\frac{1}{\tau}\iff  \left\{\begin{array}{l}
\frac{\beta'}{L_3}=\frac{L_3\beta}{\alpha^2+\beta^2}\\
 \frac{\alpha'}{L_3}= -\frac{L_3\alpha}{\alpha^2+\beta^2}
                            \end{array}
                          \right..
\end{equation}
Under the first of these transformations,
$\ln Z(\beta,\alpha)$ in \eqref{eq:48} is manifestly invariant, as required by
\eqref{eq:42} for $c=0,a=b=d=1$. For the second of these transformations, we get
\begin{equation}
  \ln Z(\beta',\alpha')=\frac{\alpha^2+\beta^2}{L_3^{2}}
    \ln Z(\beta,\alpha).\label{eq:85}
\end{equation}
When transposing and using the explicit expression in \eqref{eq:52} for the LHS,
this gives
\begin{multline}
  \label{eq:51}
  \ln Z(\beta,\alpha)=\frac{L_1L_2}{\alpha^2+\beta^2}
\Big[
  \xi(4)\frac{L_3\beta}{\alpha^2+\beta^2}
  +\xi(3)(\frac{\alpha^2+\beta^2}{L_3\beta})^{2}\\
  +2(\frac{\alpha^2+\beta^2}{L_3\beta})^{\frac{1}{2}}
   \sideset{}{'}\sum_{n_3}\sum_{l\in\mathbb N^*}(\frac{\abs{n_3}}{l})^{\frac{3}{2}}
  K_{\frac{3}{2}}(2\pi l\abs{n_3}\frac{L_3\beta}{\alpha^2+\beta^2})e^{-2\pi i l n_3
    \frac{L_3\alpha}{\alpha^2+\beta^2}}\Big].
\end{multline}
As a consequence, in the high temperature/large distance limit,
$\frac{\beta}{L_3}\ll 1$, the partition function is given by the first term of
the above equation. In the expression of the partition function in terms of the
real analytic Eisenstein \eqref{eq:40}, this term comes from the modes with
vanishing Matsubara frequencies $n_4=0$. For $\alpha=0$, this gives the standard
black body result \eqref{eq:74} when taking into account that the volume of the
system is $V=L_1L_2d$. More generally, the leading contributions are given by
the first two terms, while the others are exponentially suppressed. The high
temperature/large distance expansion of the entropy is
\begin{equation}
  \label{eq:43}
  S(\beta,\alpha)={L_1L_2}
  \Big[L_3\xi(4)\frac{4\beta^3}{(\alpha^2+\beta^2)^{3}}+
  \xi(3)\frac{3\alpha^2+\beta^2}{(L_3\beta)^{2}}+\dots\Big],
\end{equation}
up to exponentially suppressed terms. When $\alpha=0$, this gives the black body
result corrected by a temperature-independent term,
\begin{equation}
  \label{eq:53}
  S(\beta,0)=8\xi(4)
  \frac{L_1L_2d}{\beta^{3}}+
  \frac{\xi(3)}{4}\frac{L_1L_2}{d^{2}}+\dots.
\end{equation}
The leading term is directly related to the Casimir energy of the system, as
discussed in the context of the Cardy formula and modular invariance in higher
dimensions in \cite{Cappelli:1988vw,Cardy:1991kr,Shaghoulian:2015kta}.

\subsection{Complex trigonometric functions from method of images}
\label{sec:expr-terms-hyperb}

Both \eqref{eq:40} and \eqref{eq:52} can be shown to be equivalent to a single
series in terms of complex trigonometric functions,
\begin{multline}
  \label{eq:46}
  \ln\mathcal Z(\tau,\bar\tau)= \frac{\pi^2L_1L_2\mathfrak{Im}(\tau)}{4L_3^2}
  \Big[\frac{2}{45}+\sum_{l\in \mathbb N^*}\big(
  \frac{i\cot{\pi l\tau }}{(\pi l\mathfrak{Im}(\tau))^3}\\
 -\frac{1}{(\pi l\mathfrak{Im}(\tau))^2\sin^2{\pi l\tau}}
  +{\rm c.c.}\big)\Big].
\end{multline}
If the chemical potential $\alpha$ vanishes, $\tau=it=i\frac{\beta}{L_3}$,
by using $\sin(ix)=i\sinh(x)$ and $\cot(ix)=-i\coth(x)$, we get
\begin{equation}
  \label{eq:54}
  \ln{\rm Z}(t)= \frac{\pi^2L_1L_2t}{2L_3^2}
  \Big[\frac{1}{45}+\sum_{l\in \mathbb N^*}\big(
  \frac{\coth{\pi lt }}{(\pi l t)^3}
  +\frac{1}{(\pi l t)^2\sinh^2{\pi l t}}
  \big)\Big].
\end{equation}
This is the form under which the result usually appears in the Casimir
literature. It has originally been derived in \cite{Brown:1969na} by applying a
Sommerfeld-Watson transform to the Epstein zeta function \eqref{eq:59}, which in
turn had been obtained through the method of images.

\section{Discussion}
\label{sec:discussion-1}

As shown by the detailed Hamiltonian analysis, our version of perfectly
conducting boundary conditions for massless spin 2 fields are fully consistent
with gauge invariance. That the existence of such boundary conditions is
non-trivial has been mentioned for instance in \cite{Ambjorn:1981xw} (end of
section 4). It would be interesting to explore in detail how these boundary
conditions relate to those discussed recently in \cite{Witten:2018lgb}.

At the level of the free theory, an open question concerns ADM surface charges
for the set-up we are considering. This is relevant in relation to black holes
because in the Gibbons-Hawking treatment, these charges do play a crucial role.
It would be interesting to see whether they can be understood in terms of the
particles present in the spectrum, in analogy with what can be done for electric
charge in the case of the electromagnetic field between two perfectly conducting
plates. The other obvious question is to extend the analysis to more general
boundary geometries and to curved backgrounds.

With the boundary conditions that we have chosen, the Casimir force on the
``walls'' will be the same for gravitons than for photons. Trying to give a
meaning to this brings one back to the discussion in the introduction of what
type of wall that would be. More generally, in order to distinguish the result
for linearized gravity from that for electromagnetism and to begin to discuss
potential implications, one needs to consider interactions, and thus bring in
Newton's constant, in one way or another. We are looking forward to studying all
these questions in the near future.

\section*{Acknowledgements}
\label{sec:acknowledgements}

\addcontentsline{toc}{section}{Acknowledgments}

This work is supported by the F.R.S.-FNRS Belgium through conventions
FRFC PDR T.1025.14 and IISN 4.4503.15.


text:
\section{Introduction}

The rapid  neutron capture process  (r-process) is believed to  be the
mechanism for the nucleosynthesis of about half of the stable nuclei
heavier than  Iron (\cite{Burbidge,Cameron}). Explosive and neutron-rich
astrophysical environments present  ideal conditions for the r-process
to take  place. Possible candidate  sites discussed in  the literature
include  the  neutrino-driven neutron-rich wind from proto-neutron
stars (\cite{Qian96,Qian03b}), prompt explosions of collapsed stellar cores, (\cite{Sumiyoshi,Wanajo03}),  neutron  star  decompression (\cite{Meyer,Goriely}), tidal disruption in binary merger events (\cite{FRT}), outflows in Gamma-Ray Bursts (\cite{Surman}) etc.  Most  importantly,  abundance data on
r-process elements in  metal-poor stars (\cite{Sneden03}) and certain radionuclides in meteorites (\cite{QW08}) point   toward  the  distinct
possibility  of multiple r-process sites, so that more than
one of these sites may contribute to the observed abundance
pattern of the r-process elements (\cite{Truran}).

\vskip 0.2cm

In 1985, an important work by Bethe \& Wilson
brought neutrino-driven winds from Type II Supernovae (SNe) to the forefront
of the discussion about astrophysical sites for the r-process. Since then, much progress has been made in the modelling of type II SNe and neutrino winds of nascent neutron stars. Some examples of such work
are Woosley et al. (1994), Takahashi et al. (1994), Qian \& Woosley (1996), Cardall \& Fuller (1997),
Otsuki et al.  (2000), Wanajo et al. (2001), Thomson et al. (2001).
However, we have not so far arrived upon a single self-consistent
astrophysical model that reproduces the observed abundance of
r-process elements and at the same time, is consistent with spectroscopic and meteoritic data. The most natural explanation is that the observed r-process pattern follows from a superposition of
neutron capture events with differing neutron-to-seed ratios and
exposure time-scales. Producing the third peak requires
extreme values of entropy and dynamic time-scale that are a challenge for current models of Type II SNe explosions~(\cite{QW08}). Neutron star mergers provide a much larger neutron-to-seed ratio than type II SNe,
but occur so infrequently that they fail to explain the early enrichment
of r-process elements relative to Iron observed in metal-poor stars
(\cite{Qian00}).  Clearly, further study is needed to better understand the r-process conditions,
to this end we present a fairly general and flexible r-process code that is transparent and freely available to the nuclear physics and astrophysics community.

\vskip 0.2cm

\begin{figure*}[t!]
\centering
\includegraphics[width=0.65\textwidth]{figure1.pdf}
\caption{Schematic representation of r-Java's capabilities and tools (see text for meaning of symbols). More details can be found in the manual available on the quark nova project website: http://quarknova.ucalgary.ca/}
\label{rjavachart}
\end{figure*}

 The motivation behind developing our own r-process code came in part from recognizing the difficulty entailed in having access to an openly shared cross-platform program where any interested researcher can examine and modify the dynamical evolution equations and the nuclear input.  We perceive a distinct need for this, as rare-isotope experiments such as FRIB later in this decade will make available new and improved nuclear data for reaction rates relevant to the r-process. Besides, the r-process is of wide interest to the astrophysics community who might wish to explore different dynamical environments for the r-process, and not be as concerned, so to speak, about the input nuclear physics. To avoid a kind of "black-box" situation arising from limited exchange of numerical codes, we decided to build a new code that is flexible enough to incorporate different conditions and has detailed documentation. The code is capable of simulating nuclear statistical equilibrium (NSE), neutron bombardment of a static target as well as the expansion (including the relativistic case) of neutron-rich material. In this work, we discuss all of the above and also apply this code to the unique environment of an ejected and expanding neutron star crust. For the code we designed a graphical user interface (GUI) which allows for real-time interpretation and analysis of simulations through data and graphs.  We have made the r-process application freely available in Java-based format on the quark nova project website\footnote{ http://quarknova.ucalgary.ca/software/rJava/ . If using this code for calculations for a paper, we request that you please acknowledge this website and cite the online arXiv article, or the published version when it becomes available.}, where one can also find input nuclear data and supporting documentation on most aspects of our code.  We request users to test it, and welcome questions and feedback. Our aim is to make our nucleosynthesis code, which we call r-Java, and its workings as transparent as possible so that any scientific user can adapt it to their needs.

\vskip 0.2cm

This paper is organized as follows.  In section 2, we mention our sources for the nuclear input within r-Java and display results for NSE that are in agreement with previous benchmarks. In sections 3 and 4, we present the simplified reaction network (waiting-point approximation) and simulation results for heavy-element production by static and dynamic neutron capture. In section 5, the conclusions from our simulations are discussed and we look ahead to future work.

\section{The R-process Code: r-Java}

\begin{figure}[b!]
    \includegraphics[width=0.5\textwidth]{figure2.pdf}
    \caption{Abundances of nucleons and selected nuclei in the NSE state with varying temperature, at $\rho=1\times 10^{7} g.cm^{-3}$ and $Y_{e}=0.5$}
    \label{NSEvsT}

   \includegraphics[width=0.5\textwidth]{figure3.pdf}
    \caption{Abundances of selected nuclei in the NSE state with a varying electron-to-baryon fraction, at $\rho=1\times 10^{7} g.cm^{-3}$ and $T=3.5 \; T_{9}$}
    \label{NSEvsYe}

    \includegraphics[width=0.5\textwidth]{figure4.pdf}
    \caption{Abundances of nucleons and selected nuclei in the NSE state with varying density, at $Y_e=0.5$ and $T=6.5 \; T_{9}$}
    \label{NSEvsRho}
 \end{figure}

Fundamental to our code is the waiting point approximation which sets lower limits on the temperature of $2\times10^9$K and neutron density of $10^{20}$ cm$^{-3}$ for all our simulations (\cite{Cameron83}). Each of the modules of r-Java; namely NSE, static and dynamic have  their own appropriate set of input parameters and numerical solution methods (see Fig. \ref{rjavachart}).  For example, in the dynamic case, the user can choose the amount of heating from neutrinos, turn on/off nuclear fission, alpha or neutron decay and specify the density profile of the expanding material.  The nuclear input data that is provided along with r-Java uses nuclear masses from experimental sources as much as possible (\cite{ExpMasses}), and where none was available we used \cite{Bdecay}.  For the two nuclear inputs: the neutron separation energy ($S_{n}$) and the $\beta$-decay rates ($\lambda^{\beta}$), we used \cite{Bdecay}.  A freedom provided to the users of r-Java is the ability to modify any nuclear input in order to simulate the effect that such a change would have on abundance yields.  There are many different outputs that can be generated by r-Java as well as the ability to animate how a given parameter or abundance changes during a simulation run.  A unique module included in r-Java displays the periodic table and in real-time shows which elements are being created and destroyed during a simulation.  Also included is a module to study nucleosynthesis in a Quark-Nova~(\cite{Jai07}), which will be discussed in detail in subsequent papers. The following subsections are dedicated to describing in detail the physics involved in each of the modules of r-Java and presenting corresponding simulation results.

\subsection{Nuclear Statistical Equilibrium}

As a prelude to more complex situations, we begin by developing a NSE code where the determining nuclear property is the binding energy. Consider an astrophysical plasma composed of photons, free neutrons and protons, and a mix of seed nuclei such as $^{56}$Fe, with temperature and density high enough that the nuclear reactions assembling nuclei from free neutrons and protons are much faster than the expansion time scale. For large enough temperatures $T_9$=$T/(10^9)$K$\simeq 6$, the system would then reach NSE with forward and backward reaction rates equal. As long as such conditions prevail, the detailed balance equation holds

\begin{equation}
\label{detailbalance}
Z_{i}\mu_{p}+N_{i}\mu_{n}=\mu_{i} \ .
\end{equation}

where $\mu_i$ is the chemical potential of the species $i$ (with $Z_{i}, N_{i}, A_{i}$
 the corresponding atomic number, neutron number, and mass number, respectively). Assuming all the nuclei in the gas follow Maxwell-Boltzmann statistics, the particle number densities are

\begin{equation}
\label{numberdensity}
n_{i}=g_{i} \left(\frac{2\pi kT m_{i}}{h^2}\right)^{3/2} \exp{\left(\frac{\mu_{i}+B_{i}}{kT} \right)} \ .
\end{equation}

where $B$ is the binding energy and $g$ is the degeneracy factor. In addition, mass and charge conservation imply that

\begin{equation}
\label{massconservation}
\sum_{i} A_{i} Y_{i} = 1 \ ,
\end{equation}

\begin{equation}
\label{chargeconservation}
\sum_{i} Z_{i} Y_{i} = Y_{e} \ .
\end{equation}

Combining Eqs. (\ref{numberdensity}) and (\ref{detailbalance}) with the conservation equations (\ref{massconservation}) and (\ref{chargeconservation}), we obtain the following equations for the neutron and proton chemical potentials:

\begin{eqnarray}
\left\{ \begin{array}{l}
 \sum_{i} A_{i} n_{i}( \mu_{p},\mu_{n}) - \rho N_{A} = 0  \\
 \sum_{i} Z_{i} n_{i}( \mu_{p},\mu_{n}) - Y_{e} \rho N_{A} = 0
\end{array} \right.
\end{eqnarray}

For any given $Y_e$ (electron fraction), $T$ (the temperature in Kelvin) and baryon mass density $\rho$ (in g cm$^{-3}$), we then solve this system by the Newton-Raphson method for $ \mu_{n} $ and $\mu_{p}$ (Fig. \ref{rjavachart}), which we use to obtain the mass fractions of every nucleus through Eqn.(\ref{massabundance}), where $N_A$ is the Avogadro number. The nuclei that dominate the abundances distribution in NSE depend only on $\rho,T$ and $B$. High densities favor large nuclei due to the $\rho$-dependence while high temperatures favor light nuclei since the Planck distribution contains energetic photons which photo-disintegrate heavy nuclei. We found the detailed balance method (Eqn.~(\ref{detailbalance})) numerically amenable to the Newton-Raphson method. However, we should caution that the success of the Newton-Raphson method depends sensitively on a proper choice for the initial value of  $\mu_p,\mu_n$ and the starting composition.  We chose to use the NSE calculations of \cite{NSEcode} as a test of our code. We have plotted our NSE results in Fig.\ref{NSEvsT}. We also display the mass fractions as a function of $Y_e$ in Fig.\ref{NSEvsYe}. Comparing the former to the results of \cite{CliffTayler}, and the latter with Figure 3 of \cite{NSEcode}, we find very good agreement with these established results. In Fig.\ref{NSEvsRho}, we have also shown a new plot: the mass fractions obtained with our code as a function of density, for fixed $Y_e=0.5$ and $T=6.5\times 10^9$K.

\section{Reaction Network and Static Case}

\begin{figure*}[t!]
\includegraphics[width=180mm]{figure5.pdf}
\caption{Final abundances in the static simulations after one second of neutron bombardment and at different physical conditions. {\bf Top left:}  $n_{n}=10^{26}$cm$^{-3}$ and $T_9$=2, {\bf Top right:}  $n_{n}=10^{28}$cm$^{-3}$ and $T_9$=2, {\bf Bottom left:}  $n_{n}=10^{26}$cm$^{-3}$ and $T_9$=4, {\bf Bottom right:}  $n_{n}=10^{28}$cm$^{-3}$ and $T_9$=4. The solar abundance is shown in each case with vertical lines marking the location of the important peaks.}
\label{StaticResults}
\end{figure*}

\begin{figure*}[t!]
\includegraphics[width=180mm]{figure6.pdf}
\caption{Final abundances in the static simulations showing how the elements $\beta$-decay back to stability once bombardment of neutrons stops. The cases shown are one second and one hundred seconds following the end of neutron bombardment.  The left panel is the corresponding  Z versus N distribution while the right panel displays how the abundances versus charge number evolve in time.}
\label{StaticResultsBetaDecay}
\end{figure*}

In general, a whole variety of reactions can occur simultaneously, producing and destroying nucleus $i$, whose reaction induced density evolution is described by:
\begin{equation}
\label{reacnet}
\left(\frac{dn_{i}}{dt}\right)_{\rho = {\rm const}} = \sum_{j} N_{j}^{i} r_{j} + \sum_{j,k} N_{j,k}^{i} r_{j,k} \ ,
\end{equation}

where one-body reactions, such as decays and photo-disintegrations (eg. $ r_{j}=\lambda_{j} n_{j}$) and two body reactions ($r_{jk}= \frac{1}{1+\delta_{jk}} <\sigma v>_{jk} n_{j} n_{k}$ where $\delta_{jk}$ is
 the Kronecker delta) dominate. The individual $N^{i}s$ in Eqn.~(\ref{reacnet}) are positive or negative numbers specifying creation and annihilation. Then, the number abundances $Y_{i}$ and mass abundances $X_i$ are given by
\begin{eqnarray}
Y_{i} &=& \frac{n_{i}}{\rho N_{A}} \,,\\
X_{i} &=& \frac{m_{i}n_{i}}{\rho N_{A}} \,.
\label{massabundance}
\end{eqnarray}

In terms of $Y_{i}$, the nuclear reaction network equation (\ref{reacnet}) becomes

\begin{equation}
\label{abunnet}
\dot{Y}_{i} = \sum_{j} N_{j}^{i} \lambda_{j} Y_{j} + \sum_{j,k} \frac{N_{j,k}^{i}}{1+\delta_{jk}} \rho N_{A} <\sigma v>_{jk} Y_{j} Y_{k} \ .
\end{equation}

Applying this reaction network to the r-process, we will limit ourselves at the moment to the following interactions:

\begin{itemize}
\item neutron capture
\item photo-disintegration
\item $\beta$-decay
\end{itemize}

For the static simulation, we assume a static, non-expanding chunk of neutron-rich material near drip density ($\approx 4\times 10^{11}$g/cc), described by the BPS equation of state\footnote{This leaves only T, n$_n$ and Z$_0$ (the dominant nucleus at that density) as inputs since Y$_e$ is prescribed.} (\cite{Crust}), being bombarded by a large and constant flux of neutrons over a given period of time. Eqn. (\ref{abunnet}) is then

\begin{eqnarray}
\label{reducenet}
\dot{Y}(Z,A) &=& \lambda_{Z-1,A}^{\beta} Y(Z-1,A) \;+\; <\sigma v> n_{n} Y(Z,A-1)  \nonumber \\
&+& Y(Z,A+1)\lambda_{Z,A+1}^{\gamma} - Y(Z,A)(<\sigma v> n_{n} \nonumber \\
 &+& \lambda_{Z,A}^{\gamma} \;+\; \lambda_{Z,A}^{\beta}\:)  \ .
\end{eqnarray}

\noindent where $\langle\sigma v\rangle$ is the thermally averaged neutron capture cross-section. Since we have a high neutron density  ($10^{22}$-$10^{26}$cm$^{-3}$) and high temperatures ($T\gtrsim 2\times 10^{9} K$), neutron captures and photo-disintegrations occur much faster than $\beta$-decays. For such conditions we can assume that $(n,\gamma)\rightleftharpoons (\gamma,n)$ equilibrium is established within every isotopic chain. We also assume steady flow while the high density and temperature conditions remain, viz.,
\begin{equation}
\label{ngEqu}
\dot{Y}(Z,A)=0 \ .
\end{equation}
\noindent This leads to steady $\beta$-flow or the standard waiting-point approximation in which we can express the abundance ratio of two neighboring isotopes within one isotopic chain in terms of only the neutron separation energy. In terms of the reduced reaction network (\ref{reducenet}), we then have:

\begin{eqnarray}
\label{abunratio}
\lefteqn{\frac{Y(Z,A+1)}{Y(Z,A)} = } \nonumber \\
&& n_{n}\frac{g(Z,A+1)}{2g(Z,A)}\left(\frac{A+1}{A}\right)^{3/2}\left(\frac{2 \pi \hbar^2}{\mu kT}\right)^{3/2} e^{\frac{S_{n}(Z,A+1)}{kT}} \ .
\end{eqnarray}

\noindent Further, we introduce the abundance of an isotopic chain, $Y(Z)$, defined as the sum of every isotope abundance in that chain: $ Y(Z) = \sum_{A} Y(Z,A)$. With Eqn.(\ref{abunratio}), we can now express every $Y(Z,A)$ in terms of only the corresponding $Y(Z)$ and population coefficients defined through $Y(Z,A)=P(Z,A)Y(Z)$ which depend only on $n_{n}$, T and $S_{n}$. This allows us to dramatically reduce the number of equations from several thousand down to 136, the number of isotopic chains we have in our network. Now, we study the evolution of Y(Z) through $\beta$-decay:

\begin{eqnarray}
\label{eq:Yevolv1}
\dot{Y}(Z) &=& \sum \dot{Y}(Z,A) \nonumber \\
           &=& \sum(\lambda_{Z-1,A}^{\beta} Y(Z-1,A) - \lambda_{Z,A}^{\beta} Y(Z,A) )\ .
 \end{eqnarray}
 This leads to
 \begin{eqnarray}
   \dot{Y}(Z)  &=& \sum(\lambda_{Z-1,A}^{\beta} P(Z-1,A)) \times Y(Z-1) \nonumber \\
           &&  - \sum(\lambda_{Z,A}^{\beta} P(Z,A)) \times Y(Z)\ .
\end{eqnarray}

\noindent Introducing effective $\beta$-decay rates: $\lambda_{Z}^{eff}=\sum(\lambda_{Z,A}^{\beta} P(Z,A))$, we obtain a system of coupled differential equations:

\begin{equation}
\label{EvoEqu}
\dot{Y}(Z) = \lambda_{Z-1}^{eff} Y(Z-1) - \lambda_{Z}^{eff} Y(Z) \ .
\end{equation}

\noindent As the abundances $Y(Z)$ can vary by many orders of magnitude within the integration interval, equation (\ref{EvoEqu}) is a very stiff system, so we use an implicit scheme, the Backward Euler Scheme, to treat the problem.
 To present the capabilities of r-Java in the static regime we varied the neutron density number and the temperature of the material to observe their impact on the r-process abundance (see Fig. \ref{StaticResults}).

Comparing the left and right figures in the top panel of Fig.~\ref{StaticResults} shows clearly that the third r-process peak at A=195 is realized when the neutron number density is high ($n_{n}=10^{28}$cm$^{-3}$) and the temperature is close to $2\times10^9$K. In reality, in dense neutron-rich environments neutron number densities above $10^{32}$cm$^{-3}$ are possible, and we might expect large $A$ nuclei to be copiously produced, but as we have not included fission in the static simulations, we have restrained our code at present to lower neutron densities ($n_{n}=10^{26}-10^{28}$cm$^{-3}$), so that nuclei do not build up without limit.

\subsection{Static Simulation Results}

From the top panels of Fig.~\ref{StaticResults}, it can be seen that our static simulations produce a weak third r-process peak at A=195 even when the neutron number density is relatively high ($n_{n}=10^{28}$cm$^{-3}$). Higher neutron densities can produce a strong third peak ($n_{n}\sim10^{32}{\rm cm}^{-3}$), but the absence of fission presently limits our static code to moderate neutron densities. Absence of fission cycling also leads to a broad peak structure. The bottom panels of Fig. ~\ref{StaticResults} show the effects of doubling the temperature, relative to the top panel, keeping neutron density the same. As temperature increases heavy nuclei tend to disappear since photo-disintegrations of neutron-rich nuclei become highly effective freezing the nuclear flow at lower A. In this case, the peaks in the mass distribution appear at lower mass numbers compared to the corresponding usual r-process peaks. This was found to be the case even for longer periods of neutron exposure, possibly suggesting that photo-disintegrations are dominating over neutron-captures within an isotopic chain, violating the conditions of $(n,\gamma)$ equilibrium. Such a "shift" of the peaks to incorrect locations in mass number has been noticed in earlier dynamical situations, such as neutron star mergers (eg. \cite{Frei}), where it is attributed to the choice of $Y_e$. In our static simulations, $Y_e$ does not cause the shift. Rather, in our case, this effect is driven by the high temperature, which favours photo-disintegrations and takes the r-process path closer to the valley of beta-stability, thereby downplaying the role of waiting-point nuclei in determining the abundance peaks (see \cite{goriely96} for the high neutron density and low temperature, $T_{9} <2$, regime). Extremely high neutron densities are required to restore the peak to its correct position.

\subsection{The $\beta$-decay module}

Although we have implemented the waiting-point approximation, r-Java also provides the option of evolving the elements by $\beta$-decay back to stability from the time bombardment of neutrons is abruptly terminated. Thus, one can isolate key $\beta$-decay rates for  particular reaction networks and studies of freeze-out, which is known to be important, for eg., in explaining the rare-earth peak near Europium~(\cite{Sur97}). We use the Euler implicit scheme for each isotopic chain of elements with decay rate of each element exactly  the ones we used to compute the effective rates defined in Eqn.(\ref{EvoEqu}).
Shown in Figure \ref{StaticResultsBetaDecay} are final abundances from the static simulation (under same
input conditions as in Fig. \ref{StaticResults}) $\beta$-decaying one second and one hundred seconds after the bombardment of neutrons has stopped. An added convenience in r-Java is the fact that the $\beta$-decay module can be run independent of nucleosynthesis calculations by specifying the initial composition as a vector.  We added a modification of the bi-diagonal solver, so that this solver accepts a vector of variable size. Currently the $\beta$-decay module  is limited to nuclei with mass between 50 and 300.

\section{Dynamic Case}
\label{DynamicSims}

The dynamic module of r-Java allows the user to input the initial temperature ($T_0$), density ($\rho_0$), electron fraction ($Y_{\rm e,0}$) and dominantly abundant element (Z$_0$) of a chunk of material.  The material then expands at a user specified expansion time-scale ($\tau$) for a given duration.  In order for r-Java to accommodate a wide variety of ejection mechanisms we have allowed the density profile ($\rho\left(t\right)$) to be freely modified by the user\footnote{The default density profile in the dynamic module of r-Java comes from studies of non-relativistic decompression by \cite{DensityEvo}:

\begin{equation}
\label{density}
  \rho (t) = \frac{\rho_{0}}{(1+\frac{t}{2 \tau})^2} \
\end{equation}}.

 The relativistic expansion regime will be discussed in an upcoming paper, here we discuss the non-relativistic case.  Users of r-Java are given the choice to include a simple fission model for dynamic simulations, in which the last species of the network can be specified.  For our dynamic simulations we assume that the last species in the network is $Z$=$92$, $A$=$276$  which fissions into two fragments--one with $Z$=$40$ and one with $Z$=$52$. The code does not presently include other types of fission (neutron-induced, neutrino-induced, $\beta$-delayed) or $\beta$-delayed neutron emission. Further options given to the user of r-Java for dynamic simulations are the ability to specify the fraction deposited from $\beta$-decays, and allowing $\beta$-decay and neutron decay to be turned on or off.

 Once the user provides the input parameters as above,
 the initial  entropy/nucleon is determined (in $k_B$=1 units) from the initial density, electron fraction and temperature using

\begin{equation}
\label{Entropy}
s = \frac{11 \pi^2}{45}\frac{T^3}{\rho}+ \frac{\mu_e^2T}{3\rho}+\frac{5}{2m_N}+\frac{1}{m_N}{\rm ln}\left(\frac{g_0(2\pi m_Nk_BT)^{3/2}}{n_nh^3}\right)\ ,
\end{equation}

 \noindent where $\mu_e$=$(3\pi^2Y_e\,\rho)^{1/3}$ is the chemical potential of the fully degenerate relativistic electrons, $m_N$=939.1 MeV is the neutron mass and $g_0$=2 (spins) is the statistical weight of the neutron.

 Once the code begins to run, at each and  every time step, we update the abundance of nucleus $i$
if the fractional change $dY_i/Y_i <  0.1$,  else we  decrease  the time  step  until  this condition  is satisfied.  From   Eqn.~(\ref{general_equation}),  this  implies  that $|\partial   n_i/n_i-\partial\rho/\rho|<0.1$.   For   non-relativistic expansion, both $\partial  n_i/n_i$ and $\partial\rho/\rho$  are small, but for  relativistic expansion,  $\partial\rho/\rho$ can be  large at early times. We can still ensure  that the condition on the $Y_i$'s is
satisfied    by   choosing    our   time    step    small   enough that
$\partial\rho/\rho< 0.1$,  provided $\partial\rho/\rho$ dominates over
$\partial  n_i/n_i$. This  is  equivalent  to an  upper  limit on  the
entropy  which  can   be  converted to an upper limit on the initial temperature $T_{\rm up}$(K)=0.01$n_0^{1/3}Z_0$, where $n_0$ is the initial ejecta number density (cm$^{-3}$). For example, taking $Z_0=Z_{Fe}$ (iron-rich ejecta), we find that an initial density $\rho_0\sim 10^{13}$g/cc requires that the initial temperature satisfy $T_0 \lesssim 10^{11}$K.  If the temperature input by the user violates this density relation r-Java asks for a new choice of initial temperature and shows the user the upper bound.  Introducing heat  into the
system  increases   the  entropy  and  can potentially reduce   the  maximum  initial temperature    for   which    our   abundance    criterion    can   be
satisfied.  However,  heating due  to  nuclear  transmutations can  be
ignored at early  times in the case of fast relativistic  expansion, since the
expansion is much faster than typical beta-decay time-scales. Thus, the
above estimate for  the maximum initial temperature remains  valid for fast
expansions.  It is  automatically satisfied  for  the non-relativistic
case,  where $\partial\rho/\rho$  is much  smaller.  Furthermore, once
satisfied  at the  outset, the  abundance criterion  $dY_i/Y_i<0.1$ is
seen to be always satisfied for all future times. These features are naturally
built-in to our code, as detailed by the algorithm flow-chart seen in Figure \ref{Chart}.

\begin{figure}[h!]
\includegraphics[scale=0.5]{figure7.pdf}
\caption{Schematic representation of a single time step in our dynamic code.
In this chart $\rho_0$ is the initial crust density, Y$_{e,0}$ is the initial electron to baryon ratio, Y(Z$_0$) is the initial chemical abundance of the crust, T$_0$ is the initial temperature, Z$_0$ is the initial atomic number and n$_0$ is the initial particle density of the crust.  The neutron density (n$_n$) is first found using the secant method and then reaction network is solved using the Crank-Nicholson method.  After solving the reaction network evolution of the density ($\rho(t)$), heating due to $\beta$ and neutron decay and nuclear fission is calculated.  The neutron density is then re-evaluated and finally the maximal temperature change criterion is checked before moving on to the next time step.}
\label{Chart}

\end{figure}

Equation (\ref{reacnet}) for the evolution of Y(Z) must now be modified due to the density evolution, which was neglected before in the static calculations. The total derivative of $Y(Z)$ is

\begin{equation}
\label{general_equation}
 \frac {d Y_{i}}{dt} = n_{i} \left(\frac{\partial{\frac{1}{\rho}}}{\partial{t}}\right)_{n_{i}={\rm const}} + \frac{1}{\rho } \left(\frac{\partial{n_{i}}}{\partial{t}}\right)_{\rho = {\rm const}} \ .
\end{equation}

\noindent The first term, which describes changes in ejecta density, can be rewritten in term of $Y_{i}$ and then incorporated into our reaction network:

\begin{equation}
\label{yofz}
\frac{n_{i}}{N_{A}} \frac{\partial{\frac{1}{\rho}}}{\partial{t}} = -\frac{n_{i}}{N_{A}} \frac{1}{\rho ^2} \frac{\partial{\rho}}{\partial{t}} = - \frac{1}{\rho} \frac{\partial{\rho}}{\partial{t}} Y_{i} \ .
\end{equation}

\noindent For the dynamic simulation, we used a Crank-Nicholson method to solve Eqn.(\ref{yofz}). As in the static case, we were able to solve the linear system of abundance equations to find the variations in abundances within a timestep. These nuclear transmutations generate entropy and lead to heating, so to keep track of this, we inverted the entropy from these changes to compute the change in temperature. This procedure was implemented in a self-consistent manner, i.e, the new temperature was fed back to check the value of entropy. The heating of the material via $\beta$-decay is added to Eqn.(\ref{Entropy}) as $\delta s = \delta q / T $ and this whole expression is inverted to get the new temperature after each time step. We check that the temperature remains above 2$\times 10^9$K at all times during r-process, so that the waiting point approximation we employ is still valid.

\begin{figure*}[t!]
\includegraphics[width=\textwidth]{figure8.pdf}
\caption{The final abundances (solid red line) for dynamic simulations (one second duration) of nucleosynthesis from expansion of a chunk of pure Iron situated in neutron-rich matter with $Y_e=0.03$. Solar abundances are shown as black crosses.  For each panel the same initial temperature ($4\times10^9$K) and density ($1\times10^{11}$g/cc) is used.  The left column of panels use an expansion time-scale of $\tau=0.1$s and for the right column of panels the expansion time-scale is $\tau=0.001$s.  Each row considers a different density profile: {\bf top} $\rho(t)\propto \left(t/\tau\right)^{-1}$, {\bf middle} $\rho(t)\propto \left(t/\tau\right)^{-2}$, {\bf bottom} $\rho(t)\propto \left(t/\tau\right)^{-3}$.}
\label{dynsim}
\end{figure*}

\begin{figure*}
\includegraphics[width=\textwidth]{figure9.pdf}
\caption{The final Z versus N distribution is plotted for each dynamic simulation (after one second duration).  In each panel the solid green line is at $N=Z$, and the large crosses mark the location of the closed neutron shells.  For each panel the same initial temperature ($4\times10^9$K), density ($1\times10^{11}$g/cc) and electron fraction (0.03) are used.  The left column of panels use an expansion time-scale of $\tau=0.1$s and for the right column of panels, the expansion time-scale is $\tau=0.001$s.  Each row considers a different density profile: {\bf top} $\rho(t)\propto \left(t/\tau\right)^{-1}$, {\bf middle} $\rho(t)\propto \left(t/\tau\right)^{-2}$, {\bf bottom} $\rho(t)\propto \left(t/\tau\right)^{-3}$. }
\label{dynZN}
\end{figure*}

\begin{figure*}
\includegraphics[width=\textwidth]{figure10.pdf}
\caption{The evolution of temperature over the simulation run is plotted for each dynamic simulation (50\% heating from $\beta$-decays is assumed).  For each panel the same initial temperature ($4\times10^9$K), density ($1\times10^{11}$g/cc) and electron fraction (0.03) are used.  The left column of panels use a expansion time-scale of $\tau=0.1$s and for the right column of panels the expansion time-scale is $\tau=0.001$s.  Each row considers a different density profile: {\bf top} $\rho(t)\propto \left(t/\tau\right)^{-1}$, {\bf middle} $\rho(t)\propto \left(t/\tau\right)^{-2}$, {\bf bottom} $\rho(t)\propto \left(t/\tau\right)^{-3}$.  }
\label{dyntemp}
\end{figure*}

To complete the time step, we also need to explicitly calculate the neutron density, which we can do using $(n,\gamma)\rightleftharpoons(\gamma,n)$ equilibrium, as long as temperatures and densities remain high enough. This ensures equilibrium within any isotopic chain, which will instantly redistribute $Y(Z)$ between all isotopes in accordance with the neutron separation energy. Once again, we use the population coefficients which depend on the neutron density, so to get the equilibrium value of neutron density at each timestep, we solve the implicit equation of mass conservation:

\begin{equation}
\label{nn_equation}
 n_{n} m_{n} = \rho N_{a} \left(1 - \sum A(n_{n})\right) \ ,
\end{equation}

\noindent where $\sum A$ is the mass of all nuclei {\it except} neutrons. This quantity depends on the neutron density through the population coefficients, which give the effective mass of every isotopic chain. We used a secant method to solve this equation and get the equilibrium value of $n_{n}$. Then the system is completely described and the code can move to the next time step. We can run the code as long as the waiting point approximation is valid, for $n_{n} > 10^{20} cm^{-3}$ and $T > 2\times 10^{9} K$ (see algorithm in Figure \ref{Chart}).

\subsection{Dynamic Simulation Results}

In Fig. \ref{dynsim}, we explore the effect on final abundance of changing the density profile of the expanding material.  To do this we considered the following general density profile:

\begin{equation}
\rho \left(t \right) = \frac{\rho_0}{\left(1+\frac{t}{\tau}\right)^{\alpha}} \, .
\end{equation}

\noindent For Fig. \ref{dynsim}, Fig. \ref{dynZN} and Fig. \ref{dyntemp} in the top panels $\alpha=1$, the middle panels $\alpha=2$ and the bottom panels $\alpha=3$, where increasing values of $\alpha$ roughly simulate increasing degrees of freedom for expansion.  The left column of these figures use an expansion time-scale of $\tau=0.001$s while for the right column $\tau=0.1$s.  We left the remaining input parameters the same for all panels ($T_0=4\times10^9$K, Y$_{e,0}$=0.03, Z$_0$=26 and a simulation duration of one second).  From the left column of Fig. \ref{dynsim} it can be seen that the density drops off too rapidly to produce elements heavier than A=120 in all but the $\alpha=1$ case.  This is because the expansion is too fast for neutron captures to occur for a sufficiently long period of time to make heavy elements.  An opposite trend can be seen in the $\tau=0.1$s cases where the heavy element production increases with the degrees of freedom for expansion.  This is because for slow expansions, the density, particularly for low degrees of  freedom ($\alpha$=1) remains high for an extended period of time, blocking beta-decays that are needed to move the flow up to higher mass numbers.
Consequently, the third peak is absent or weak. For $\alpha$=2 or 3, the density fall off is rapid enough to compensate this effect, so a strong 3rd peak is seen.  The effect of the beta-decay blocking can be clearly seen in the Z vs. N plots shown in Fig. \ref{dynZN}.  In the top panel, right column of Fig. \ref{dynZN}, it is clear that the slower expansion case allows elements to build up beyond the magic numbers, while for faster expansion (top-left panel of Fig. \ref{dynZN}) or multiple degrees of freedom (right column of 2nd and 3rd panel) there exists clear steps at each magic number, revealing the imprint of the extra stability. It is interesting to note that the final abundances seen in the top left ($\tau=0.001$, $\alpha=1$) and bottom right ($\tau=0.1$, $\alpha=3$) panels of both Fig. \ref{dynsim} and \ref{dynZN} are remarkably similar.  This similarity is due to the fact that after one second of expansion the effects of the different $\alpha$ and $\tau$ values will directly offset one another leaving the final density of the expanding chunk the same.  The slight difference in final abundances of the two simulations are due to temperature effects.  Finally, Fig. \ref{dyntemp} displays the temperature evolution for each of our dynamical simulations.  The left column of Fig. \ref{dyntemp} exemplifies how in the short expansion time-scale adiabatic cooling dominates the heating from beta-decays, while in the longer expansion time-scale cases (the right column of Fig. \ref{dyntemp}) beta decays have a strong heating effect early on, specially for $\alpha$=2,3.

\section{Discussion and Conclusion}

We have constructed and tested a nucleosynthesis code for the r-process under static and dynamic conditions. Called r-Java, it is made available online through an easy-to-use interface at http://quarknova.ucalgary.ca/ . The user is given a number of input options as well as the ability to easily modify the nuclear data used in the code. In this paper, we have studied NSE, as well as both static and dynamic simulations.

Exploring static simulations for a range of neutron irradiation and temperatures, we found that neutron density $\sim 10^{28}$cm$^{-3}$ and above is required for a strong third r-process peak. With increasing temperature beyond $T_9$=2, nuclei photo-disintegrate extremely fast, before significant pile-up can occur at the true waiting-point nuclei. This resulted in abundance peaks at smaller mass numbers. This feature is compensated by increasing the neutron density; however, the effect of fission will then play a role in determining the final abundances. We intend to include fission-rates in our static simulations in the future, so that higher neutron densities can be reliably addressed. We have included a $\beta$-decay module that allows the user to track the flow of $\beta$-decays once neutron bombardment stops. This should be useful in studying freeze-out conditions and quasi-equilibrium reaction networks.

Using the dynamic module of r-Java we compared different density profiles for two expansion time-scales.  We found that at the shorter expansion time-scale ($\tau=0.001$s) only the slowest evolving density profile ($\rho(t)\propto \rho_0\left(t/\tau\right)^{-1}$) produced a significant abundance of elements above A=120.  For the longer expansion time-scale ($\tau=0.1$s) the heavy abundance yield increased with faster evolving density profiles, where the strongest A=195 peak occurred in the simulation that considered $\rho(t)\propto \rho_0\left(t/\tau\right)^{-3}$.  At present, all our calculations are carried out in the waiting-point approximation.

The large qualitative differences in the r-process pattern we observe as a function of dynamical expansion parameters are important for alternate sites of the r-process, such as neutron star mergers or spherically symmetric decompressing neutron star matter. In the case of mergers, tidal forces can lead to a "tube of toothpaste" quasi-1D ejection, whereas spherical decompressions resemble the 3D case. From Fig.\ref{dynsim}, we see that fast expansion is required if heavy elements are made in mergers, while less energy (slower expansion) is required for spherical decompression. This is the case, for eg., in Quark-Novae, which
will be discussed in a follow-up paper. We also note the curious fact that the peak at A = 3, corresponding to Tritium, is formed from the sequence of reactions comprised of neutron decay followed by formation of Deuterium and then Tritium. Since the half-life of Tritium is relatively short (12.32 years), a large amount of its daughter element, He-3, would be formed and dispelled into the surrounding space. In the case of a dual-shock Quark-Nova (\cite{LO,ouyed2010}), this can have important consequences, through spallation on supernova ejecta, for cosmic rays and anomalous abundance ratios in some supernovae (eg., Cas A).

Our future efforts are aimed in three directions: (1) applications of r-Java to Neutron Star Mergers and the Quark-Nova scenario. The Quark-Nova is an ejection mechanism for neutron star material powered by the phase conversion from hadronic to deconfined quark matter occurring at the core of a neutron star~(\cite{ODD,KOJ}). The extreme energetics and high neutron-to-seed ratio makes the Quark-Nova an ideal candidate for an astrophysical r-process site~(\cite{Jai07}); (2) including the effects of freeze-out through a set of quasi-equilibrium reaction networks and other processes neglected in Eq.(\ref{eq:Yevolv1}) such as beta-delayed neutron emission, to obtain a more realistic pattern of final abundances, which can then be used for schematic chemical evolution studies; (3) based on~\cite{Nieb}, further numerical study of the dynamics of the nuclear-quark conversion front inside the neutron star to better constrain the Lorentz factor of the ejecta, and possibly develop a semi-analytic approach to constrain the parameter space, similar to what has been done for neutrino-driven winds from supernovae.

\section{Acknowledgments}  RO, MK and NK are supported by the Natural Sciences and Engineering Research Council of Canada.  NK acknowledges support from Alberta Ingenuity and the Killam Trusts. PJ acknowledges start-up funding from California State University Long Beach.

redpajama_set_name: RedPajamaArXiv