# Nonequilibrium phase transitions and absorbing states in a model for the
dynamics of religious affiliation

Nuno Crokidakis<EMAIL_ADDRESS>Instituto de Física, Universidade
Federal Fluminense, Niterói, Rio de Janeiro, Brazil

###### Abstract

We propose a simple model to describe the dynamics of religious affiliation.
For such purpose, we built a compartmental model with three distinct
subpopulations, namely religious committed individuals, religious noncommitted
individuals and not religious affiliated individuals. The transitions among
the compartments are governed by probabilities, modeling social interactions
among the groups and also spontaneous transitions among the compartments.
First of all, we consider the model on a fully-connected network. Thus, we
write a set of ordinary differential equations to study the evolution of the
subpopulations. Our analytical and numerical results show that there is an
absorbing state in the model where only one of the subpopulations survive in
the long-time limit. There are also regions of parameters where some of the
subpopulations coexist (two or three). We also verified the occurrence of two
distinct critical points. In addition, we also present Monte Carlo simulations
of the model on two-dimensional square lattices, in order to analyze the
impact of the presence of a lattice structure on the critical behavior of the
model. Comparison of the models’ results with data for religious affiliation
in Northern Ireland shows a good qualitative agreement. Finally, we considered
the presence of inflexible individuals in the population, i.e., individuals
that never change their states. The impact of such special agents on the
critical behavior of the model is also discussed.

###### keywords:

Dynamics of social systems , Collective phenomena , Nonequilibrium phase
transitions , Monte Carlo simulation

††journal: Elsevier

##  1 Introduction

The decline of religious affiliation has been observed in some countries
recently. For example, in USA the decline in church membership is primarily a
function of the increasing number of Americans who express no religious
preference. A 2021 poll shown that over the past two decades, the percentage
of Americans who do not identify with any religion has grown from $8\%$ in
1998-2000 to $13\%$ in 2008-2010 and $21\%$ in 2019-2021 [1]. In Europe, it
has been observed an accelerated decline in religious belief since 1980 [2].

The study of the dynamics of religion evolution was attracted the attention of
researchers through the world [3, 4, 5, 6, 7, 8]. Considering religions as
population dynamics, the authors in [3] considered Weibull and lognormal
functions to fit the distribution of number of adherents of religions. The
same authors also considered a statistical physics like approach to the
problem [4]. Another study considered the growth dynamics of religious
activities. It was found that the distribution of annual logarithmic growth
rates exhibits a non-Gaussian behavior, being well described by a Laplace
distribution, for distinct size scales, and that the standard deviation of the
growth rates scales with size as a power law [5]. In order to explain the
decline of religious affiliation in some countries, some works considered
models based on coupled differential equations [6, 7]. In opposition to such
mentioned studies, a recent work considered a discrete modeling of the
propagation of religious ideas in order to maximize the number of religious
individuals and minimize the individuals who renounce religion and do not
practise religious rites, considering an optimal control strategy [8].

Compartmental models are simple and important tools concerning the modeling of
a diverse set of social dynamics, including the addoption of social habits [7,
9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. Following the ideas of references
[6, 7], we consider in this work a three-state compartmental model to study
the dynamics of religious affiliation. The model takes into account
subpopulations $X$ (religious committed individuals), $Y$ (not religious
affiliated individuals) and $Z$ (religious noncommitted individuals). Our
target is to study the critical behavior of the model. Considering the model
defined on fully-connected networks and on 2d square lattices, we found that
the system presents two distinct critical points. These critical points
separates three distinct regions, with distinct long-time behaviors. One of
the regions represents an absorbing state where the dynamics becomes frozen.
After the analysis of such model, we consider the presence of
inflexible/zealot individuals in the subpopulations of religious committed
individuals (X) and not religious ones (Y). Inflexible, or zealot, individuals
are the ones that never change their status - in this case, either X or Y
individuals will remain in their states during all the dynamics. This
realistic situation can model the presence of radical individuals concerning
religions, and the presence of inflexible individuals can dramastically change
the dynamics and the critical behavior of social models [20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31].

This work is organized as follows. In section 2 we present the initial
formulation of the model, in the absence of inflexible individuals in the
population. We divide section 3 in two subsections, one to discuss analytical
and numerical results for such initial model defined on a fully-connected
network, and other to discuss results of Monte Carlo simulations of the model
defined on the top of two-dimensional square lattices. In section 4, we
consider the presence of zealot individuals in the population, and study the
impact of such presence in the critical behavior of the initial model in
section 5. Finally, in section 6 we present our final remarks.

##  2 Model I: Three-state model with no inflexible agents

Based on the works [6, 7], we consider a population of $N$ individuals or
agents, where each agent $i$ ($i=1,2,...,N$) can be in one of three possible
states, namely the subpopulations $X$ (religious committed individuals), $Y$
(not religious affiliated individuals) and $Z$ (religious noncommitted
individuals). It can be argued that the consideration of such religious
noncommitted group ($Z$) instead of only two groups ($X$ and $Y$) is a
realistic feature for built models to study dynamics of religious
affiliations. Indeed, the majority of Europe’s Christians are non-practicing,
but they differ from religiously unaffiliated people in their views on God,
attitudes toward Muslims and immigrants, and opinions about religion’s role in
society [32].

Since we are considering a social contagion model, the probabilities related
to changes in agents’ compartments represent the possible contagions. The
transitions among compartments are as following:

  * 1.

$X\stackrel{{\scriptstyle p}}{{\rightarrow}}Z$: a religious committed agent
(X) becomes a religious noncommitted individual (Z) with probability $p$;

  * 2.

$Z\stackrel{{\scriptstyle q}}{{\rightarrow}}X$: a religious noncommitted agent
(Z) becomes a religious committed individual (X) with probability $q$ if
he/she is in contact with religious committed individuals (X);

  * 3.

$Z\stackrel{{\scriptstyle u}}{{\rightarrow}}Y$: a religious noncommitted agent
(Z) becomes a not religious affiliated individual (Y) with probability $u$ if
he/she is in contact with not religious affiliated individuals (Y);

  * 4.

$Y\stackrel{{\scriptstyle w}}{{\rightarrow}}X$: a not religious affiliated
agent (Y) becomes a religious committed individual (X) with probability $w$ if
he/she is in contact with religious committed individuals (X);

In the above rules, $p$ represents a spontaneous transitions, for which
religious committed agents can lose interest in in attending churches or
temples, but if they are questioned, they affirm they belong to a given
religion (any religion). Thus, religious committed agents (X) can become
religious noncommitted individuals (Z). The probability $q$ represents an
“infection” probability, i.e., the probability that a religious noncommitted
agent (Z) becomes a religious committed individual (X). We argue that the
$\textbf{Z}\to\textbf{X}$ transition is more difficult to occur in comparison
with $\textbf{X}\to\textbf{Z}$ transition, and for such reason we considered
that the $\textbf{Z}\to\textbf{X}$ transition is due to to social pressure of
religious committed agents (X) over religious noncommitted individuals (Z). In
addition, religious noncommitted individuals (Z) can be influenced by not
religious affiliated individuals (Y) to become not religious affiliated
individuals (Y). This transition is ruled by a social pressure probability
$u$. Finally, not religious affiliated individuals (Y) can become religious
committed individuals (X) due to social pressure of religious committed
individuals (X). This pressure is modeled by a probability $w$.

In the next section we will discuss our results.

##  3 Results: Model I

###  3.1 Fully-connected population

Considering the above-mentioned rules and a population defined on a fully-
connected network, we can write the following mean-field equations for the
compartments’ evolution:

$\displaystyle\frac{dx(t)}{dt}$ $\displaystyle=$
$\displaystyle-p\,x(t)+q\,x(t)\,z(t)+w\,x(t)\,y(t)$ (1)
$\displaystyle\frac{dy(t)}{dt}$ $\displaystyle=$
$\displaystyle-w\,x(t)\,y(t)+u\,y(t)\,z(t)$ (2)
$\displaystyle\frac{dz(t)}{dt}$ $\displaystyle=$
$\displaystyle-u\,y(t)\,z(t)-q\,x(t)\,z(t)+p\,x(t)$ (3)

For simplicity, we consider a fixed population, i.e., at each time step $t$ we
have the normalization condition $x(t)+y(t)+z(t)=1$, where we defined the
population densities $x(t)=X(t)/N$, $y(t)=Y(t)/N$ and $z(t)=Z(t)/N$.

Let us start analyzing the time evolution of the previous defined quantities
$x,y$ and $z$. We numerically integrated the Eqs. (1) - (3), and the results
are exhibited in Figure 1. For simplicity we fixed $w=0.05$, $q=0.12$ and
$u=0.10$, and vary the parameter $p$. One can see that the increase of $p$
causes the increase of $y$ and the decrease of $x$ and $z$. It is also
interesting to observe that, for sufficient high values of $p$ like $p=0.08$
(panel (c)) and $p=0.15$ (panel (d)), the densities evolve in time, and in the
steady states we observe the survival of only the not religious affiliated
subpopulation $y$, i.e., for $t\to\infty$ we have $x=z=0$ and $y=1$. This last
result will be discussed in more details analytically in the following.

Figure 1: (Color online) Time evolution of the three densities of agents
$x(t)$, $y(t)$ and $z(t)$ for the mean-field formulation of the model, based
on the numerical integration of Eqs. (1) - (3). The fixed parameters are
$w=0.05$, $q=0.12$ and $u=0.10$, and we varied the parameter $p$: (a)
$p=0.030$, (b) $p=0.045$, (c) $p=0.080$, (d) $p=0.150$.

As we observed in Figure 1, the quantities $x(t),y(t)$ and $z(t)$ evolve in
time, and after some time they stabilize. In such stationary states, the time
derivatives of Eqs. (1) - (3) are zero. Denoting the stationary densities as
$x=x(t\to\infty)$, $y=y(t\to\infty)$ and $z=z(t\to\infty)$, one can obtain
three distinct solutions for the stationary states, namely (see Appendix),

$\displaystyle x=\begin{cases}1-\frac{p}{q}&if\ p<p_{c}^{(1)}\\\
\frac{u\,(w-p)}{w\,(u+w-q)}&if\ p_{c}^{(1)}<p<p_{c}^{(2)}\\\ 0&if\
p>p_{c}^{(2)}\end{cases}$ (4)

$\displaystyle y=\begin{cases}0&if\ p<p_{c}^{(1)}\\\
1-\frac{(u+w)\,(w-p)}{w\,(u+w-q)}&if\ p_{c}^{(1)}<p<p_{c}^{(2)}\\\ 1&if\
p>p_{c}^{(2)}\end{cases}$ (5)

$\displaystyle z=\begin{cases}\frac{p}{q}&if\ p<p_{c}^{(1)}\\\
\frac{w-p}{u+w-q}&if\ p_{c}^{(1)}<p<p_{c}^{(2)}\\\ 0&if\
p>p_{c}^{(2)}\end{cases}$ (6)

where the critical points are given by 111The stability analysis of the three
solutions (4) - (6) is discussed in the Appendix.

$\displaystyle p_{c}^{(1)}$ $\displaystyle=$ $\displaystyle\frac{w\,q}{u+w}$
(7) $\displaystyle p_{c}^{(2)}$ $\displaystyle=$ $\displaystyle w$ (8)

Taking $y$ as the order parameter of the model (see the Appendix), the
solution for $p>p_{c}^{(2)}$, i.e., $(x,y,z)=(0,1,0)$, defines an absorbing
state where there are only not religious affiliated individuals in the
population for long times. Looking for Fig. 1 again, this absorbing state can
be seen in panels (c) and (d) for long times. For the values of the considered
parameters, namely $w=0.05$, $q=0.12$ and $u=0.10$, we obtain from Eqs. (7)
and (8) the values $p_{c}^{(1)}=0.04$ and $p_{c}^{(2)}=0.05$. Thus for the
considered values ($p=0.08$ in panel (c) and $p=0.15$ in panel (d)), the
system evolves to steady states where $x=z=0$ and $y=1$, as predicted
analytically in Eqs. (4) - (6).

Figure 2: (Color online) Stationary densities $x$, $y$ and $z$ as functions of
the probability $p$ for $w=0.05$, $q=0.12$ and $u=0.10$ for the mean-field
formulation of the model. The lines were obtained from Eqs. (4) - (6). For the
considered parameters, the critical points are $p_{c}^{(1)}=0.04$ and
$p_{c}^{(2)}=0.05$. The absorbing state where $x=z=0$ and $y=1$ is observed
for $p>p_{c}^{(2)}=0.05$. It is important to note that these behaviors are
present also for other parameter values, and what is shown here works as a
pattern.

Based on the results above, we plot in Fig. 2 the steady state values of the
three subpopulation densities $x$, $y$ and $z$ as functions of $p$. For this
graphic, again we fixed the parameters $w=0.05$, $q=0.12$ and $u=0.10$. For
such values, we have $p_{c}^{(1)}=0.04$ and $p_{c}^{(2)}=0.05$. As discussed
above, for $p>p_{c}^{(2)}$ the system becomes frozen in the absorbing phase,
where there are only not religious affiliated agents $y$ in the population
after a long time, i.e., $y=1$ and $x=z=0$. For $p<p_{c}^{(1)}$ we have a
region where the populations of religious individuals coexist, i.e., we have
the religious committed and religious noncommitted individuals coexisting in
the stationary states, whereas the population of not religious agents
disappears. The realistic region where the three populations $x,y$ and $z$
coexist in the steady states is obtained for $p_{c}^{(1)}<p<p_{c}^{(2)}$.

Figure 3: (Color online) Phase diagram of the mean-field model in the plane
$p$ vs $w$ for $q=0.12$ and $u=0.10$. The full line represents the critical
point $p_{c}^{(1)}$ given by Eq. (7) and the dotted line the other critical
point $p_{c}^{(2)}$ given by Eq. (8). I denotes the region where the the
system falls in the absorbing phase $y=1,x=z=0$, II denotes the region where
the three densities $x$, $y$ and $z$ coexist and III denotes the region where
the subpopulation $y$ disappears and we observe the coexistence of
subpopulations $x$ and $z$.

The competition among the contagions cause the occurrence of such three
regions in the model [20]. From one side we have religious committed
individuals influencing religious noncommitted ones to increase their beliefs
about the religion, with probability $q$. On the other hand, we have the
social pressure of not religious affiliated individuals over religious
noncommitted ones, with probability $u$, in order to make such religious
noncommitted agents to became not religious affiliated. We also have another
social pressure of the religious committed individuals over not religious
affiliated ones, with probability $w$, that intend to convert such not
religiuos affiliated agents into religious commited ones. Finally, it is
important to mention the parameter $p$, that drives the only transition of the
model that does not depend on a direct social interaction. That parameter
models the spontaneous decrease of interest in the religion of the
individuals. Notice also that the parameter $w$ is present in both expressions
for the nonequilibrium critical points $p_{c}^{(1)}$ and $p_{c}^{(2)}$,
suggesting the importance of such parameter for the dynamics of the model.

For clarity, we exhibit in Fig. 3 the phase diagram of the model in the plane
$p$ versus $w$, separating the three above discussed regions. In Fig. 3, the
absorbing phase with $y=1$ and $x=z=0$ is located in region I for
$p<p_{c}^{(1)}$, the coexistence phase is denoted by II for
$p_{c}^{(1)}<p<p_{c}^{(2)}$ (where the three densities $x,y$ and $z$ coexist)
and region III denotes the region where the subpopulation $y$ disappears and
we observe the coexistence of subpopulations $x$ and $z$.

Figure 4: (Color online) Comparison between data for the fraction of Northern
Irish population who are deemed to be religiously committed (group X, black
squares) over time and the model (full red line). The line is the time
evolution of $x(t)$, that was obtained from the numerical integration of Eqs.
(1) - (3). Parameters are $w=0.047$, $q=0.078$, $u=0.32$ and $p=0.043$. Data
was obtained from reference [7].

Finally, in order to verify if this simple three-compartment model is able to
qualitatively reproduce real data, we compare the model’s results with data
for religious affiliation in Northern Ireland, available in ref. [7]. In such
reference [7] we found data of the total number of communicant members of the
Presbyterian Church as recorded in church reports, which results in a density
of religious committed individuals. Data were collected from 1951 to 2011 in
time steps of 10 years, thus in Fig. 4 the initial time $t=0$ represents the
fraction of religious committed individuals for 1951, $t=1$ represents the
fraction of religious committed individuals for 1961, and so on. Since the
data is only for the fraction of religious committed individuals, we plot the
density $x(t)$ obtained from the model, together with the data. The curve for
$x(t)$ was obtained from the numerical integration of Eqs. (1) - (3). For the
numerical results, we considered the parameters are $w=0.047$, $q=0.078$,
$u=0.32$ and $p=0.043$, obtained from a least-square fit. One can see in Fig.
4 that the numerical results qualitatively correspond to the data. It is
important to notice that, for such parameters, the critical points are
$p_{c}^{(1)}=0.047$ and $p_{c}^{(2)}\approx 0.0099$. Thus, the fitted
parameters represents the population in the region where the three
subpopulations $X,Y$ and $Z$ coexist ($p_{c}^{(1)}<p<p_{c}^{(2)}$), which is a
realistic region of the model (region II of fig. 3).

###  3.2 Two-dimensional square lattice

In this subsection we present some numerical results of Monte Carlo
simulations of the model on the top of 2d square lattices. Thus, we built an
agent-based formulation of the compartmental model proposed in the last
subsection. For this purpose, the algorithm to simulate the model is defined
as follows:

Figure 5: (Color online) Time evolution of the three densities of agents
$x(t)$, $y(t)$ and $z(t)$ for the model defined on a 2d square lattice with
linear size $L=100$. The fixed parameters are $w=0.05$, $q=0.12$ and $u=0.10$,
and we varied the parameter $p$: (a) $p=0.030$, (b) $p=0.038$, (c) $p=0.060$,
(d) $p=0.120$. Results are averaged over $100$ independent simulations.

  * 1.

we generate a $L$ x $L$ square lattice with a population size $N=L^{2}$ and
periodic boundary conditions;

  * 2.

given an initial condition $X(0),Y(0)$ and $Z(0)$, we randomly distribute
these agents in the lattice sites;

  * 3.

at each time step, every lattice site is visited in a sequential order;

  * 4.

if a given agent $i$ is in $X$ state, we generate a random number $r$ in the
range $[0,1]$. If $r<p$, the agent $i$ changes to state $Z$;

  * 5.

if a given agent $i$ is in $Z$ state, we choose at random one of his/her
nearest neighbors, say $j$. If such neighbor $j$ is in $X$ state, we generate
a random number $r$ in the range $[0,1]$. If $r<q$, the agent $i$ changes to
state $X$. On the other hand, if the neighbor $j$ is in $Y$ state, we generate
a random number $r$ in the range $[0,1]$. If $r<u$, the agent $i$ changes to
state $Y$;

  * 6.

if a given agent $i$ is in $Y$ state, we choose at random one of his/her
nearest neighbors, say $k$. If such neighbor $k$ is in $X$ state, we generate
a random number $r$ in the range $[0,1]$. If $r<w$, the agent $i$ changes to
state $X$.

One time step is defined by the visit of all lattice sites. The agents’ states
were updated synchronously, i.e., we considered parallel updating, as it is
standard in probabilistic cellular automata in order to avoid correlations
between consecutive steps [34, 35]. In addition, all results are averaged over
$100$ independent simulations.

Considering the same parameters of the previous subsection, namely $w=0.05$,
$q=0.12$ and $u=0.10$, and varying the parameter $p$, we exhibit in Fig. 5 the
time evolution of the system on a square lattice of linear size $L=100$. We
can see a similar behavior observed in the fully-connected case, i.e., we can
have the coexistence of $x$ and $z$ populations in the steady states (panel
(a)), the coexistence of the three subpopulations (panel (b)) and absorbing
states with $y=1$, $x=z=0$ (panels (c) and (d)). However, the time to reach
the steady states is longer than in the mean-field case, as it is usual in
agent-based models defined on lattices [36].

We can also compute the stationary densities $x,y$ and $z$ to compare with the
mean-field case. In Fig. 6 we exhibit the stationary values of the three
subpopulations as function of the probability $p$ for $w=0.05$, $q=0.12$ and
$u=0.10$, considering a square lattice with linear size $L=100$. We observe
similar behaviors for such quantities in comparison with the mean-field
calculations presented in section 3.1. However, as it is usual in lattice
models, we observe a distinct location of the critical points $p_{c}^{(1)}$
and $p_{c}^{(2)}$ in comparison with the mean-field predictions [33]. For the
same set of parameters, looking for the data for the 2d case we have
$p_{c}^{(1)}\approx 0.0364$ and $p_{c}^{(2)}\approx 0.0402$, whereas in the
previous subsection we observed $p_{c}^{(1)}=0.04$ and $p_{c}^{(2)}=0.05$. The
differences can be pointed to the correlations generated by the lattice
structure, that are absent in the fully-connected case.

Figure 6: (Color online) Stationary densities $x$, $y$ and $z$ as functions of
the probability $p$ for $w=0.05$, $q=0.12$ and $u=0.10$ for the model defined
on a 2d square lattice with linear size $L=100$. The absorbing state where
$x=z=0$ and $y=1$ is observed for $p>\approx 0.04$. Results are averaged over
$100$ independent simulations.

For better visualization of the lattice and the subpopulations, we plot in
Fig. 7 some snapshots of the stationary states of the model. We considered the
same lattice size used before, namely $L=100$, the same fixed parameters
$w=0.05$, $q=0.12$ and $u=0.10$, and we plot four distinct values of $p$. The
distinct colors represent the subpopulations $x$ (black), $y$ (red) and $z$
(blue). For low values of $p$ like $p=0.010$, we observe the dominance of the
$x$ population (upper panel, left side), and no $y$ population. If we increase
$p$ for $p=0.030$ we observe an increase of the $z$ population and the
decrease of the $x$ population, and the $y$ population remains absent (upper
panel, right side). For $p=0.038$ we are in the region where the three
subpopulations coexist (lower panel, left side). Finally, for sufficient high
values of $p$ like $p=0.050$ we observe that only the $y$ population survives
in the long-time limit (lower panel, right side), which represents the
previous mentioned absorbing state of the model.

Figure 7: (Color online) Snapshots of the stationary states of the model on a
square lattice of linear size $L=100$. The distinct colors represent the
subpopulations $x$ (black), $y$ (red) and $z$ (blue). The fixed parameters are
$w=0.05$, $q=0.12$ and $u=0.10$, and we varied the parameter $p$: (a)
$p=0.010$, (b) $p=0.030$, (c) $p=0.038$, (d) $p=0.050$.

##  4 Model II: Presence of inflexible individuals

In this section we consider the presence of inflexible individuals in the
population. Such individuals will be included in the subpopulations of
religious committed individuals (X) and not religious ones (Y). Inflexible, or
zealot, individuals are the ones that never change their status - in this
case, either X or Y individuals will remain in their states during all the
dynamics. This is analogous to introducing quenched disorder in the system,
considerably altering its critical behavior.

In this case, as considered recently for the dynamics of alcohol consumption
[20], we consider that we have a certain number $X_{I}$ of inflexible agents
related to religious commited individuals, i.e., individuals that are
committed with a given religion and will always devoted to such reiligion,
i.e., they will keep this behavior during all their lives. In addition, we
also considered a certain number $Y_{I}$ of inflexible agents related to not
religious individuals, i.e., individuals that do not follow any religion and
do not intend to do it. Thus, the total number of religious committed
individuals are given by $X_{total}=X+X_{I}$ and the total number of not
religious individuals are given by $Y_{total}=Y+Y_{I}$, where $X$ and $Y$
denote the noninflexible individuals of each class. Since we are not
considering inflexible individuals for the religious noncommitted individuals
$Z$, we have simply $Z_{total}=Z$.

For this model, the transitions among compartments are the following:

  * 1.

$X\stackrel{{\scriptstyle p}}{{\rightarrow}}Z$: a noninflexible religious
committed agent (X) becomes a religious noncommitted individual (Z) with
probability $p$;

  * 2.

$Z\stackrel{{\scriptstyle q}}{{\rightarrow}}X$: a religious noncommitted agent
(Z) becomes a religious committed individual (X) with probability $q$ if
he/she is in contact with religious committed individuals (X and $X_{I}$);

  * 3.

$Z\stackrel{{\scriptstyle u}}{{\rightarrow}}Y$: a religious noncommitted agent
(Z) becomes a not religious affiliated individual (Y) with probability $u$ if
he/she is in contact with not religious affiliated individuals (Y and
$Y_{I}$);

  * 4.

$Y\stackrel{{\scriptstyle w}}{{\rightarrow}}X$: a noninflexible not religious
affiliated agent (Y) becomes a religious committed individual (X) with
probability $w$ if he/she is in contact with religious committed individuals
(X and $X_{I}$);

For such a case, the normalization condition becomes

$x(t)+y(t)+z(t)+x_{I}+y_{I}=1~{},$ (9)

where $x_{I}=X_{I}/N$ and $y_{I}=Y_{I}/N$ are the inflexible/zealot
subpopulations in the compartments $X$ and $Y$, respectively, that are
constant in time.

In the next section we discuss the impact of such zealots in the critical
behavior of the model.

##  5 Results: Model II

###  5.1 Fully-connected population

Considering the above-mentioned rules for the model in the presence of
inflexible individuals and a population defined on a fully-connected network,
we can write the following mean-field equations for the compartments’
evolution:

$\displaystyle\frac{dx(t)}{dt}$ $\displaystyle=$
$\displaystyle-p\,x(t)+q\,[x(t)+x_{I}]\,z(t)+w\,[x(t)+x_{I}]\,y(t)$ (10)
$\displaystyle\frac{dy(t)}{dt}$ $\displaystyle=$
$\displaystyle-w\,[x(t)+x_{I}]\,y(t)+u\,[y(t)+y_{I}]\,z(t)$ (11)
$\displaystyle\frac{dz(t)}{dt}$ $\displaystyle=$
$\displaystyle-u\,[y(t)+y_{I}]\,z(t)-q\,[x(t)+x_{I}]\,z(t)+p\,x(t)$ (12)

Figure 8: (Color online) Stationary densities of the subpopulations as
functions of the probability $p$ for $w=0.05$, $q=0.12$ and $u=0.10$, for the
mean-field formulation of the model in the presence of inflexible individuals.
The lines were obtained from the numerical integration of Eqs. (10) - (12).
The densities of inflexible/zealots are: (a) $x_{I}=0.1,y_{I}=0.0$, (b)
$x_{I}=0.0,y_{I}=0.1$ and (c) $x_{I}=0.1,y_{I}=0.1$.

Since we are interested in the critical behavior of the model, we present in
Fig. 8 only the stationary densities of the subpopulations as functions of
$p$. The fixed parameters are the same used in the remaining of the work,
i.e., $w=0.05$, $q=0.12$ and $u=0.10$. To analyze the impact of the distinct
classes of inflexible individuals $x_{I}$ and $y_{I}$, we plot in Fig. 8 (a)
the stationary values only in the presence of inflexible agents in the class
of religious committed individuals. In this case, we exhibit results for
$x_{I}=0.1$ and $y_{I}=0.0$. We can see the phase where $y=0$ in the steady
states persists for low values of $p$, and we observe there is a critical
point associated with the transition from such phase where $y=0$ to $y\neq 0$.
The second phase transition (active-absorbing transition) is eliminated in the
presence of inflexible agents in the class of religious committed individuals,
i.e., we can see that the subpopulations of religious committed individuals
and not religious ones will not disappear of the population for long times,
for all values of $p$ (see the inset of Fig. 8 (a)). In other words, the
absorbing state is suppressed in the presence of $X_{I}$ inflexible
individuals.

In Fig. 8 (b) we exhibit the stationary densities only in the presence of
inflexible agents in the class of not religious affiliated individuals. In
this case, we exhibit results for $x_{I}=0.0$ and $y_{I}=0.1$. We can see the
collective state where the class of not religious affiliated individuals
disappear of the population is eliminated. However, we observe the persistence
of the active-absorbing phase transition. However, such phase is composed now
by $Y$ and $Y_{I}$ individuals, and the remaining subpopulations $X$ and $Z$
disappear.

Finally, the join effect of the inflexible individuals in the population are
shown in Fig. 8 (c), where we show results for $x_{I}=0.1$ and $y_{I}=0.1$. We
see that both phase transitions are now suppressed, and we have the
coexistence of all the subpopulations $X,X_{I},Y,Y_{I}$ and $Z$ for all values
of $p>0$. Since the coexistence of the distinct classes of individuals is a
realistic scenario, we can argue that the presence of inflexible individuals
in both classes (religious committed agents and not religious affiliated ones)
makes the model more realistic.

###  5.2 Two-dimensional square lattice

In this subsection we present some numerical results of Monte Carlo
simulations of the model with inflexible agents on the top of 2d square
lattices. Thus, we built an agent-based formulation of the compartmental model
proposed in the last subsection. We will keep the nomenclature of the last
subsection, i.e., the noninflexible population in the three classes is given
by $X,Y$ and $Z$, and the inflexible population is given by $X_{I}$ and
$Y_{I}$. Thus, the algorithm to simulate the model is defined as follows:

  * 1.

we generate a $L$ x $L$ square lattice with a population size $N=L^{2}$ and
periodic boundary conditions;

  * 2.

we choose at random $X_{I}$ ($Y_{I}$) lattice sites and freeze their opinions
in the religious committed individuals (not religious affiliated individuals),
in such a way that such individuals will never change opinion during all the
dynamics;

  * 3.

the remaining of individuals in the population are distributed in the
noninflexbile states $X,Y$ and $Z$ with a given initial condition $X(0),Y(0)$
and $Z(0)$, and the are randomly distributed over the lattice;

  * 4.

at each time step, every lattice site is visited in a sequential order;

  * 5.

if a given agent $i$ is in $X$ state, we generate a random number $r$ in the
range $[0,1]$. If $r<p$, the agent $i$ changes to state $Z$;

  * 6.

if a given agent $i$ is in $Z$ state, we choose at random one of his/her
nearest neighbors, say $j$. If such neighbor $j$ is in $X$ or $X_{I}$ states,
we generate a random number $r$ in the range $[0,1]$. If $r<q$, the agent $i$
changes to state $X$. On the other hand, if the neighbor $j$ is in $Y$ or
$Y_{I}$ states, we generate a random number $r$ in the range $[0,1]$. If
$r<u$, the agent $i$ changes to state $Y$;

  * 7.

if a given agent $i$ is in $Y$ state, we choose at random one of his/her
nearest neighbors, say $k$. If such neighbor $k$ is in $X$ or $X_{I}$ states,
we generate a random number $r$ in the range $[0,1]$. If $r<w$, the agent $i$
changes to state $X$.

As in section 3.2, one time step is defined by the visit of all lattice sites,
and the agents’ states were updated syncronously. In addition, all results are
averaged over $100$ independent simulations.

Considering the same parameters of the previous subsection, namely $w=0.05$,
$q=0.12$ and $u=0.10$, and varying the parameter $p$, we exhibit in Fig. 9 the
stationary densities of the subpopulations for: $x_{I}=0.1$ and $y_{I}=0.0$
(panel (a)), $x_{I}=0.0$ and $y_{I}=0.1$ (panel (b)) and $x_{I}=0.1$ and
$y_{I}=0.1$ (panel (c)). The behavior is similar to the observed in the mean-
field case, but again the critical points, when they exist, are located in
smaller values of $p$ in comparison with the mean-field case. The suppression
of the absorbing phase is also observed for the model in the presence of only
$X_{I}$ inflexible agents. In addition, the coexistence of two subpopulations
does not occur anymore, as in the model defined in a fully-connected
population. Finally, the model in the presence of both $X_{I}$ and $Y_{I}$
inflexible individuals shows the coexistence of all subpopulations for all
values of $p$, denoting a more realistic situation.

Figure 9: (Color online) Stationary densities of the subpopulations as
functions of the probability $p$ for $w=0.05$, $q=0.12$ and $u=0.10$ for the
model in the presence of inflexible individuals in a square lattice of linear
size $L=100$. The densities of inflexible/zealots are: (a)
$x_{I}=0.1,y_{I}=0.0$, (b) $x_{I}=0.0,y_{I}=0.1$ and (c)
$x_{I}=0.1,y_{I}=0.1$. Results are averaged over $100$ independent
simulations.

##  6 Final Remarks

In this work we study a simple model for the dynamics of religious
affiliation. Fot this purpose, we built a three-state compartmental model, and
the transitions among such compartments are ruled by probabilities. The three
states are: religious committed individuals (X), religious noncommitted
individuals (Z) and not religious affiliated individuals (Y). The transitions
among the states are ruled by probabilities.

First of all, we study the model on a fully-connected population. Thus, we
write the mean-field equations for the evolutions of the three subpopulations.
Some analytical results were obtained for the stationary states of the model.
Such solutions, together with the numerical integration of the model’s
equations, reveal that there are three distinct collective states or phases in
the stationary states of the model, depending on the range of parameters: (I)
a phase where only the Y subpopulation survives; (II) a region where the
subpopulations X and Z survive and the Y population disappers; (III) a region
where the three subpopulations coexist. Region (I) represents a region of
theoretical interest for nonequilibrium Statistical Physics, since it
represents an absorbing state, where the system is frozen and no transitions
among compartments occur anymore. Region (II) represents a realistic
representation of a society, where there are religious individuals (some of
them committed with the religion, and some others no committed) and not
religious ones. The transition from region II to region I represents an
active-absorbing nonequilibrium phase transition that appears to be in the
directed percolation universality class. The model is also capable to
qualitatively reproduce data for religious affiliation in Northern Ireland. In
addition, we also presented in this work results of agent-based Monte Carlo
simulations of the model on two-dimensional square lattices. The results show
a similar behavior in comparison with the mean-field case, but the critical
points are localized for lower values of the control parameter of the phase
transition. In addition, the times to the system achieves stationary states
are higher in comparison with the mean-field case.

After such initial analysis, we considered the presence of inflexible
individuals in the population. For this purpose, we included inflexible/zealot
individuals in the subpopulations of religious committed individuals (fixed
density $x_{I}$) and not religious ones (fixed density $y_{I}$). Inflexible,
or zealot, individuals are the ones that never change their status - in this
case, either X or Y individuals will remain in their states during all the
dynamics. The densities of such inflexible agents are fixed in time,
representing quenched disorder in the language of magnetic systems [38]. For
the practical point of view it is a realistic situation that can model the
presence of radical individuals concerning religions, and the presence of
inflexible individuals can dramastically change the dynamics and the critical
behavior of social models [39]. We verified that the presence of inflexible
agents only in the religious committed individuals can destroy the absorbing
phase where only $Y$ individuals survive in the long run. On the other hand,
in the presence of zealots only in the not religious affiliated subpopulation
we observe the persistence of such absorbing state, but the phase where the
$Y$ population disappears is not present anymore. Finally, in the presence of
inflexible agents in both religious committed individuals and not religious
affiliated ones the only collective macroscopic state is the coexistence of
all the subpopulations, i.e., we have no phase transitions anymore.

## Acknowledgments

The author acknowledges financial support from the Brazilian scientific
funding agencies Conselho Nacional de Desenvolvimento Científico e Tecnológico
(CNPq, Grant 308643/2023-2) and Fundação de Amparo à Pesquisa do Estado do Rio
de Janeiro (FAPERJ, Grant 203.217/2017).

##  Appendix A Analytical considerations

In this appendix we will detail some of the analytical calculations
considering the mean-field formulation of the model (fully-connected
population, section $3.1$).

Let us start considering the long-time limit, $t\to\infty$, in Eq. (1). Taking
$dx/dt=0$, we obtain two solutions, namely

$\displaystyle x$ $\displaystyle=$ $\displaystyle 0$ (13) $\displaystyle z$
$\displaystyle=$ $\displaystyle\frac{p-w\,y}{q}$ (14)

For Eq. (2), taking $dy/dt=0$ we also obtain two solutions,

$\displaystyle y$ $\displaystyle=$ $\displaystyle 0$ (15) $\displaystyle x$
$\displaystyle=$ $\displaystyle\frac{u\,z}{w}$ (16)

For Eq. (3), we obtain the following relation,

$(p-q\,z+u\,z)\,x=u\,z\,(1-z)$ (17)

If Eq. (16) is valid, Eq. (17) gives us two other solutions,

$\displaystyle z$ $\displaystyle=$ $\displaystyle 0$ (18) $\displaystyle z$
$\displaystyle=$ $\displaystyle\frac{w-p}{u+w-q}$ (19)

Observe that, for the region of parameters where Eq. (18) is valid ($z=0$),
Eq. (16) gives us $x=0$, leading to $y=1$ from the normalization condition.
This solution $(x,y,z)=(0,1,0)$ represents the absorbing state of the model,
where there are only not religious affiliated individuals in the population
for long times, as discussed in section $3.1$.

We considered $z=0$ as a valid solution to obtain the absorbing state. If now
we consider Eq. (19) valid, Eq. (16) gives us another solution for $x$, and
after we can consider the the normalization condition to obtain $y$. These
solutions are given by

$\displaystyle x$ $\displaystyle=$ $\displaystyle\frac{u\,(w-p)}{w\,(u+w-q)}$
(20) $\displaystyle y$ $\displaystyle=$ $\displaystyle
1-\frac{(u+w)\,(w-p)}{w\,(u+w-q)}$ (21)

Now we can consider the solution given by Eq. (15). Using this result in Eq.
(14), we obtain $z=p/q$. Since this solution is valid for $y=0$, the
normalization condition gives us $x=1-p/q$. These results defines another
solution for the model, namely $(x,y,z)=(1-p/q,0,p/q)$.

From the previous three solutions, and considering $y$ as the order parameter
of the model, we see that we have the solutions $y=0$, $y=1$ and $y$ given by
Eq. (21). From Eq. (21) we can obtain the two critical points of the model.
This solution can be rewriten as $y\sim(p-p_{c}^{(1)})^{\beta}$, that
approaches zero for

$p_{c}^{(1)}=\frac{w\,q}{u+w}$ (22)

and present a critical exponent $\beta=1$, as in the directed percolation
problem [33].

The other critical point can be obtained looking again for Eq. (21). This
solution achieves the absorbing state where $y=1$ for $p=w$, defining the
second critical point,

$p_{c}^{(2)}=w$ (23)

The above equations define the three distinct solutions discussed in section
$3.1$, as well the two critical points of the model.

The local stability of the above equilibrium points can be inferred from the
eigenvalues of the Jacobian matrix $J$ obtained from the system of
differential equations (1) - (3). An equilibrium point is locally
asymptotically stable if all eigenvalues of $J$ have negative real parts [37].
The eigenvalues can be obtained from $det(J-\lambda\,I)=0$, where $I$ is the
identity matrix. The Jacobian matrix is given by taking the partial
derivatives of Eqs. (1) - (3), namely

$J=\begin{bmatrix}\frac{\partial{\dot{x}}}{\partial{x}}&\frac{\partial{\dot{x}}}{\partial{y}}&\frac{\partial{\dot{x}}}{\partial{z}}\\\
\frac{\partial{\dot{y}}}{\partial{x}}&\frac{\partial{\dot{y}}}{\partial{y}}&\frac{\partial{\dot{y}}}{\partial{z}}\\\
\frac{\partial{\dot{z}}}{\partial{x}}&\frac{\partial{\dot{z}}}{\partial{y}}&\frac{\partial{\dot{z}}}{\partial{z}}\end{bmatrix}$

Taking the derivatives, the matrix becomes

$J=\begin{bmatrix}-p+w\,y+q\,z&w\,x&q\,x\\\ -w\,y&-w\,x+u\,z&u\,y\\\
p-q\,z&-u\,z&-u\,y-q\,x\end{bmatrix}$

For the absorbing state solution $(x,y,z)=(0,1,0)$, the eigenvelues of the
Jacobian matrix are $\lambda_{1}=0,\lambda_{2}=-u$ and $\lambda_{3}=w-p$. We
have $\lambda_{2}<0$ since $u>0$, and $\lambda_{3}<0$ if we have $p>w$. This
is exactly the condition we obtained above, i.e., the absorbing state solution
is valid for $p>p_{c}^{(2)}=w$.

Another equilibrium is given by $(x,y,z)=(1-p/q,0,p/q)$. The eigenvelues of
the Jacobian matrix are $\lambda_{1}=0,\lambda_{2}=p-q$ and
$\lambda_{3}=(u+w)(p/q)-w$. We have $\lambda_{2}<0$ for $p<q$ which leads to
$x=1-p/q<1$ and $z=p/q<1$, which is expected since $x$ and $z$ are
subpopulation densities. We have $\lambda_{3}<0$ for $p<wq/(u+w)$. This is the
condition we obtained above, i.e., the solution $(x,y,z)=(1-p/q,0,p/q)$ is
valid for $p<p_{c}^{(1)}=wq/(u+w)$. The third solution (coexistence of the
three subpopulations $x,y$ and $z$) is then stable for
$p_{c}^{(1)}<p<p_{c}^{(2)}$.

## References

  * [1] U.S. Church Membership Falls Below Majority for First Time, $https://news.gallup.com/poll/341963/church-membership-falls-below-majority-first-time.aspx$
  * [2] M. Dogan, Accelerated Decline of Religious Beliefs in Europe, Comparative Sociology 1(2), 127-149 (2002).
  * [3] M. Ausloos, F. Petroni, Statistical dynamics of religions and adherents, Europhys. Lett. 77, 38002 (2007).
  * [4] M. Ausloos, F. Petroni, Statistical dynamics of religion evolutions, Physica A 388, 4438 (2009).
  * [5] S. Picoli, R. S. Mendes, Universal features in the growth dynamics of religious activities, Phys. Rev. E 77, 036105 (2008).
  * [6] D. M. Abrams, H. A. Yaple, Dynamics of social group competition: Modeling the decline of religious affiliation, Phys. Rev. Lett. 107, 088701 (2011).
  * [7] M. McCartney, D. H. Glass, A three-state dynamical model for religious affiliation, Physica A 419, 145 (2015).
  * [8] A. Labzai, B. Khajji, O. Balatif, M. Rachik, Discrete mathematical modeling of the propagation of religious ideas: Optimal control approach,Commun. Math. Biol. Neurosci. 2020:21 (2020).
  * [9] A. L. Oestereich, M. A. Pires, N. Crokidakis, D. O. Cajueiro, Optimal rewiring in adaptive networks in multi-coupled vaccination, epidemic and opinion dynamics, Chaos, Solitions $\&$ Fractals 176, 114125 (2023).
  * [10] N. Crokidakis, Modeling the impact of civilian firearm ownership in the evolution of violent crimes, App. Math. Comput. 429, 127256 (2022).
  * [11] N. Crokidakis, A simple mechanism leading to first-order phase transitions in a model of tax evasion, Int. J. Mod. Phys. C 33, 2250075 (2022).
  * [12] N. Crokidakis, Radicalization phenomena: Phase transitions, extinction processes and control of violent activities, Int. J. Mod. Phys. C 34, 2350100 (2023).
  * [13] N. Crokidakis, Recent violent political extremist events in Brazil and epidemic modeling: The role of a SIS-like model on the understanding of spreading and control of radicalism, Int. J. Mod. Phys. C 35, 2450015 (2024).
  * [14] S. Galam, M. A. Javarone, Modeling Radicalization Phenomena in Heterogeneous Populations, PLoS ONE 11(5): e0155407 (2016).
  * [15] A. K. Fantaye, Z. K. Birhanu, Mathematical Model and Analysis of Corruption Dynamics with Optimal Control, Journal of Applied Mathematics 2022: 1-16 (2022).
  * [16] D. K. Mamo, Modeling the spread dynamics of racism in cyberspace, Journal of Mathematical Modeling 8(2), 105-122 (2020).
  * [17] M. ur Rahman, M. Arfan, Z. Shah, E. Alzahrani, Evolution of fractional mathematical model for drinking under Atangana-Baleanu Caputo derivatives, Phys. Scr. 96, 115203 (2021).
  * [18] J. Sooknanan, T. A. R. Seemungal, Criminals and their models - a review of epidemiological models describing criminal behaviour, App. Math. Comput. 458, 128212 (2023).
  * [19] N. Crokidakis, A three-state kinetic agent-based model to analyze tax evasion dynamics, Physica A 415, 321-328 (2014).
  * [20] N. Crokidakis, L. Sigaud, Role of inflexible minorities in the evolution of alcohol consumption, J. Stat. Mech. 093403 (2022).
  * [21] S. Galam, The role of inflexible minorities in the breaking of democratic opinion dynamics, Physica A 381, 366-376 (2007).
  * [22] N. Crokidakis, J. S. Sá Martins, Can honesty survive in a corrupt parliament?, Int. J. Mod. Phys. C 29, 1850094 (2018).
  * [23] S. Galam, S. Moscovici, Towards a theory of collective phenomena: Consensus and attitude changes in groups, Eur. J. of Social Psychology 21, 49 (1991)
  * [24] A. C. R. Martins, S. Galam, Building up of individual inflexibility in opinion dynamics, Phys. Rev. E 87, 042807 (2013)
  * [25] N. Crokidakis, V. H. Blanco, C. Anteneodo, Impact of contrarians and intransigents in a kinetic model of opinion dynamics, Phys. Rev. E 89, 013310 (2014).
  * [26] M. Mobilia, A. Petersen, S. Redner, On the Role of Zealotry in the Voter Model, J. Stat. Mech. P08029 (2007).
  * [27] L. C. F. Latoski, W. G. Dantas, J. J. Arenzon, Curvature-driven growth and interfacial noise in the voter model with self-induced zealots, Phys. Rev. E 106, 014121 (2022).
  * [28] N. Khalil, T. Galla, Zealots in multistate noisy voter models, Phys. Rev. E 103, 012311 (2021).
  * [29] G. Verma, A. Swami, K. Chan, The impact of competing zealots on opinion dynamics, Physica A 395, 310-331 (2014).
  * [30] P. P. Klamser, M. Wiedermann, J. F. Donges, R. V. Donner, Zealotry effects on opinion dynamics in the adaptive voter model, Phys. Rev. E 96, 052315 (2017).
  * [31] M. Mobilia, Nonlinear q-voter model with inflexible zealots, Phys. Rev. E 92, 012803 (2015).
  * [32] Being Christian in Western Europe, $https://www.pewresearch.org/religion/2018/05/29/being-christian-in-western-europe/$
  * [33] J. Marro, R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, England, 1999).
  * [34] S. Biswas, Mean-field solutions of kinetic-exchange opinion models, Phys. Rev. E 84, 056106 (2011).
  * [35] B. Schonfisch, A. de Roos, Synchronous and asynchronous updating in cellular automata, Biosystems 51, 123-143 (1999).
  * [36] D. Stauffer, S. M. Moss de Oliveira, P. M. C. de Oliveira, J. S. Sá Martins, Biology, Sociology, Geology by Computational Physicists (Elsevier, 2006).
  * [37] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (Springer, New York, 2002).
  * [38] D. P. Belanger, Spin Glasses and Random Fields, edited by A. P. Young (World Scientific, Singapore, 1998).
  * [39] C. Castellano, S. Fortunato, V. Loreto, Statistical physics of social dynamics, Rev. Mod. Phys. 81, 591 (2009).