$\displaystyle=p\left(\bm{\Phi},\bm{\sigma},\bm{z},\bm{\Psi},\bm{\Omega},\bm{\tau}^{2}_{\phi},\bm{\xi}_{\phi},\bar{\bm{\lambda}}^{2}_{\phi},\bm{\upsilon}_{\phi},\bm{\tau}^{2}_{\psi},\bm{\xi}_{\psi},\bar{\bm{\lambda}}^{2}_{\psi},\bm{\upsilon}_{\psi}\right)$
(42)
$\displaystyle=\prod_{q=1}^{nq}\prod_{j=1}^{m}\prod_{i=1}^{n}\Bigg{[}p(\bm{\phi}_{ij(q)}|\bar{\bm{\lambda}}_{ij(q),\phi}^{2},\tau_{ij(q),\phi}^{2},\sigma_{ij(q)})p(\sigma_{ij(q)})p(\tau_{ij(q),\phi}^{2}|\xi_{ij(q),\phi})p(\xi_{ij(q),\phi})\Bigg{.}$
$\displaystyle\>\>\>\>\>\>\cdot\Bigg{.}p(\bm{z}_{ij(q)})\prod_{k=1}^{l_{\phi}}p(\bar{\lambda}_{ij(q),k,\phi}^{2}|\upsilon_{ij(q),k,\phi})p(\upsilon_{ij(q),k,\phi})\Bigg{]}$
$\displaystyle\>\>\>\>\>\>\cdot\prod_{r=1}^{R}\Bigg{[}\Bigg{.}p(\bm{\psi}_{r}|\bar{\bm{\lambda}}_{r,\psi}^{2},\tau_{r,\psi}^{2},\omega_{r})p(\omega_{r})p(\bm{a}_{r})p(\tau_{r,\psi}^{2}|\xi_{r,\psi})p(\xi_{r,\psi})\prod_{s=1}^{l_{\psi}}p(\bar{\lambda}_{r,s,\psi}^{2}|\upsilon_{r,s,\psi})p(\upsilon_{r,s,\psi})\Bigg{]}.$

Given this partitioning, it can be shown that the solution to the optimization
problem can be obtained by sequentially iterating over the densities

$q(\bm{\theta}_{l}|\bm{x})\propto\text{exp}\>\>\mathbb{E}_{q(\bm{\theta}_{(-l)}|\bm{x})}\left(\text{log}\>\>p(\bm{\theta}_{l}|\bm{\theta}_{(-l)},\bm{x})\right)$
(43)

where $\bm{\theta}_{(-l)}$ denotes all elements of $\bm{\theta}$, excluding
those in the $l^{th}$ group, $l=1,\cdots,L$. Consequently, the variational
posterior can be obtained by calculating the variational expectation of the
conditional posterior. Generally, the accuracy of the variational
approximation depends on how well the chosen partitioning matches the
independence structure of the parameters in the target posterior. In essence,
this is where we trade off computational tractability and speed with accuracy.
For a general discussion of this issue see e.g. Ormerod et al. (2017). In
order to arrive at the final variational densities, we now need to insert 41
and 42 into 43. Given that this expression relies on expectations of
logarithms of standard densities it is rather easy to derive, however, we
suppress it here due to its very lengthy nature. The expressions for
$q(\bm{\theta}_{l}|\bm{x})$ are presented in algorithm 2 together with the
full variational algorithm.

Define $\kappa_{1(q)}=\frac{1-2q}{q(1-q)}$ and
$\kappa_{2(q)}^{2}=\frac{2}{q(1-q)}$ for each quantile level $q$.

Step 1: Estimate the factors $\hat{F}$ using the variational VBQFA algorithm
proposed in Korobilis and Schröder (2022) for each quantile level, $q$, and
for the variable blockings in 4.

Step 2: Conditional on the factors, estimate the remaining parameters of the
model.

begin

$[\bm{1}]$ Update the parameters in the measurement equation:

for  _$q=1:nq$_ do

for  _$i=1:m$_ do

for  _$j=1:n$_ do

$[\text{i}]$ Update $\bm{\phi}_{ij(q)}$ from
$q(\bm{\phi}_{ij(q)}|\bm{x})=N(\bm{\mu}_{ij(q)}^{\phi},\bm{\Sigma}_{ij(q)}^{\phi})$,
with
$\bm{\mu}^{\phi}_{ij(q)}=\bm{\Sigma}_{ij(q)}^{\phi}\left\\{\hat{\bm{F}}^{\prime}diag(\tilde{\bm{\nu}}_{ij(q)}^{-1})\tilde{\bm{y}}_{ij(q)}\right\\}$,
$\bm{\Sigma}_{ij(q)}^{\phi}=\left\\{\hat{\bm{F}}^{\prime}diag(\tilde{\bm{\nu}}_{ij(q)})^{-1}\hat{\bm{F}}+\tilde{\bm{\Lambda}}_{ij(q)}^{-1}\right\\}^{-1}$,
where
$\tilde{\bm{\nu}}_{ij(q)}^{-1}=\kappa_{2(q)}^{-2}\mathbb{E}\left(\sigma_{ij(q)}^{-1}\right)\mathbb{E}\left(\bm{z}_{ij,\bullet(q)}^{-1}\right)$,
$\tilde{\bm{y}}_{ij(q)}=\bm{y}_{ij,\bullet(q)}-\kappa_{1(q)}\mathbb{E}\left(\bm{z}_{ij,\bullet(q)}\right)$,
and
$\tilde{\bm{\Lambda}}_{ij(q)}^{-1}=diag\left(\mathbb{E}\left(\bar{\bm{\lambda}}_{ij(q),\bullet,\phi}^{-2}\right)\mathbb{E}\left(\tau_{ij(q),\phi}^{-2}\right)\right)$.

Assign: $\mathbb{E}\left(\bm{\phi}_{ij(q)}\right)=\bm{\mu}_{ij(q)}^{\phi}$ and
$\mathbb{E}\left(\bm{\phi}_{ij(q)}^{2}\right)=\bm{\mu}_{ij(q)}^{2,\phi}+diag\left(\Sigma_{ij(q)}^{\phi}\right)$

$[\text{ii}]$ Update the parameters of the horseshoe prior from

$q(\bm{\bar{}}{\lambda}_{ij(q),k,\phi}^{2}|\bm{x})=IG(a_{\bar{\lambda},\phi},b_{\bar{\lambda},\phi})=IG\left(1,\frac{\mathbb{E}(\bm{\phi}_{ij,k(q)}^{2})}{2}+\mathbb{E}\left(\frac{1}{\upsilon_{ij(q),k,\phi}}\right)\right)$,

$q({\upsilon}_{ij(q),k,\phi}|\bm{x})=IG(a_{\upsilon,\phi},b_{\upsilon,\phi})=IG\left(1,\mathbb{E}\left(\frac{1}{\bar{\bm{\lambda}}_{ij(q),k,\phi}^{2}}\right)+b_{\phi}^{-2}\mathbb{E}\left(\frac{1}{\tau^{2}_{ij,\phi}}\right)\right)$,

$q(\tau_{ij(q),k,\phi}^{2}|\bm{x})=IG(a_{\tau,\phi},b_{\tau,\phi})=IG\left(1,b_{\phi}^{-2}\mathbb{E}\left(\frac{1}{\upsilon_{ij(q),k,\phi}}\right)+\mathbb{E}\left(\frac{1}{\xi_{ij(q),\phi}}\right)\right)$,

$q({\xi}_{ij(q),\phi}|\bm{x})=IG(a_{\xi,\phi},b_{\xi,\phi})=IG\left(\frac{l_{\phi}+1}{2},1+\sum_{k=1}^{l_{\phi}}\mathbb{E}\left(\frac{1}{\tau_{ij(q),k,\phi}^{2}}\right)\right)$,

for $k=1,\cdots,l_{\phi}$.

Assign:
$\mathbb{E}\left(\frac{a_{\bar{\lambda},\phi}}{b_{\bar{\lambda},\phi}}\right)$,
$\mathbb{E}\left(\frac{a_{\upsilon,\phi}}{b_{\upsilon,\phi}}\right)$,
$\mathbb{E}\left(\frac{a_{\tau,\phi}}{b_{\tau,\phi}}\right)$, and
$\mathbb{E}\left(\frac{a_{\xi,\phi}}{b_{\xi,\phi}}\right)$.

$[\text{iii}]$ Update the latent quantile indicators, $\bm{z}_{ij,t(q)}$, from
$q(\bm{z}_{ij,t(q)}|\bm{x})=GIG(1/2,\delta_{z},\rho_{t,z})$ with

$\delta_{z}$ =
$\mathbb{E}\left(\frac{1}{\sigma_{ij(q)}}\right)\frac{\kappa_{1(q)}^{2}}{\kappa_{2(q)}^{2}}$,
$\rho_{t,z}=\mathbb{E}\left(\frac{1}{\sigma_{ij(q)}}\right)\frac{\left\\{y_{ij,t(q)}-\mathbb{E}(\bm{\phi}_{ij(q)})^{\prime}\bm{\hat{}}{\bm{f}_{t}}\right\\}^{2}+\hat{\bm{f}_{t}}^{\prime}\bm{\Sigma}_{ij(q)}^{\phi}\hat{\bm{f}_{t}}}{\kappa_{2(q)}^{2}}$,

for $t=1,\cdots,T$

Assign:
$\mathbb{E}(z_{ij,t(q)})=\frac{\sqrt{\rho_{t,z}}K_{3/2}\left(\sqrt{\delta_{z}\rho_{t,z}}\right)}{\sqrt{\delta_{z}}K_{1/2}\left(\sqrt{\delta_{z}\rho_{t,z}}\right)}$,
$\mathbb{E}\left(\frac{1}{z_{ij,t(q)}}\right)=\frac{\sqrt{\delta_{z}}K_{3/2}\left(\sqrt{\delta_{z}\rho_{t,z}}\right)}{\sqrt{\rho_{t,z}}K_{1/2}\left(\sqrt{\delta_{z}\rho_{t,z}}\right)}-\frac{1}{\rho_{t,z}}$

$[\text{iv}]$ Update the factor regression variances, $\sigma_{ij(q)}$, from
$q(\sigma_{ij(q)}|\bm{x})=IG(r_{\sigma},s_{\sigma})$ with

$r_{\sigma}=r_{0}+3T$ and
$s_{\sigma}=s_{0}+\sum_{t=1}^{T}\left[\mathbb{E}\left(\frac{1}{z_{ij,t(q)}}\right)\frac{M}{2\kappa_{2(q)}^{2}}-\kappa_{1(q)}^{2}\frac{y_{ij,t(q)}-\mathbb{E}(\bm{\phi}_{ij(q)})^{\prime}\bm{\hat{}}{\bm{f}_{t}}}{\kappa_{2(q)}^{2}}+\left(1+\frac{\kappa_{1(q)}^{2}}{2\kappa_{2(q)}^{2}}\right)\mathbb{E}(z_{ij,t(q)})\right]$
and
$M=\left\\{y_{ij,t(q)}-\mathbb{E}(\bm{\phi}_{ij(q)})^{\prime}\bm{\hat{}}{\bm{f}_{t}}\right\\}^{2}+\hat{\bm{f}_{t}}^{\prime}\bm{\Sigma}_{ij(q)}^{\phi}\hat{\bm{f}_{t}}$.

Assign:
$\mathbb{E}\left(\frac{1}{\sigma_{ij(q)}}\right)=\frac{r_{\sigma}}{s_{\sigma}}$.

end for

end for

end for

$[\bm{2}]$ Update the diagonal elements of $\Omega$ from
$q(\omega_{r}|\bm{x})=IG(r_{\omega},s_{\omega})$ with

$r_{\omega}=r_{h}+T/2$ and
$s_{\omega}=s_{h}+\sum_{t=1}^{T}(\hat{\bm{f}}_{r,t}-\hat{\bm{F}}_{t-1}\bm{\psi}_{r}^{\prime})^{2})$,
for all $r=1,\cdots,R$.

Assign:
$\mathbb{E}\left(\frac{1}{\omega_{r}}\right)=\frac{r_{\omega}}{s_{\omega}}$.

$[\text{3}]$ Sample the off-diagonal elements of $\bm{\Omega}$, $\bm{A}$, and
the VAR coefficients $\bm{\Psi}$

for  _$r=1:R$_ do

$[\text{i}]$ Sample $\bm{\psi}_{r}$ and $\bm{a}_{r}$ jointly from
$q(\bm{\psi}_{r},\bm{a}_{r}|\bm{x})=N(\bm{\mu}^{\psi}_{r},\bm{\Sigma}_{r}^{\psi})$,
with
$\bm{\mu}^{\psi}_{r}=\bm{\Sigma}_{r}^{\psi}\left\\{\bm{X}^{\prime}\mathbb{E}(\sigma_{\omega}^{-1})\hat{\bm{f}}_{r}\right\\}$,
$\bm{\Sigma}_{r}^{\psi}=\left\\{\bm{X}^{\prime}\mathbb{E}(\sigma_{\omega}^{-1})\bm{X}+\tilde{\bm{\Lambda}}_{r}^{-1}\right\\}^{-1}$,
where
$\bm{X}=[L\bm{F}^{\prime},L\bm{G}^{\prime},\bm{E}_{1:r-1}^{\prime}]^{\prime}$,
$\bm{e}_{r}=\hat{\bm{f}}_{r}-\bm{\hat{}}{\bm{F}}\mathbb{E}(\bm{\psi}_{r})^{\prime}-\bm{E}_{1:r-1}^{\prime}\mathbb{E}(\bm{a}_{r})^{\prime}$,
and
$\tilde{\bm{\Lambda}}_{r}^{-1}=diag\left([\mathbb{E}(\bar{\bm{\lambda}}_{r,\bullet,\psi}^{-2}),\bm{\mu}_{r,a}/\bm{\Sigma}_{r,a}]\right)$.

$[\text{ii}]$ Update the parameters of the horseshoe prior from
$q(\bar{\lambda}_{r,s,\psi}^{2}|\bm{x})=IG(a_{\bar{\lambda},\psi},b_{\bar{\lambda},\psi})=IG\left(1,\frac{\mathbb{E}(\bm{\psi}_{r,s}^{2})}{2}+\frac{1}{\mathbb{E}(\upsilon_{r,s,\psi})}\right)$,

$q({\upsilon}_{r,s,\psi}|\bm{x})=IG(a_{{\upsilon},\psi},b_{{\upsilon},\psi})=IG\left(1,\mathbb{E}\left(\frac{1}{\bar{\lambda}_{r,s,\psi}^{2}}\right)+b_{\psi}^{-2}\mathbb{E}\left(\frac{1}{\tau_{r,\psi}^{2}}\right)\right)$,

$q(\tau_{r,s,\psi}^{2}|\bm{x})=IG(a_{{\tau},\psi},b_{{\tau},\psi})=IG\left(1,b_{\psi}^{-2}\mathbb{E}\left(\frac{1}{{\upsilon}_{r,s,\psi}}\right)+\mathbb{E}\left(\frac{1}{{\xi}_{r,\psi}}\right)\right)$,

$q({\xi}_{r,\psi}|\bm{x})=IG(a_{{\xi},\psi},b_{{\xi},\psi})=IG\left(\frac{l_{\psi}+1}{2},1+\sum_{s=1}^{l_{\psi}}\mathbb{E}\left(\frac{1}{\tau_{r,s,\psi}^{2}}\right)\right)$,

for $s=1,\cdots,l_{\psi}$.

Assign:
$\mathbb{E}\left(\frac{a_{\bar{\lambda},\psi}}{b_{\bar{\lambda},\psi}}\right)$,
$\mathbb{E}\left(\frac{a_{\upsilon,\psi}}{b_{\upsilon,\psi}}\right)$,
$\mathbb{E}\left(\frac{a_{\tau,\psi}}{b_{\tau,\psi}}\right)$, and
$\mathbb{E}\left(\frac{a_{\xi,\psi}}{b_{\xi,\psi}}\right)$.

end for

end

Note: We use the rate parametrization of the $IG$ distribution throughout.
$K_{p}(\bullet)$ denotes the Bessel function of order $p$ and $L$ denotes the
Lag operator.

Algorithm 2 2-step Variational Bayes QFAVAR algorithm

Finally, note that the horseshoe priors take a slightly different form in our
VB algorithm with

$\displaystyle\bm{\phi}_{ij(q)}|\left\\{\bar{\lambda}_{ij(q),k,\phi},\upsilon_{ij(q),k,\phi}\right\\}_{k=1}^{l_{\phi}},\tau_{ij(q),\phi},\xi_{ij(q),\phi}$
$\displaystyle\sim N(\bm{0},\bar{\bm{\Lambda}}_{ij(q),\phi}),$ (44)
$\displaystyle\bar{\lambda}^{2}_{ij(q),k,\phi}|\upsilon_{ij(q),k,\phi}$
$\displaystyle\sim
G^{-1}\left(\frac{1}{2},\frac{1}{\upsilon_{ij(q),k,\phi}}\right),\>\>\>\text{for
}k=1,\cdots,l_{\phi},$
$\displaystyle\upsilon_{ij(q),k,\phi}|\tau^{2}_{ij(q),\phi}$
$\displaystyle\sim
G^{-1}\left(\frac{1}{2},\frac{1}{b_{\phi}^{2}\tau^{2}_{ij(q),\phi}}\right),\>\>\>\text{for
}k=1,\cdots,l_{\phi},$ $\displaystyle\tau^{2}_{ij(q),k,\phi}|\xi_{ij(q),\phi}$
$\displaystyle\sim G^{-1}\left(\frac{1}{2},\frac{1}{\xi_{ij(q),\phi}}\right),$
$\displaystyle\xi_{ij(q),\phi}$ $\displaystyle\sim
G^{-1}\left(\frac{1}{2},1\right),$

and

$\displaystyle\bm{\psi}_{r}|\left\\{\bar{\lambda}_{r,s,\psi},\upsilon_{r,s,\psi}\right\\}_{s=1}^{l_{\psi}},\tau_{r,\psi},\xi_{r,\psi}$
$\displaystyle\sim N(\bm{0},\bar{\bm{\Lambda}}_{r,\psi}),$ (45)
$\displaystyle\bar{\lambda}^{2}_{r,s,\psi}|\upsilon_{r,s,\psi}$
$\displaystyle\sim
G^{-1}\left(\frac{1}{2},\frac{1}{\upsilon_{r,s,\psi}}\right),\>\>\>\text{for
}s=1,\cdots,l_{\psi},$ $\displaystyle\upsilon_{r,s,\psi}|\tau^{2}_{r,\psi}$
$\displaystyle\sim
G^{-1}\left(\frac{1}{2},\frac{1}{b_{\psi}^{2}\tau^{2}_{r,\psi}}\right),\>\>\>\text{for
}s=1,\cdots,l_{\psi},$ $\displaystyle\tau^{2}_{r,s,\psi}|\xi_{r,\psi}$
$\displaystyle\sim G^{-1}\left(\frac{1}{2},\frac{1}{\xi_{r,\psi}}\right),$
$\displaystyle\xi_{r,\psi}$ $\displaystyle\sim
G^{-1}\left(\frac{1}{2},1\right),$

where $b_{\phi}$ and $b_{\psi}$ are hyperparameters, which we set to $0.0001$
in application.

##  Appendix B Additional empirical results

###  B.1 Forecast evaluation

Figure 1: Cumulative quantile score (QS) loss for forecast horizon $h=6$. Four
models are compared, QAR (yellow circled line), QAR-X (purple dashed line),
QDFM (red dotted line) and QFAVAR (blue solid line). The out-of-sample
evaluation period shown on the x-axis is 2011Jan to 2022Dec-$h$. First two
rows show quantile scores for the 10th and 90th percentiles of inflation and
third and fourth rows show quantile scores for the 10th and 90th percentiles
of industrial production. Figure 2: Cumulative quantile score (QS) loss for
forecast horizon $h=12$. Four models are compared, QAR (yellow circled line),
QAR-X (purple dashed line), QDFM (red dotted line) and QFAVAR (blue solid
line). The out-of-sample evaluation period shown on the x-axis is 2011Jan to
2022Dec-$h$. First two rows show quantile scores for the 10th and 90th
percentiles of inflation and third and fourth rows show quantile scores for
the 10th and 90th percentiles of industrial production. Figure 3: Cumulative
quantile score (QS) loss for forecast horizon $h=24$. Four models are
compared, QAR (yellow circled line), QAR-X (purple dashed line), QDFM (red
dotted line) and QFAVAR (blue solid line). The out-of-sample evaluation period
shown on the x-axis is 2011Jan to 2022Dec-$h$. First two rows show quantile
scores for the 10th and 90th percentiles of inflation and third and fourth
rows show quantile scores for the 10th and 90th percentiles of industrial
production.

###  B.2 Detailed quantile IRFs

Figure 4: Detailed quantile impulse response functions, shock to global
inflation (GINF). Figure 5: Detailed quantile impulse response functions,
shock to global supply chain pressure index (GSCPI). Figure 6: Detailed
quantile impulse response functions, shock to global (US) financial conditions
(FCI). Figure 7: Detailed quantile impulse response functions, shock to global
economic policy uncertainty (GEPU).


# Gamilaraay kinship revisited: incidence of recessive disease is dynamically
traded-off against benefits of cooperative behaviours

Jared M. Field<EMAIL_ADDRESS>School of Mathematics and
Statistics, University of Melbourne, Melbourne, Australia

###### Abstract

Traditional Indigenous marriage rules have been studied extensively since the
mid 1800s. Despite this, they have historically been cast aside as having very
little utility. This is, in large part, due to a focus on trying to understand
broad-stroke marriage restrictions or how they may evolve. Here, taking the
Gamilaraay system as a case study, we instead ask how relatedness may be
distributed under such a system. We show, remarkably, that this system
dynamically trades off kin avoidance to minimise incidence of recessive
diseases against expected levels of cooperation, as understood formally
through Hamilton’s rule.

## Introduction

In his hugely influential work Ancient Society Lewis Morgan wrote that the
Gamilaraay kinship system represents ‘the most primitive form of society
hitherto discovered’morgan1877ancient . Such sentiment persisted for decades,
with many academics reluctant to throw out old ideas of Social Evolutionism,
whereby all societies were perceived to be on a sliding scale of savagery to
civilisation (see Chapter 2 of mcconvell2018skin  for a survey of such ideas
in the Australian context). More rigorous (and slightly less colonial) work
arrived with the seminal study of Levi-Strauss levi2020structures . Here, for
example, the mathematician André Weil used group theory to understand
precisely who can marry whom in kinship systems, such as that of the
Gamilaraay, where descent operates in cycles. Since then most work has
followed Weil to attempt to understand broad-stroke marriage restrictions or
indeed how they might come about, with studies in this area more generally
remaining hotly debated (see for example itao2020evolution ;
read2020simulation ; itao2020reply ). Here, we take a different tact entirely
and instead ask how relatedness may be distributed in such societies. More
than this, we ask how might this distribution change over time. In this way,
we aim to better grasp not just the form of the Gamilaraay kinship system but
also its function.

To this end, we worked out the exact distribution of common ancestors that may
appear in a general Gamilaraay pedigree. Coupled with kinship coefficients,
this allowed for the derivation of an average kinship coefficient for a couple
permitted to reproduce in this system. Under the assumption that inbreeding
coefficients are negligible at the great grandparent stage, we found that the
entire Gamilaraay nation would need to reduce to just 24 individuals so that
couples may (on average) be as closely related as full first cousins. While
striking, inbreeding coefficients may of course not always be negligible. For
this reason we dropped this assumption, and instead asked how average kinship
coefficients change over time. We found that there is, in fact, an equilibrium
point towards which average kinship is pushed. This point is a function of
reproductive population size, and so dynamically re-calibrates as that
quantity changes. In this way, average kinship may both decrease and increase
in reaction to historical values. Remarkably, the Gamilaraay kinship system
dynamically trades-off the incidence of recessive diseases against the
benefits of cooperative behaviours, formally understood through Hamilton’s
Rule frank2013naturalhist . Otherwise put, the Gamilaraay kinship system does
not simply strive to create kin restrictions, but instead has an underlying
utility function that also accounts for social benefits. Contrary to the
beliefs of Morgan, Gamilaraay culture appears to be anything but primitive.

The rest of this paper is organised as follows: we start by providing
background on the Gamilaraay system. Following this, we map a general pedigree
by sub-types, or marriageable groups, to show where uncertainty may first
enter a family tree. Next, we use combinatoric arguments to derive explicit
distributions of common ancestors first on a particular branch, and then
across all branches at the great grandparent stage. In turn, we use this to
derive average kinship coefficients for a general permitted couple and
consider their dynamics. Finally, we contextualise our findings in the current
theoretical biology literature, discuss their consequences and consider future
directions of work.

## Gamilaraay System

In the Gamilaraay kinship system, as with many others of this region (see for
example mcconvell2018skin  for a more linguistic treatment), all individuals
are divided into two groups, namely Guwaymadhaan (Dark or Heavy Blood) and
Guwaygaliyarr (Light Blood). Each child inherits belonging to a particular
group directly from their mother, so that, for example, a Guwaymadhaan mother
produces a Guwaymadhaan child. These two groups are split further in two so
that, for the purpose of marriage, there are four subgroups in total. These
subgroupings are such that a child belongs to the subgroup opposite to that of
their mother. In the Gamilaraay system, for example, a Yibadhaa woman produces
Gambuu sons and Buudhaa daughters (while the language is different, due to
gender, the types are the same.). If a Buudhaa woman reproduces, her children
then once again become Yibaay/Yibadhaa, which will not however be the case for
her Gambuu brother. Marriage rules are then such that one can only marry an
individual that is of opposite blood type, and opposite subgrouping to their
father. This is presented more clearly in Fig. 1. Here general labels are
adopted both to ease mathematical analysis later and because this system
extends beyond the Gamilaraay nation. See fison1880kamilaroi ;
parker1905euahlayi  for older anthropological treatments of the Gamilaraay
system and Spearim2014a ; Spearim2014b  for a modern primary source account. A
few obvious kin marriage restrictions, from Fig. 1, become clear. For example,
an individual may not marry siblings, any maternal cousins related through
females, nor any aunties or uncles. Paternal cousins are possible (likely
half, due to marriage norms), but as we will show, the probability is low.

Figure 1: Kinship dynamics of the four group Gamilaraay system. Solid lines
represent matrilineal inheritance of subgrouping whereas dashed lines link
subgroups that are permitted to marry. Individuals are only permitted to marry
someone of the opposite bloody type (CD type people, if an individual is of
type AB, for example) and opposite subgrouping to their father (C type people,
if an individual is B type). Note that this ensures that an individual can
only marry someone who is of type different to themselves, and both their
parents. By virtue of inheritance rules it also ensures individuals do not
marry siblings, aunties, uncles and maternal cousins of any order related
through females.

## General Pedigree

To work out the average relatedness of a couple, first note that by virtue of
this system, uncertainty enters pedigrees only at the great-grandparent stage.
This is most clear in Fig. 2. Here, we take a general individual of
subgrouping type B, and trace their pedigree. This individual will have a
mother of type A and father of type D (denoting movement to the left and right
as female and male ancestors, respectively). Tracing back one generation
further, a type A individual will have a type B mother and C father whereas a
type D individual will have a type C mother and B father. While there are two
Bs and two Cs at this stage, they must be distinct individuals owing to their
different sexes. Repeating this process, we see that at the great-grandparent
stage there will be four types As and four type Bs. However, they will be
split across the sexes. In other words, there will be precisely two male As,
two females As, two male Ds, and two female Ds. At this stage, then, it is now
possible that one individual may occupy two distinct spots in this
genealogical tree. In the next section, we work out the precise probability of
one of those spots being filled by the same or distinct individuals.

Figure 2: General pedigree of an individual in the Gamilaraay system following
the subgrouping notation as in Fig. 1. Sex in this figure is encoded by left
(female) and right (male) branches. While this pedigree is presented for an
individual of type B, due to symmetry, this general form applies to all
individuals. From Fig. 1. we see that a type B individual will have a type A
mother and type D father. In turn, a type A individual will have a type B
mother and type C father. A type D individual, however, will have a type C
mother and type B father. Note that uncertainty only enters at the great-
grandparent stage, where one individual may occupy two distinct spots in this
potential geneological tree. For example, one type A female (left branches)
may potentially fill both type A female positions (first and seventh). This
will make her the same mother to the type B female and male in the generation
below, and hence increase the relatedness of the AD couple one generation
later.

## Distribution of common ancestors on one branch

Suppose there are $n$ reproductive individuals in the population belonging to
our subgroup and sex of interest. Further, suppose we assign each of these
individuals a unique label between $1,...,n$. With this notation, we may
enumerate every possible combination of individuals (to fill the two permitted
genealogical spots) by multisets of the form $\\{1,1\\}$, $\\{1,2\\}$,
$\\{1,3\\}$, … , $\\{n-1,n\\}$, $\\{n,n\\}$. $\\{1,1\\}$, for example,
represents the instance where the first individual occupies both spots.
$\\{6,1\\}$ represents the instance where the sixth and first individual
occupy a spot each, and so forth. The total number of different ways this
section of the pedigree may be realised will then be given by the total number
of multisets of the above form. This, in turn, is well known to be given by
the binomial coefficient:

${n+1\choose 2}.$ (1)

The number of instances where an individual occupies both spots (which is to
say, of the form $\\{3,3\\}$, for example) is simply the total number of
individuals $n$. Assuming each instance is equally likely, this in turn means
that the probability that both spots are occupied by the same individual in
general will be:

$\displaystyle p_{common}$ $\displaystyle=\frac{n}{{n+1\choose 2}},$ (2)
$\displaystyle=\frac{2}{n+1}.$ (3)

If the two spots are occupied by different individuals, the first may be
filled by $n$ individuals, while the second may be filled by the remaining
$(n-1)$. As we do not care about order, we need to factor out repeats by $2!$
(the ordered couple $(1,2)$ for our purposes is equivalent to $(2,1)$). Hence,
the number of instances where individuals are distinct will be given by:

$\frac{n\cdot(n-1)}{2!},$ (4)

so that the probability the two spots are filled by distinct individuals is

$\displaystyle p_{distinct}$
$\displaystyle=\frac{\frac{n\cdot(n-1)}{2!}}{{n+1\choose 2}},$ (5)
$\displaystyle=\frac{n-1}{n+1}.$ (6)

## Distribution of common ancestors across all branches

Recall that at the great-grandparent stage there are eight individuals in
total, spread across the two sexes and four subgroups. In the previous
section, we worked out the probability that two permitted spots of a
particular subtype are filled by the same or different individuals. We now
need to combine these probabilities to work out the distribution of any number
of common ancestors at this generation. In other words, the probability that
there are $4,6,7$ or $8$ distinct individuals in this generation in total.

To do so, note that there are four subtypes in this generation: male and
female type As, and male and female type Ds, for this pedigree. If we encode
the event that one individual occupies both genealogical spots of a subtype by
$1$ and the other case by $2$ then we can use 4-tuples to count ancestors. For
example, the 4-tuple $(1,2,2,2)$ represents the instance where the first
subtype has one individual in both spots, while the other subtypes have
distinct individuals. In other words, of the total possible eight ancestors at
this stage there are only seven. In turn this would mean that our couple of
interest in Fig. 2, individuals A and D in generation one, will be related.
With this notation, it is clear that the total number of combinations where
there are seven ancestors will be given by the binomial coefficient $4\choose
1$, so that the probability of having seven ancestors at the great grandparent
stage is

$p_{7}={4\choose 1}\cdot(p_{common})\cdot(p_{distinct})^{3}.$ (7)

With similar reasoning the other probabilities can be shown to be given by

$\displaystyle p_{4}$ $\displaystyle={4\choose 4}\cdot(p_{common})^{4},$ (8)
$\displaystyle p_{5}$ $\displaystyle={4\choose
3}\cdot(p_{common})^{3}\cdot(p_{distinct}),$ (9) $\displaystyle p_{6}$
$\displaystyle={4\choose 2}\cdot(p_{common})^{2}\cdot(p_{distinct})^{2},$ (10)
$\displaystyle p_{8}$ $\displaystyle={4\choose 0}\cdot(p_{distinct})^{4}.$
(11)

## Average kinship coefficient

If there were no uncertainty in our problem, and we knew exactly how many
common ancestors were at a given stage, calculating coefficients of kinship, a
measure of relatedness useful in the study of disease, would be simple. For
example, suppose there were one common ancestor at the great grandparent
stage, a type D female. This would mean that one copy of the tree as per Fig.
3. would appear in this particular pedigree. To work out the coefficient of
kinship of the AD couple, we would simply need to count the number of
individuals in the path between them, and raise $1/2$ to that number. In this
case, assuming that their common ancestor D is not the result of related
individuals, their coefficient would be $1/2^{5}$. If instead their common
ancestor D is in turn the offspring of related individuals, their coefficient
of kinship would be accordingly higher. In particular, we can write their
coefficient of kinship as $1/2^{5}\cdot(1+f_{D})$, where $f_{D}$ is the
coefficient of inbreeding of D wright1922coefficients ; crow1986basic .

More generally, if there are $k$ common ancestors at this stage then the
coefficient of kinship of the $AD$ couple can be written as

$f_{AD}=\sum^{k}_{i=1}\frac{1}{2^{5}}(1+f_{i}),$ (12)

where $f_{i}$ is the coefficient of inbreeding of each of the $k$ common
ancestors. Averaging over the distribution that the $f_{i}$ are drawn from we
can write

$\displaystyle\bar{f}_{AD}$
$\displaystyle=\sum^{k}_{i=1}\frac{1}{2^{5}}(1+\bar{f}),$ (13)
$\displaystyle=\frac{1}{2^{5}}(k+k\bar{f}).$ (14)

There is also however uncertainty in the total number of common ancestors $k$.
Averaging over $k$, using the probabilities from the previous section, we
obtain the average coefficient of kinship:

$\displaystyle\bar{\bar{f}}_{AD}$
$\displaystyle=\sum^{4}_{k=0}\frac{1}{2^{5}}(k+k\bar{f})p_{8-k},$ (15)
$\displaystyle=\frac{1}{2^{5}}(\bar{k}+\bar{k}\bar{f}),$ (16)

where bars indicate averages.

In Fig. 4, we plot $\bar{\bar{f}}_{AD}$ as a function of $n$, one eighth of
the total number of reproductive individuals in the population. Here we assume
$\bar{f}=0$, an assumption we drop later. References lines are added to show
relatedness of full first cousins (top grey) and full second cousins (bottom
grey). Observe that the total population would need to reduce to $24$
individuals $(n=3)$ for average relatedness of couples to be the same as full
first cousins.

Figure 3: Section of a pedigree assuming one common ancestor at the great
grandparent stage. To calculate the kinship coefficient of the couple AD, we
need to count the number of individuals in the paths between them (5) and
raise $1/2$ to that number. Here, $f=1/2^{5}$.  Figure 4: Average coefficient
of kinship $\bar{\bar{f}}$ of a randomly sampled couple operating under the
Gamilaraay system as a function of $n$, one eighth of the total reproductive
population size. Here we assume $\bar{f}$, inbreeding at the great grandparent
stage, is zero (an assumption we drop later). Grey reference lines indicate
$\bar{\bar{f}}$ for full first (top) and full second (bottom) cousins. Note
that, for a couple to be related as closely (on average) as full first
cousins, the total population size needs to reduce to $24$ individuals
($n=3$).

## Dynamics

While in the previous section we calculated the average kinship coefficient
for an $AD$ couple, our equation holds for any permitted mating. To derive the
equivalent equation for a $BC$ couple, we would simply need to reseed the
original pedigree with either an $A$ or $D$ child. In light of this, we can
drop the $AD$ from equation (16) to write

$\bar{\bar{f}}=\frac{1}{2^{5}}(\bar{k}+\bar{k}\bar{f}),$ (17)

which is the average kinship coefficient of any couple operating under the
Gamilaraay system. To understand how kinship coefficients may change over
time, note that the coefficient of inbreeding of an individual $f_{i}$, is in
fact the coefficient of kinship of their parents (see crow1986basic ). In
turn, this means the average coefficient of inbreeding is the average
coefficient of kinship in the generation prior. With this definition, we can
introduce a new index $j$ that tracks generations and replace the inbreeding
coefficient with the kinship coefficient one generation earlier to rewrite
(17) as:

$\bar{\bar{f}}_{j}=\frac{1}{2^{5}}(\bar{k}+\bar{k}\bar{\bar{f}}_{j-1}),$ (18)

where $\bar{\bar{f}}_{j}$ and $\bar{\bar{f}}_{j-1}$ are the average kinship
coefficients at the current and great great grandparent stage.

The change over generations can then be written as

$\displaystyle\bar{\bar{f_{j}}}-\bar{\bar{f}}_{j-1}$
$\displaystyle=\frac{1}{2^{5}}(\bar{k}+(\bar{k}-2^{5})\cdot\bar{\bar{f}}_{j-1})$
(19)
$\displaystyle=\frac{\bar{k}}{2^{5}}(1+(1-\frac{2^{5}}{\bar{k}})\cdot\bar{\bar{f}}_{j-1}),$
(20)

where (20) simply factors out positive terms (note that $\bar{k}$ approaches
zero only as $n$ approaches infinity). In turn, this means the sign of
$\bar{\bar{f_{j}}}-\bar{\bar{f}}_{j-1}$ is fully determined by the other
factor. In particular, the inequality

$\bar{\bar{f_{j}}}-\bar{\bar{f}}_{j-1}<0$ (21)

will hold precisely when

$1+(1-\frac{2^{5}}{\bar{k}})\cdot\bar{\bar{f}}_{j-1}<0.$ (22)

In turn, noting that $2^{5}/\bar{k}$ is always greater than unity, this is
equivalent to

$\bar{\bar{f}}_{j-1}\geq\frac{-1}{1-\frac{2^{5}}{\bar{k}}}:=f_{critical}.$
(23)

In other words, if in the past the average kinship coefficient is greater than
this $f_{critical}$ then the future average kinship will be smaller. The
converse is also true. Note that $f_{critical}$, due to $\bar{k}$, is a
function of $n$. In this way, if the population is neither growing nor
decreasing then $f_{critical}$ is an equilibrium. If $n$ is changing (perhaps,
for example, due to drought), then so too will $f_{critical}$. In this way,
the point towards which average kinship is pushed is in fact dynamic. This is
shown in Fig. 5 where the average kinship coefficient $\bar{\bar{f}}$ is
plotted alongside $f_{critical}$, both as functions of $n$. Note that, for
large $n$, $f_{critical}$ and $\bar{\bar{f}}$ coincide. This is not true when
$n$ is small. This means that for small $n$ the population is pushed toward a
state that resembles one in which ancestors were slightly related. Most
importantly, however, observe that average kinship coefficients may both
increase and decrease in reaction to historical values, and changes in $n$.

## Costs and benefits

High kinship coefficients, and so too high inbreeding coefficients (a zygote
may be viewed as the sampling of genes from parents), are well known to incur
costs arising from increased incidence of recessive diseases crow1986basic .
In this way, the value of decreasing $\bar{\bar{f}}$ is immediately clear. The
reason why it may be desirable for a system to sometimes tend to increase
$\bar{\bar{f}}$ is less obvious. Note, however, that the coefficient of
kinship we have calculated is two times the coefficient of relatedness
crow1986basic , a quantity that is critical to our understanding of the
evolution of altruism gardner2011genetical . More specifically, Hamilton’s
Rule tells us that we may expect a costly altruistic act to spread through a
population if the benefit to others, weighted by this relatedness, outweighs
the costs to the self frank2013naturalhist . An increase in $\bar{\bar{f}}$,
then, we should expect to also lead to an increase in cooperation, as
understood through social evolution theory gardner2011genetical . In turn this
means that the Gamilaraay system dynamically trades off the costs of being too
closely related against the benefits of cooperative behaviours.

Figure 5: Average coefficient of kinship $\bar{\bar{f}}$ of a randomly sampled
couple operating under the Gamilaraay system as a function of $n$, one eighth
of the total reproductive population size. Overlayed in pink is
$f_{critical}$, above (below) which future average kinship will decrease
(increase). If the reproductive population size is stable then $f_{critical}$
is a true equilibrium, otherwise it recalibrates with changes in $n$.

## Discussion

Colonial expansion of what was to become Australia, from around the mid-1800s,
led to a flurry of anthropological studies of its original inhabitants
taplin1878narrinyeri ; fison1880kamilaroi ; mathews1898victorian ;
mathews1900divisions ; mathews1901ethnological ; parker1905euahlayi ;
brown1918notes . These works, if nothing else, provided a broad stroke
understanding of kinship on the continent. A more rigorous study came only
decades later with the seminal work of Levi-Strauss levi2020structures .
Indeed, in the appendix to this work, the mathematician André Weil of the
Bourbaki set used group theory to provide deeper insight into precisely who
can marry whom in these sorts of systems levi2020structures . Mathematical
interest, no doubt, was sparked by the oddity of cycles present in rules of
descent typical to the continent. Work steadily increased throughout the
following century (see for example hiatt1967authority ; maddock1969alliance ;
hiatt1996arguments ), with studies on the resulting incest taboos remaining
hotly debated itao2020evolution ; read2020simulation ; itao2020reply . All of
these studies, however, attempt to determine broad kin restrictions. In other
words, following Weil, who can reproduce with whom? Here, instead, we ask how
relatedness will be distributed in a society operating under such kinship
systems. In particular, we take the Gamilaraay kinship system as a case study
and ask: how related on average will a permitted couple be? Further, how might
this average relatedness change over time?

To this end, we started by observing that the Gamilaraay kinship system
categorises all individuals into four distinct subtypes. We then noted that,
coupled with descent cycles, uncertainty in pedigrees enters the Gamilaraay
system only at the great grandparent stage, as in Fig. 2. A consequence of
this is that it is only then that two genealogical spots in a potential
Gamilaraay family tree may be occupied by the same individual. Following this,
we worked out the precise probability that such an event is realised, as a
function of $n$, the number of reproductive individuals in a subtype of
interest. With this, it was then possible to write down the exact distribution
of common or distinct ancestors for a general pedigree at this stage. We then
linked this distribution with kinship coefficients (employed often in the
study of disease wright1922coefficients ; crow1986basic ) to derive an average
kinship coefficient for a general couple operating under this system.

Remarkably, we found that the entire Gamilaraay nation would need to reduce to
just 24 reproductive individuals for couples to be as closely related on
average as full first cousins. This result assumes, however, that inbreeding
coefficients at the great grandparent stage are negligible. While reasonable,
this will of course not always be true. For this reason, we then dropped this
assumption and asked instead how might kinship coefficients change over time.
We found that kinship coefficients may both increase and decrease _in
reaction_ to historical values, and changes in $n$. In particular, future
values will decrease if in the past they were larger than a critical value
which we label $f_{critical}$. Conversely, future values will increase if
instead in the the past they were smaller than this $f_{critical}$. In this
way, the Gamilaraay kinship trades-off the costs of being too closely related
(such as incidence of recessive diseases wright1922coefficients ;
crow1986basic ) against the benefits of cooperative behaviours, accrued
through increasing relatedness, as understood formally through Hamilton’s Rule
frank2012natural ; gardner2011genetical .

K Langloh Parker, in an early colonial account, remarked on this very system
‘the Blacks were early scientists in some of their ideas, being before Darwin
with the evolution theory … I rather think the Central Australians have the
key to it. One old man here was quite an Ibsen with his ghastly version of
heredity.’ parker1905euahlayi . In this way, there was some early recognition
(missed almost entirely by early anthropologists) of the scientific
underpinnings of these systems. For future work, it will prove fruitful to
reanalyse other kinship systems whose utility has not yet been fully
understood. In particular, taking a more probabilistic approach to gauge how
relatedness is distributed in societies operating under a given system,
instead of searching for broad kin restriction rules, will be insightful.

## Competing interests

We have no competing interests.

## Authors’ contributions

JMF carried out the research and wrote the manuscript.

## Acknowledgements

Gratitude is given to Garruu Winanga-li-Gii for his patience in learning about
the Gamilaraay system. The author is also grateful for the hospitality and
support of the Sydney Mathematical Research Institute (SMRI) during which many
revisions were made.

## Funding

JMF is funded by a McKenzie Fellowship at the University of Melbourne.

## References

  * (1) Morgan LH.  Ancient Society, Chicago: Charles H.  Kerr And Company. 1877;.
  * (2) McConvell P, Kelly P, Lacrampe S.  Skin, kin and clan: The dynamics of social categories in Indigenous Australia.  ANU Press; 2018.
  * (3) Lévi-Strauss C.  Les structures élémentaires de la parenté. vol. 2.  Walter de Gruyter GmbH & Co KG; 2020.
  * (4) Itao K, Kaneko K.  Evolution of kinship structures driven by marriage tie and competition.  Proceedings of the National Academy of Sciences. 2020;117(5):2378–2384.
  * (5) Read DW, Parkin R.  Simulation model for kinship structures is ethnographically invalid.  Proceedings of the National Academy of Sciences of the United States of America. 2020;117(17):9166.
  * (6) Itao K, Kaneko K.  Reply to Read and Parkin: Our model correctly expresses the ethnographic nature of the cultural incest taboo and kinship structures.  Proceedings of the National Academy of Sciences. 2020;117(17):9167–9168.
  * (7) Frank SA.  Natural selection. VII. History and interpretation of kin selection theory.  Journal of Evolutionary Biology. 2013;26(6):1151–1184.
  * (8) Fison L, Howitt AW.  Kamilaroi and Kurnai: Group Marriage and Relationship, and Marriage by Elopement… Also the Kurnai Tribe: Their Customs in Peace and War.  George Robertson; 1880.
  * (9) Parker KL.  The Euahlayi Tribe: a study of Aboriginal life in Australia. vol. 1.  Library of Alexandria; 1905.
  * (10) Spearim P.  Garruu (Uncle) Paul Spearim (Winanga-li-Gii) of the Gamilaraay Nation Tells the Story of Guwaymadhan, The Dark Blood Group.  The National Indigenous Times. 2014;13(344):18.
  * (11) Spearim P.  Garruu (Uncle) Paul Spearim (Winanga-li-Gii) of the Gamilaraay Nation tells the story of BUWADJARR and the Creation Story Part 2.  The National Indigenous Times. 2014;13(342):18.
  * (12) Wright S.  Coefficients of inbreeding and relationship.  The American Naturalist. 1922;56(645):330–338.
  * (13) Crow JF, et al.  Basic concepts in population, quantitative, and evolutionary genetics.  WH Freeman and Company; 1986.
  * (14) Gardner A, West SA, Wild G.  The genetical theory of kin selection.  Journal of evolutionary biology. 2011;24(5):1020–1043.
  * (15) Taplin G.  The Narrinyeri: An account of the tribes of South Australian aborigines inhabiting the country around the lakes Alexandrina, Albert, and Coorong, and the lower part of the River Murray: their manners and customs; Also, an account of the mission at Point Macleay.  Wigg; 1878.
  * (16) Mathews RH.  The Victorian Aborigines: their initiation ceremonies and divisional systems.  American Anthropologist. 1898;11(11):325–343.
  * (17) Mathews RH.  Divisions of the South Australian Aborigines.  Proceedings of the American Philosophical Society. 1900;39(161):78–93.
  * (18) Mathews RH.  Ethnological notes on the Aboriginal tribes of the Northern Territory.  Queensland Geographical Journal. 1901;16:69–90.
  * (19) Brown AR.  Notes on the Social Organization of Australian Tribes.  The Journal of the Royal Anthropological Institute of Great Britain and Ireland. 1918;48:222–253.
  * (20) Hiatt LR.  Authority and reciprocity in Australian aboriginal marriage arrangements.  The Australian Journal of Anthropology. 1967;6(10):468.
  * (21) Maddock K.  Alliance and entailment in Australian marriage.  The Australian Journal of Anthropology. 1969;7(1):19.
  * (22) Hiatt I, Hiatt L.  Arguments about Aborigines: Australia and the evolution of social anthropology.  Cambridge University Press; 1996.
  * (23) Frank SA.  Natural selection. IV. The price equation.  Journal of evolutionary biology. 2012;25(6):1002–1019.


# The effect of thermal radiation on singularities in the dark universe

I. Brevik,1 A. V. Timoshkin,2,3 Tanmoy Paul4,5
iver.h.brevik@ntnu.noalex.timosh<EMAIL_ADDRESS>1) Department
of Energy and Process Engineering, Norwegian University of Science and
Technology, N-7491 Trondheim, Norway
2) Tomsk State Pedagogical University, Kievskaja Street, 60, 634061 Tomsk,
Russia
3) International Laboratory of Theoretical Cosmology, Tomsk State University
of Control Systems and Radio Electronics, Lenin Avenue, 36, 634050 Tomsk,
Russia
4) Department of Physics, Chandernagore College, Hooghly - 712 136.
(5) Department of Theoretical Physics, Indian Association for the Cultivation
of Science, 2A $\&$ 2B Raja S.C. Mullick Road, Kolkata - 700 032, India

###### Abstract

Cosmological models with an inhomogeneous viscous dark fluid, coupled with
dark matter in the Friedmann- Robertson-Walker (FRW) flat universe, are
considered. The influence of thermal effects caused by Hawking radiation on
the visible horizon is studied, in connection with the classified type I and
type III singularities which are known to occur within a finite amount of
time. Allowance of thermal effects implies that a transition to a type II
singularity can take place, in a finite time. We take into account a bulk
viscosity of the dark fluid, observing the equation of state in the case of
radiation, and find that there is a qualitative change in the singular
universe of type I: it may pass into a singularity of type III, or it may
avoid the singularity at all.

##  I Introduction

After it was observed that the universe is exposed to an accelerated
expansion, it is of great interest to study the nature of dark energy which is
responsible for this acceleration Nojiri:2017ncd ; Shi:2011sa ;
Perivolaropoulos:2006ce ; Arefeva:2006ido . In the era of dark energy the
universe can be qualitatively described using an exotic effective fluid with
negative pressure (more appropriately called a positive tensile stress), that
satisfies an unusual equation of state Nojiri:2005sx ; Nojiri:2006zh ;
Nojiri:2004pf ; Setare:2007bx ; Capozziello:2005mj ; Cognola:2006eg ;
Nojiri:2005sr ; Nojiri:2020wmh ; Odintsov:2020zct . Dark energy can be
characterized by an equation of state parameter
$\omega_{\mathrm{eff}}=p_{\mathrm{eff}}/\rho_{\mathrm{eff}}$ . Depending on
the values of $\omega_{\mathrm{eff}}$ several possibilities of the universe
are possible: for example, phantom behaviour occurs for
$\omega_{\mathrm{eff}}<-1$, de-Sitter evolution occurs for
$\omega_{\mathrm{eff}}=-1$, and quintessential behaviour occurs for
$-1<\omega_{\mathrm{eff}}<-1/3$. The experimental value of the parameter of
the equation of state is however determined with insufficient accuracy to be
able to unambiguously determine the phase in which our Universe is located.
Today the value of this parameter lies within the following limits
$\omega_{\mathrm{eff}}=1.04^{+0.09}_{-0.10}$ Nakamura:2010zzi .

One of the properties of phantom dark energy is the prediction of a Big Rip
singularity (type I) in the future. For this kind of singularity the scale
factor and the Hubble function go to infinity at a finite time, called $t_{s}$
or $t_{\mathrm{rip}}$. This is the most destructive type of singularity
Caldwell:1999ew . There are several less drastic future singularities, namely
soft singularities classified as type II, III and IV. All these singularities
are purely classic in nature.

From a physical viewpoint, an increase in the Hubble function should be
expected to lead to an increase in temperature. At high temperatures,
especially near the singularity, thermal radiation should appear. Thermal
radiation is associated with Hawking radiation, which effectively should be
generated at the apparent horizon of the FRW universe Cai:2009ph . Hawking’s
radiation is manifested in black holes and is associated with the existence of
the visible black hole horizon, as well as the visible horizon of cosmic
events in de-Sitter space. Hawking’s thermal spectrum radiation should be
taken to appear in the late universe at high temperatures shortly before its
rupture. Accounting for thermal radiation will allow a qualitative change in
the classical description and give a more realistic picture of the future of
the universe.

Recently, the effect of thermal radiation on future singularities of types I,
II, III and IV was studied in Ref. Nojiri:2020sti . It was shown that with
singular universes of types I and III, as well as for the Little Rip universe,
there occurs a qualitative change in the singularity due to thermal effects.
The singularities end up as type II singularities. In universes of types II
and IV there is no qualitative change in the final state.

The model of a non-viscous fluid in cosmology is an idealized case. Dark
energy universe with a viscous fluid was studied in Refs. Gron:1990ew ;
Brevik:2017msy ; Cataldo:2005qh ; Brevik:2005gx ; Brevik:2001ed ;
Brevik:2010jv ; Brevik:2005bj ; Li:2009mf ; Brevik:2017juz . The cosmic
viscosity property was taken into account in connection with the Big Rip type
singularity Nojiri:2003vn , and also in connection with the singularities of
type II, III and IV Barrow:2013ria ; Nojiri:2004ip . Models of dark energy
interacting with dark matter, in which the singularities of one of these four
types are formed, were considered in Ref. Myrzakul:2013qka ;
Capozziello:2009hc .

The article Elizalde:2014ova  studied a combined phantom/fluid model
consisting of a viscous dark fluid (dark energy) with a linear inhomogeneous
equation of state and dark matter with a linear homogeneous equation of state
in a spatially flat FRW metric.

The purpose of the present article is to study the influence of thermal
radiation, taking into account the interaction of viscous dark fluid with dark
matter. In particular, we will focus on the change of singular behavior of the
late-time universe.

##  II The effect of thermal radiation on the formation of singularities in
the late universe

Let us consider a spatially-flat FRW universe

$\displaystyle ds^{2}=-dt^{2}+a^{2}(t)\delta_{ij}dx^{i}dx^{j},$ (1)

where $a(t)$ is the scale factor. We will be interested in the case when the
effective parameter in the equation of state takes values in the vicinity of
$-1$. Then the following kinds of evolution of the accelerating universe are
possible : phantom, quintessence, or a de Sitter expansion. The question
arises about how this evolution will end up in the future. The answer depends
on the behavior of the time-dependent parameters of the effective equation of
state. We will be interested in dark energy universes, in which there occur
future singularities within a finite, or an infinite, time. For such universes
the Nojiri-Odintsov-Tsujikawa classification was given in Ref. Nojiri:2005sx
(see also Odintsov:2018uaw ). Singularities arise when one or more of the
central cosmological parameters diverge: the scale factor $a(t)$, the
effective (total) energy density $\rho_{\mathrm{eff}}$, the effective (total)
pressure $p_{\mathrm{eff}}$, or higher derivatives of the Hubble function.

In the limit $t\rightarrow t_{s}$, the following classes of singularities can
be distinguished:

  * •

Type I (Big Rip): $a\rightarrow\infty$, $\rho_{\mathrm{eff}}\rightarrow\infty$
and $p_{\mathrm{eff}}\rightarrow\infty$. This class of singularities includes
the case when $\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$ are finite at
$t=t_{s}$. A Big Rip leads to the decay of gravitationally bound objects large
on a cosmological scale.

  * •

Type II (“sudden” singularity): $a\rightarrow a_{s}$ ,
$\rho_{\mathrm{eff}}\rightarrow\rho_{s}$ and
$\big{|}p_{\mathrm{eff}}\big{|}\rightarrow\infty$, where $a_{s}\neq 0$ and
$\rho_{s}$ are constant. That is a pressure singularity.

  * •

Type III: $a\rightarrow a_{s}$, $\rho_{\mathrm{eff}}\rightarrow\infty$ and
$p_{\mathrm{eff}}\rightarrow\infty$. This type singularity is milder than Type
I but stronger than Type II.

  * •

Type IV: $a\rightarrow a_{s}$, $\rho_{\mathrm{eff}}\rightarrow 0$ and
$p_{\mathrm{eff}}\rightarrow 0$, but the higher derivatives of the Hubble
function $H$ diverge. This type also includes the case where
$\rho_{\mathrm{eff}}$ and/or $p_{\mathrm{eff}}$ are finite for $t=t_{s}$.

Here $\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$ can be calculated by the
following expressions,

$\displaystyle\rho_{\mathrm{eff}}=\frac{3}{\kappa^{2}}H^{2}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}p_{\mathrm{eff}}=-\frac{2}{\kappa^{2}}\big{(}2\dot{H}+3H^{2}\big{)}$
(2)

where $\kappa^{2}=8\pi G$ and $H=\frac{\dot{a}}{a}$ is the Hubble parameter.
It may be mentioned that the effective energy density $\rho_{\mathrm{eff}}$
and the effective pressure $p_{\mathrm{eff}}$ may include the contribution
from the modified gravity.

However, the singularity is not the only possible outcome of the evolution of
our universe in the phantom phase. It was shown in Odintsov:2018uaw ;
Frampton:2011aa ; Frampton:2011sp ; Frampton:2011rh  that if the cosmic energy
density remains constant or monotonically increases, then, depending on the
asymptotic behavior of the Hubble parameter $H$ Odintsov:2018uaw , all
possible types of evolution of our universe can be divided into four
categories:

  * •

Big Rip: $H(t)\rightarrow\infty$ when $t=t_{s}<\infty$.

  * •

Little Rip: $H(t)\rightarrow\infty$ when $t\rightarrow\infty$.

  * •

Cosmological constant : $H(t)=\mathrm{constant}$.

  * •

Pseudo Rip: $H(t)\rightarrow H_{\mathrm{}\infty}$ when $t\rightarrow\infty$,
where $H_{\mathrm{\infty}}$ is a constant.

Here we would like to mention that both the Little Rip and Pseudo Rip models
are nonsingular.

We will study the cosmological models induced by the inhomogeneous viscous
dark fluids coupled with dark matter, in terms of the parameters appearing in
the equation of state (EoS). Let us consider the following formulation of the
EoS of an inhomogeneous viscous fluid in flat FRW space-time
Capozziello:2005pa , namely

$\displaystyle p=\omega(\rho,t)\rho-3H\xi(H,t)~{}~{}~{},$ (3)

where $\xi(H,t)$ is the bulk viscosity, which depends on the Hubble parameter
$H$ and on the cosmic time $t$. According to the thermodynamic set up, we
naturally assume that $\xi(H,t)>0$.

We will take the following form for the thermodynamic (EoS) parameter $\omega$
Capozziello:2005pa ,

$\displaystyle\omega(\rho,t)=\omega_{1}(t)\big{(}A_{0}\rho^{\alpha-1}-1\big{)}$
(4)

where $A_{0}\neq 0$ and $\alpha\geq 1$ are constants. (Note that $A_{0}$ is
nondimensional only if $\alpha=1$. If $\alpha=3/2$ for instance, the dimension
will be cm2 in geometric units.)

We choose the bulk viscosity as Capozziello:2005pa ,

$\displaystyle\xi(H,t)=\xi_{1}(t)\big{(}3H\big{)}^{n}$ (5)

with $n>0$.

Let us consider the influence of thermal effects on the change in
singularities of types I, taking into account the viscosity property of a dark
fluid and its interaction with dark matter. Since the temperature of the
universe increases near the singularity, thermal radiation is generated, as
mentioned above. From statistical physics, the energy density of thermal
radiation is proportional to the fourth power of the absolute temperature.
Therefore, near the future singularity where the Hubble parameter becomes very
high, we assume that the thermal energy density has the form Nojiri:2020sti

$\displaystyle\rho_{\mathrm{rad}}=\lambda H^{4}~{}~{},$ (6)

where $\lambda$ is a positive constant.

Taking thermal radiation into account, the FRW equation is modified as follows
Nojiri:2020sti ,

$\displaystyle\frac{3}{\kappa^{2}}H^{2}=\rho_{\mathrm{eff}}+\lambda H^{4}$ (7)

From Eq.(7) it follows that when the evolution time of the late universe is
much less than the singularity time, the first term of the equation makes the
greatest contribution. While near the singularity time, the second term makes
the largest contribution. We analyze equation (7) further, by solving it with
respect to the square of the Hubble parameter $H^{2}$

$\displaystyle
H^{2}=\frac{1}{2\lambda}\bigg{[}\frac{3}{\kappa^{2}}\pm\sqrt{\frac{9}{\kappa^{4}}-4\lambda\rho_{\mathrm{eff}}}\bigg{]}~{}~{}.$
(8)

In the following, we will apply the cosmological models of a viscous fluid
from the article Elizalde:2014ova .

##  III Singular behavior of late-time Universe taking into account the
viscosity of a fluid and its interaction with dark matter

We start from the simplest, constant case, namely $\omega(\rho,t)=\omega_{0}$
and will consider different forms for the bulk viscosity.

###  III.1 Constant viscosity

Let us consider the case of constant bulk viscosity $\xi(H,t)=\xi_{0}>0$. The
Hubble function has the form Elizalde:2014ova

$\displaystyle
H(t)=\frac{\xi_{0}\kappa^{2}}{\big{(}1+\omega_{0}\big{)}\bigg{(}1-\sqrt{C_{1}}~{}\exp{\big{[}3\xi_{0}\kappa^{2}t/2\big{]}}\bigg{)}},$
(9)

with $C_{1}$ a nondimensional constant. This model does not so far take into
account interactions with dark matter. Observe that $H$ diverges for
$t\rightarrow
t_{s}=-\frac{2}{3\xi_{0}\kappa^{2}}\ln{\big{(}\sqrt{C_{1}}\big{)}}$, thus a
Big Rip singularity appears. Let us now see if the type of singularity will
change if we take into account thermal radiation near the singularity time.

We calculate the scale factor,

$\displaystyle a(t)=e^{\int
Hdt}=a_{0}\bigg{(}1-\frac{1}{\sqrt{C_{1}}}~{}\exp{\big{[}-3\xi_{0}\kappa^{2}t/2\big{]}}\bigg{)}^{\frac{2}{3(1+\omega_{0})}}~{}~{},$
(10)

where $a_{0}$ is an integration constant. Then the effective energy density in
terms of the scale factor is

$\displaystyle\rho_{\mathrm{eff}}=\bigg{(}\sqrt{\frac{3}{\rho_{0}}}~{}\frac{\xi_{0}\kappa}{1+\omega_{0}}\bigg{)}^{2}\bigg{[}\frac{a(t)}{a_{0}}\bigg{]}^{3(1+\omega_{0})}~{}e^{-3\xi_{0}\kappa^{2}t}~{}~{}.$
(11)

Let us return to Eq.(8). Since $H^{2}$ is a real number, we must have

$\displaystyle\frac{9}{\kappa^{4}}-4\lambda\rho_{\mathrm{eff}}\geq 0~{}~{}.$
(12)

We get in this model

$\displaystyle\frac{9}{\kappa^{4}}-4\lambda
Ae^{-Bt}\bigg{[}\frac{a(t)}{a_{0}}\bigg{]}^{C}\geq 0~{}~{},$ (13)

where
$A=\bigg{(}\sqrt{\frac{3}{\rho_{0}}}\frac{\xi_{0}\kappa}{1+\omega_{0}}\bigg{)}^{2}$,
$B=3\xi_{0}\kappa^{2}$ and $C=3(1+\omega_{0})$. The above inequality puts a
restriction on the scale factor,

$\displaystyle a(t)\leq a_{0}\bigg{(}\frac{9}{4\lambda
A\kappa^{4}}\bigg{)}^{\frac{1}{C}}~{}e^{\frac{B}{C}t}$ (14)

under the condition $C<0$ ($\omega_{0}<-1$), which corresponds to the case of
phantom dark energy.

Taking into account the thermal radiation $\rho_{\rm rad}$, we obtain from Eq.
(8) that there exists another upper limit $a_{\rm max}$ for the scale
parameter,

$\displaystyle a(t)\leq a_{\mathrm{max}}=a_{0}\bigg{(}\frac{4\lambda
A\kappa^{4}}{9}\bigg{)}^{\frac{1}{3(1+\omega_{0})}}~{}~{},$ (15)

which corresponds to the instant $t_{\mathrm{max}}$,

$\displaystyle
t_{\mathrm{max}}=-\frac{2}{3\xi_{0}\kappa^{2}}\ln{\bigg{\\{}\sqrt{C_{1}}\bigg{(}1-\frac{2\kappa^{2}\sqrt{A\lambda}}{3}\bigg{)}\bigg{\\}}}~{}~{}.$
(16)

This is thus a singularity of another type than that arising from the Hawking
radiation.

We calculate the difference between $t_{\mathrm{max}}$ and $t_{s}$

$\displaystyle
t_{\mathrm{max}}-t_{s}=-\frac{2}{3\xi_{0}\kappa^{2}}\ln{\bigg{\\{}1-\frac{2\kappa^{2}\sqrt{A\lambda}}{3}\bigg{\\}}}>0~{}~{},$
(17)

which shows that $t_{\mathrm{max}}$ is larger than $t_{s}$. In the limit
$t\rightarrow t_{\mathrm{max}}$, $a\rightarrow a_{\mathrm{max}}$. From Eqs.(2)
and (3), one can calculate the effective energy density $\rho_{\mathrm{eff}}$
and the effective pressure $p_{\mathrm{eff}}$. Then in the limit $t\rightarrow
t_{\mathrm{max}}$, the effective energy density and effective pressure become

$\displaystyle\rho_{\mathrm{max}}=\rho_{\mathrm{eff}}(t_{\mathrm{max}})=\frac{9\xi_{0}^{2}\kappa^{2}}{(1+\omega_{0})^{2}}\bigg{[}\frac{1}{1+C_{1}\bigg{(}\frac{2\kappa^{2}\sqrt{A\lambda}}{3}-1\bigg{)}}\bigg{]}^{2}~{}~{},$
(18)

and

$\displaystyle
p_{\mathrm{max}}=p_{\mathrm{eff}}(t_{\mathrm{max}})=\bigg{|}\omega_{0}\rho_{\mathrm{max}}-3\xi_{0}H_{\mathrm{max}}\bigg{|}=\frac{3\xi_{0}^{2}\kappa^{2}}{(1+\omega_{0})}\bigg{\\{}\frac{\bigg{|}\frac{2\omega_{0}-1}{1+\omega_{0}}-C_{1}\bigg{(}\frac{2\kappa^{2}\sqrt{A\lambda}}{3}-1\bigg{)}\bigg{|}}{\bigg{[}1+C_{1}\bigg{(}\frac{2\kappa^{2}\sqrt{A\lambda}}{3}-1\bigg{)}\bigg{]}^{2}}\bigg{\\}}~{}~{}.$
(19)

respectively. In the general case the values of the scale factor, energy
density and effective pressure turn out to be finite, but higher derivatives
of H diverge. Thus, a cosmological finite-time future singularity is not
formed. This behavior is due to the influence from the viscosity of the dark
fluid, which compensates for the effect of thermal radiation. However, if the
radiation parameter $\lambda$ goes to $\lambda_{0}$, with

$\displaystyle\lambda_{0}=\frac{3}{C_{1}}\bigg{(}\frac{(C_{1}-1)(1+\omega_{0})}{2\xi_{0}\kappa^{3}}\bigg{)}^{2}~{}~{},$
(20)

then $\rho_{\mathrm{max}}\rightarrow\infty$ and
$\big{|}p_{\mathrm{max}}\big{|}\rightarrow\infty$. This is a type III
singularity. It is milder than type I, but stronger than type II.

If $\zeta_{0}\rightarrow\infty$, then $\lambda_{0}\rightarrow 0$. Thus, the
viscosity weakens the effect of thermal radiation.

###  III.2 Viscosity proportional to the Hubble parameter

Let us consider the case where the viscosity is proportional to the Hubble
parameter, $\xi(H,t)=3\tau H$, the constant $\tau$ being positive. The Hubble
parameter becomes Elizalde:2014ova

$\displaystyle
H(t)=\frac{\kappa}{\sqrt{3}}\bigg{\\{}\frac{\delta{\gamma}~{}\sqrt{C_{1}}}{3\theta\sqrt{C_{1}}+\exp{\big{[}-\tilde{\eta}t/2\big{]}}}\bigg{\\}}~{}~{}.$
(21)

Here the following designations are introduced:

$\eta=\delta\gamma^{2},\quad\gamma=\frac{\kappa}{\sqrt{3}}\sqrt{1+\frac{1}{r}},\quad\frac{\tilde{\eta}}{\eta}=r,\quad\theta=1+\omega_{0}-9\tau{\gamma}^{2},$

where the constant $r$ is associated with the influence of dark matter and is
equal to the ratio of the energy density of dark matter to the density of dark
energy. The dimensions are $[\gamma]=~{}$cm, $[\eta]=[\tilde{\eta}]=$ cm-1,
$[\delta]=$cm-2, $[\tau]=~{}$cm-2. Further, the constant parameter $\delta$ is
responsible for the interaction with dark matter. If
$\omega_{0}<-1+9\tau{\gamma}^{2}$, then $\theta<0$ and consequently $H$
diverges at $t\rightarrow
t_{s}=-\frac{2}{\tilde{\eta}}\ln{\bigg{(}-3\theta\sqrt{C_{1}}\bigg{)}}$, which
leads to the appearance of a singularity of the type Big Rip. We consider
again the behavior of the late-time universe near the singularity, taking into
account the effect of thermal radiation.

Let us calculate the scale factor,

$\displaystyle a(t)=a_{0}~{}\exp{\big{[}\frac{\delta{\gamma}~{}\kappa
t}{3\sqrt{3}~{}\theta}\big{]}}\bigg{(}3\theta\sqrt{C_{1}}+\exp{\big{[}-\tilde{\eta}t/2\big{]}}\bigg{)}^{\frac{2\delta{\gamma}~{}\kappa}{3\sqrt{3}~{}\theta\tilde{\eta}}}$
(22)

and express the effective energy density in terms of the scale factor as

$\displaystyle\rho_{\mathrm{eff}}=C_{1}~{}\big{(}\delta{\gamma}\big{)}^{2}\bigg{[}\frac{a(t)}{a_{0}}\bigg{]}^{-\frac{2}{\tilde{\delta}}}~{}\exp{\big{[}\frac{2\alpha^{\prime}}{\tilde{\delta}}t\big{]}}~{}~{},$
(23)

where $\alpha^{\prime}=\frac{\delta\gamma\kappa}{3\sqrt{3}\theta}$ and
$\tilde{\delta}=\frac{2\delta\gamma\kappa}{3\sqrt{3}\theta}$. The inequality
in Eq.(12) implies the following restriction on the scale factor,

$\displaystyle a(t)\leq
a_{0}~{}\bigg{(}\frac{2\delta{\gamma}~{}\kappa^{2}\sqrt{\lambda
C_{1}}}{3}\bigg{)}^{\tilde{\delta}}~{}e^{\alpha^{\prime}t}~{}~{}.$ (24)

Since $\theta<0$, then $\alpha^{\prime}<0$ and scale factor values are limited
by the maximum number $a_{\mathrm{max}}$ given by

$\displaystyle a(t)\leq
a_{\mathrm{max}}=a_{0}~{}\bigg{(}\frac{2\delta{\gamma}~{}\kappa^{2}\sqrt{\lambda
C_{1}}}{3}\bigg{)}^{\tilde{\delta}}~{}~{},$ (25)

which corresponds to the instant

$\displaystyle
t_{\mathrm{max}}=\frac{2}{\tilde{\eta}}~{}\ln{\bigg{\\{}\frac{2\delta{\gamma}~{}\kappa^{2}\sqrt{\lambda
C_{1}}-3}{9\theta\sqrt{C_{1}}}\bigg{\\}}}~{}~{},$ (26)

Let us find the difference between $t_{s}$ and $t_{\mathrm{max}}$,

$\displaystyle
t_{s}-t_{\mathrm{max}}=-\frac{2}{\tilde{\eta}}~{}\ln{\bigg{\\{}1-\frac{2}{3}\delta{\gamma}~{}\kappa^{2}\sqrt{\lambda
C_{1}}\bigg{\\}}}>0~{}~{}.$ (27)

Hence it follows that $t_{s}$ is larger than $t_{\mathrm{max}}$. Thus in the
limit $t\rightarrow t_{\mathrm{max}}$, the effective energy density and the
effective pressure are given by

$\displaystyle\rho_{\mathrm{max}}=\rho_{\mathrm{eff}}(t_{\mathrm{max}})=\frac{1}{\theta^{2}}\bigg{(}\delta{\gamma}-\frac{1}{\kappa^{2}\sqrt{\lambda
C_{1}}}\bigg{)}^{2}~{}~{},$ (28)

and

$\displaystyle\big{|}p_{\mathrm{max}}\big{|}=\big{|}p(t_{\mathrm{max}})\big{|}=\big{|}\omega_{0}\rho_{\mathrm{max}}-9\tau
H_{\mathrm{max}}^{2}\big{|}=\big{|}\omega_{0}-3\tau\kappa^{2}\big{|}\rho_{\mathrm{max}}~{}~{},$
(29)

respectively. In the limit when the radiation parameter
$\lambda\rightarrow\lambda_{0}$ where $\lambda_{0}$ is given by

$\displaystyle\lambda_{0}=\bigg{(}\frac{3}{\sqrt{\rho_{0}}\,\delta\,\tilde{\gamma}~{}\kappa^{2}}\bigg{)}^{2}~{}~{},$
(30)

$\rho_{\mathrm{max}}\rightarrow 0$ and
$\big{|}p_{\mathrm{max}}\big{|}\rightarrow 0$. However, the higher derivatives
of the Hubble function do not diverge. Thus, the formation of a singularity of
IV type does not occur in this model.

We see that viscosity softens the singularity or avoids it altogether. This is
due to a decrease in the pressure of a viscous fluid due to the viscosity term
in the equation of state.

###  III.3 Inhomogeneous fluid with variable parameter $\omega$

In this section we will assume that the thermodynamic parameter
$\omega(\rho,t)$ is a function of the energy density of the fluid. Let us
choose it to have the form

$\displaystyle\omega(\rho,t)=A_{0}\rho^{\alpha-1}-1~{}~{},$ (31)

where $A_{0}\neq 0$ is a dimensional constant. The bulk viscosity is taken to
be proportional to $H^{n}$,

$\displaystyle\xi(H,t)=\tau\big{(}3H\big{)}^{n}$ (32)

with $\tau$ and $n$ positive. In the case $n=2\alpha-1$, the energy density
becomes Elizalde:2014ova

$\displaystyle H(t)=H_{0}\bigg{\\{}\rho_{0}\exp{[(\alpha-\frac{1}{2})\eta
t]}+\frac{\mu}{\tilde{\eta}}\bigg{\\}}^{\frac{2}{1-2\alpha}}~{}~{},\quad
H_{0}=\frac{\kappa}{\sqrt{3}}\sqrt{\rho_{0}},$ (33)

where $\mu$ is a dimensionless constant and $\alpha\neq\frac{1}{2}$. Here
$t_{0}$ is the present time, and $\rho_{0}=\rho(t_{0})$, The dimension is
$[\rho_{0}]=~{}$cm-4.

For $\alpha>\frac{1}{2}$, then at $t\rightarrow
t_{s}=\frac{1}{\tilde{\eta}\big{(}\alpha-1/2\big{)}}\ln{\bigg{(}-\frac{\mu}{C_{1}\tilde{\eta}}\bigg{)}}$,
the Hubble parameter diverges and we obtain again the Big Rip singularity.
Further, if we consider the case $\alpha=3/2$ the scale factor turns out to be

$\displaystyle a(t)=a_{0}\bigg{(}1+\frac{\mu}{C_{1}\tilde{\eta}}e^{-\eta
t}\bigg{)}^{\frac{rH_{0}}{\mu}}$ (34)

and consequently the energy density in terms of the scale factor becomes

$\displaystyle\rho_{\mathrm{eff}}=\rho_{0}(C_{1}e^{\eta
t})^{-2}\bigg{[}\frac{a(t)}{a_{0}}\bigg{]}^{\frac{2\mu}{rH_{0}}}$ (35)

The inequality in Eq.(12) implies the following restriction on the scale
factor,

$\displaystyle a(t)\leq
a_{0}\bigg{(}\frac{2\kappa^{2}\sqrt{\lambda\rho_{0}}}{3C_{1}}\bigg{)}^{\frac{rH_{0}}{\mu}}\exp[\frac{rH_{0}\eta
t}{\mu}]$ (36)

Scale factor values are limited by the maximum number $a_{\mathrm{max}}$,

$\displaystyle a(t)\leq
a_{\mathrm{max}}=a_{0}\bigg{(}\frac{2\kappa^{2}\sqrt{\lambda\rho_{0}}}{3C_{1}}\bigg{)}^{\frac{rH_{0}}{\mu}}$
(37)

which corresponds to the time

$\displaystyle
t_{\mathrm{max}}=\frac{1}{\eta}~{}\ln\frac{3\mu}{\tilde{\eta}\big{(}2\kappa^{2}\sqrt{\lambda\rho_{0}}-3C_{1})}~{}~{}.$
(38)

The difference between $t_{s}$ and $t_{\mathrm{max}}$ takes the following
form,

$\displaystyle
t_{s}-t_{\mathrm{max}}=\frac{1}{\eta}~{}\ln\big{(}1-\frac{2\kappa^{2}}{3C_{1}}\sqrt{\lambda\rho_{0}}\big{)}<0~{}~{}.$
(39)

It turns out that $t_{\mathrm{max}}$ is larger than $t_{s}$. Thus, it due to
thermal radiation, the time of formation of the cosmological singularity
changes qualitatively. In the limit $t\rightarrow t_{\mathrm{max}}$, the
effective energy density and the effective pressure are given by

$\displaystyle\rho_{\mathrm{max}}=\rho_{\mathrm{eff}}(t_{\mathrm{max}})=\rho_{0}\left[\frac{\tilde{\eta}}{\mu}\left(1-\frac{3C_{1}}{2\kappa^{2}\sqrt{\lambda\rho_{0}}}\right)\right],$
(40)

and

$\displaystyle\big{|}p_{\mathrm{max}}\big{|}=\big{|}p(t_{\mathrm{max}})\big{|}=\bigg{|}\bigg{(}A_{0}\rho_{\mathrm{max}}^{\frac{1}{2}}-1\bigg{)}\rho_{\mathrm{max}}-\tau\big{(}3H_{\mathrm{max}}\big{)}^{3}\bigg{|}=\bigg{|}\bigg{(}A_{0}-3\sqrt{3}~{}\tau\kappa^{3}\bigg{)}\rho_{\mathrm{max}}^{\frac{1}{2}}-1\bigg{|}\rho_{\mathrm{max}}~{}~{},$
(41)

respectively. Thereby, the value of the energy density and effective pressure
becomes finite, while higher derivatives of the Hubble function diverge.
Consequently, a cosmological singularity is not formed.

##  IV Conclusion

We have investigated the singular behavior of the dark universe, taking into
account the thermal effects caused by Hawking radiation, the viscosity
properties of the dark fluid, and its interaction with dark matter on the
visible horizon of the FRW universe. According to the study carried out in
Nojiri:2020sti  for an ideal fluid, near the singularity it is necessary to
take into account the Hawking thermal radiation, which leads to a change in
the type of the singularity. In a dark universe with singularities of types I
and III with a finite formation time, a transition to a type II singularity
occurs.

Models with an inhomogeneous viscous dark fluid interacting with dark matter
were considered in Elizalde:2014ova . It is shown that, in our case, a
transition from singularity of type I to singularity III type is possible due
to the influence of thermal radiation. Singularities may be absent, due to the
viscosity of the dark fluid and its interaction with dark matter. The absence
of a singularity in some models is explained by the fact that both the
presence of a bulk viscosity in the equation of state of a dark fluid, and the
presence of a thermal radiation term in the Friedmann equation, are
proportional to a power of the Hubble parameter. As a result, the effect of
thermal radiation near the singularity may be neutralized by a viscous fluid.

As was shown in Ref. Astashenok:2021 , taking thermal radiation into account
does not weaken the agreement of cosmological models with astronomical
obeservations.

### Acknowledgments

This work was supported in part by Ministry of Education of Russian
Federation, Project No FEWF-2020-0003 (A.V.T.).

## References

  * (1) S. Nojiri, S. D. Odintsov and V. K. Oikonomou, Modified gravity theories on a nutshell: Inflation, bounce and late-time evolution, Phys. Rept. 692 (2017), 1-104 doi:10.1016/j.physrep.2017.06.001 [arXiv:1705.11098 [gr-qc]].
  * (2) Y. Shi, A cyclic cosmological model based on the $f(\rho)$ modified theory of gravity, [arXiv:1106.0341 [physics.gen-ph]].
  * (3) L. Perivolaropoulos, Accelerating universe: observational status and theoretical implications, AIP Conf. Proc. 848 (2006) no.1, 698-712 doi:10.1063/1.2348048 [arXiv:astro-ph/0601014 [astro-ph]].
  * (4) I. Y. Aref’eva and I. V. Volovich, On the null energy condition and cosmology, Theor. Math. Phys. 155 (2008), 503-511 doi:10.1007/s11232-008-0041-8 [arXiv:hep-th/0612098 [hep-th]].
  * (5) S. Nojiri, S. D. Odintsov and S. Tsujikawa, Properties of singularities in (phantom) dark energy universe, Phys. Rev. D 71 (2005), 063004 doi:10.1103/PhysRevD.71.063004 [arXiv:hep-th/0501025 [hep-th]].
  * (6) S. Nojiri and S. D. Odintsov, The new form of the equation of state for dark energy fluid and accelerating universe, Phys. Lett. B 639 (2006), 144-150 doi:10.1016/j.physletb.2006.06.065 [arXiv:hep-th/0606025 [hep-th]].
  * (7) S. Nojiri and S. D. Odintsov, The final state and thermodynamics of dark energy universe, Phys. Rev. D 70 (2004), 103522 doi:10.1103/PhysRevD.70.103522 [arXiv:hep-th/0408170 [hep-th]].
  * (8) M. R. Setare, Interacting generalized Chaplygin gas model in non-flat universe, Eur. Phys. J. C 52 (2007), 689-692 doi:10.1140/epjc/s10052-007-0405-5 [arXiv:0711.0524 [gr-qc]].
  * (9) S. Capozziello, S. Nojiri and S. D. Odintsov, Dark energy: The equation of state description versus scalar-tensor or modified gravity, Phys. Lett. B 634 (2006), 93-100 doi:10.1016/j.physletb.2006.01.065 [arXiv:hep-th/0512118 [hep-th]].
  * (10) G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov and S. Zerbini, Dark energy in modified Gauss-Bonnet gravity: Late-time acceleration and the hierarchy problem, Phys. Rev. D 73 (2006), 084007 doi:10.1103/PhysRevD.73.084007 [arXiv:hep-th/0601008 [hep-th]].
  * (11) S. Nojiri and S. D. Odintsov, Inhomogeneous equation of state of the universe: Phantom era, future singularity and crossing the phantom barrier, Phys. Rev. D 72 (2005), 023003 doi:10.1103/PhysRevD.72.023003 [arXiv:hep-th/0505215 [hep-th]].
  * (12) S. Nojiri, S. D. Odintsov, V. K. Oikonomou and T. Paul, Unifying holographic inflation with holographic dark energy: a covariant approach, Phys. Rev. D 102 (2020) no.2, 023540 doi:10.1103/PhysRevD.102.023540 [arXiv:2007.06829 [gr-qc]].
  * (13) S. D. Odintsov, V. K. Oikonomou and T. Paul, From a bounce to the dark energy era with $F(R)$ gravity, Class. Quant. Grav. 37 (2020) no.23, 235005 doi:10.1088/1361-6382/abbc47 [arXiv:2009.09947 [gr-qc]].
  * (14) K. Nakamura et al. [Particle Data Group], Review of particle physics, J. Phys. G 37 (2010), 075021 doi:10.1088/0954-3899/37/7A/075021
  * (15) R. R. Caldwell, A phantom menace?, Phys. Lett. B 545 (2002), 23-29 doi:10.1016/S0370-2693(02)02589-3 [arXiv:astro-ph/9908168 [astro-ph]].
  * (16) R. G. Cai and N. Ohta, Horizon thermodynamics and gravitational field equations in Horava-Lifshitz Gravity, Phys. Rev. D 81 (2010), 084061 doi:10.1103/PhysRevD.81.084061 [arXiv:0910.2307 [hep-th]].
  * (17) S. Nojiri and S. D. Odintsov, The dark universe future and singularities: the account of thermal and quantum effects, Phys. Dark Univ. 30 (2020), 100695 doi:10.1016/j.dark.2020.100695 [arXiv:2006.03946 [gr-qc]]; A. Astashenok, S. D. Odintsov and V. K. Oikonomou, submitted.
  * (18) Ø. Grøn, Viscous inflationary universe models, Astrophys. Space Sci. 173 (1990), 191-225 doi:10.1007/BF00643930
  * (19) I. Brevik, Ø. Grøn, J. de Haro, S. D. Odintsov and E. N. Saridakis, Viscous cosmology for early- and late-time Universe, Int. J. Mod. Phys. D 26 (2017) no.14, 1730024 doi:10.1142/S0218271817300245 [arXiv:1706.02543 [gr-qc]].
  * (20) M. Cataldo, N. Cruz and S. Lepe, Viscous dark energy and phantom evolution, Phys. Lett. B 619 (2005), 5-10 doi:10.1016/j.physletb.2005.05.029 [arXiv:hep-th/0506153 [hep-th]].
  * (21) I. Brevik, J. M. Børven and S. Ng, Viscous brane cosmology with a brane-bulk energy interchange term, Gen. Rel. Grav. 38 (2006), 907-915 doi:10.1007/s10714-006-0271-8 [arXiv:gr-qc/0512026 [gr-qc]].
  * (22) I. Brevik and S. D. Odintsov, On the Cardy-Verlinde entropy formula in viscous cosmology, Phys. Rev. D 65 (2002), 067302 doi:10.1103/PhysRevD.65.067302 [arXiv:gr-qc/0110105 [gr-qc]].
  * (23) I. Brevik, S. Nojiri, S. D. Odintsov and D. Saez-Gomez, Cardy-Verlinde formula in FRW Universe with inhomogeneous generalized fluid and dynamical entropy bounds near the future singularity, Eur. Phys. J. C 69 (2010), 563-574 doi:10.1140/epjc/s10052-010-1425-0 [arXiv:1002.1942 [hep-th]].
  * (24) I. Brevik and O. Gorbunova, Dark energy and viscous cosmology, Gen. Rel. Grav. 37 (2005), 2039-2045 doi:10.1007/s10714-005-0178-9 [arXiv:gr-qc/0504001 [gr-qc]].
  * (25) B. Li and J. D. Barrow, Does bulk viscosity create a viable unified dark matter model?, Phys. Rev. D 79 (2009), 103521 doi:10.1103/PhysRevD.79.103521 [arXiv:0902.3163 [gr-qc]].
  * (26) I. Brevik, E. Elizalde, S. D. Odintsov and A. V. Timoshkin, Inflationary universe in terms of a van der Waals viscous fluid, Int. J. Geom. Meth. Mod. Phys. 14 (2017) no.12, 1750185 doi:10.1142/S0219887817501857 [arXiv:1708.06244 [gr-qc]].
  * (27) S. Nojiri and S. D. Odintsov, Quantum de Sitter cosmology and phantom matter, Phys. Lett. B 562 (2003), 147-152 doi:10.1016/S0370-2693(03)00594-X [arXiv:hep-th/0303117 [hep-th]].
  * (28) J. D. Barrow and S. Cotsakis, Geodesics at sudden singularities, Phys. Rev. D 88 (2013), 067301 doi:10.1103/PhysRevD.88.067301 [arXiv:1307.5005 [gr-qc]].
  * (29) S. Nojiri and S. D. Odintsov, Quantum escape of sudden future singularity, Phys. Lett. B 595 (2004), 1-8 doi:10.1016/j.physletb.2004.06.060 [arXiv:hep-th/0405078 [hep-th]].
  * (30) S. Myrzakul, R. Myrzakulov and L. Sebastiani, Inhomogeneous viscous fluids in FRW universe and finite-future time singularities, Astrophys. Space Sci. 350 (2014), 845-853 doi:10.1007/s10509-014-1799-9 [arXiv:1311.6939 [gr-qc]].
  * (31) S. Capozziello, M. De Laurentis, S. Nojiri and S. D. Odintsov, Classifying and avoiding singularities in the alternative gravity dark energy models, Phys. Rev. D 79 (2009), 124007 doi:10.1103/PhysRevD.79.124007 [arXiv:0903.2753 [hep-th]].
  * (32) E. Elizalde, V. V. Obukhov and A. V. Timoshkin, Inhomogeneous viscous dark fluid coupled with dark matter in the FRW universe, Mod. Phys. Lett. A 29 (2014) no.25, 1450132 doi:10.1142/S0217732314501326 [arXiv:1406.7653 [gr-qc]].
  * (33) S. D. Odintsov and V. K. Oikonomou, Dynamical Systems Perspective of Cosmological Finite-time Singularities in $f(R)$ Gravity and Interacting Multifluid Cosmology, Phys. Rev. D 98 (2018) no.2, 024013 doi:10.1103/PhysRevD.98.024013 [arXiv:1806.07295 [gr-qc]].
  * (34) P. H. Frampton, K. J. Ludwick and R. J. Scherrer, Pseudo-rip: Cosmological models intermediate between the cosmological constant and the little rip, Phys. Rev. D 85 (2012), 083001 doi:10.1103/PhysRevD.85.083001 [arXiv:1112.2964 [astro-ph.CO]].
  * (35) P. H. Frampton, K. J. Ludwick and R. J. Scherrer, The Little Rip, Phys. Rev. D 84 (2011), 063003 doi:10.1103/PhysRevD.84.063003 [arXiv:1106.4996 [astro-ph.CO]].
  * (36) P. H. Frampton, K. J. Ludwick, S. Nojiri, S. D. Odintsov and R. J. Scherrer, Models for Little Rip dark energy, Phys. Lett. B 708 (2012), 204-211 doi:10.1016/j.physletb.2012.01.048 [arXiv:1108.0067 [hep-th]].
  * (37) S. Capozziello, V. F. Cardone, E. Elizalde, S. Nojiri and S. D. Odintsov, Observational constraints on dark energy with generalized equations of state, Phys. Rev. D 73 (2006), 043512 doi:10.1103/PhysRevD.73.043512 [arXiv:astro-ph/0508350 [astro-ph]].
  * (38) V. Astashenok, S. D. Odintsov and V. K. Oikonomou, Dark energy and cosmological horizon thermal effects, Phys. Rev. D 103 (2021), 043514 doi:10.1103/PhysRevD.103.043514.