\section{Introduction}

Recent advancements in observational astronomy have now made  long-term
continuous observation of astrophysical sources possible. These
observations are crucial for obtaining consistent and conclusive
interpretation of flux distributions.  The study of long-term flux
distributions is an important tool for analysing and characterising
the stochastic variations in  astrophysical sources.
Characterisation of flux distribution can provide insight into the
 underlying physical process which drives the
variability in the source. For example, the central limit theorem
states that the addition of a larger number of independent random variables will
result in a Gaussian distribution. Therefore, a Gaussian distribution of
flux would imply additive models, where a linear summation
of components contributes to the emission. On the other hand, a log-normal
distribution (which is a Gaussian distribution when the flux is taken in log scale) would suggest that the
mechanism driving the emission is of a multiplicative nature \citep{Lyubarskii1997, Uttley2005}.

Over the past decade,  log-normal flux distributions have been
observed  in  X-ray light curves of compact accreting system like
X-ray binaries and Active Galactic Nuclei (AGN)
\citep{Uttley2001,Vaughan2003,Gaskell2004,Uttley2005}. The emission in
these sources is mainly produced by accretion discs.  This possibly connects
the log-normal behaviour of observed flux distributions to the
accretion disc.
The log-normal flux variations in the  accretion disc have been
explained by the propagation fluctuation model
\citep{Lyubarskii1997}. In this model,  fluctuations in the mass
accretion rate are produced on timescale corresponding to the local viscous time
scales. They propagate inward and couple together
to produce a multiplicative behaviour in the inner part of disk
\citep{Lyubarskii1997,Uttley2001,Uttley2005,Arevalo2006,McHardy2010}.
Log-normal behaviour of flux distribution has also been observed in
sources which are dominated by the non-thermal jet emission
\citep{Giebels2009,Ackermann2015,Sinha2016, Shah2018, Meyer2019, Khatoon2020}. One
possible realisation for such observations is that  disk fluctuations
are possibly imprinted on the jet emission
\citep{Giebels2009,McHardy2010,Shah2018} i.e. the fluctuations from the disk
are efficiently transmitted to the jet and thereby modulating the jet
emission accordingly.
However, the minute time scale variation observed in high energy light
curves \citep{Gaidos1996,Aharonian2007,Paliya2015}
reflect that the jet emission should be independent of accretion disc
fluctuations \citep{Narayan2012}. In such cases, the log-normal
distribution in flux can be produced by the linear Gaussian
perturbations in the intrinsic particle acceleration time-scales in
the acceleration region \citep{Sinha2018}. Moreover, the fluctuations
in the diffuse escape time scales of the emitting electrons would
produce flux distribution shapes other than Gaussian or
log-normal. Alternatively, \citet{Biteau2012} have shown that under
specific conditions, the additive shot noise-model can also produce a
log-normal flux distribution. For example the Doppler boosting of
emission from a large number of randomly oriented mini-jets
results in a flux distribution with features similar to that of a
log-normal distribution.

Two of the common tests used to differentiate between a Gaussian and
a log-normal distribution are the skewness \citep{Zwillinger2000} and the Anderson-Darling (AD) \citep{Stephens1974}.
These tests can be used to determine at some confidence level (a) that the
distribution is not consistent with a Gaussian one and (b) that it
is consistent with a log-normal one. Thus, they are first step prelude
for more detailed investigations such as direct fitting of the binned
distributions. Also, they can be implemented on data with a relatively
small number of points. However, the reliability and effectiveness of
the tests depend on several aspects of the flux light curve such as
the amplitude of the variability, the length of the data and the percentage
of measurement error that the data has. Further, the slope of the power spectral density which characterizes the variability
also has an important impact on the reliability of these tests. In this work, using simulated
data, we study these effects on the effectiveness and reliability of the
tests. Our aim is to quantitatively understand these effects and provide
guidelines to a user as to the minimum number of points (as a function
of the variability strength and measurement errors)  required
for an effective implementation of AD and skewness test. The parameter values corresponding to different
confidence levels for the  AD test, have been numerically
computed in the literature for the case when the lightcurve has an underlying
power spectrum independent of frequency i.e. for white noise. However, lightcurves
often have power spectrum with a power-law frequency dependence and we tabulate
the corresponding parameter values for different confidence levels as
a function of the power-law index. Thus we facilitate the correct use of the
test to realistic lightcurves. Moreover, there are subtle and important
issues regarding the effect of time-binning of the data on the nature
of the resultant flux distribution. In general the binning changes nature of non-Gaussian distribution such as log-normal one, i.e. an intrinsic log-normal distribution may not remain a log-normal when the light curve is binned. Since the observed lightcurves are always binned to some time-bin, it is not clear why several astrophysical systems display log-normal distributions? In this work, we address this paradoxical issue by analysing the simulated lightcurves with the skewness and AD test.

The framework of this paper is as following:-   In section (\S\ref{sec:skew_AD_white}), we carry a detailed investigation of the lightcurves derived from white noise, for the effective and reliable use of skewness and AD test to check the Gaussianity/log-normality of flux distribution. In section (\S\ref{sec:PL}), we explore the efficiency of the skewness and AD tests for light curves with power-law power spectra. We applied the derived results of our simulation to observed flux distribution of one of bright  \emph{Fermi} blazar  in section \S\ref{sec:sim_app}.  In section \S\ref{sec:conc}, we provide summary and discussion of the work with important guidelines for the effective use of skewness and AD test.

\section{Systems with white noise}\label{sec:skew_AD_white}

Skewness is measure of the asymmetry of a distribution and its value  can be zero, positive or negative. For a flux distribution of $N$ points, the sample skewness, which is third standardised moment is computed as the Fisher-Pearson coefficient i.e.,
\begin{equation}
\kappa_F=\frac{\frac{1}{N}\sum_{i\rightarrow 1}^{N}(x_i-\bar{x})^3}{s^3}
\end{equation}
where $\bar{x}$ and $s^2$ are mean and variance of the sample, these two quantities are first and second statistical moments, respectively.
In case of symmetric distribution
  such as a Gaussian one, the skewness value will be zero. Thus, if the
  skewness of a flux distribution is shown to be inconsistent with zero,
  while the skewness of its log flux distribution is shown to be
  consistent with zero, then one can ascertain that the underlying
  distribution is consistent with a log-normal one and inconsistent
  with a Gaussian one.\\
  In addition to skewness test, Anderson Darling (AD) test is also used to check the Gaussianity or log-normality of the flux distribution. AD test is a statistical test which checks the null hypothesis that a sample is drawn from a particular distribution (e.g, normal distribution). The AD test returns the test statistic ($TS$-value) and critical value ($cv$). The $cv$  are function of the number of points in the lightcurve and are obtained numerically at significance levels of 15\%, 10\%, 5\%, 2.5\%, 1\%. The null hypothesis that the sample is drawn from normal distribution is rejected at a significance level, if the returned TS values is greater than the critical value at that significance level. We have used $cv$ at 5\% significance level in this work, which implies that we accept the null hypothesis at 95\% confidence level.  In this section, we consider the case when the distribution is derived from white noise i.e. when the power spectrum is independent of frequency.

  For a finite length lightcurve from a Gaussian flux distribution,
  the measured skewness, $\kappa_F$  will deviate
  from zero with a standard deviation of $\sqrt(6/n)$ \citep{Tabachnick1996}, where $n$ is
  the number of bins in the lightcurve.  A user would then
  check if  the skewness difference of a given lightcurve
\begin{equation}\label{eq:skew_diff}
\Delta\kappa_F=|\kappa_F|-L\sqrt{\frac{6}{n}}
\end{equation}
exceeds zero for
a certain confidence, $L$ then its true  $\kappa_F$ is nonzero. Since it is sensible to use 95\% confidence level, therefore we carried the analysis for $L=2$ (i.e. 2-sigma confidence). For a
lightcurve whose intrinsic distribution is log-normal, $\Delta\kappa_F$
would exceed zero only for a sufficiently large number of points, $n$.
To quantify the minimum number of points required, we simulated
lightcurves corresponding to a log-normal distribution.
In practice, we generated a time-series of the logarithm flux from
  a Gaussian distribution  with centroid $\mu_{LF} = 0$ and some values
  of width  $\sigma_{LF}$ and length $n$, and then converted it to a
  flux lightcurve by taking the exponential of each value. Given the values of $\mu_{LF}$ and $\sigma_{LF}$, one can obtain the theoretical estimate of mean of log-normal distribution $\mu_F$, using $\mu_{F}=\exp(\mu_{LF}+\frac{\sigma_{LF}}{2})$.
  For different sets of  $\sigma_{LF}$ and $n$, we simulated 3000 such
  flux lightcurves, where $\sigma_{LF}$ ranges from 0.1--1. For each of the 3000 simulated
  lightcurves corresponding to a fixed value of  $\sigma_{LF}$ and length $n$, we
  computed the mean $\mu_{F}$, square root of the variance
  $\sigma_{F}$  and the skewness difference $\Delta\kappa_{F}$. The fraction,
  $f$ of the 3000 lightcurves for which  $\Delta\kappa_{F} > 0$ was noted,
  which increased as the number of bins in the lightcurve, $n$ increased.
  Figure \ref{fig:skew_ad_number} shows the minimum number of bins required,  $N_{50}$ and $N_{90}$,
  for 50\% and 90\%  of the lightcurves to have $\Delta\kappa_{F} > 0$,
  as a function of the average normalised standard deviation, defined as
  \begin{equation}\label{eq:sigma_norm}
\sigma_{F,norm}=\frac{\bar{\sigma_F}}{\bar{\mu_F}}
  \end{equation}

  The Figure \ref{fig:skew_ad_number} should be interpreted as
  giving an estimate of the
  minimum number of bins a lightcurve should have in order to
  use skewness to reject the hypothesis that its flux distribution is
  Gaussian. The minimum number of bins required for 90\%  of the lightcurves to have $\Delta\kappa_{F} > 0$ ranges from 25 for $\sigma_{F,norm} \sim 1.1$ to 730 for $\sigma_{F,norm} \sim 0.1$.

  While,  $\Delta\kappa_{F}$ can be used to reject the hypothesis that the
  lightcurve is from a Gaussian distribution, in principle one can go further
  and compare the measured skewness $\kappa_{F}$ with the theoretical value
  for a log-normal distribution,
\begin{equation}\label{eq:pred_skew}
\kappa_T=\frac{\sigma_F^3}{\mu_F^3}+\frac{3\sigma_F}{\mu_F}
\end{equation}
However, from our simulation we found that this is not reliable. The
measured value of skewness deviates significantly from the theoretical value
for  finite number of
data points and large $\sigma_F$. The comparison between the
observed and theoretical skewness values are given in appendix
\ref{apx:SS}.

\begin{figure*}
\centering
\includegraphics[width=0.48\textwidth,angle=270]{ad_skew_sigma_number_2SE.eps}

\caption{Left Panel: The minimum number of bins required,  $N_{50}$ and $N_{90}$,
  for 50\% and 90\%  of the lightcurves to have $\Delta\kappa_{F} > 0$
  as a function of normalised standard deviation $\sigma_{F,norm}$. Right Panel: Same as the
  left one except that the condition is for the AD test, i.e. $TS_{diff} > 0$.
}
\label{fig:skew_ad_number}
\end{figure*}

Measurement errors of the flux can have a significant impact on the skewness
test. If the errors are significant, it would
require a larger number of flux points for a distribution obtained from an intrinsically log-normal
one, to reject the hypothesis that the distribution is Gaussian. To quantify
this effect, we introduce measurement errors which are normally distributed
to the simulated lightcurves. We  define a parameter, error ratio, $R$ as the
ratio of the sigma of the error distribution $\sigma_{err}$ to that of
the intrinsic flux $\sigma_{F}$.   The measurement error is generated as random normal distribution with standard deviation, $\sigma_{err}=R*\sigma_F$. We have added this random normal distribution to the intrinsic simulated log-normal light curves. The AD and skewness tests are then applied to the modified log-normal distribution (combination of random noise and intrinsic log-normal distribution). Figure \ref{fig:skew_ad_frac_error} shows
the variation of $N_{90}$ (left panel)  as a function
of average normalised standard deviation $\sigma_{F.norm}$ for $R = 0.1, 0.5$
and $1.0$. The Figures shows that if the error fraction is larger than 0.5,
the number of points in the lightcurve required to reject the Gaussianity of
the distribution is significantly larger than the case when there are no
measurement errors.

Moreover, in some circumstances the measurement error on the flux depends on the
flux value, $F_i$ i.e. $\sigma_{err}\propto {F_i}^{\alpha}$. We have simulated
lightcurves for such cases and show the result for $\alpha = 0.1$ and $1.0$
in Figure \ref{fig:skew_sigma_number_error_slope}. Here $N_{90}$ is plotted
against $\sigma_{F.norm}$ for $R = 0.1$ and $1.0$. For $\alpha = 1.0$, the
number of points required to reject the Gaussianity is less than for $\alpha = 0$.

Similarly,  we repeated the above analysis for the AD test. For a lightcurve, we define a measure
\begin{equation}\label{eq:ts_diff}
TS_{diff.}=TS-cv_{5\%}
\end{equation}
where $cv_{5\%}$ is the critical value such that the probability that TS is
larger for a Gaussian distribution is 5\%. Analogous to the
skewness test, we use simulations to estimate $N_{90}$ ($N_{50}$),
the minimum number of bins for which $90\%$ ($50\%$) of the simulated
lightcurves have a $TS_{diff} > 0$. The results are presented in the right
panels of Figures \ref{fig:skew_ad_number},  \ref{fig:skew_ad_frac_error} and \ref{fig:skew_sigma_number_error_slope}.

For a log-normal distribution, the skewness and AD test will accept
the normality of the logarithm of the flux at more than 95\% confidence level by definition. However, if there is
measurement error in the flux, then this may not be so and for sufficiently
large error and large number of data points, the test would indicate that
the log of the flux distribution is not Gaussian. We simulated lightcurves
with measurement errors and find that when the number of points exceeds
$10^{5}$ ($10^4$), the skewness test fails when the error ratio is 0.2 (0.5).
Similar results are obtained for the AD test.

As a thumb rule, our results show that both  skewness and AD tests can be
used to determine that the
distribution is not Gaussian, for lightcurves having more than 100 points and with
a normalised standard deviation (fractional r.m.s of variability) greater
than 30\%. The ratio of the sigma of the measurement error to that of
the intrinsic variability should be less $0.2$. The Figures presented can
be used to gauge the reliability of the tests for other values.

The measurement errors for a lightcurve can be decreased by binning the
lightcurve in time. Moreover, most measurement processes involves summing
the fluxes for a minimum time which is the time resolution of the instrument.
While a Gaussian distribution remains a Gaussian one with binning in time, this
is not the case for an intrinsically log-normal one. We have performed simulations to illustrate this point quantitatively. We binned  simulated lightcurves from a log-normal distributions such that the resultant binned lightcurve has
400 data points and a $\sigma_L \sim 0.3$. The logarithm of the fluxes of these
binned lightcurves were then subjected to the skewness and AD test, and the
fraction of the lightcurves which failed the tests were noted. This fraction
for the skewness and AD test are plotted as a function of the number points
averaged in Figure \ref{fig:flux_bin_skew_ad}. It is clear from the Figure that binning changes the
intrinsic log distribution to some other distribution.

Thus, an intrinsically ``true'' log-normal distribution which has a white noise
power spectrum cannot be determined to be such, if the lightcurve has been binned. If there is a cutoff in the power spectrum such that there is no intrinsic
variability on time-scales shorter than a particular value and if the time
resolution is smaller than this value, then the distribution maybe determined to
be log-normal. This is further discussed in the last section.

\begin{figure*}
\centering
\subfigure{\includegraphics[scale=0.49,angle=270]{err_skew_ad_sigma_number_2SE.eps}}
\caption{Left panel: Minimum number of flux points, $N_{90}$, required for 90\% of light curves to have $\Delta \kappa_F>0$, is plotted as function of  normalize standard deviation  $\sigma_{F,norm}$. The three curves with points shown by colours red (filled square), green (open circle) and blue (filled circle) are for error fraction, R=0.1, 0.5 and 1 respectively. Right Panel: Same as left one except the condition is for the AD test i.e. $TS_{diff}>0$.}
\label{fig:skew_ad_frac_error}
\end{figure*}

\begin{figure*}
\subfigure{\includegraphics[scale=0.48,angle=270]{slope_err_skew_ad_sigma_number_2SE.eps}}\quad \hspace{-0.5cm}
\caption{Minimum number of flux bins, $N_{90}$,  as function of normalised standard deviation of flux
  distribution $\sigma_{F,norm}$,  when the measurement error depend on the flux value.  The left and right panel are obtained using skewness and AD tests respectively. In both panels, black and red curves correspond to $\alpha=$ 0.1 and 1.0 respectively. The dashed and solid representation of curves are for R= 0.1 and 1.0 respectively.}
\label{fig:skew_sigma_number_error_slope}
\end{figure*}

\begin{figure}
\subfigure{\includegraphics[scale=0.33,angle=270]{ad_skew_bin_lognor_frac.eps}}\quad \hspace{-0.5cm}
\caption{Top panel: The fraction for which log-normality is rejected in the binned flux distribution $f(N)$ as function of the number of flux points over which averaging is done.
The solid red curve with filled square points corresponds to skewness test, while the dotted blue curve with open circles corresponds to AD test.}
\label{fig:flux_bin_skew_ad}
\end{figure}

\section{Systems with power-law noise}\label{sec:PL}

For many physical systems, the power spectra generated from the lightcurves
can be described as a power-law i.e. $P(f) \propto f^{-\beta}$. In this section,
we explore the efficiency of the skewness and AD tests for such
systems. Random time series representing power-law
power spectra and whose distribution are Gaussian can be generated
using the algorithm proposed by \citet{Timmer1995}.  Using this algorithm, we generated a large sample of light curves for different sets of $n$ and $\beta$, the minimum and maximum frequency of the corresponding power law power spectra are chosen as $f_{min}=\frac{1}{n}$ and $f_{max}=0.5$ , respectively. For $\beta\leq1$, the time-series generated  with \citet{Timmer1995} algorithm will have Gaussian PDF for all values of $n$  \citep{Morris2019} . However, for steeper power spectral index ($\beta>1$),
the light curves becomes more non-stationary and generated time-series are less likely to be described by Gaussian PDF \citep{Alston2019, Morris2019}. Therefore, in this section we mainly focus on the results for $\beta\lesssim1$.

In the previous section, the analysis was based on systems with white
noise i.e. for the special case of $\beta = 0$. It should be noted that
 the standard deviation of
skewness as well as the critical values for the AD tests
are valid only for white noise  and cannot be used
for the general  non-zero $\beta$ cases. Since we did not find these
values in the Literature, we first undertook a simulation study to
estimate them.

In order to estimate standard error of skewness for nonzero $\beta$-cases, we simulated lightcurves having power law PSD  for different sets of $n$ and $\beta$. For each pair of $n$ and $\beta$, we generated 2000 random light curves and measured the skewness value in each trial. We plotted the standard deviation of skewness, $\sigma_k$ (standard deviation obtained from 2000 skewness
values) with $n$ in Figure \ref{fig:std_skew_err}.  For $\beta = 0$, the variation in $\sigma_k$ with
$n$ is as expected i.e.  $\sigma_\kappa = \sqrt(6/n)$ and this is shown as
a solid line. For $\beta < 0.5$ the standard deviation of skewness is close to that
for $\beta = 0$, but for larger $\beta$ there is significant deviation.
We empirically fit the variation with $n$ using a generic function
\begin{equation}\label{eq:skew_err}
\sigma_\kappa (n)=(a/n)^b
\end{equation}
and tabulate the best fit values of $a$ and $b$ for different values
of $\beta$ in Table \ref{table:skew_err}.
The fits are shown by dotted lines in Figure \ref{fig:std_skew_err}.
For a stationary time-series, the error (i.e. the standard deviation
  of the skewness) should decrease with the number of elements, such that for
  a long enough series , the skewness is well defined and convergent.
  This is indeed confirmed here for white noise ($\beta = 0$) which is expected to be
  stationary.   For $\beta\lesssim1$, the slope of standard deviation of skewness with $n$ gets flattened (i.e. its dependence on $n$ decreases) showing that
  the system is approaching non-stationarity.
Further, for $\beta>1$, \citet{Morris2019} have shown that the light curves generated with the \citet{Timmer1995} algorithm are less likely to be Gaussian due non-stationarity. This is further established here as in such cases the standard deviation in skewness remains almost constant with $n$. For example, as shown in Figure \ref{fig:std_skew_err}, for $\beta\gtrsim 1.5$, $\sigma_{\kappa}$ is independent of $n$. Thus for such cases the tests are not reliable. The  divergence of tests for $\beta>1$ can thus be mainly attributed to increase in non-stationarity, where the variance of individual lightcurve segments is not constant in time but instead varies about a mean value determined by the underlying PSD \citep{Uttley2005}.

\begin{figure}
\centering
\subfigure{\includegraphics[scale=0.32,angle=270]{std_skew_err.eps}}\quad \hspace{-0.5cm}
\caption{The standard deviation of skewness if plotted against the number of flux points for different values of the power spectrum
index $\beta$. The solid line is the expected standard deviation $\sigma_\kappa = \sqrt{6/n}$ valid for the case when $\beta = 0$.
The dotted lines are the fitted lines using Equation \ref{eq:skew_err}, with best fit values tabulated in Table \ref{table:skew_err}.}
\label{fig:std_skew_err}
\end{figure}

\begin{table}
\centering
\begin{tabular}{lcr}
\hline
$\beta$ &  a  &  b \\
\hline
\hline
0.0 &  5.96$\pm$0.17 & 0.499$\pm$0.003\\
0.1 &  5.63$\pm$0.21 & 0.494$\pm$0.003\\
0.2 &  6.86$\pm$0.25 & 0.513$\pm$0.004\\
0.3 &  5.89$\pm$0.27 & 0.497$\pm$0.004\\
0.4 &  6.04$\pm$0.33 & 0.496$\pm$0.005\\
0.5 &  5.58$\pm$0.24 & 0.482$\pm$0.004\\
0.6 &  4.54$\pm$0.23 & 0.451$\pm$0.004\\
0.7 &  3.92$\pm$0.19 & 0.418$\pm$0.004\\
0.8 &  1.81$\pm$0.18 & 0.345$\pm$0.005\\
0.9 &  0.58$\pm$0.07 & 0.269$\pm$0.004\\
1.0 &  0.06$\pm$0.02 & 0.189$\pm$0.005\\

\hline
\hline
\end{tabular}
\caption{Skewness error fit parameters `a' and `b' obtained by fitting the equation (\ref{eq:skew_err}) to the variation of skewness standard deviation with $n$ for different values of $\beta$.}
\label{table:skew_err}
\end{table}

Similarly using 2000 simulated lightcurves,  we estimated the critical
values for the AD test for different values of $\beta$ and
significance levels as a function of the number of bins and the results
are tabulated in Table \ref{table:critical_val}. As a check we note that the critical values obtained
for $\beta=0$ are similar to those given in Literature (written in parenthesis
in Table  \ref{table:critical_val}).

Having obtained the standard deviation for skewness and the critical values
of AD test applicable for lightcurves having power-law power spectra, we
perform similar analysis as was done in the previous section for white noise.
In the end of the previous section it was shown that if a lightcurve
having a log-normal distribution and white noise, is binned, the resultant
distribution is no longer a log-normal one and both the skewness and AD test
will reject the hypotheses that it is log-normal. We first test the effect of
binning on lightcurves generated from a power-law distribution. We simulated a sample of Gaussian time-series using \citet{Timmer1995} algorithm, which were exponentiated in order to obtain lightcurves having log-normal distribution. The exponential model reproduces the observed properties for both stationary and non-stationary light curves \citep{Alstonet2019, Alston2019}. It should be noted that the shape of power spectra may not be preserved by exponentiating a time-series. However, \citet{Uttley2005} showed that the effect is negligible for PSD's observed in astrophysical system and it is acceptable to consider the shape of power spectra of log-normal distribution same as that of input linear time-series. The generated log-normal lightcurves are binned in order to have lightcurves with 400 points and a $\sigma_L \sim 0.3$. The distribution of the logarithm of the fluxes were tested using
skewness and AD tests. In Figure \ref{fig:flux_pl_bin_skew_ad}, the fraction of the lightcurves which failed the tests are plotted as a function of the binning factor. As expected
for no binning ($N_{average} = 1$), the fraction $f_{\kappa}$ and $f_{AD}$ is 0.05
since a 95\% confidence has been taken for the tests. As shown earlier,
for white noise ($\beta = 0$), the fractions increase with binning showing that
the distribution is no longer consistent with being a log-normal. However,
for $\beta = 0.5$ and $\beta = 1.0$, the fractions are more or less independent
of binning indicating that the log-normal nature is maintained after binning
for systems where the power spectrum has a power-law shape with index greater than 0.5. A plausible reason that binning preserves the log-normality could be related to low power in high frequency components for lightcurve having steeper PSD.  The maximum frequency in case of binned lightcurve can be written as $f_{max}=\frac{1}{2N_{averaged}\Delta t}$, where $\Delta t$ is the time bin of unbinned (or original) light curve, the maximum frequency decreases as the binning (or $N_{averaged}$) increases.  In case of steeper PSD, most power is contained in low frequency components and less power will be at higher frequencies. In such cases, the effect of truncation of high frequencies will be less, hence the properties of the original light curve (e.g, log-normal) are preserved. However, in case of white noise, the variability at higher frequencies strongly influence the behaviour of the time-series, hence removal of high frequencies from the PSD will change the time series property effectively.

\begin{figure*}
\subfigure{\includegraphics[scale=0.48,angle=270]{skew_ad_pl_bin_lognor_frac.eps}}
\caption{Left panel:The fraction $f_{\kappa}(N)$ for which log-normality of the binned light curves is rejected ($\Delta \kappa>0$), plotted against the  number of averaged bins. The length of each binned light curve is 400.
The three curves viz. violet curve with open circle points, red curve with filled square points and blue curve with filled triangle points corresponds to the light curves which are generated from  power law noise with $\beta=0.0$, 0.5 and 1.0 respectively. Right panel: Same as the left one except the condition used is for AD test.}
\label{fig:flux_pl_bin_skew_ad}
\end{figure*}

We continue the analysis by estimating the minimum number of bins required to
reject the hypothesis that the flux distribution of such lightcurves are Gaussian.
The left panel of Figure \ref{fig:skew_frac}, displays the minimum number of points required for a log-normal lightcurve, such that for 90\% of them $\Delta \kappa_F > 0$, $N_{90}$, as a
function of the normalised standard deviation.
The plots are for $\beta = 1.0$, $0.5$ and $0.0$, where $\beta = 0.0$ is
a repetition of the results obtained in the previous section for white noise.
The right panel of Figure \ref{fig:skew_frac}, is the same as the left except that it is
for the case when $TS_{diff} >0$ for the AD test. We note that the results are
similar for $\beta = 0.5$ and white noise, but there is marked difference for
$\beta = 1.0$, where a substantially larger number of bins is required to
reject Gaussianity for normalised standard deviations less than $0.5$ as
compared to white noise. The effect of measurement error is shown in
Figure \ref{fig:skew_frac_error}, where $N_{90}$ as a function of the normalised standard deviation
is shown for different values of the error ratio, $R$ for $\beta = 1.0$.
As expected increasing the ratio of measurement errors to the intrinsic
variation, $R = \sigma_{err}/\sigma_F$), requires a larger number of data points
to reject Gaussianity.

Thus, the skewness and AD tests to identify a system has a log-normal
distribution can only be effectively used if the power spectra  has
an index, $\beta \gtrsim 0.5$, otherwise flux binning would modify the
distribution. On the other hand, for large $\beta \sim  1$, the number
of bins required to reject the hypothesis that the distribution is Gaussian,
is significantly larger ($\gtrsim1000$), than for lower values
of $\beta$. The most effective analysis can be undertaken when
$\beta \sim 0.5$, error ration $R<0.2$, the normalised standard deviation of the lightcurve
is $\gtrsim 30\%$ and the number of bins is $> 100$.

\begin{figure*}
\centering
\includegraphics[width=0.46\textwidth,angle=270]{ad_skew_sigma_number_pl_noise_beta.eps}

\caption{The minimum number of flux points required,  $N_{90}$ such that
  90\% of the lightcurve generated will  have $\Delta\kappa_{F} > 0$ (left panel) and $TS_{diff}>0$ (right panel)
  as a function of normalised $\sigma_{F,norm}$. In both panels, the three curves viz. green curve with filled squares, blue curve with open circles and red curve with filled circles corresponds to $\beta$=0.0, 0.5 and 1.0 respectively.
}
\label{fig:skew_frac}
\end{figure*}

\begin{figure*}
\centering
\subfigure{\includegraphics[scale=0.495,angle=270]{err_skew_ad_sigma_number_plnoise_beta_1.eps}}
\caption{Left panel: Minimum number of flux points, $N_{90}$, required for 90\% of the light curves generated from power law noise with $\beta=1.0$, to have $\Delta \kappa_F>0$, is plotted as function of  normalise standard deviation  $\sigma_{F,norm}$. The three curve's with points shown by colours red (filled square), green (open circle) and blue (filled circle) are for error fraction, R=0.1, 0.5 and 1 respectively. Right Panel: Same as left one except the condition  $TS_{diff}>0$ used is for the case AD test.}
\label{fig:skew_frac_error}
\end{figure*}

\begin{table*}
\caption{critical values of AD test for power law noise with  $\beta=0.0, 0.5, 1.0$ obtained from simulations. The critical values are obtained at significance level  of 15\%, 10\%, 5\%, 2.5\% and 1\% .}
\label{table:critical_val}
\begin{adjustbox}{width=0.4\textwidth,center=\textwidth}
\begin{tabular}{lccc}
\hline
 Slope & number of & significance  & CV calculated  (CV given)  \\
($\beta$) &	 data points (n)	&  level	\%	& 	\\
\hline
\hline
0.0 & 4096   & 	15	&	0.56 (0.575)	\\
0.0 & 2048	&	15	&	0.57 (0.575)	\\
0.0 & 1024 	&	15	&	0.56 (0.574)	\\
0.0 & 512  	&       15	&	0.56 (0.572) \\
0.0 & 100	&	15	&	0.55 (0.555) \\
\hline
0.0 & 4096 	&	10	&	0.62 (0.655)	\\
0.0 & 2048   &	10	&	0.63 (0.655)\\
0.0 & 1024 	&	10	&	0.63 (0.653)\\
0.0 & 512 	&       10	&	0.63 (0.651)\\
0.0 & 100	&	10	&	0.68 (0.632)\\
\hline
0.0  & 4096 	&	5	&	0.75 (0.786)	\\
0.0  & 2048 	&	5	&	0.74 (0.785)\\
0.0  & 1024 	&	5	&	0.75 (0.784)\\
0.0  & 512  	&       5	&	0.75 (0.781)\\
0.0 & 100 & 5 & 0.75 (0.759) \\

\hline
0.0  & 4096 	&	2.5	&	0.88 (0.917)	\\
0.0  & 2048 	&	2.5	&	0.87 (0.916)\\
0.0  & 1024	&	2.5	&	0.89 (0.914)\\
0.0  & 512  	&       2.5	&	0.87 (0.911)	\\
0.0  & 100  	&       2.5	&	0.86 (0.885) \\
\hline
0.0  & 4096 &	1.0	&	1.02 (1.091)	\\
0.0  & 2048 	&	1.0	&	1.04 (1.09)	\\
0.0  & 1024 &	1.0	&	1.06 (1.088)\\
0.0  & 512  	&       1.0	&	1.03 (1.084)\\
0.0  & 100  	&       1.0	&	1.02 (1.053)\\
\hline
\hline
0.5 & 4096  &	15	&	0.58\\
0.5 & 2048  &	15	&	0.58	\\
0.5 & 1024  &	15	&	0.58	\\
0.5 & 512  &  15	&	0.57	\\
0.5 & 100	&	15	&	0.56 \\
\hline
0.5 & 4096 &  10	&	0.66	\\
0.5 & 2048 &  10	&	0.65	\\
0.5 & 1024 &  10	&	0.65	\\
0.5 & 512  &   10	&	0.65	\\
0.5 & 100  &   10	&	0.68 \\
\hline
0.5 & 4096 & 5	&	0.80	\\
0.5 & 2048 & 	5	&	0.79	\\
0.5 & 1024 &	5	&	0.79	\\
0.5 & 512  &  5	&	0.78	\\
0.5 & 100 & 5 & 0.76 \\
\hline
0.5 & 4096 & 2.5	&	0.94	\\
0.5 & 2048 & 	2.5	&	0.93	\\
0.5 & 1024 & 	2.5	&	0.92	\\
0.5 & 512  &   2.5	&	0.91	\\
0.5 & 100  &   2.5	&	0.88 \\
\hline
0.5 & 4096 & 	1.0	&	1.11	\\
0.5 & 2048 & 	1.0	&	1.12	\\
0.5 & 1024 &	1.0	&	1.09	\\
0.5 & 512  &  1.0	&	1.04	\\
0.5 & 100  &  1.0	&	1.05 \\
\hline
\hline
1.0 & 4096 & 	15	&	2.62	\\
1.0 & 2048 & 	15	&	2.09	\\
1.0 & 1024 & 	15	&	1.57	\\
1.0 & 512  &       15	&	1.15	\\
1.0	& 100	&	15	&	0.69 \\
\hline
1.0 & 4096 & 	10	&	3.30	\\
1.0 & 2048 & 	10	&	2.65	\\
1.0 & 1024 & 	10	&	1.95	\\
1.0 & 512  &       10	&	1.40	\\
1.0 & 100  &       10	&	0.67 \\
\hline
1.0 & 4096 & 	5	&	4.70	\\
1.0 & 2048 & 	5	&	3.72	\\
1.0 & 1024 &   5	&	2.68	\\
1.0 & 512  &        5	&	1.64	\\
1.0 & 100 & 5 & 0.98 \\
\hline
1.0 & 4096 & 	2.5	&	5.85	\\
1.0 & 2048 & 	2.5	&	4.83	\\
1.0 & 1024 & 	2.5	&	3.47	\\
1.0 & 512  &       2.5	&	2.42	\\
1.0 & 100  &       2.5	&	1.17	\\
\hline
1.0 & 4096 & 	1.0	&	7.76	\\
1.0 & 2048 &    1.0	&	6.67	\\
1.0 & 1024 & 	1.0	&	4.90	\\
1.0 & 512  &       1.0	&	3.26	\\
1.0 & 100  &       1.0	&	1.43   \\
\hline
\hline
\\
\end{tabular}
\end{adjustbox}
\end{table*}

\section {Application of simulation results on blazar $\gamma$-ray observations}\label{sec:sim_app}
The  $\gamma$-ray flux distribution for a sample of bright \emph{Fermi}-blazars has been studied in detail using skewness, AD test and histogram fitting  by \cite{Shah2018}. To check the application of our results derived in section \S 3,  we  consider one of the bright \emph{Fermi} blazar viz. 3FGL\,J0730.2-1141 from the sample studied by \cite{Shah2018}. Using the skewness and AD test, \cite{Shah2018} showed that the monthly binned  $\gamma$-ray flux distribution  of this source is log-normal distribution.
The number of flux points (after applying quality cuts) in the considered distribution are 99 with normalised standard deviation, $\sigma_{F,norm}\sim0.68$. The skewness, $\kappa$ and AD TS values of the logarithm of flux distribution are obtained as -0.35 and 0.62 respectively. In order to check the reliability of the AD and skewness test for confirming the log-normality of  flux distribution, we computed the error ratio `R' and power spectra of the monthly binned $\gamma$-ray light curve. The R-value for the observed light curve is obtained as 0.0417. We computed the averaged power spectrum by dividing the monthly binned light curve into four segments, each segment of duration 24 months with a binning time resolution of one month. The power law fit ($P(f)\propto f^{-\beta}$) to the averaged power spectrum, shown in Figure \ref{fig:psd_fit}, results in $\beta=1.02\pm0.22$.  Corresponding to $\beta\sim1.0$ and n=99, the standard deviation on skewness $\sigma_{\kappa}$  and $cv_{5\%}$ are obtained as $\sim$0.25 and $\sim$0.98 respectively (see Table \ref{table:skew_err} and \ref{table:critical_val}). On using the $\kappa$, $\sigma_\kappa$, TS and $cv_{5\%}$, the values of $\Delta \kappa$ and $TS_{diff}$ of the logarithm of flux distribution are obtained as -0.15 and -0.36 respectively, both values suggesting the log-normal behaviour of the flux distribution.  Further, as shown in Figure \ref{fig:skew_frac_error}, for $\beta\sim 1$, $\sigma_{F,norm}\sim0.68$ and $R<0.1$, the available number of flux points, n=99 are sufficient to satisfy the condition required to reject the Gaussianity of flux distribution at 90\% confidence.
\begin{figure}
\centering
\includegraphics[width=0.34\textwidth,angle=270]{psd_fsrq.eps}
\caption{Power spectrum of the monthly binned $\gamma$-ray light curve of 3FGL\,J0730.2-1141. The red solid line is the power law fit with $\beta=1.02\pm0.22$.}
\label{fig:psd_fit}
\end{figure}

\section{Summary and Discussion}\label{sec:conc}

In this work, using simulations, we have studied the efficiency and reliability
of the skewness and AD tests to study whether a
lightcurve corresponds to a log-normal distribution. The tests are
used to ascertain whether the lightcurve's flux distribution is inconsistent
with a Gaussian one, while the distribution of the log of flux is
consistent with being a normal distribution. Our motivation is to
provide a prescription or guidelines for the effective use of these tests for
this endeavour. The recommended prescription is described below, while
an example of its implementation is given in the previous section.

$\bullet$ For an intrinsically log-normal distribution, the
    efficiency of the skewness and AD tests to reject that the distribution
    is a Gaussian one, is similar. The minimum number of bins required in a
    lightcurve to reject Gaussianity is nearly equal for both the tests.

$\bullet$ It is important to ascertain the nature of the power spectrum of the
lightcurve by estimating the index $\beta$ of the power-law describing the
spectrum. If $\beta \sim 0$ (i.e. white noise) then binning changes the
distribution and since all lightcurves are binned to some extent, estimating
its flux distribution is meaningless. If $\beta \gtrsim 1.5$, then the
variance of the skewness is large and independent of n, and again the skewness and AD tests are not
meaningful.

$\bullet$ If $\beta \sim 0.5$ or $\sim 1$, then we have provided empirical
fits to the standard deviations of the skewness and tabulated critical values for the AD test.
These are different from the ones for white noise and we computed them using
simulations since we could not find these values in the Literature. The
skewness standard deviation and critical values can be used to determine significance by
which the tests reject that the distribution is normal.

$\bullet$ The minimum number of bins required to show that a log-normal distribution is not consistent with a Gaussian one, $N_{90}$ has been computed
as a function of normalised standard deviation and shown in Figures \ref{fig:skew_frac} and  \ref{fig:skew_frac_error} for
different power spectrum index $\beta$ and levels of measurement error.
These can be used {\it a priori} to indicate whether a particular
observation plan will produce lightcurves where these tests can be applied.
Moreover, if a lightcurve having larger number of points shows that its
flux distribution is consistent with a Gaussian one, then one can say with
some confidence that its distribution is not a log-normal one.

$\bullet$ The presence of Gaussian measurement errors increases $N_{90}$ and
we recommend that the error ratio, R  which is the ratio of
the standard deviation of error to the intrinsic standard deviation of flux, to be $ < 0.2$. We suggest that
the data should be binned in time till $R < 0.2$. One needs to take care while constructing a histogram for characterising a PDF shape. The histogram bin width should be chosen such that the width of bin is larger than the measurement error in that bin. Smaller histogram bin-width will introduce bias in the shape of PDF, as the flux point with large error can not be located in single bin. Larger bin-width in the histogram will be required for the larger measurement errors, which some times produces histogram bins with empty counts. Therefore, in order to reduce the measurement error, one should bin the data till the error ratio, $R$ is optimised.  Our work suggest $R< 0.2$.

$\bullet$ For a system having an intrinsic log-normal distribution, the skewness and AD test will accept the normality of the logarithm of the flux at more than 95\% confidence level. However if there is significant Gaussian measurement error in the flux, the skewness and AD test may not accept the normality of the logarithm of the flux at 95\% confidence level. We show that for an error fraction of $R\lesssim0.2$, this would require more than $10^5$ data points. Thus for most practical considerations this effect would be negligible.

$\bullet$ The average skewness of the flux of a log normal distribution
differs from its theoretical value for finite number of bins. Moreover, its
variance is large. Hence, We do not recommend that the value of the skewness
of the flux to be used to show that the distribution is a log-normal one.

This analysis is pertinent to the case when one has a lightcurve with
moderate number of bins say of the order of hundreds. For a larger data set,
one can directly fit the histogram of the flux (or the log of the flux) with
different distributions to ascertain its nature. Naturally, for a data
set with moderate number of bins, the power spectra estimate will not be
precise and the power spectra index $\beta$ would only be known
approximately. That is the reason why we have only considered
$\beta \sim 0.5$
and $1.0$, and not undertaken numerically expensive analysis for finer values
of $\beta$. This implies that any result obtained from the skewness and AD
tests should come with the caveat that they depend on the uncertainty of
estimating the power spectrum. For similar reasons, we have not considered
the more generic situation when the power spectrum can be described as a
broken power-law. We note that the variance of the skewness and the
AD critical values listed in this work are for power-law power spectra
and are not applicable for a broken power-law. For typical cases, where
a flat power spectrum ( e.g. $\beta \sim 0.5$) breaks into a
steeper one (e.g. $\beta \sim 1.5$) beyond a break frequency, we recommend
that the lightcurve be rebinned to a bin size corresponding to the
break frequency and then apply the skewness and AD tests. Since the
power-law index is typically steep beyond the break frequency, rebinning will
not lead to significant loss of information and in general will decrease
the measurement noise level.

As an important aside, we have addressed a conceptually paradoxical
issue regarding non-Gaussian flux distributions.  While Gaussian distributions
remain Gaussian on addition, i.e. Gaussian flux distribution remain Gaussian
when the lightcurve is rebinned, this is in general not true for other
 distributions such as a log-normal one. Thus, it is not clear
why several astrophysical systems display log-normal distributions when
the lightcurves they are estimated from, are always binned to some
time-bin? We resolve this issue by showing that if the power spectra
of the lightcurves can be described by a power law with index $\beta \gtrsim 0.5$,
then the nature of the distribution is invariant to rebinning. Thus for
a system to be scale free in time and have a log-normal flux distribution its
power spectra must be steeper than 0.5. This insight may have important
consequences for the model development and understanding of systems with
log-normal distributions.

Understanding the nature of flux distribution using first analysis tools
like skewness and AD that can be applied to moderate size data sets, is an
important step laying the foundation for more elaborate studies. The results
presented here provide a prescription and guidelines for the effective
and reliable use.

\section{ACKNOWLEDGEMENTS}
We thank the anonymous referee for his critical assessment of our work. The referee comments and suggestions have helped us in improving this manuscript. We thank a UGC-UKIERI Thematic Partnership for support. We acknowledge the use of Fermi-LAT
data provided by Fermi Science Support Center (FSSC).

\section{Data availability}
The Python codes and data used in this article will be shared on reasonable request to the corresponding author, Zahir Shah (email: zahir@iucaa.in or shahzahir4@gmail.com).

\bibliographystyle{mnras}


\section{Introduction}
The recent detections of gravitational waves by the LIGO-Virgo collaboration, in particular, the multimessenger binary-neutron star (BNS) merger event GW170817~\cite{Abbott2017}, motivate studies of the transport properties of dense nuclear matter at temperatures and densities relevant to BNS
mergers~\cite{Alford2018a,Harutyunyan:2016a,Harutyunyan:2016b,Harutyunyan2018b,Alford2019a,Alford2019b}. The mass of the post-merger object typically would exceed the maximum mass of a neutron star and, as a consequence, it would collapse to a black hole on the timescales ranging from tens of milliseconds up to seconds depending on the mass of the post-merger object~\cite{Faber2012:lrr,Baiotti2017,Baiotti2019}. While gravitational waves in the post-merger phase have not been observed in the GW170817 event due to lack of detector sensitivity at high frequencies,
numerical relativity studies of BNS mergers in their highly non-linear regime predict intense emission of gravitational waves in the kHz frequency range during the initial phase of post-merger phase lasting typically 10 ms~(for recent simulations see, for example, \cite{Endrizzi2018,Most2019,Ciolfi2019,Tsokaros2019}). Improvements at least by a factor of 2 compared to advanced LIGO  design sensitivity are necessary to measure the dominant frequency component of the signal of GW170817-like event and by factors of 4-5 to observe sub-dominant features of post-merger signal~\cite{Torres-Rivas2019}.
It is expected that dissipation or damping of matter flows in the merged stars could influence the gravitational waves emitted during the post-merger phase. Recent estimates of the role of the thermal conduction and shear and bulk viscosities indicate that, of these, damping of density oscillations via bulk viscosity has the strongest
influence~\cite{Alford2018a,Alford2019a,Alford2019b}. After a brief introduction to the problem of computation of the bulk viscosity, we present an extension of our recent work~\cite{Alford2019b} on the bulk viscosity of nucleonic matter which includes the estimates of timescales of the damping of the oscillations by the bulk viscosity. The relevance of the bulk viscosity will be assessed by comparing the damping timescales of density oscillations to the characteristic timescales of the initial phase of post-merger
 $\sim 10$ ms (over which the post-merger object is expected to emit intense gravitational waves) as well as the longer-term phase  $\sim 1$ s.

The bulk viscosity of nuclear matter at temperatures up to about 1\,MeV has been studied extensively~\cite{Sawyer1979ApJ,Sawyer1980ApJ,Sawyer1989,Haensel1992PhRvD,Haensel:2000vz,Haensel:2001mw,Haensel:2001em,Dong2007,Alford2010JPhG,Alford:2010jf,Kolomeitsev2015} in the context of oscillations of neutron stars and, in particular, as a source of damping of (unstable) $r$-mode oscillations; for reviews see~\cite{Schmitt2017,Kokkotas2016EPJA}. More recently, interest in neutron star mergers has motivated studies of bulk viscosity of nucleonic matter at temperatures up to several tens of MeV, covering both the regime where neutrinos escape from neutron stars and the regime where they are trapped~\cite{Alford2019a,Alford2019b}. Such high temperatures significantly affect the phase space occupation of the fermions and, therefore, the rates of the weak-interaction processes \cite{Alford2018b}. Once the temperature becomes high enough for neutrinos to be trapped, they affect the composition of matter and ensure that direct Urca processes are always kinematically possible, and modified Urca is a subleading correction.

In this contribution, we extend recent work~\cite{Alford2019b} to compute the timescales associated with bulk-viscous damping of oscillations. We use two models of the equation of state (EoS) and associated composition of dense matter which are based on the relativistic density functional theory of nuclear matter. The overall picture is that, at densities from around $n_0$ (nuclear saturation density) to around $3 n_0$, the bulk viscosity reaches a maximum at temperature $T\simeq 2\div 6$\,MeV, which is in the regime where neutrinos are not (or not completely) trapped. At these temperatures, therefore, the  damping timescale is at a minimum, with values ranging from about 10 ms at low densities $n_B\leq n_0$ down to milliseconds (or even tenths of milliseconds depending on the EoS)
 at $n_B\simeq 3n_0$. This means that bulk viscous damping can have noticeable effects during the $\sim 10$\,ms of initial  (gravitational-wave-emission) phase of the post-merger. At higher temperatures where the neutrino-trapping occurs the bulk viscosity falls by orders of magnitude, which implies much longer damping timescales which are larger than the characteristic timescales involved.

This paper is organized as follows. In~Section~\ref{sec:bulk} we discuss the
formalism for computing the bulk viscosity and the approximations involved in such a computation. Our focus is on the beta equilibration processes of neutron decay and electron capture and the microscopic relaxation rates associated with these processes. Section~\ref{sec:results} starts with a brief discussion of the thermodynamic properties of nuclear matter derived from density functional theory in Section~\ref{sec:DFT}. Bulk viscosity and the oscillation damping timescale are discussed in Sections~\ref{sec:bulk_visc} and \ref{sec:damping}. Our main results are summarized in Section~\ref{sec:summary}. We use natural units ($\hbar = c = k_B = 1$) and the metric signature $g_{\mu\nu} = \textrm{diag}(1, -1, -1, -1)$.

\section{Urca processes and bulk viscosity}
\label{sec:bulk}

We start with a brief reminder of the bulk viscosity of nuclear matter
composed of neutrons, protons, electrons, muons; for more details see~\cite{Alford2019a,Alford2019b}. For~simplicity, we will neglect the muonic contribution to equilibration rates, although~we include their contribution to static thermodynamic quantities such as~susceptibilities.

Above the trapping temperature $T_{\rm tr}\simeq 5$\,MeV, the neutrino mean-free-path is smaller than the size of a neutron star, so neutrinos are trapped in the merger region.  Under these conditions beta equilibrium is established via neutron decay and electron capture, and their inverse processes
\begin{eqnarray}\label{eq:Urca_1}
n\rightleftarrows p + e^-+\bar{\nu}_e \ ,\\
\label{eq:Urca_2}
p + e^- \rightleftarrows n+\nu_e \ .
\end{eqnarray}

In $\beta$-equilibrium the chemical potentials of particles obey the relation
\begin{eqnarray}\label{eq:eq_condition}
\mu_p+\mu_e = \mu_n+\mu_{\nu} \ .
\end{eqnarray}
The particle fractions of baryonic matter for any given temperature $T$, baryon number density
$n_B=n_n+n_p$  and lepton densities $n_{L_l}=n_l+n_{\nu_l}=Y_{L_l}n_B$
(since we ignore muon reactions, the lepton fractions $Y_{L_l}$ need to be fixed for each flavor separately) are found by imposing the beta-equilibrium
condition~\eqref{eq:eq_condition}, a similar condition for muons $\mu_p+\mu_\mu = \mu_n+\mu_{\nu_\mu}$ and the charge neutrality condition $n_p=n_e+n_\mu$.

At lower temperatures $T\lesssim T_{\rm tr}$ the neutrino mean free path is larger than the
size of a neutron star so the system is neutrino-transparent. Neutrinos cannot occur in initial
states, therefore  the reactions~\eqref{eq:Urca_1} and \eqref{eq:Urca_2} proceed only in one direction (from left to right).
To determine the composition of matter in this case we use
the ordinary zero-temperature $\beta$-equilibrium conditions
$\mu_n=\mu_p+\mu_e$ and $\mu_\mu=\mu_e$.  Reference~\cite{Alford2018b}
found that there are significant corrections to these conditions
at $T\gtrsim 1$\,MeV; nevertheless the bulk viscosity in the neutrino transparent regime is not affected significantly~\cite{Alford2019a}.

If the matter is driven out of $\beta$-equilibrium, for~example by
compression and rarefaction, the~left and right-hand sides of
Equation~\eqref{eq:eq_condition} do not balance anymore. The~
deviation from $\beta$-equilibrium is then measured by a quantity
\begin{eqnarray}\label{eq:mu_Delta}
\mu_\Delta = \mu_n+\mu_{\nu}-\mu_p-\mu_e.
\end{eqnarray}
As a result, the Urca processes~\eqref{eq:Urca_1} and
\eqref{eq:Urca_2} will go faster in one direction than in the other
until the beta equilibrium of matter is restored.

Consider now small-amplitude density oscillations
in nuclear matter with a frequency $\omega$.
 The baryon and lepton conservation implies for periodic
 perturbations $\delta n_B(t),\delta n_L(t) \sim e^{i\omega t}$
\begin{eqnarray}\label{eq:cont_i}
\delta n_i(t) =-\frac{\theta}{i\omega}\, n_{i0},
\quad i=\{B,L\},
\end{eqnarray}
where $n_{B0}=n_{n0}+n_{p0}$ and $n_{L0}=n_{e0}+n_{\nu 0}$
are the unperturbed background densities of baryons and leptons,
and $\theta$ is the divergence of fluid velocity.
The compression and rarefaction of matter implies
perturbations in particle densities which can be
separated into instantaneous equilibrium and
non-equilibrium parts
\begin{eqnarray}\label{eq:dens_j}
n_j(t)=n_{j0}+\delta n_j(t), \quad \delta n_j(t)=
\delta n^{\rm eq}_j(t)+\delta n'_j(t),\quad j=\{n,p,e,\nu\},
\end{eqnarray}
where $n_{j0}$ are the static values of particle densities. The
variations $\delta n^{\rm eq}_j(t)$ stand for the shifts of the
equilibrium state for the instantaneous values of $n_{B}(t)$ and
$n_{L}(t)$, whereas $\delta n'_j(t)$ are the deviations of the
particle densities from those equilibrium values. There exist
two choices of the instantaneous equilibrium state. Below, we
follow our recent work~\cite{Alford2019b}, for the alternative
see Ref.~\cite{Huang2010}. We compare below these two approaches
and explain why they give the same result for the bulk viscosity.

The non-equilibrium perturbations $\delta n'_j(t)$ drag
matter out of chemical equilibrium by leading to a small
chemical potential shift~\eqref{eq:mu_Delta} which can
be written in terms of particle densities as
\begin{eqnarray}\label{eq:delta_mu}
\mu_\Delta(t)=A_{n}\delta n_n(t)+A_{\nu} \delta n_\nu(t) -
A_{p} \delta n_p(t) -A_{e}\delta n_e(t),
\end{eqnarray}
where $A_n=A_{nn}-A_{pn}$, $A_p=A_{pp}-A_{np}$,
and $A_e=A_{ee}$, $A_\nu=A_{\nu\nu}$ with
\begin{eqnarray}\label{eq:A_j}
A_{ij} = \frac{\partial \mu_i}{\partial n_j}.
\end{eqnarray}
The off-diagonal elements $A_{np}$ and $A_{pn}$ are non-zero
because of the cross-species strong interaction between neutrons
and protons. For small amplitude density oscillations we only need to evaluate the derivatives in Equation~\eqref{eq:A_j} at $\mu_\Delta=0$.

If there were no flavor-changing weak processes,
then the particle densities would just oscillate
around their static equilibrium values according to
\begin{eqnarray}\label{eq:cont_j}
\frac{\partial}{\partial t} \delta {n}_j^0(t)=
-\theta n_{j0}\quad\Rightarrow\quad
\delta {n}_j^0(t) = -\frac{\theta}{i\omega}\, n_{j0}.
\end{eqnarray}
The weak interactions lead to an imbalance between
the rates of direct and inverse Urca processes which
in the ``subthermal regime'' ($\mu_\Delta\ll T$) can be written
as~\cite{Sawyer1989,Haensel1992PhRvD,Huang2010}
\begin{eqnarray}\label{eq:lambda_def}
\Gamma_\Delta \equiv \Gamma_p-\Gamma_n =\lambda\mu_\Delta,\quad\lambda > 0,
\end{eqnarray}
with $\Gamma_p$ and $\Gamma_n$ being the production rates of protons
and neutrons, respectively. The~production rate~\eqref{eq:lambda_def}
should be added to the right hand sides of Equation~\eqref{eq:cont_j} with
a plus sign for $p,e$ and a minus sign for $n,\nu$, for~example,
\begin{eqnarray}\label{eq:cont_n_weak}
\frac{\partial}{\partial t}\delta n_n(t) =
-\theta  n_{n0} -\lambda\mu_\Delta(t),\quad etc.
\end{eqnarray}

Substituting here Equation~\eqref{eq:delta_mu}, exploiting the
conditions $\delta n_B =\delta n_n +\delta n_p$,
$\delta n_p=\delta n_e+\delta n_\mu$ and
$\delta n_L =\delta n_e+\delta n_\nu$, and using Equation~\eqref{eq:cont_i}
for $n_B$, $n_L$ and an analogous equation for $n_\mu$
(as muons are assumed not to participate in any reactions, their
fraction is conserved) one finds
\begin{eqnarray}\label{eq:delta_n}
\delta n_n &=& -\frac{i\omega n_{n0} +\lambda
  (A_{p}+A_e+A_\nu)n_{B0}
-\lambda A_\nu n_{L0}-\lambda(A_e+A_\nu)n_{\mu 0}}
{i\omega+\lambda A}\frac{\theta}{i\omega},\\
\label{eq:delta_p}
\delta n_p &=& -\frac{i\omega n_{p0}
+\lambda A_{n}n_{B0}+\lambda A_\nu n_{L0}+\lambda(A_e+A_\nu)n_{\mu 0}}
{i\omega+\lambda A}\frac{\theta}{i\omega},\\
\label{eq:delta_e}
\delta n_e &=& -\frac{i\omega n_{e0}
+\lambda A_{n}n_{B0}+\lambda A_\nu n_{L0}-\lambda(A_n+A_p)n_{\mu 0}}
{i\omega+\lambda A}\frac{\theta}{i\omega},\\
\label{eq:delta_nu}
\delta n_\nu &=& -\frac{i\omega n_{\nu 0}+\lambda
  (A_{n}+A_{p}+A_e)n_{L0} -\lambda A_{n} n_{B0}
  +\lambda(A_n+A_p)n_{\mu 0}}{i\omega+\lambda A}\frac{\theta}{i\omega},
\end{eqnarray}
with a ``beta-disequilibrium--proton-fraction'' susceptibility given by
\begin{eqnarray}\label{eq:A_def}
A=\sum_{i} A_{i}
=\frac{\partial\mu_n}{\partial n_n}
+\frac{\partial\mu_p}{\partial n_p}
-\frac{\partial\mu_n}{\partial n_p}
-\frac{\partial\mu_p}{\partial n_n}
+\frac{\partial\mu_e}{\partial n_e}
+\frac{\partial\mu_\nu}{\partial n_\nu}
= \biggl(\dfrac{\partial \mu_\Delta}
{\partial n_n} \biggr)_{\!n_B}.
\end{eqnarray}
Equations~\eqref{eq:delta_n}--\eqref{eq:delta_nu}
are the extensions of Eqs.~(37)--(39) of Ref.~\cite{Alford2019b}
as we included non-zero muon density $n_\mu$ here, which
was previously neglected. The final formula for the bulk
viscosity, however, remains the same after this addition.

In the next step we find $\delta n_j^{\rm eq}$ using the
definition of the instantaneous $\beta$-equilibrium state:
$A_{n}\delta n^{\rm eq}_n+A_\nu\delta n^{\rm eq}_\nu-
A_{p}\delta n^{\rm eq}_p-A_e\delta n^{\rm eq}_e=0$, which gives
\begin{eqnarray}\label{eq:delta_n_eq}
\delta n^{\rm eq}_n &=& \frac{-(A_{p}+A_e+A_\nu)n_{B0}
+A_\nu n_{L0} +(A_e+A_\nu)n_{\mu 0}}{A} \frac{\theta}{i\omega},\\
\label{eq:delta_p_eq}
\delta n^{\rm eq}_p &=& -\frac{A_{n} n_{B0}
+ A_\nu n_{L0} +(A_e+A_\nu)n_{\mu 0}}{A}\frac{\theta}{i\omega},\\
\label{eq:delta_e_eq}
\delta n^{\rm eq}_e &=& -\frac{A_{n} n_{B0}
+ A_\nu n_{L0} -(A_n+A_p)n_{\mu 0}}{A}\frac{\theta}{i\omega},\\
\label{eq:delta_nu_eq}
\delta n^{\rm eq}_\nu &=& \frac{-(A_{n}+A_{p}+A_e)n_{L0}
+A_{n} n_{B0}-(A_n+A_p)n_{\mu 0}}{A} \frac{\theta}{i\omega}.
\end{eqnarray}
Note that these expressions are the solutions of the balance
equations~\eqref{eq:cont_n_weak} in the limit of infinite
relaxation rate $\lambda\to\infty$, which implies necessarily
$\mu_\Delta\to 0$. One can check this also by pushing
$\lambda$ to infinity directly in general solutions
\eqref{eq:delta_n}--\eqref{eq:delta_nu} which will reduce
then to Equations~\eqref{eq:delta_n_eq}--\eqref{eq:delta_nu_eq}.

Now the non-equilibrium density perturbations can be
found according to Equation~\eqref{eq:dens_j}
\begin{eqnarray}\label{eq:delta_j'}
\delta n'_p = \delta n'_e =-\delta n'_n = -\delta n'_\nu
  =\frac{ C}{A(i\omega+\gamma)}\theta,
\end{eqnarray}
where $\gamma=\lambda A$ has a dimension of frequency and
measures the relaxation rate of particle densities to their
equilibrium values, and
\begin{eqnarray}\label{eq:C_def}
C = n_{n0} A_{n}+n_{\nu 0} A_\nu - n_{p0}A_{p} - n_{e0}A_e
= n_B\biggl(\dfrac{\partial \mu_\Delta}{\partial n_B}
\biggr)_{\!Y_n}
\end{eqnarray}
is the ``beta-disequilibrium--baryon-density'' susceptibility,
with $Y_n=n_n/n_B$ being the neutron fraction.
The non-equilibrium part of the pressure -- the so-called
bulk viscous pressure, can be now computed as
\begin{eqnarray}\label{eq:Pi}
\Pi =\sum_j
\frac{\partial p}{\partial n_j}\delta n'_j
=\sum_{lj} n_{l0}A_{lj}\delta n'_j,
\end{eqnarray}
where we used the Gibbs-Duhem relation $dp=sdT+\sum_l n_l d \mu_l$.
Substituting the solutions~\eqref{eq:delta_j'} in Equation~\eqref{eq:Pi}
one finds the bulk viscosity from the definition $\Pi=-\zeta\theta$
\begin{eqnarray}\label{eq:zeta}
\zeta = \frac{C^2}{A}\frac{\gamma}{\omega^2+\gamma^2}.
\end{eqnarray}
The susceptibility prefactor $C^2/A$ is a pure
thermodynamic quantity and depends only on the EoS,
whereas the relaxation rate $\gamma=\lambda A$
[Equation~\eqref{eq:lambda_def}] depends on the
microscopic scattering amplitudes of weak interactions. It is seen from Equations~\eqref{eq:lambda_def} and
\eqref{eq:A_def} that $\gamma$ is actually the derivative
\begin{eqnarray}\label{eq:gamma}
\gamma =\biggl(\dfrac{\partial \Gamma_\Delta}
{\partial \mu_\Delta} \biggr)
\biggl(\dfrac{\partial \mu_\Delta}
{\partial n_n} \biggr)_{\!n_B}
= \biggl(\dfrac{\partial \Gamma_\Delta}
{\partial n_n} \biggr)_{\!n_B} \ .
\end{eqnarray}
Thus, the quantity $\gamma$ measures how the proton net
production rate increases when the neutron fraction
increases at fixed baryon density, {\it i.e.}, $\gamma$
measures how fast the system reacts to a change in the
chemical composition of matter. The quantity $\gamma^{-1}$
has a dimension of time and can be interpreted as a
relaxation time of the system to its beta-equilibrium state.

The bulk viscosity in Equation~\eqref{eq:zeta} has the
classic resonant form and for density oscillations of a
given frequency, $\omega$ attains a resonant maximum at
the temperature where the relaxation rate matches the oscillation
frequency, $\gamma(T)=\omega$.

The value of the bulk viscosity at that maximum is
\begin{eqnarray}\label{eq:zeta_max}
\zeta_{\rm max}=\frac{C^2}{2A\omega} \ ,
\end{eqnarray}
which is independent of the microscopic interaction rates.

In the regime of slow equilibration, where $\gamma\ll\omega$,
the bulk viscosity takes the form
\begin{eqnarray}\label{eq:zeta_slow}
\zeta_{\rm slow} = \frac{C^2}{A}
\frac{\gamma}{\omega^2} {= \dfrac{2 \gamma}{\omega}\, \zeta_{\rm max}}\, .
\end{eqnarray}
In the fast equilibration regime $\gamma\gg\omega$,
and the bulk viscosity, in this case, reduces to
\begin{eqnarray}\label{eq:zeta_fast}
\zeta_{\rm fast} = \frac{C^2}{A\gamma}
{= \dfrac{2\omega}{\gamma} \zeta_{\rm max} }\, .
\end{eqnarray}
The physical reason for this resonant maximum is easy to understand.
In the limit where the relaxation rate is much smaller than the oscillation frequency
there would effectively be an additional conserved quantity since both proton number and neutron number would be conserved. The proton fraction would be independent of density, so density oscillations would not drive the system out of chemical equilibrium. There would then be no bulk viscosity.

In the opposite limit of fast equilibration or slow density
oscillations, $\gamma\gg\omega$, weak interactions are able
to restore the chemical equilibrium of matter on timescales much
smaller than the oscillation period. This means that the matter is
practically always beta-equilibrated
while undergoing compression and rarefaction and,
therefore, will not experience any bulk viscosity.

Above we computed the bulk viscous pressure using its standard definition
\begin{eqnarray}\label{eq:Pi_inf}
\Pi =\delta P -\delta P_{\rm eq}
=\delta P(\lambda) -\delta P(\lambda\to\infty),
\end{eqnarray}
where $\delta P$ is the shift of the pressure from its
static equilibrium value $P_0=P(n_{j0})$ for arbitrary
perturbations $\delta n_{j}$, and $\delta P_{\rm eq}$
is the instantaneous shift of the equilibrium pressure,
which depends on $\delta n_{j}^{\rm eq}$ and does not
contribute to the bulk viscous pressure.

According to the two limiting cases of vanishing bulk
viscosity discussed above the bulk viscous pressure
can be defined also in an alternative way
\begin{eqnarray}\label{eq:Pi_0}
\Pi =\delta P(\lambda) -\delta P(\lambda\to 0)
=\delta P -\delta P_{0},
\end{eqnarray}
where $\delta P_{0}$ is the instantaneous shift of the pressure
in a certain non-equilibrium state with conserved particle fractions
which corresponds to the limit $\lambda\to 0$.
According to Equation~\eqref{eq:Pi_0} we can take as alternatives
of the beta-equilibrium shifts $\delta n^{\rm eq}_j(t)$ the shifts $\delta n^0_j(t)$ given by Equation~\eqref{eq:cont_j}, as these are the solutions
of exact balance equations~\eqref{eq:cont_n_weak} in the limit
$\lambda\to 0$. This was just the choice of the equilibrium
state in Ref.~\cite{Huang2010}. Then instead of
Equation~\eqref{eq:delta_j'} we will have
\begin{eqnarray}\label{eq:delta_j'_alt}
\delta n'_n(t) = \delta n'_\nu(t)=
-\delta n'_p(t) = -\delta n'_e(t)=
\frac{C}{A(i\omega+\gamma)}\frac{\gamma\theta}{i\omega}.
\end{eqnarray}
Computing the bulk viscous pressure from Equation~\eqref{eq:Pi}
we obtain the same result for $\zeta$, as expected.

\subsection{Urca process rates}
\label{sec:urca_rates}

As it was shown in Ref.~\cite{Alford2019b} the equilibration rate
$\lambda$ of the neutrino-trapped matter is dominated by the electron
capture process because the neutron decay rate is exponentially
damped at low temperatures as $\sim\exp(-\mu{_\nu}/T)$.
The microscopic $\beta$-equilibration rate of the electron
capture process is given by~\cite{Alford2019b}
\begin{eqnarray} \label{eq:Gammap_def}
\Gamma_{n}(\mu_\Delta) = 2\tilde G^2
\prod_j \int\!\! \frac{d^3p_j}{(2\pi)^3}
\bar{f}(p_1){f}(p_2){f}(p_3)\bar{f}(p_4)\,
(2\pi)^4\delta(p_1-p_2-p_3+p_4),
\end{eqnarray}
where $\tilde{G}^2\equiv G_F^2\cos^2 \theta_c (1+3g_A^2)$,
$G_F=1.166\cdot 10 ^{-5}$ GeV$^{-2}$ is the Fermi coupling
constant, $\theta_c$ is the Cabibbo angle with
$\cos\theta_c=0.974$, and $g_A=1.26$ is the axial-vector
coupling constant; the index $j$ runs over the four
participating particles $j=\{n,p,e,\nu\}$, $f(p_j)$ are the
Fermi distributions of particles and $\bar{f}(p_j)=1-f(p_j)$.
The inverse process $\Gamma_p$ will be given by an analogous
expression by replacing all $f_j$ functions with $\bar{f_j}$.

For small departures from $\beta$-equilibrium $\mu_\Delta
\ll T$ the imbalance between the direct and inverse
rates can be linearized in $\mu_\Delta$ with the
coefficient of the linear expansion given by~\cite{Alford2019b}
\begin{eqnarray}\label{eq:lambda_final}
\lambda = \frac{ m^{*2}\tilde G^2 T^5}{8\pi^5}
\!\int_{-\infty}^{\infty} dy~g(y)\!
\int_{z_0}^{\infty}  dz \, {\cal L}(y,z)
\int_{x_0}^{\infty}  dx\, (x+\alpha_{\nu})
(y+\alpha_e+x)f(x)\bar f(x+y),
\end{eqnarray}
where $m^\ast$ is the effective nucleon mass, $\alpha_j = \mu_j^*/T$,
$\mu_j^*$ are the effective chemical potentials of particles,
$f(x)=(e^x+1)^{-1}$ and $g(x)=(e^x-1)^{-1}$ are the Fermi
and the Bose distributions, respectively,
$x_0=(z-y-\alpha_e-\alpha_{\nu})/2$,
$z_0=|y+\alpha_e-\alpha_{\nu}|$, and
\begin{eqnarray}\label{eq:log_y_0}
 {\cal L}(y,z) = \ln \Bigg\vert \frac{1+\exp(-y_0)}
 {1+\exp\left(-y_0-y\right)}\Bigg\vert,
\quad
y_0 = \frac{m^*}{2Tz^2}\bigg(\alpha_n-\alpha_p+y
-z^2\frac{T}{2m^*}\bigg)^2-\alpha_p.
\end{eqnarray}

In the case of neutrino-transparent matter the equilibration rate
$\lambda$ is given by a similar expression, see Refs.~\cite{Alford2019a,Alford2019b}
for details. The low-temperature
limit of $\lambda$ in neutrino-trapped matter is given by~\cite{Alford2019b}
\begin{eqnarray} \label{eq:lambda_deg_trap}
\lambda = \frac{1}{12\pi^3}m^{*2}\tilde G^2 T^2
\, p_{Fe}p_{F\nu} (p_{Fp}+p_{Fe}+p_{F\nu}-p_{Fn}),
\end{eqnarray}
and in the neutrino-transparent matter~\cite{Haensel:2000vz,Alford2019a,Alford2019b}
\begin{eqnarray} \label{eq:lambda_deg_trans}
\lambda = \frac{17}{240\pi} m^{*2} \tilde{G}^2
T^4 \, p_{Fe}\,\theta( p_{Fp}+p_{Fe} -p_{Fn}),
\end{eqnarray}
where $p_{Fj}$ are the Fermi-momenta of the particles.
The $\theta$-function in Equation~\eqref{eq:lambda_deg_trans}
blocks the direct Urca processes at low densities
where the proton and electron Fermi momenta are not
sufficiently large to guarantee the momentum conservation.
In the case where neutrinos are trapped in matter the
momentum conservation can be always satisfied for
certain particle momenta and the rate is always finite.

\section{Numerical results}
\label{sec:results}

We start with the thermodynamics of nuclear matter, which is derived
from the relativistic density functional theory based on
phenomenological baryon-meson Lagrangians, for reviews
see~\cite{Weber_book,Sedrakian2007PrPNP}. We use the density-dependent
baryon-meson coupling model DDME2~\cite{Lalazissis2005} as applied to
finite-temperature nucleonic matter, see for details
Ref.~\cite{Colucci2013,Li_PLB_2018,Li2019ApJ}.

\subsection{Thermodynamics of nuclear matter}
\label{sec:DFT}

The Lagrangian density of the model considered reads
\begin{eqnarray}\label{eq:lagrangian}
{\cal L} & = &
\sum_N\bar\psi_N\bigg[\gamma^\mu \left(i\partial_\mu-g_{\omega}
\omega_\mu - \frac1 2 g_{\rho}{\bm \tau}\cdot{\bm \rho}_\mu\right)
- (m_N - g_{\sigma}\sigma)\bigg]\psi_N
+ \sum_{\lambda}\bar\psi_\lambda(i\gamma^\mu\partial_\mu -
m_\lambda)\psi_\lambda\nonumber\\
 & + & \frac1 2 \partial^\mu\sigma\partial_\mu\sigma-\frac1 2
m_\sigma^2\sigma^2 - \frac{1}{4}\omega^{\mu\nu}\omega_{\mu\nu} + \frac1 2
m_\omega^2\omega^\mu\omega_\mu -
\frac{1}{4}{\bm \rho}^{\mu\nu}{\bm \rho}_{\mu\nu} + \frac1 2
m_\rho^2{\bm \rho}^\mu\cdot{\bm \rho}_\mu,
\end{eqnarray}
where $N$ sums over nucleons, $\lambda$ -- over the leptons,
and $\psi_i$ are the fermionic Dirac fields with masses $m_i$.
The meson fields $\sigma,\omega_\mu$, and ${\bm \rho}_\mu$ are the
effective mediators of strong interaction between baryons,
$\omega_{\mu\nu}=\partial_\mu\omega_\nu-\partial_\nu\omega_\mu$
and ${\bm \rho}_{\mu\nu}=\partial_\mu{\bm \rho}_\nu-\partial_\nu{\bm \rho}_\mu$
are the field strength tensors of $\omega_\mu$ and ${\bm \rho}_\mu$
mesons, respectively, $m_{\sigma}$, $m_{\omega}$, and $m_{\rho}$ are
meson masses, and $g_{i}$ are the baryon-meson coupling constants
with $i=\sigma,\omega,\rho$.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.45\columnwidth, keepaspectratio]{fractions_trap0.pdf}
\hskip 0.5cm
\includegraphics[width=0.45\columnwidth, keepaspectratio]{fractions_trap1.pdf}
\caption{ Particle fractions as functions of the baryon
density for (a) neutrino-transparent matter ($Y_\nu=0$) at
$T=1$\,MeV; (b) neutrino-trapped matter at $Y_L=0.1$ and $T=50$\,MeV. }
\label{fig:fractions}
\end{center}
\end{figure}
\begin{figure}[!]
\begin{center}
\includegraphics[width=0.45\columnwidth,keepaspectratio]{C2_to_A_vs_dens.pdf}
\caption{ The susceptibility prefactor $C^2/A$ as a function
  of the baryon density for three values of the temperature
  for the neutrino-transparent matter (upper curves labelled as $\nu$-trans.) and neutrino-trapped matter (lower curves labelled as $\nu$-trap.).
  The upper and lower triangles show the results of Ref.~\cite{Alford2019a} for models IUFSU and DD2, respectively.}
\label{fig:C2A_dens}
\end{center}
\end{figure}
The density dependence of the particle fractions $Y_j = n_j/n_B$ is
shown in Fig.~\ref{fig:fractions}. The left panel is for a low
temperature $T=1$\,MeV where the system is neutrino-transparent,
and the right panel shows the results for neutrino-trapped matter at
temperature $T=50$\,MeV for lepton fraction fixed at $Y_L=0.1$.  In the
first case muons appear only above a certain baryon density
$n_B\gtrsim n_0$, where the condition $\mu_e\geq m_\mu\simeq 106$\,MeV
is satisfied, whereas in the neutrino-trapped case the muons' threshold
disappears.

Within the framework of the model above the susceptibilities $A$ and $C$
in the non-relativistic limit for nucleons are given by (see Ref.~\cite{Alford2019b} for details)
\begin{eqnarray}\label{eq:A_fianl}
A &=& \frac{1}{{I}_{2n}}+\frac{1}{{I}_{2p}}+
\frac{1}{{I}_{2e}}+\frac{2}{{I}_{2\nu}}
+\left(\frac{g_\rho}{m_\rho}\right)^2,\\
\label{eq:C_fianl}
C &=& \frac{n_n}{{I}_{2n}}-\frac{n_p}{{I}_{2p}}-\frac{n_e}{{I}_{2e}}
+2\frac{n_\nu}{{I}_{2\nu}}-g_\rho \rho_{03}
+\frac{g_\sigma \sigma}{2m^{*2}}
\left(\frac{{I}_{4n}}{{I}_{2n}} - \frac{{I}_{4p}}{{I}_{2p}}\right),
\end{eqnarray}
where
\begin{eqnarray}\label{eq:I_nonrel}
{I}_{q i}= \frac{1}{\pi^2 T}\int_0^\infty p^q dp \,f_i \bar f_i,
\end{eqnarray}
$\sigma$ is the $\sigma$ meson mean field and $\rho_{03}$
is the $\rho$ meson mean field which is non-zero
in asymmetric nuclear matter. These expressions for susceptibilities
are derived for isothermal density perturbations, which is the case
only if the thermal conduction is fast enough to smoothen the
temperature gradients during one period of
oscillation~\cite{Alford2018a,Alford2019b}.
This might happen, {\it e.g.}, in the presence of turbulent flows in the merger region which could generate temperature and density variations on distance scales
of the order of a few hundred meters.
An order-of-magnitude estimate of
thermal relaxation timescale is given
in Ref.~\cite{Alford2018a}, which is
1 sec $\times (z_{\rm typ}/{\rm km})^2 (T/10\, {\rm MeV})^2$, where
$z_{\rm typ}$ is the typical length
scale of thermal gradients. Assuming
$z_{\rm typ}\simeq 100$ m and temperatures $1\div 10$ MeV (this is the temperature range where the bulk viscosity is relevant to mergers, see below), the thermal relaxation time will lie in the interval 0.1-10~ms, which is below the characteristic timescale of binary neutron star mergers. Thus, for thermal gradients on this distance scale the assumption of isothermal matter is the relevant one.
On the scales over which thermal conduction is inefficient the
matter should be treated as iso-entropic. The isothermal and
adiabatic susceptibilities differ at most by a factor of 2 in
the relevant density and temperature range, see
Ref.~\cite{Alford2019a} for further details.

Figure~\ref{fig:C2A_dens} shows the ratio $C^2/A$ of
susceptibilities as a function of density for three values of
the temperature.  The susceptibility $C$ is an increasing function
of density. At sufficiently high temperatures $T\gtrsim 30$\,MeV,
it is negative at low density and crosses zero
at a temperature-dependent critical density
where the proton fraction reaches a minimum as a
function of the density. At that point, the system becomes
scale-invariant, as the compression does not drive the matter
out of beta equilibrium, and the bulk viscosity vanishes at
that critical point.

The ratio $C^2/A$ grows rapidly with the density in both cases of
neutrino-transparent and neutrino-trapped matter and is sensitive to
the temperature only close to the point where $C=0$. We see that
$C^2/A$ is approximately an order of magnitude smaller in the
neutrino-trapped matter the reason being much larger values of $A$
dominated by the contribution of neutrinos as compared to neutrino-transparent matter.

\subsection{Beta relaxation rates and bulk viscosity}
\label{sec:bulk_visc}

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.435\columnwidth,keepaspectratio]{gamma_vs_temp_trap0.pdf}
\hskip 0.5cm
\includegraphics[width=0.45\columnwidth,keepaspectratio]{gamma_vs_temp_trap1.pdf}
\caption{ The beta-relaxation rate $\gamma$ as a function of the
  temperature for (a) neutrino-transparent and (b) neutrino-trapped
  regimes for several values of baryon density. Note the very different $y$-axis scales.
  The horizontal lines in panel (a) correspond to the fixed values of oscillation frequency
  $f=1$ kHz (solid line) and $f=10$ kHz (dashed line).}
\label{fig:gamma_temp}
\end{center}
\end{figure}

Figure~\ref{fig:gamma_temp} shows the relaxation rate $\gamma=\lambda A$
[Equations~\eqref{eq:lambda_def} and \eqref{eq:gamma}] as a
function of temperature for various densities. The equilibration rate
$\lambda$ of the neutrino-trapped matter is dominated by the electron capture
process because the neutron decay rate is exponentially damped at low
temperatures as $\lambda_1\sim\exp(-\mu{_\nu}/T)$, whereas $\lambda_2$
has approximately a quadratic increase with the temperature as
suggested by Equation~\eqref{eq:lambda_deg_trap}. It is
seen from Fig.~\ref{fig:gamma_temp}(b) that the relaxation
rate $\gamma$ of the neutrino-trapped matter is several
orders of magnitude larger than the oscillation frequencies
$f=\omega/2\pi\simeq 1$ kHz which are typical to neutron star
mergers. This means that the neutrino-trapped matter is always in the
fast equilibration regime where the bulk viscosity is independent of
the oscillation frequency and is given by Equation~\eqref{eq:zeta_fast}.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.45\columnwidth,keepaspectratio]{zeta_vs_dens_trap0_DDME2.pdf}
\hskip 0.5cm
\includegraphics[width=0.45\columnwidth,keepaspectratio]{zeta_vs_dens_trap1_DDME2.pdf}
\caption{ The density dependence of the bulk viscosity for
  (a) neutrino-transparent matter; (b) neutrino-trapped matter at
  $Y_L=0.1$ for various values of the temperature. The oscillation frequency is $f=1$ kHz
  in panel (a); the bulk viscosity in panel (b) is frequency independent.
  Note the very different $y$-axis scales.}
\label{fig:zeta_dens}
\end{center}
\end{figure}
\begin{figure}[!]
\begin{center}
\includegraphics[width=0.45\columnwidth, keepaspectratio]{zeta_vs_temp_trap0_DDME2.pdf}
\hskip 0.5cm
\includegraphics[width=0.465\columnwidth, keepaspectratio]{zeta_vs_temp_trap1_DDME2.pdf}
\caption{ The temperature dependence of the bulk viscosity
for several values of the baryon density marked on the plots
for (a) neutrino-transparent matter; (b) neutrino-trapped matter
at $Y_L=0.1$. The oscillation frequency is $f=1$ kHz in panel (a); the bulk viscosity
 in panel (b) is frequency independent.
Note the very different $y$-axis scales.  The shaded region in panel (b)
indicates that the results therein are extrapolations to the regime where they are not applicable.
}
\label{fig:zeta_temp}
\end{center}
\end{figure}

In the neutrino-transparent regime, in contrast, the relaxation rate
is much slower, reaching values in the kHz range where it
can resonate with typical density oscillations in mergers.
We see this in Figure~\ref{fig:gamma_temp}(a) where the
relaxation rate $\gamma$ crosses the $\omega = 2\pi$\,kHz (corresponding to $f=1$\,kHz) line
at temperatures $4\div 5$\,MeV indicating that the neutrino-transparent
matter possesses a resonant maximum at those temperatures, as it was
found also in Refs.~\cite{Alford2018a,Alford2019a}.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.45\columnwidth, keepaspectratio]{zeta_vs_temp_DDME2.pdf}
\hskip 0.5 cm
\includegraphics[width=0.45\columnwidth, keepaspectratio]{zeta_vs_temp_NL3.pdf}
\caption{ The bulk viscosity as a function of temperature for
 several values of the baryon density marked on the plots for
 (a) model DDME2; (b) model NL3. In the temperature range
 $5\leq T\leq 10$\,MeV the bulk viscosity is obtained by
 interpolating between the results of neutrino-transparent
 and neutrino-trapped matter. The oscillation frequency is
 fixed at $f=1$ kHz, and the lepton fraction is fixed at $Y_L=0.1$. }
\label{fig:zeta_temp_final}
\end{center}
\end{figure}

The density dependence of the bulk viscosity is shown in
Fig.~\ref{fig:zeta_dens}. The oscillation frequency is fixed at
$f=1$ kHz in the case of neutrino-transparent matter whereas
the neutrino-trapped matter features a frequency-independent
bulk viscosity, as discussed above.
The bulk viscosity of neutrino-transparent matter mainly increases
with the density, the increase being faster at low temperatures
$T\lesssim 3$\,MeV where $\omega\gg \gamma$, as seen from
Fig.~\ref{fig:gamma_temp}.

The density dependence of the bulk viscosity of neutrino-trapped matter mainly follows that of the susceptibility $C^2/A$ because $\gamma$ depends weakly on the baryon density in this case. At sufficiently high temperatures $T\gtrsim 30$\,MeV there are sharp drops of the bulk viscosity to zero related to the fact that the matter becomes scale-invariant at certain critical values of density. Note that for temperature range $5\div 10$\,MeV we show the results for both neutrino-transparent and neutrino-trapped cases to account for the uncertainty in the exact value of neutrino trapping temperature $T_{\rm tr}$ which is supposed to lie in that range.

Figure~\ref{fig:zeta_temp} plots the dependence of the bulk viscosity
on the temperature. As already discussed above, the bulk viscosity of
neutrino-transparent matter attains its maximum at the temperature
where {$\gamma(T)=\omega$}, whereas the temperature dependence of the bulk viscosity
in the neutrino-trapped matter is mainly decreasing (up to point where
the matter becomes scale-invariant) because the relaxation rate is
already too fast: the resonant maximum would be at lower temperatures.
Since the relaxation rate $\gamma$ rises as $T^2$ [Equation~\eqref{eq:lambda_deg_trap}]
we expect [from  Equation~\eqref{eq:zeta_fast}] that
$\zeta\propto T^{-2}$ in this regime.

In Fig.~\ref{fig:zeta_temp_final}, panel (a) we combine the results
obtained for neutrino-transparent and neutrino-trapped matter by
interpolating the results between these two regimes in the temperature
range $5\leq T\leq 10$\,MeV which is regarded as the transition
region. Close to the transition temperature, the bulk viscosity is much
larger in the neutrino transparent regime because of much lower beta
relaxation rate and larger susceptibility $C^2/A$ as well. As a
result, the resonant peak of the bulk viscosity occurs always at or
below the neutrino-trapping temperature. Hence we can already
anticipate that the bulk viscosity is going to play an important role
in the dynamics of neutron star mergers in the regime of
neutrino-transparent rather than neutrino-trapped matter.

For comparison, we show also the bulk viscosity of nuclear matter
for an alternative  model NL3~\cite{Lalazissis1997} in panel (b) of
Fig.~\ref{fig:zeta_temp_final}. This model has density-independent
meson-nucleon couplings but contains non-linear terms in the $\sigma$-meson
field. The results obtained within the two models differ mainly
in the low-temperature regime, where the  model NL3 features much higher
viscosities than the model DDME2 above the saturation density. The reason for this is that NL3
model has a threshold of direct Urca opening at around $n_B\simeq 2.5 n_0$
whereas the model DDME2 does not have a threshold up to densities
$n_B=5n_0$. Because of the threshold, the model NL3 has much faster relaxation
rates at $n_B\gtrsim 2n_0$ than the model DDME2, and, as
a consequence, the maximum of bulk viscosity is shifted to lower
values of the temperature as compared to the DDME2 model.

\subsection{Estimating damping timescales}
\label{sec:damping}

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.45\columnwidth, keepaspectratio]{compress.pdf}
\caption{ The isothermal incompressibility of nuclear matter as a function of density for several temperatures for models DDME2 and NL3. }
\label{fig:compress}
\end{center}
\end{figure}

In this subsection we examine the characteristic timescales
of damping of density oscillations by the bulk viscosity in neutron
star mergers. The energy density stored in baryonic oscillations
with amplitude $\delta n_B$ is given by
\begin{eqnarray}\label{eq:osc_energy}
\epsilon =\frac{1}{2} \frac{\partial^2\varepsilon}{\partial n_B^2}
(\delta n_B)^2=\frac{K}{18}\frac{(\delta n_B)^2}{n_B},
\end{eqnarray}
where
\begin{eqnarray}\label{eq:compress}
K=9n_B\frac{\partial^2\varepsilon}{\partial n_B^2}
\end{eqnarray}
is the isothermal incompressibility of nuclear matter.
The energy dissipation rate by the bulk viscosity per unit volume is
\begin{eqnarray}\label{eq:energy_diss}
\frac{d\epsilon}{dt}=\frac{\omega^2\zeta}{2}
\left(\frac{\delta n_B}{n_B}\right)^2.
\end{eqnarray}
The characteristic timescale required for damping of oscillations
$\tau_\zeta =\epsilon/(d\epsilon/dt)$ is then given by
\begin{eqnarray}\label{eq:damping_time}
\tau_{\zeta} =\frac{1}{9}\frac{Kn_B}{\omega^2\zeta}.
\end{eqnarray}

The incompressibility of nuclear matter is shown in Fig.~\ref{fig:compress} for the two parametrizations discussed above. It is an increasing function of the density and at low densities decreases from its value at zero temperature as the temperature is increased.

Figures~\ref{fig:tau_damp1}a and ~\ref{fig:tau_damp1}b show, for two EoSs, the damping timescales of oscillations with frequency $f=1$\,kHz. We use the interpolation of bulk viscosity between neutrino-transparent and neutrino-trapped regimes shown in Fig.~\ref{fig:zeta_temp_final}.  As nuclear
incompressibility depends weakly on the temperature, the temperature
dependence of the damping timescale closely follows that of the
bulk viscosity. As a result, the damping timescale attains a minimum
at the temperature where the bulk viscosity has a maximum for a fixed
value of the density. The minimal value of the damping timescale is
\begin{eqnarray}\label{eq:damping_time_min}
\tau_{\zeta}^{\rm min} =\frac{2}{9\omega}\frac{Kn_B}{C^2/A}.
\end{eqnarray}
In the limits of slow and fast equilibration the damping
timescale is given by
\begin{eqnarray}\label{eq:damping_time_limits}
\tau^{\rm slow}_\zeta = \frac{1}{9\gamma}\frac{Kn_B}{C^2/A},\qquad
\tau^{\rm fast}_\zeta =\frac{\gamma }{9\omega^2}
\frac{Kn_B}{C^2/A}.
\end{eqnarray}
Thus, in contrast to the bulk viscosity, the damping timescale
becomes frequency-independent in the low-temperature (slow equilibration)
regime, and decreases with the frequency at high temperature
(fast equilibration) regime.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.45\columnwidth, keepaspectratio]{tau_damp1_DDME2.pdf}
\hskip 0.5cm
\includegraphics[width=0.45\columnwidth, keepaspectratio]{tau_damp1_NL3.pdf}
\caption{ The oscillation damping timescale as a function of temperature
for various densities and for frequency fixed at $f=1$ kHz. We interpolate between the neutrino-transparent and neutrino-trapped regime, as in Fig.~\ref{fig:zeta_temp_final}.
Panel (a): model DDME2; (b): model NL3.
}
\label{fig:tau_damp1}
\end{center}
\end{figure}

The density dependence of $\tau_{\zeta}$ reflects the density
dependence of the ratio of nuclear incompressibility and the
bulk viscosity. The density dependence of these two quantities almost compensates each other in the neutrino-transparent regime of the DDME2 model. In the case of NL3 model the density dependence of the bulk viscosity dominates and the damping
timescale mainly decreases with density. The reason for this is the increase of the reaction rates with density as a result of fast opening of phase space for
direct Urca reactions in the NL3 model.

\begin{figure}[!]
\begin{center}
\includegraphics[width=0.45\columnwidth, keepaspectratio]{tau_damp10_DDME2.pdf}
\hskip 0.5cm
\includegraphics[width=0.45\columnwidth, keepaspectratio]{tau_damp10_NL3.pdf}
\caption{The oscillation damping timescale as a function of temperature for various densities, interpolating between
the neutrino-transparent and neutrino-trapped
regimes as in Fig.~\ref{fig:tau_damp1} but for $f=10$ kHz;
(a) model DDME2; (b) model NL3.}
\label{fig:tau_damp10}
\end{center}
\end{figure}

The gray shaded areas in Fig.~\ref{fig:tau_damp1} show the temperature regions  where the damping timescale is smaller than the characteristic timescales for the early  ($\simeq 10$\,ms, dark shaded areas) and long-term ($\simeq 1$ s, light shaded areas) post-merger evolution, respectively. It is seen that in massive neutron star mergers, bulk viscosity strongly damps density oscillations at densities $n_B\gtrsim n_0$ in the temperature range 3 $\lesssim T\lesssim $ 6 MeV for model DDME2 and 1 $\lesssim T\lesssim$ 6 MeV for model NL3. On the timescale of long-term evolution
the damping is efficient also at lower densities, and the whole range of temperatures where the bulk viscosity  is relevant is 2 $\lesssim T\lesssim$ 7 MeV for
model DDME2 and 1 $\leq T\leq$ 7 MeV for model NL3.

\begin{figure}[t]
\begin{center}
\includegraphics[width=\columnwidth, keepaspectratio]{timescales1.pdf}
\caption{
Bulk viscous damping timescale as a function of
the density and temperature for model DDME2 (left panel) and model
NL3 (right panel). The density oscillation frequency is fixed at
$f=1$ kHz. The white solid lines correspond to the characteristic timescale
$\tau =10$\,ms, the blue dashed lines -
to the timescale $\tau =1$\,s. The areas colored in  red are the
regions where the bulk viscous dissipation becomes important in
the dynamics of neutron star mergers.
}
\label{fig:tau1_color}
\end{center}
\end{figure}
\begin{figure}[!]
\begin{center}
\includegraphics[width=\columnwidth, keepaspectratio]{timescales10.pdf}
\caption{ Same as Fig.~\ref{fig:tau1_color} but for frequency $f=10$ kHz.
}
\label{fig:tau10_color}
\end{center}
\end{figure}

Comparing the two panels of Fig.~\ref{fig:tau_damp1}  we observe that at low temperatures $T\leq 4$ MeV and at densities $n_B\geq 2 n_0$ the damping timescales are much shorter for the EoS model NL3 which has a direct Urca threshold. Thus, the dynamics and observational signatures of post-merger objects potentially  contain information on whether the direct Urca processes are operative in the high-density domain of the neutron stars.

The temperature where the damping timescale reaches its minimum for
$f=1$ kHz is located around $T\simeq 5$\,MeV at densities below the direct
Urca threshold and around $T\simeq 3$\,MeV above the threshold.
These results agree well with the results of Ref.~\cite{Alford2019a}
obtained within the Fermi-surface approximation. The exact computation,
however, obtained using exact beta-equilibrium condition for
the neutrino-transparent matter
at finite temperatures~\cite{Alford2018b} suggests that the
minimum of $\tau_\zeta$ is shifted to lower temperatures at
densities which are below the direct Urca threshold, and,
as a result, the minimum always appears around $T\simeq 3$\,MeV~\cite{Alford2019a}
(note that the authors of Ref.~\cite{Alford2019a} included also
the modified Urca processes in their calculations, which, however,
do not change the location of the maximum of bulk viscosity as
their contribution is subdominant above $T=2 $\,MeV).

The density-dependence of the damping timescale found here differs from that of Ref.~\cite{Alford2019a} where $\tau_\zeta$ was found to reach its minimum at low densities $n_B\lesssim n_0$. Apart from this, we find much lower values for the damping timescale at the minimum. This discrepancy arises because of the non-relativistic approximation for nucleon susceptibilities used in this work. This approximation works well at low densities, but strongly overestimates the susceptibility $C^2/A$ at higher densities $n_B\geq 2n_0$.

The triangles in Fig.~\ref{fig:C2A_dens} show the values of the susceptibility $C^2/A$ obtained in
Ref.~\cite{Alford2019a} for the models DD2 and IUFSU. We see that,
although the relativistic corrections to the spectrum of nucleonic excitations
are about $20 \%$ at $n_B=2n_0$, they need to be included in the susceptibilities.
This will require also a fully relativistic study of the
beta-equilibration rates which is relegated to a future work.

In Fig.~\ref{fig:tau_damp10}  we show the damping timescale for
10 kHz density oscillations. The minimum value of $\tau_\zeta$, in this case, is smaller than in the case of $f=1$ kHz by factors between 2 and 10, and the values of $\tau_\zeta$ in the neutrino-trapped regime are smaller by two orders of magnitude. However, the damping timescales of neutrino-trapped matter always remain larger than a second, since the bulk viscosity is not high enough to affect the evolution of mergers in this regime.

Figures~\ref{fig:tau1_color} and
\ref{fig:tau10_color} show the dependence of the bulk viscous damping timescale on the density and temperature colormap for the oscillation frequency fixed at $f=1$ kHz and $f=10$ kHz, respectively. The white solid and blue dashed lines show where the damping timescale becomes equal to the  timescales of 10 ms and 1 s,
respectively. In the areas shaded in dark red the bulk viscous damping timescale is $\tau_\zeta\leq 10$ ms, therefore, the  damping of density oscillations by the bulk viscosity is very efficient in those regimes. In the regions shaded in blue the
role of the bulk viscosity in damping of oscillations is negligible, as the damping timescale $\tau_\zeta\geq 1$ s there.

For completeness we comment also on how our results will change if larger lepton fractions are considered. The case $Y_L=0.4$ was studied in our previous work Ref.~\cite{Alford2019b}, where the bulk viscosity was shown to be reduced
by factors from 1 to 3 as compared to the $Y_L=0.1$ case. The pressure, and, therefore, also the nuclear incompressibility is only slightly sensitive to the lepton fraction. As a consequence, the damping timescales in
the $Y_L=0.4$ case will be larger than in the $Y_L=0.1$ case by factors of a few, but the overall quantitative picture will remain the same.

In closing, we stress again that at densities $n_B\geq 2n_0$ the relativistic corrections to the spectrum of nucleonic excitations
become important for the computation of the bulk viscosity, and our results at high densities need to be improved accordingly. Also the appearance of hyperons and  other heavy baryons needs to be taken into account. Finally, we note that in the case of hybrid stars with quark cores, the bulk viscosity of quark matter can be important for damping of density oscillations (for computations in the case of cold compact stars see \cite{
Madsen1992,Drago2005,Alford:2006gy,Alford:2007rw,Manuel:2007pz,SaD2007b,SaD2007a,Alford:2008pb,Huang2010,Wang2010PhRvD}).

\section{Summary}
\label{sec:summary}

We have reviewed the computation, ingredients, and approximations involved in computations of bulk viscosity of nuclear matter at finite temperatures relevant to binary neutron star mergers. The bulk viscosity arises from the direct Urca $\beta$-equilibration reactions. A novel ingredient relative to the studies of cold neutron stars is the trapped neutrino component coexisting with the nuclear matter at temperatures $T\gtrsim 5$\,MeV.  The concrete computations were carried out with the relativistic density functional approach to the EoS of nuclear matter with two different parametrizations.

At a given value of oscillation frequency $\omega\equiv 2\pi f$ the bulk viscosity
shows the standard resonant form \eqref{eq:zeta}, with a maximum where the beta relaxation rate $\gamma$ matches $\omega$.  This resonant maximum is achieved in the temperature range where neutrinos escape from the merger region, since the relaxation rate at temperatures of a few MeV is sufficiently low to match the density oscillation frequency. The reason for lower relaxation rates as compared to the neutrino-trapped case is the suppression of the direct Urca processes at the relevant temperatures and densities.

When the temperature rises to the threshold for neutrino trapping  ($T\sim 5\,{\rm MeV}$) the bulk viscosity experiences a sharp fall by several orders of magnitude as the relaxation rate rises and the material enters the fast beta-equilibration regime with $\gamma{\gg} \omega$. In this regime the bulk viscosity is independent of the frequency and decreases with the temperature approximately as $\zeta\propto T^{-2}$. At temperatures of about $30$\,MeV a new feature appears: the bulk viscosity drops to zero at the temperature where the beta-disequilibrium--baryon-density susceptibility $C$ vanishes, and then rises again at higher temperatures.  The susceptibility vanishes because
the particle fractions become independent of the density and the material becomes scale-invariant.

The main new result of this work concerns the timescales of damping of density oscillations in neutron star mergers by the bulk viscous dissipation. As an input we used the results for the bulk viscosity in Ref.~\cite{Alford2019b}.  We find that the damping timescale has a minimum as a function of temperature, which is located at temperatures in the range 3$\div$6 MeV for various densities. Assuming oscillation frequency of 1\,kHz we find that the damping timescale at its minimum is of the order of ms, {\it i.e.}, much shorter in the entire density range considered than the characteristic timescales of initial ($\sim 10$ms) and long-term ($\sim 1$s) post-merger evolution. We further find that the timescales of damping of density oscillations are shorter at the higher densities.  If the temperature is above the neutrino trapping temperature, the damping timescales are much longer
as the bulk viscosity is strongly suppressed. Finally, we note that bulk viscous dissipation could be of interest in the context of hydrodynamics simulations of supernovas, where electron capture rates on protons and nuclei could be out of equilibrium (for recent numerical simulations, see \cite{Mezzacappa2015,Fischer2017,Connor2018ApJ,Burrows2020MNRAS}).

\vspace{6pt}

\funding{The research of M. A. was funded by the U.S. Department of
  Energy, Office of Science, Office of Nuclear Physics under Award
  Number $\#$DE-FG02-05ER41375.  The research of A. H. and A. S. was
  funded by the Volkswagen Foundation (Hannover, Germany) grant
  No. 97029 and the European COST Action ``PHAROS'' (CA16214).  The
  research of A.S.\ was funded by Deutsche Forschungsgemeinschaft
  Grant No.  SE 1836/5-1. }

\acknowledgments{We thank Steven Harris, Kai Schwenzer, Alex Haber for
  discussions. M. A. and A. H. acknowledge the hospitality of Frankfurt Institute for Advanced Studies.}

\conflictsofinterest{The authors declare no conflict of interest.}

\reftitle{References}

\externalbibliography{yes}