# Catalytic evolution of cooperation in a population with behavioural
bimodality

Anhui Sheng  School of Physics and Information Technology, Shaanxi Normal
University, Xi’an 710061, P. R. China  Jing Zhang  College of Information
Science and Technology, Donghua University, Shanghai 201620, P. R. China
School of Physics and Information Technology, Shaanxi Normal University, Xi’an
710061, P. R. China  Guozhong Zheng  School of Physics and Information
Technology, Shaanxi Normal University, Xi’an 710061, P. R. China  Jiqiang
Zhang  School of Physics, Ningxia University, Yinchuan 750021, P. R. China
Weiran Cai  School of Computer Science, Soochow University, Suzhou 215006, P.
R. China  Li Chen<EMAIL_ADDRESS>School of Physics and Information
Technology, Shaanxi Normal University, Xi’an 710061, P. R. China

###### Abstract

The remarkable adaptability of humans in response to complex environments is
often demonstrated by the context-dependent adoption of different behavioral
modes. However, the existing game-theoretic studies mostly focus on the
single-mode assumption, and the impact of this behavioral multimodality on the
evolution of cooperation remains largely unknown. Here, we study how
cooperation evolves in a population with two behavioral modes. Specifically,
we incorporate Q-learning and Tit-for-Tat (TFT) rules into our toy model,
where prisoner’s dilemma game is played and we investigate the impact of the
mode mixture on the evolution of cooperation. While players in Q-learning mode
aim to maximize their accumulated payoffs, players within TFT mode repeat what
their neighbors have done to them. In a structured mixing implementation where
the updating rule is fixed for each individual, we find that the mode mixture
greatly promotes the overall cooperation prevalence. The promotion is even
more significant in the probabilistic mixing, where players randomly select
one of the two rules at each step. Finally, this promotion is robust when
players are allowed to adaptively choose the two modes by real-time
comparison. In all three scenarios, players within the Q-learning mode act as
catalyzer that turns the TFT players to be more cooperative, and as a result
drive the whole population to be highly cooperative. The analysis of Q-tables
explains the underlying mechanism of cooperation promotion, which captures the
“psychologic evolution” in the players’ mind. Our study indicates that the
variety of behavioral modes is non-negligible, and could be crucial to clarify
the emergence of cooperation in the real world.

##  I 1\. Introduction

Cooperation exists in a wide range of natural and human societies Maynard
Smith and Szathmary (1997), from intra- and inter-species cooperation in
ecology to the labor division of individuals in socio-economic activities
Zelenski  _et al._ (2015); Cheney (2011); Dawkins (2006). While the
cooperation is widely observed in human beings and other animals, it’s
actually counterintuitive since defection in many scenarios is the rational
choice in terms of higher payoff, especially in some social dilemma situations
Rapoport and Chammah (1965). The key question is: how does cooperation evolve?
which is also listed as one of the twenty-five grand questions in this century
Pennisi (2005). Understanding the underlying mechanism of cooperation is of
fundamental importance in theory and has practical implications for solving
many pressing challenges, such as the climate negotiation, the trade war, and
regional conflicts, etc Hardin (1968); Kollock (1998); Milinski  _et al._
(2006).

Within the framework of evolutionary game theory, plenty of works based on the
prototypical model such as prisoner’s dilemma game (PDG) Poundstone (1993);
Doebeli and Hauert ; Axelrod (1980a); Perc and Szolnoki (2008); Axelrod
(1980b); Milinski and Wedekind (1998) have been done to address the mechanisms
for the emergence of cooperation and its maintenance. Important progresses
have been made, where several mechanisms are revealed Nowak (2006); Perc  _et
al._ (2017), such as kin selection Hamilton (1964), direct Trivers (1971) and
indirect reciprocity Nowak and Sigmund (1998), network reciprocity Nowak and
May (1992); Wang  _et al._ (2013), dynamical reciprocity Liang  _et al._
(2022), group selection Keller (1999); Maynard Smith (1964), punishment and
reward Sigmund  _et al._ (2001); Xia  _et al._ (2023), and social diversity
Perc and Szolnoki (2008); Santos  _et al._ (2008); Liang  _et al._ (2021) etc.

Recently, the reinforcement learning (RL) Watkins and Dayan (1992a); Sutton
and Barto (2018) as a new paradigm has been applied to the game-theoretic
studies, which has garnered increasing attention. Distinct from the previous
framework of social learning Bandura and Walters (1977), where individuals
imitate the strategies of those who have higher payoffs, the decision-making
in RL is based upon past experience to maximize cumulative rewards through
interaction with the environment Sutton and Barto (2018); Silver  _et al._
(2018). The rationale of RL is deeply rooted within human cognition and has
solid experimental evidence from neuroscience Lee  _et al._ (2012); Rangel
_et al._ (2008). New insights have been obtained by RL for the emergence of
cooperation Zhang  _et al._ (2020); Wang  _et al._ (2022, 2023); He  _et al._
(2022); Ding  _et al._ (2023); Geng  _et al._ (2022); Zhang  _et al._ (2024a)
and trust Zheng  _et al._ (2024), resource allocation Andrecut and Ali (2001);
Zhang  _et al._ (2019); Zheng  _et al._ (2023); Zhang  _et al._ (2024b), and
other collective human behaviors Zhang  _et al._ (2020); Tomov  _et al._
(2021); Shi and Rong (2022); Ding  _et al._ (2024). Till now, Q-learning
Watkins (1989); Watkins and Dayan (1992a) is the commonly adopted RL algorithm
Watkins and Dayan (1992b); Van Hasselt  _et al._ (2016); Wunder  _et al._
(2010); Shi and Rong (2022), where there is a Q-table scoring different
actions within different states, which provides the explainability of the
algorithm. By monitoring the evolution of Q-table, we are able to study the
preference change, which captures the psychological evolution in the players’
mind. This is critical for understanding the logic of human behaviors, and is
lacking in the traditional framework.

Although huge progresses are made in the paradigm of both social learning and
reinforcement learning, most of the game-theoretic work assumes single-mode
scenario. Namely, individuals make their decision according to a single
behavior mode, e.g. the imitation rule Szabó and Fath (2007); Perc  _et al._
(2017) (e.g. Fermi updating rule, Moran rule etc.), the reinforcement learning
Sutton and Barto (2018) (Q-learning, SASA, the Bush-Mosteller model etc.), and
reactive strategies Axelrod and Hamilton (1981); Nowak and Sigmund (1992,
1993) (e.g. tit-for-tat and generous tit-for-tat). This assumption is,
however, in contrast to observations in the real world, where individuals
behave within one mode in a given scenario, but switch to a different mode
once the scenario is changed. Behavioral experiments Traulsen  _et al._ (2010)
indicates that we human actually adopt multiple behavioral modes and the
proportions of different modes are also time-varying. In fact, such behavioral
multimodality is crucial to the adaptability to complex surroundings, and is
the key to our survival.

Some recent game-theoretic work has started to account for the diversity of
behavioral modes, in order to understand their impact on the emergence of
cooperation. Previous studies show the cooperation prevalence could either be
promoted by mixed mode of payoff-based imitation and conformity-based rule
Szolnoki and Perc (2015); Szolnoki and Chen (2018), or be deteriorated by mode
mixture of imitation and innovation Amaral and Javarone (2018). Interestingly,
when the population is composed of normal players and zealots Masuda (2012) or
“good samaritans” Zheng  _et al._ (2022), this mixture could unexpectedly lift
the prevalence of cooperation or fairness. In Ref. Ma  _et al._ (2023), Ma et
al. study the evolution of cooperation within bimodality of Fermi rule and
TFT, and they show the cooperation level is significantly promoted in both
quenched and annealed manners, and this promotion is robust to the underlying
topology. Han et. al propose a hybrid learning combining social learning and
the Bush-Mosteller model as self-learning Han  _et al._ (2022), they reveal
cooperation prevalence is promoted within some parameter regions. As the new
paradigm of RL gradually releases its potential in explaining cooperation, we
are interested in the following questions: what is impact of the behavioral
multimodality on the evolution of cooperation when the reinforcement learning
is engaged?

In this work, we focus on the impact of behavioral multimodality in the
evolution of cooperation. Specifically, we investigate the bimodal mixture of
two fundamentally different update rules — TFT and Q-learning. As a new
paradigm, players within Q-learning mode resort to continuously policy
modification as introspection, while players within TFT mode just repeat what
their opponents have done to them Nowak and Sigmund (1993). The mode mixture
is implemented in three manners: structural mixing, probabilistic mixing, and
adaptive mixing. In all three implementations, the prevalence of cooperation
is significantly improved compared to either single-mode scenario. By
monitoring the Q-table, we clarify the mechanism of the cooperation promotion
and the impact of the bimodality, which provides valuable insights into the
psychologic evolution of population.

The rest of this paper is structured as follows. Sec. 2 describe our model
with three implementations. Sec. 3 shows the results of the cooperation
evolution. Sec. 4 discusses the underlying mechanism behind the cooperation
promotion. Finally, Sec. 5 provides conclusions and discussions of this work.

State Action | C ($a_{1}$) | D ($a_{2}$)
---|---|---
$s_{1}=(a_{1},0)$ | $Q_{s_{1},a_{1}}$ | $Q_{s_{1},a_{2}}$
$s_{2}=(a_{2},0)$ | $Q_{s_{2},a_{1}}$ | $Q_{s_{2},a_{2}}$
$s_{3}=(a_{1},1)$ | $Q_{s_{3},a_{1}}$ | $Q_{s_{3},a_{2}}$
$s_{4}=(a_{2},1)$ | $Q_{s_{4},a_{1}}$ | $Q_{s_{4},a_{2}}$
$s_{5}=(a_{1},2)$ | $Q_{s_{5},a_{1}}$ | $Q_{s_{5},a_{2}}$
$s_{6}=(a_{2},2)$ | $Q_{s_{6},a_{1}}$ | $Q_{s_{6},a_{2}}$
$s_{7}=(a_{1},3)$ | $Q_{s_{7},a_{1}}$ | $Q_{s_{7},a_{2}}$
$s_{8}=(a_{2},3)$ | $Q_{s_{8},a_{1}}$ | $Q_{s_{8},a_{2}}$
$s_{9}=(a_{1},4)$ | $Q_{s_{9},a_{1}}$ | $Q_{s_{9},a_{2}}$
$s_{10}=(a_{2},4)$ | $Q_{s_{10},a_{1}}$ | $Q_{s_{10},a_{2}}$
Table 1: Q-table for every player within the Q-learning mode. The state of a
given player $s_{1,...,10}$ corresponds to the 10 combination of player’s own
action ($a_{1}$ or $a_{2}$) and the number of cooperators in his four nearest
neighbours.

##  II 2\. Methods

We consider a population of $N$ individuals who play the Prisoner’s Dilemma
Game (PDG) on a square lattice $L\times L=N$ with periodic boundary
conditions, and $L=100$ if not stated otherwise. For each paired players, each
can choose to be either cooperator (C) or defector (D), both players get a
reward $R$ for mutual cooperation, mutual defection yields the punishment $P$
to each, and mixed encounter gives the cooperator the sucker’s payoff $S$ yet
the temptation $T$ for the defector. The order of these four payoffs $T>R>P>S$
and $2R>T+S$ are required for the PDG. For clarity, the payoff matrix of the
PDG is as follows:

$\Pi=\left(\begin{array}[]{ll}\Pi_{CC}&\Pi_{CD}\\\
\Pi_{DC}&\Pi_{DD}\end{array}\right)=\left(\begin{array}[]{cc}R&S\\\
T&P\end{array}\right),$ (1)

where the values of these four payoffs are reconfigured as $R=1$, $S=-r$,
$T=1+r$, $P=0$ in our study, making a strong version of PDG. The dilemma lies
in the fact that although mutual cooperation is more beneficial for the
collective payoff, they still defect driven by the self-interest and mutual
defection is the only Nash equilibrium of the game Nowak and May (1992). The
value of $r$ $(0<r<1.0)$ controls the strength of the social dilemma, the
larger value of $r$, the stronger the social dilemma.

\begin{overpic}[width=130.08731pt]{SM1-eps-converted-
to.pdf}\put(2.5,65.0){(a)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM2-eps-
converted-
to.pdf}\put(2.5,65.0){(b)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM3-eps-
converted-to.pdf}\put(2.5,65.0){(c)}\end{overpic}

Figure 1: The evolution of cooperation with SM. $\omega$ is the proportion of
players within the TFT mode. (a) Time series of cooperation level $f_{C}$ for
different proportions of players within the TFT mode. (b) The prevalence of
cooperation $f_{C}$ as a function of mode proportion $\omega$, separately
computed for people within TFT mode, Q-learning mode, and for the whole
population; the inset provides three time series of $f_{C}$ for the whole
population for three cases of $\omega\rightarrow 1$. (c) Heat map for the
cooperation prevalence $f_{C}$ in the $r-\omega$ parameter space. The game
parameter $r=0.1$ in (a) and (b). Other parameters: $\varepsilon=0.01$,
$\alpha=0.1$, and $\gamma=0.9$.

In our model, individuals play the game within two possible behavioral modes,
i.e. the TFT rule and Q-learning. Within TFT mode, players just repeat what
their neighbors play against it. Within the Q-learning mode Watkins and Dayan
(1992b), each player is equipped with a Q-table that can be regarded as a
policy and guides one’s decision making, see Table 1. The idea behind Q-table
is to score the different action $a\in A$ within the different state $s\in S$,
reflected in the value of the corresponding item $Q_{s,a}$. In our model, the
action set $A=\\{C,D\\}$, and the state set $S=\\{s_{1},...,s_{10}\\}$, which
includes the number of cooperators in its four nearest neighbors and its own
action in the last round. The items in the table $Q_{s,a}$ are the action-
value function that estimates the value of action $a$ within the given state
$s$; if $Q_{s,a}>Q_{s,\hat{a}}$, this means that the action $a$ is of higher
value than action $\hat{a}$ within state $s$, and is suggested to be chosen by
Q-learning. Within Q-learning mode, players revise their Q-tables to maximize
their accumulated payoffs.

In our bimodal behaviors model, we consider three mixing implementations –
_structure mixing_ (SM), _probabilistic mixing_ (PM), and _adaptive mixing_
(AM). In the SM scenario, each player is assigned with one of the two
behavioral modes at the very beginning, specifically the TFT mode is chosen
with probability $\omega$ and Q-learning with $1-\omega$, and they behave in
the same mode throughout the whole evolution. Instead, players randomly choose
one of the two modes at each round in the PM, with the same probability
denotation. Different from these two implementations where the probabilities
are fixed, the AM procedure follows an adaptive manner: when the payoff of
player $i$ is lower than the value in previous round i.e.
$\Pi_{i}(t)<\Pi_{i}(t-1)$, this unsatisfactory outcome makes him switch from
the adopted mode to the other.

To be specific, we first detail the procedure of SM as follows:

(1) At the beginning ($t=0$), every player chooses the TFT and Q-learning mode
at random, respectively with the probability $\omega$ and $1-\omega$. Once the
mode is chosen, the player will stick to it throughout the whole evolution.
For initialization, each player is endowed an action $a\in A$ with an equal
probability, and for those within the Q-learning mode, each has a Q-table with
all items $Q_{s,a}$ being initialized from 0 to 1 independently.

(2) At round $t>0$, if player $i$ is within TFT mode, one of its four
neighbors is randomly chosen (say player $j$), and then $i$ adopts the action
chosen by player $j$ in the last round, i.e. $a_{i}(t)=a_{j}(t-1)$.

(3) Otherwise, if player $i$ is within the Q-learning mode, he makes his moves
according to the Q-learning algorithm. With the probability $\epsilon$, player
$i$ randomly choose an action $a\in A$ to make the trial-and-error
exploration. Otherwise, player $i$ selects the action based upon

$\displaystyle
a_{i}(t)\leftarrow\arg\max_{a}\left\\{Q_{s_{i},a_{i}}(t)\right\\},a_{i}\in{A},$
(2)

where $\arg\max_{a}\left\\{Q_{s_{i},a}(t)\right\\}$ is the action
corresponding to the maximum Q-value in the row of state $s_{i}$, where
$s_{i}=(a_{i}(t-1),n_{i}^{C}(t-1))$ and $n_{i}^{C}(t-1)$ is the number of
cooperators within player $i$’s four neighbors at the end of round $t-1$.

(4) After all players make their moves [i.e. complete step (2) or (3)], those
who are within the Q-learning mode proceed to revise their Q-tables as
follows:

$\displaystyle
Q_{s_{i},a_{i}}(t\\!+\\!1)\\!=\\!(1\\!-\\!\alpha)Q_{s_{i},a_{i}}(t)\\!+\\!\alpha\left[\Pi_{i}(t)\\!+\\!\gamma\max_{a_{i}^{\prime}}Q_{s_{i}^{\prime},a_{i}^{\prime}}(t)\right],$
(3)

where $s_{i},a_{i}$ represent the current state and action of player $i$, and
$s_{i}^{\prime}=(a_{i}(t),n_{i}^{C}(t))$ is the new state and $a_{i}^{\prime}$
is the new possible action at $t+1$. $\alpha\in(0,1]$ is the learning rate,
which reflects the intensity of memory effect; the larger value of $\alpha$,
the more amount of Q-value is renewed, meaning more forgetful of the player.
$\Pi_{i}(t)$ is the player $i$’s total payoff by playing against its four
neighbors at the current round $t$, which follows the payoff matrix Eq. (1).
$\gamma\in[0,1]$ is the discount factor, determining the weight of the future
reward since $\max_{a_{i}^{\prime}}Q_{s_{i}^{\prime},a_{i}^{\prime}}(t)$ is
the maximum Q-value one can expect within its new state $s_{i}^{\prime}$. This
completes the Q-table and state update, and one single round is done when
there are players within the Q-learning mode.

(5) Repeat steps(2)-(4) until the system reaches statistically stable state.

\begin{overpic}[width=130.08731pt]{PM1-eps-converted-
to.pdf}\put(2.5,65.0){(a)}\end{overpic}\begin{overpic}[width=130.08731pt]{PM2-eps-
converted-
to.pdf}\put(2.5,65.0){(b)}\end{overpic}\begin{overpic}[width=130.08731pt]{PM3-eps-
converted-to.pdf}\put(2.5,65.0){(c)}\end{overpic}

Figure 2: The evolution of cooperation with PM. $\omega$ is the probability
that player evolves according to the TFT rule within each round. (a) Time
series of cooperation level $f_{C}$ for several typical $\omega$. (b) The
prevalence of cooperation $f_{C}$ versus $\omega$, of all player groups,
separately computed for people within TFT mode, Q-learning mode, and for the
whole population; the inset provides three time series of $f_{C}$ for the
whole population for three cases of $\omega\rightarrow 1$. (c) Heat map for
the cooperation prevalence $f_{C}$ in the $r-\omega$ parameter space. The game
parameter $r=0.1$ in (a) and (b). Other parameters: $\varepsilon=0.01$,
$\alpha=0.1$, and $\gamma=0.9$.

Note that, the above procedure follows a typical synchronous update process,
and each round [steps (2)-(4)] constitutes a Monte Carlo step, where every
player updates exactly once.

The difference in the procedure for PM is within in step (1), where no
behavioral mode needs to fixed for the player. Instead, at the beginning of
each round $t$, each player randomly selects the TFT mode with the probability
$\omega$, otherwise the Q-learning mode is chosen. Notice that, since all
players could behave within two modes from time to time, every player will be
assigned with a Q-table at the very beginning ($t=0$) in step (1).

Compared to PM, the procedure for AM differs in the mode selection. Within
step (1), one of the two modes is randomly chosen for initialization, but at
the beginning of each new round $t+1$, player $i$ compares the payoff obtained
at present round $\Pi_{i}(t)$ with the one at the last round $\Pi_{i}(t-1)$;
if $\Pi_{i}(t)\geq\Pi_{i}(t-1)$, player $i$ keeps within the old behavioral
mode, otherwise player $i$ switches to the other mode. Notice that in both PM
and AM, players revise their Q-tables at each time step even though they may
not behave in the Q-learning mode to mimic the observed reality where
individuals continuously draw lesson from their experiences. The pseudocode
for the details of three implementations are given in Appendix A.

As can be seen, the probability $\omega$ determines the mixture ratio of the
two behavioral modes in SM and PM. As $\omega\rightarrow$ 0 or 1, our mode is
recovered to the evolution with pure Q-learning algorithm, or with pure TFT
rule Nowak and Sigmund (1992, 1993), and we are interested in the mixture case
where $0<\omega<1$. By contrast, this mixture ratio is adaptively adjusted in
AM implication, where the evolution within either pure mode is rarely seen. In
our practice, we fix $\varepsilon=0.01$ to allow some degree of random
exploration. We compute the cooperation fraction $f_{C}$ defined as the number
of action $C$ divided by $N$ as the primary order parameter to measure the
overall preference of the population.

##  III 3\. Results

###  III.1 3.1. Structural mixing

Let’s first examine the impact of mode mixture on the evolution of cooperation
with the SM implementation, where $\omega$ is the proportion of players within
TFT mode, and they remain in the same mode throughout the evolution. Fig. 1
reports the results on the 2d square lattice with SM. As the benchmarks, let’s
check the two extreme cases (i.e. $\omega=0,1$), which corresponds to the
single-mode scenarios of Q-learning and TFT mode, respectively. Specifically,
when the population is within pure Q-learning mode, the resulting
$f_{C}\approx 0.29$, meaning that players learn to cooperate to some degree
but the prevalence is not high. When the population is within pure TFT mode,
the resulting cooperation prevalence depends on the initial condition, here in
our study $f_{C}\approx f_{C}(t=0)\approx 0.5$.

When the two modes are mixed ($0<\omega<1$), their cooperation prevalence is
promoted in the given cases (i.e. $\omega=0.25,0.5,0.75$) compared to the two
single-mode scenarios, see Fig. 1(a). In all three mixing cases, $f_{C}>0.5$
and the prevalence is highest for $\omega=0.75$. Given the pure TFT scenario,
this indicates there is a non-monotonic dependence of cooperation prevalence
on $\omega$. This is confirmed by plotting the cooperation prevalence versus
the mode proportion $\omega$, shown in Fig. 1(b). As seen, the overall
prevalence $f_{C}$ increases with escalating $\omega$, almost full cooperation
is reached when $0.8<\omega<1$. Detailed statistical analysis shows the
cooperation prevalence is slightly higher for those within the TFT mode than
those Q-learners. The trend of $\omega\rightarrow 1$ indicates that the
addition of very few Q-learners into the TFT population is sufficient to drive
the whole system into full-cooperation state, and this is indeed confirmed by
our simulations [the inset within Fig. 1(b)]. This means that, individuals in
Q-learning mode play the role of cooperation catalysts, very few of them are
able to promote the overall cooperation.

Fig. 1(c) provides the heat map for the cooperation prevalence in $r-\omega$
parameter domain to further investigate the impact of the game parameter $r$.
As expected, $f_{C}$ is reduced as the dilemma strength $r$ is increased, and
the transitions from high cooperation to low cooperation state are all of
continuous type.

\begin{overpic}[width=173.44534pt]{AM1-eps-converted-
to.pdf}\put(2.5,65.0){(a)}\end{overpic}\begin{overpic}[width=173.44534pt]{AM2-eps-
converted-to.pdf}\put(2.5,65.0){(b)}\end{overpic}

Figure 3: Time series of cooperation prevalence with AM. (a) The prevalence of
cooperation $f_{C}$ computed separately for players with TFT and Q-learning,
together with the whole population for comparison. (b) The time evolution of
probabilities of adopting the two modes by the population. Other parameters:
$\varepsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$, and $r=0.1$.

###  III.2 3.2. Probabilistic mixing

Fig. 2 shows the impact of the mode mixture on the evolution cooperation with
the PM implementation. Bear in mind, here the behavioral mode for each players
is switching from time to time, and $\omega$ is interpreted as the probability
to behave within the mode of TFT rule at each round.

Compared to Fig. 1, the two single-mode scenarios ($\omega=0,1$), shows
exactly the same prevalence in cooperation, as expected. The difference,
however, lies in the observations in the mixture cases, where the cooperation
prevalences in all three given cases ($\omega=0.25,0.5,0.75$) jump to be very
high $f_{C}\rightarrow 1$ in the long term. Fig. 2(b) clearly shows this
difference compared to Fig. 1(b). We observe that there is a wide of the
parameter $\omega$, where the cooperation prevalence $f_{C}>0.8$. Also, the
prevalences computed separately for people within TFT and Q-learning mode show
that the two are now statistically the same, no quantitative difference is
detected. Once again, when $\omega\rightarrow 1$, the presence of very few
Q-learning mode acts as the cooperation catalysts, driving the population into
almost full cooperation [see the inset within Fig. 2(b)]. Though, the time
series of $f_{C}$ indicates that the evolution with too few of Q-learning mode
requires longer transients to have $f_{C}\rightarrow 1$. The heat map shown in
Fig. 2(c) shows a widened region for cooperation compared to Fig. 1(c), and
the cooperation-defection transition becomes more abrupt.

Note that, the long term evolution (e.g. $t>10^{5}$) in Fig. 2(a) exhibits
some large fluctuations in the cooperation prevalence, though a decent
cooperation prevalence is still maintained. These fluctuations are due to the
Q-learning players, who turn to defection to obtain higher payoffs thus reduce
$f_{C}$, but the presence of players within TFT mode inhibits this trend. This
forces the Q-learning players turn back to cooperation again and the
prevalence increase again. This process repeats again and again, and an
oscillation is formed. Detailed discussion is provided in Appendix B.

###  III.3 3.3. Adaptive mixing

Different from the above two implementations, where the two mode fractions are
predetermined, the mode adoption in AM is adjusted from time to time according
to the payoff changes. Fig. 3 shows a typical evolution where both the
cooperation preference and the probabilities of mode adoption are monitored.
Fig. 3(a) shows that after a slight decay in $f_{C}$ for $t<10^{2}$, the
cooperation prevalence continuously increases till $f_{C}\rightarrow 1$ is
reached around $t>10^{4}$. Detailed examination shows that the prevalence
obtained within Q-learning and TFT are almost identical. Further evolution
($t>10^{5}$), however, some fluctuations in cooperation appear, with the same
observations made in PM [see Fig. 2(a)].

Fig. 3(b) shows that the adoption probability for the two modes remains
largely identical. Though a detailed examination shows that the probability of
adopting Q-learning mode is slightly larger than TFT mode, and later there is
a crossover at around $t=10^{3}$.

##  IV 4\. Mechanism analysis

_Structure mixing (SM)_ — To reveal the mechanism why the mode mixing promotes
the cooperation prevalence, let’s move back the evolution within the SM
implementation. To gain some intuition, let’s first examine the spatial
pattern of states with different probabilities $\omega$, shown in Fig. 4. Fig.
4(a) illustrates that when all players behave within the Q-learning mode
(i.e., $\omega=0$ ), the majority of players in the group choose strategy $D$,
and no obvious clusters are seen for those cooperators. Instead, when the
population behave within only TFT mode, the evolution experiences a coarse-
graining-like process that forms several large clusters [see Fig. 4(c)]. The
configurations of these clusters are time-varying, but the cooperation
prevalence $f_{C}$ remains at one half on average. The typical mixture case is
shown in Fig. 4(b), where half of the population adopts the TFT mode while the
other half behaves with Q-learning ($\omega=0.5$). The resulting cooperation
prevalence $f_{C}\approx 0.86$, much higher than the value in either single-
mode scenario. Different from Fig. 4(c), no obvious cluster configuration is
observed, the pattern is rather similar to Fig. 4(a). This means that the
emergence of cooperation could not be attributed to the network reciprocity, a
well-known mechanism for structured population in the imitation learning
framework Szabó and Fath (2007); Perc  _et al._ (2017).

\begin{overpic}[width=130.08731pt]{SM_pattern1-eps-converted-
to.pdf}\put(2.5,65.0){(a)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_pattern2-eps-
converted-
to.pdf}\put(2.5,65.0){(b)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_pattern3-eps-
converted-to.pdf}\put(2.5,65.0){(c)}\end{overpic}

Figure 4: Typical spatial patterns on $2d$ square lattice with SM at
$t=2\times 10^{6}$ for three $\omega$: (a) $\omega=0$, (b) $\omega=0.5$, (c)
$\omega=1$. Both the strategy and mode information are color coded, e.g. blue
pixels ($D_{Q}$) denote the players are within the Q-learning mode and adopt
defection as their strategies. $f_{C}\approx 0.29$, 0.86, and 0.51 in (a-c),
respectively. Other parameters: $\varepsilon=0.01$, $\alpha=0.1$,
$\gamma=0.9$, and $r=0.1$.

\begin{overpic}[width=173.44534pt]{SM_DeltaQ1-eps-converted-
to.pdf}\put(2.5,65.0){(a)}\end{overpic}\begin{overpic}[width=173.44534pt]{SM_DeltaQ2-eps-
converted-to.pdf}\put(2.5,65.0){(b)}\end{overpic}

Figure 5: The temporal evolution of the average Q-value difference
$\overline{\Delta Q}_{s_{j}}$ within all 10 states in SM for (a) $\omega=0$
and (b) $\omega=0.5$. The inset in (a) is the probability density function
(PDF) of the Q-learning players with $\Delta Q^{i}_{s_{j}}>0$ in different
states at $t=2\times 10^{6}$, where $\Delta
Q^{i}_{s_{j}}=Q_{s_{j},C}^{i}-Q_{s_{j},D}^{i}$. The black dashed line in both
plots corresponds to $\overline{\Delta Q}_{s_{j}}=0$ for reference. Other
parameters: $\varepsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$, $r=0.1$.

To grasp the mechanism, let’s turn to the evolution of the Q-tables, which
captures the preference evolution of the players within the Q-learning mode.
Specifically, we define the average Q-value difference within a given state
$s_{j}\in S$ as

$\displaystyle\overline{\Delta
Q}_{s_{j}}=\frac{1}{N}\sum_{i=1}^{N}\left(Q_{s_{j},C}^{i}-Q_{s_{j},D}^{i}\right).$
(4)

When $\overline{\Delta Q}_{s_{j}}>0$, players are more inclined to choose
cooperation within state $s_{j}$; conversely, they prefer to defect when
$\overline{\Delta Q}_{s_{j}}<0$. As the benchmark, let’s first look at the
scenario of pure Q-learning mode ($\omega=0$), see Fig. 5(a). It shows that
the $\overline{\Delta Q}_{s_{j}}$ of all states go up and down, but
eventually, all of them become negative $\overline{\Delta Q}_{s_{j}}<0$
($i=1,2,...,10$). This means that on average players prefer to defect in all
possible states, though the probability density function (PDF) shown [inset of
Fig. 5(a)] indicates that there is still some probability for $\Delta
Q^{i}_{s_{j}}>0$ at the individual level, which explains the low but non-
vanishing cooperation preference in the pure mode of Q-learning. Notice that,
the preference for those cooperators (within $s_{1,3,5,7,9}$) and defectors
(within $s_{2,4,6,8,10}$) have opposite trends by evolution initially, though
all become negative, in particularly for those who are surrounded by defectors
are determined to defect where $\overline{\Delta Q}_{s_{1,2}}\rightarrow-4$.

The cooperation preference in the typical mixing case is dramatically
different, as shown in Fig. 5(b). As can be seen, the average Q-value
differences $\overline{\Delta Q}_{s_{j}}$ initially are close to zero, meaning
that no obvious preference is developed when they are unfamiliar with the
surroundings. After the transient ($t\gtrsim 10^{4}$), three kinds of
preferences are formed as clearly shown in $\overline{\Delta Q}_{s_{j}}$. For
those surrounded by four cooperators ($s_{9,10}$), they reciprocate with C as
$\overline{\Delta Q}_{s_{9,10}}>15$ in the end; instead for those surrounded
by four defectors ($s_{1,2}$), they reciprocate with D as $\overline{\Delta
Q}_{s_{1,2}}<-15$. For the rest cases, $\overline{\Delta Q}_{s_{3,...,8}}>0$
means that they all prefer cooperation on average, and detailed examination
shows that the more cooperators in the neighborhood, the stronger preference
in cooperation (with a larger value of $\overline{\Delta Q}_{s_{j}}$).

To further unveil why mixing promotes cooperation, we show the time series of
the cooperation prevalence $f_{C}$ of players in TFT and Q-learning mode
separately for three typical probabilities ($\omega=0.25,0.5,0.75$), and
compare them to the overall prevalence, shown in Fig. 6(a-c). The results
indicate that these prevalences are actually different. In particular, when
the players within the Q-learning mode dominate (e.g. $\omega=0.25$), the
value of $f_{C}$ is obviously larger than the one of Q-learning, though the
difference is vanishing as $\omega$ becomes larger. This means the distinct
role of two behavioral modes in the evolution.

\begin{overpic}[width=130.08731pt]{SM_ts1-eps-converted-
to.pdf}\put(2.5,65.0){(a)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_ts2-eps-
converted-
to.pdf}\put(2.5,65.0){(b)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_ts3-eps-
converted-to.pdf}\put(2.5,65.0){(c)}\end{overpic}

\begin{overpic}[width=130.08731pt]{SM_ts4-eps-converted-
to.pdf}\put(2.5,65.0){(d)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_ts5-eps-
converted-
to.pdf}\put(2.5,65.0){(e)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_ts6-eps-
converted-to.pdf}\put(2.5,65.0){(f)}\end{overpic}

Figure 6: Time series of the two kinds of cooperation prevalences with SM for
different $\omega$: (a, d) $\omega=0.25$, (b, e) $\omega=0.5$, (c, f)
$\omega=0.75$. The first row (a-c) are the cooperation prevalence computed
separately for players within the TFT, Q-learning mode, and the whole
population. The second row (d-f) are the cooperation propensity of the players
who are within TFT versus TFT, TFT versus Q-learning, Q-learning versus
Q-learning. Other parameters: $\varepsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$,
$r=0.1$.  Figure 7: The temporal evolution of the average Q-value difference
$\overline{\Delta Q}_{s_{j}}$ within all 10 states in PM for $\omega=0.25$.
The black dashed line corresponds to $\overline{\Delta Q}_{s_{j}}=0$ for
reference. Other parameters: $\varepsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$,
$r=0.1$.

Since two behavioral modes are present in the population, there are three
types of interaction by combination:

(i) _Q-learning versus Q-learning_ (or Q-Q in short) interactions: when only
Q-learners play against Q-learners, they prefer to defect by learning on
average, and only a low level of cooperation is achieved ($f_{C}\approx 0.29$)
according to the analysis of Fig. 5(a).

(ii) _TFT versus TFT_ (T-T) interactions: since TFT players only repeat what
their neighbors did to them, the population experiences a coarse-graining-like
process and forms several clusters with the cooperation prevalence being
approximately equal to the initial value [see Fig. 4(c)]. In our case,
$f_{C}\approx 0.5$ if there are only TFT players.

(iii) _Q-learning versus TFT_ (Q-T) interactions: when Q-learners encounter
TFT players, they are able to learn that C is a better choice to maximize the
rewards, since the action choice of D immediately leads to revenge from their
TFT neighbors. This is almost guaranteed if there are two or more TFT
neighbors. The confirmation is given by a detailed analysis of
$\overline{\Delta Q}_{s_{j}}$ for Q-learners surrounded by different TFT
players in Appendix C.

Put together, the mechanism behind the promotion of cooperation in mode-mixing
population is as follows: initially, for the domains where Q-Q or T-T
interactions are dense the cooperation prevalence $f_{C}$ keeps at the value
$f_{C}(t=0)\approx 0.5$; the promotion comes from the Q-T interactions, where
Q-learners learn to be cooperators, and as a result the TFT counterparts then
also reciprocate by being cooperators. The high cooperation prevalence within
Q-T interactions then triggers a series of propagation of C through T-T
interactions (ii), and the same level of $f_{C}$ is expected for type (ii) and
(iii) is expected. This is confirmed in Fig. 6(d-f). However, the prevalence
in Q-Q pair is qualitatively lower, and this is because the cooperation
prevalence in type (i) is lower and their cooperation prevalence is to
difficult to improve.

The nature difference in type (i-iii) then also explains why the cooperation
prevalence for players in TFT mode is higher than those with Q-learning. This
is because the promotion originated from T-Q interactions fail to propagate to
the Q-Q interactions, resulting in a relatively lower $f_{C}$ in Q-learning
players. The overall picture also explains the dependence of $f_{C}$ on the
parameter $\omega$ seen in Fig. 1: the larger proportion of TFT players
$\omega$ leads to less Q-Q interactions, whereby a higher cooperation
prevalence $f_{C}$ is expected. Notice that, since Q-T interaction is the
source of cooperation promotion, a small amount of players need to play within
the Q-learning mode, playing the role of catalyst to trigger the cooperation
promotion for the whole population.

_Probabilistic mixing (PM)_ — The above mechanism for SM is still working for
explaining the cooperation promotion in the PM implementation, where Q-T
interactions still account for the source of the promotion. The difference in
PM is that, Q-Q interactions are now hard to persist for any given paired
players, as the Q-Q interactions are destroyed from time to time due to the
probabilistic mode switching in PM. As a consequence, there is less chance for
Q-Q interactions to persist, and the cooperation prevalence is thus driven by
either Q-T or propagates along T-T interactions, finally a high level of
$f_{C}$ is seen [see Fig. 2]. And since $\omega$ only controls the probability
for the two modes, the specific value for $0<\omega<1$ does not influence the
resulting $f_{C}$ as no Q-Q interactions in SM are present for the given $r$.
Though, Fig. 2(c) indicates that for the case of larger $r$, a higher $\omega$
is required for a decent level of cooperation, as this increases the
interaction of type (ii) but reduces the interactions of type (i), and keeps
some interactions of type (iii).

In fact, we can also monitor the evolution of cooperation preference by
computing the average Q-value difference $\overline{\Delta Q}_{s_{j}}$, as
shown in Fig. 7. It shows that the preference in six states $\overline{\Delta
Q}_{s_{j}}>0$ ($j=5,6,7,8,9,10$), while the rest four $\overline{\Delta
Q}_{s_{j}}<0$ ($j=1,2,3,4$), this division becomes most prominent at round
$10^{4}$. Compared to Fig. 5(b), both show the similar reciprocity property:
cooperate when all neighbors cooperate, while defect when all neighbors
defect. Afterwards, however, there are declines all absolute values of
$\overline{\Delta Q}_{s_{j}}$, this is especially true for $\overline{\Delta
Q}_{s_{9,10}}$. This means players within Q-learning mode may attempt to
defect to achieve higher rewards, but the presence of TFT mode immediately
feedback defection to their action. This inhibits the decline of
$\overline{\Delta Q}_{s_{j}}$, and finally both values of $\overline{\Delta
Q}_{s_{9,10}}$ remain positive to guarantee that the choice of C is
statistically more preferred.

Figure 8: The PDF of players’ cooperation preference (i.e. the probability for
$\Delta Q^{i}_{s_{j}}=Q_{s_{j},C}^{i}-Q_{s_{j},D}^{i}>0$) within all 10 state,
and the density of these states $\rho_{s}$ per se. Data is sampled from the
last $10^{4}$ before the evolution of $10^{6}$ steps is completed. Other
parameters: $\varepsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$, $r=0.1$.

_Adaptive mixing (AM)_ — Compared to the PM implementation, the logic of mode
adoption differs only in the details, i.e. the mode is not adopted
probabilistically but according to a predetermined rule given in model part.
From the interaction point of view, the mechanism behind the cooperation
prevalence is essentially the same as PM.

To precisely characterize the cooperation prevalence, we sample the value of
$\Delta Q^{i}_{s_{j}}$ at the individual level defined as $\Delta
Q^{i}_{s_{j}}=Q_{s_{j},C}^{i}-Q_{s_{j},D}^{i}$, where $\Delta Q^{i}_{s_{j}}>0$
means that player $i$ prefers $C$ as its action within state $s_{j}$,
otherwise $D$ is favored. Fig. 8 shows the percentage of $\Delta
Q^{i}_{s_{j}}>0$ in all 10 states. We see that very low probabilities are seen
in $s_{1,2}$, meaning players choose to defect when they are surrounded by
defectors. By contrast, when the surroundings become cooperative, e.g.
$s_{9,10}$, the probability to cooperate becomes high. For most states (e.g.
$s_{5-10}$), the cooperation probability is larger than 0.6. In fact, Fig. 8
also shows that the state $s_{10}$ is dominating, where the density
$\rho_{s_{10}}$ is nearly as high as $90\%$, meaning that players are mostly
located in a very cooperative surrounding, whereby a decent level of
cooperation is expected since the probability of $\Delta Q^{i}_{s_{10}}>0$
goes to be 1.

##  V 5\. Conclusion and discussion

Motivated by the behavior multimodality in the real world, here we develop a
toy model of mixing Q-learning and TFT rules. While Q-learning is a typical RL
algorithm, the decision-making is based upon the learned experience, TFT is a
well-known reactive strategy. The combination of these two is not
contradictory, but complementary. In the first implementation (SM), where each
player randomly chooses one mode and sticks to it, we reveal that the mode
mixture is able to promote and a higher proportion of players in TFT mode
generally yields a higher level of cooperation. In the second implementation
(PM), players behave probabilistically in the either mode in each round, the
cooperation promotion is even more pronounced than the case of SM, where
nearly full cooperation is expected for a wide range of mixing probability. In
the third implementation (AM), players are allowed to adopt both modes
adaptively, where also nearly full cooperation is seen and the two modes are
adopted by the population with approximately equal chance.

Mechanism analysis shows that in all three implementations, players in
Q-learning modes play the role of cooperation catalyzer, where they choose to
cooperate with their TFT neighbors and ultimately drive the whole population
to a high cooperation level. Since Q-tables reflect the behavioral preference
of individuals and exactly captures the thoughts in players’ mind, monitoring
the statistical properties of Q-tables thus provides important information
about the psychologic evolution of the population. This is crucial for
understand the logic of human behaviors, and is the great advantage of
Q-learning algorithm, which is lacking in the previous framework of social
learning Szabó and Fath (2007); Perc  _et al._ (2017).

Given the nontrivial findings obtained in our work, we hope that efforts could
be made to infer the different types of learning from the behavioral
experiments. This is an important step towards understanding human behaviors,
including cooperation, fairness, honesty, and so on.

##  VI Acknowledgments

This work was supported by the National Natural Science Foundation of China
[Grants Nos. 12075144,12165014], and Fundamental Research Funds for the
Central Universities (GK202401002).

##  Appendix A Appendix A: The pseudocode for the three implementations

To provide detailed procedure of how the three implementations are realized,
the pseudocode for SM, PM, and AM is shown below [see Fig. 9], respectively.

\begin{overpic}[width=130.08731pt]{Algorithm_SM.png} \end{overpic}

\begin{overpic}[width=130.08731pt]{Algorithm_PM.png} \end{overpic}

\begin{overpic}[width=130.08731pt]{Algorithm_AM.png} \end{overpic}

Figure 9: The procedure for our three implementations, where Algorithm 1, 2, 3
are the pseudocodes for SM, PM, and AM, respectively.

##  Appendix B Appendix B: The phenomena of cooperation oscillation

To examine the large fluctuations present in the long term evolution, Fig. 10
provides a more clear illustration, where oscillatory fluctuations are seen in
the inset. To develop an intuition of how the oscillation is formed, some
typical snapshots are shown in Fig. 11. As seen, the fluctuations all come
from the players in Q-learning mode, who turn from the choice of C to D. This
is triggered by the exploration events which indeed bring them higher payoffs
in short term, and thus the number of defective Q-learning players increases
[Fig. 11(a) to Fig. 11(b)]. As the surroundings become more defective, other
Q-learners tend also to defect, and thus defectors appear in the clustered
form. However, as defective Q-learning players increases, their payoffs are
continually decreased, and this decreasing trend is strengthened by their TFT
neighbors as them immediately also choose D. This lowered reward trend reduces
the corresponding Q-values of action C, and therefore those defective players
in Q-learning mode turn back to C to maximize their payoffs. But as the
surrounding goes back to a high cooperation level [e.g. Fig. 11(c)], the
exploration events again cause the decline of $f_{C}$ and the process repeats
again and again.

Figure 10:  Time series of cooperation level $f_{C}$ at $\omega=0.25$. The
inset is a blowup of the time series in which the oscillatory phenomenon
starts. Other parameters: $\varepsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$,
$r=0.1$.

\begin{overpic}[width=130.08731pt]{PM_pattern_t=43000-eps-converted-
to.pdf}\put(2.5,67.0){(a)}\end{overpic}\begin{overpic}[width=130.08731pt]{PM_pattern_t=50000-eps-
converted-
to.pdf}\put(2.5,67.0){(b)}\end{overpic}\begin{overpic}[width=130.08731pt]{PM_pattern_t=53000-eps-
converted-
to.pdf}\put(2.5,67.0){(c)}\end{overpic}\begin{overpic}[width=130.08731pt]{PM_pattern_t=67000-eps-
converted-
to.pdf}\put(2.5,67.0){(d)}\end{overpic}\begin{overpic}[width=130.08731pt]{PM_pattern_t=71000-eps-
converted-
to.pdf}\put(2.5,67.0){(e)}\end{overpic}\begin{overpic}[width=130.08731pt]{PM_pattern_t=78000-eps-
converted-to.pdf}\put(2.5,67.0){(f)}\end{overpic}

Figure 11: Typical spatial patterns on $2d$ square lattice with PM for several
time steps: (a) $t=43000$, (b) $t=50000$, (c) $t=53000$, (d) $t=67000$, (e)
$t=71000$, (f) $t=78000$. Both the strategy and mode information are color
coded, e.g. blue pixels ($D_{Q}$) denote the players are within the Q-learning
mode and adopt defection as their strategies. Other parameters:
$\varepsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$, and $r=0.1$.

##  Appendix C Appendix C: The interaction of Q-learning versus TFT

To examine the interaction between Q-learning and TFT players in details, we
classify all Q-learning players into five categories according to the number
of TFT players in their neighborhood (i.e. 0, 1, 2, 3, 4). Fig. 12 shows the
evolution of cooperation preference in terms of $\overline{\Delta Q}_{s_{j}}$
for Q-learners for all ten states in SM. When there are no TFT players in the
Q-learning players, they prefer defection in any state [Fig. 12(a)]. Fig.
12(b) shows that Q-learning players start to choose cooperation when there are
one TFT player in their neighborhood. Remarkably, when there are two or more
TFT players in their neighborhood, Q-learning players always prefer
cooperation within any state, shown Fig. 12(c-e). This means that in Q-T
interactions, Q-learning players successfully learn to cooperate to maximize
their rewards, and this cooperation propensity increases with the number of
their TFT neighbors.

\begin{overpic}[width=130.08731pt]{SM_DeltaQ_0TFT-eps-converted-
to.pdf}\put(2.5,67.0){(a)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_DeltaQ_1TFT-
eps-converted-
to.pdf}\put(2.5,67.0){(b)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_DeltaQ_2TFT-
eps-converted-
to.pdf}\put(2.5,67.0){(c)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_DeltaQ_3TFT-
eps-converted-
to.pdf}\put(2.5,67.0){(d)}\end{overpic}\begin{overpic}[width=130.08731pt]{SM_DeltaQ_4TFT-
eps-converted-to.pdf}\put(2.5,67.0){(e)}\end{overpic}

Figure 12: The temporal evolution of the average Q-value difference
$\overline{\Delta Q}_{j}$ within all 10 states in SM for $\omega=0.25$. The
cases where the focal player in Q-learning mode are with 0, 1, 2, 3, 4 TFT
players in the neighborhood are computed separately in (a)-(e). The black
dashed line in both plots corresponds to $\overline{\Delta Q}_{s_{j}}=0$ for
reference. Other parameters: $\varepsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$,
$r=0.1$.

## References

  * Maynard Smith and Szathmary (1997) J. Maynard Smith and E. Szathmary,  _The Major Transitions in Evolution_ (OUP Oxford, 1997).
  * Zelenski  _et al._ (2015) J. M. Zelenski, R. L. Dopko, and C. A. Capaldi, Journal of Environmental Psychology 42, 24 (2015).
  * Cheney (2011) D. L. Cheney, Proceedings of the National Academy of Sciences 108, 10902 (2011).
  * Dawkins (2006) R. Dawkins,  _The Selfish Gene_ (Oxford University Press, New York, US, 2006).
  * Rapoport and Chammah (1965) A. Rapoport and A. M. Chammah,  _Prisoner’s dilemma: A study in conflict and cooperation_ , Vol. 165 (University of Michigan press, 1965).
  * Pennisi (2005) E. Pennisi, Science 309, 93 (2005).
  * Hardin (1968) G. Hardin, Science 162, 1243 (1968).
  * Kollock (1998) P. Kollock, Annual Review of Sociology 24, 183 (1998).
  * Milinski  _et al._ (2006) M. Milinski, D. Semmann, H.-J. Krambeck, and J. Marotzke, Proceedings of the National Academy of Sciences 103, 3994 (2006).
  * Poundstone (1993) W. Poundstone,  _Prisoner’s dilemma: John von Neumann, game theory, and the puzzle of the bomb_ (Anchor, 1993).
  * (11) M. Doebeli and C. Hauert, Ecology Letters 8, 748.
  * Axelrod (1980a) R. Axelrod, Journal of Conflict Resolution 24, 3 (1980a).
  * Perc and Szolnoki (2008) M. Perc and A. Szolnoki, Physical Review E 77, 011904 (2008).
  * Axelrod (1980b) R. Axelrod, Journal of Conflict Resolution 24, 379 (1980b).
  * Milinski and Wedekind (1998) M. Milinski and C. Wedekind, Proceedings of the National Academy of Sciences 95, 13755 (1998).
  * Nowak (2006) M. A. Nowak, Science 314, 1560 (2006).
  * Perc  _et al._ (2017) M. Perc, J. J. Jordan, D. G. Rand, Z. Wang, S. Boccaletti, and A. Szolnoki, Physics Reports 687, 1 (2017).
  * Hamilton (1964) W. Hamilton, Journal of Theoretical Biology 7, 17 (1964).
  * Trivers (1971) R. L. Trivers, The Quarterly review of biology 46, 35 (1971).
  * Nowak and Sigmund (1998) M. A. Nowak and K. Sigmund, Nature 393, 573 (1998).
  * Nowak and May (1992) M. A. Nowak and R. M. May, Nature 359, 826 (1992).
  * Wang  _et al._ (2013) Z. Wang, A. Szolnoki, and M. Perc, Scientific Reports 3, 1183 (2013).
  * Liang  _et al._ (2022) R. Liang, Q. Wang, J. Zhang, G. Zheng, L. Ma, and L. Chen, Physical Review E 105, 054302 (2022).
  * Keller (1999) L. Keller,  _Levels of selection in evolution_ , Vol. 66 (Princeton University Press, 1999).
  * Maynard Smith (1964) J. Maynard Smith, Nature 201, 1145 (1964).
  * Sigmund  _et al._ (2001) K. Sigmund, C. Hauert, and M. A. Nowak, Proceedings of the National Academy of Sciences 98, 10757 (2001).
  * Xia  _et al._ (2023) C. Xia, J. Wang, M. Perc, and Z. Wang, Physics of Life Reviews 46, 8 (2023).
  * Santos  _et al._ (2008) F. C. Santos, M. D. Santos, and J. M. Pacheco, Nature 454, 213 (2008).
  * Liang  _et al._ (2021) R. Liang, J. Zhang, G. Zheng, and L. Chen, Physica A 567, 125726 (2021).
  * Watkins and Dayan (1992a) C. Watkins and P. Dayan, Machine Learning 8, 279 (1992a).
  * Sutton and Barto (2018) R. S. Sutton and A. G. Barto,  _Reinforcement Learning: An Introduction_ (MIT press, 2018).
  * Bandura and Walters (1977) A. Bandura and R. H. Walters,  _Social Learning Theory_ (Englewood cliffs Prentice Hall, 1977).
  * Silver  _et al._ (2018) D. Silver, T. Hubert, J. Schrittwieser, I. Antonoglou, M. Lai, A. Guez, M. Lanctot, L. Sifre, D. Kumaran, T. Graepel, T. Lillicrap, K. Simonyan, and D. Hassabis, Science 362, 1140 (2018).
  * Lee  _et al._ (2012) D. Lee, H. Seo, and M. W. Jung, Annual Review of Neuroscience 35, 287 (2012).
  * Rangel  _et al._ (2008) A. Rangel, C. Camerer, and P. R. Montague, Nature Reviews Neuroscience 9, 545 (2008).
  * Zhang  _et al._ (2020) S. Zhang, J. Zhang, L. Chen, and X. Liu, Nonlinear Dynamics 99, 3301 (2020).
  * Wang  _et al._ (2022) L. Wang, D. Jia, L. Zhang, P. Zhu, M. Perc, L. Shi, and Z. Wang, Nonlinear Dynamics 108, 1837 (2022).
  * Wang  _et al._ (2023) L. Wang, L. Fan, L. Zhang, R. Zou, and Z. Wang, New Journal of Physics 25, 073008 (2023).
  * He  _et al._ (2022) Z. He, Y. Geng, C. Du, L. Shi, and Z. Wang, New Journal of Physics 24, 123038 (2022).
  * Ding  _et al._ (2023) Z. Ding, G. Zheng, C. Cai, W. Cai, L. Chen, J. Zhang, and X. Wang, Chaos, Solitons & Fractals 175, 114032 (2023).
  * Geng  _et al._ (2022) Y. Geng, Y. Liu, Y. Lu, C. Shen, and L. Shi, Applied Mathematics and Computation 427, 127182 (2022).
  * Zhang  _et al._ (2024a) J. Zhang, Z. Rong, G. Zheng, J. Zhang, and L. Chen, Journal of Physics: Complexity 5, 025006 (2024a).
  * Zheng  _et al._ (2024) G. Zheng, J. Zhang, J. Zhang, W. Cai, and L. Chen, New Journal of Physics 26, 053041 (2024).
  * Andrecut and Ali (2001) M. Andrecut and M. Ali, Physical Review E 64, 067103 (2001).
  * Zhang  _et al._ (2019) S. Zhang, J. Dong, L. Liu, Z. Huang, L. Huang, and Y. Lai, Physical Review E 99, 032302 (2019).
  * Zheng  _et al._ (2023) G. Zheng, W. Cai, G. Qi, J. Zhang, and L. Chen, arXiv preprint arXiv:2312.14970 (2023).
  * Zhang  _et al._ (2024b) S. Zhang, J. Dong, H. Zhang, Y. Lü, J. Wang, and Z. Huang, Frontiers of Physics 19, 1 (2024b).
  * Tomov  _et al._ (2021) M. S. Tomov, E. Schulz, and S. J. Gershman, Nature Human Behaviour 5, 764 (2021).
  * Shi and Rong (2022) Y. Shi and Z. Rong, IEEE Transactions on Circuits and Systems II: Express Briefs 69, 2463 (2022).
  * Ding  _et al._ (2024) Z. Ding, J. Zhang, G. Zheng, W. Cai, C. Cai, L. Chen, and X. Wang, Chaos, Solitons & Fractals 184, 114971 (2024).
  * Watkins (1989) C. J. C. H. Watkins,  _Learning from delayed rewards (Ph.D. thesis)_ , Ph.D. thesis (1989).
  * Watkins and Dayan (1992b) C. J. Watkins and P. Dayan, Machine Learning 8, 279 (1992b).
  * Van Hasselt  _et al._ (2016) H. Van Hasselt, A. Guez, and D. Silver, in  _Proceedings of the AAAI conference on artificial intelligence_, Vol. 30 (2016).
  * Wunder  _et al._ (2010) M. Wunder, M. Littman, and M. Babes, in  _Proceedings of the 27th International Conference on International Conference on Machine Learning_, ICML’10 (Omnipress, Madison, WI, USA, 2010) p. 1167–1174.
  * Szabó and Fath (2007) G. Szabó and G. Fath, Physics Reports 446, 97 (2007).
  * Axelrod and Hamilton (1981) R. Axelrod and W. D. Hamilton, Science 211, 1390 (1981).
  * Nowak and Sigmund (1992) M. A. Nowak and K. Sigmund, Nature 355, 250 (1992).
  * Nowak and Sigmund (1993) M. Nowak and K. Sigmund, Nature 364, 56 (1993).
  * Traulsen  _et al._ (2010) A. Traulsen, D. Semmann, R. D. Sommerfeld, H.-J. Krambeck, and M. Milinski, Proc. Natl. Acad. Sci. U.S.A. 107, 2962 (2010).
  * Szolnoki and Perc (2015) A. Szolnoki and M. Perc, Journal of The Royal Society Interface 12, 20141299 (2015).
  * Szolnoki and Chen (2018) A. Szolnoki and X. Chen, New Journal of Physics 20, 093008 (2018).
  * Amaral and Javarone (2018) M. A. Amaral and M. A. Javarone, Phys. Rev. E 97, 042305 (2018).
  * Masuda (2012) N. Masuda, Scientific Reports 2, 1 (2012).
  * Zheng  _et al._ (2022) G. Zheng, J. Zhang, R. Liang, L. Ma, and L. Chen, Journal of Physics: Complexity 3, 035004 (2022).
  * Ma  _et al._ (2023) L. Ma, J. Zhang, G. Zheng, R. Liang, and L. Chen, Chaos, Solitons & Fractals 171, 113452 (2023).
  * Han  _et al._ (2022) X. Han, X. Zhao, and H. Xia, Chaos, Solitons & Fractals 164, 112684 (2022).


# Asymptotic Normality of the Coefficients of the Morgan-Voyce Polynomials

Moussa Benoumhani  Department of Mathematics, University of M’Sila, M’Sila,
Algeria<EMAIL_ADDRESS>,  Bernhard Heim  Lehrstuhl A für Mathematik,
RWTH Aachen University, 52056 Aachen, Germany<EMAIL_ADDRESS>and  Markus Neuhauser  Kutaisi International University, 5/7, Youth Avenue,
Kutaisi, 4600 Georgia<EMAIL_ADDRESS>
###### Abstract.

We study arithmetic and asymptotic properties of polynomials provided by
$Q_{n}(x):=x\sum_{k=1}^{n}k\,Q_{n-k}(x)$ with initial value $Q_{0}(x)=1$. The
coefficients satisfy a central limit theorem and a local limit theorem
involving Fibonacci numbers. We apply methods of Berry and Esseen, Harper,
Bender, and Canfield.

###### Key words and phrases:

Central Limit Theorem, Fibonacci Numbers, Local Limit Theorem, Singularity
Analysis.

###### 2010 Mathematics Subject Classification:

Primary 60F05, 11B39; Secondary 05A16

##  1\. Introduction and main results

Let $g$ be a normalized arithmetic function. Let
$G(t):=\sum_{n=1}^{\infty}g(n)\,t^{n}$ be regular at $t=0$. We are interested
in asymptotic properties of the double sequence $a(n,k)$ of coefficients of
the polynomials $p_{n}(x)$ defined by

(1.1) $\sum_{n=0}^{\infty}p_{n}(x)\,t^{n}=\frac{1}{1-x\,G(t)}.$

Let $g(n)=\frac{1}{n!}$. Then the coefficients are asymptotically normal
[Ha67, Be73]. They calculate the number of partitions of an $n$-set having
exactly $k$ labeled blocks. Thus, are equal to $k!\,S_{n,k}$, where $S_{n,k}$
are the Stirling numbers of the second kind ([Be73], section 3 applications).
It is obvious that the case $g(n)=1$ leads to the Central Limit Theorem by de
Moivre–Laplace, since

$a(n,k)=\binom{n-1}{k-1}$

(cf. [HNT20]).

In this paper we study the case $g(n)=n$ in detail and show that the
underlying expected values and variances involve Fibonacci numbers.

We denote the associated polynomials and their coefficients by

$Q_{n}(x)=\sum_{k=1}^{n}A_{n,k}\,x^{k}.$

Let $B_{n}(x)$ denote the Morgan-Voyce polynomials ([Ko01], Chapter 41). Then
$Q_{n+1}(x)=x\,B_{n}(x)$. In 1959, Morgan-Voyce [MV59] discovered importance
of $B_{n}(x)$ in the study of electric ladder networks of resistors. Further,
in 1967 and 1968, several papers appeared by Swamy, Basin, Hoggatt, Jr., and
Bicknell. These polynomials are closely related to Fibonacci polynomials. From
[HNT20], we know that

$A_{n,k}:=\binom{n+k-1}{2k-1},\qquad n\in\mathbb{N},\,0\leq k\leq n.$

The definition (1.1) is equivalent to $Q_{n}(x)=x\sum_{k=1}^{n}kQ_{n-k}(x)$
with initial value $Q_{0}(x)=1$. This hereditary recurrence relation can be
reduced to a three term recurrence relation

(1.2) $Q_{n+2}(x)-(2+x)Q_{n+1}(x)+Q_{n}(x)=0,\qquad n\geq 0,$

with initial values $Q_{1}(x)=x$ and $Q_{2}=\left(x+2\right)x$. The double
sequence $A_{n,k}$ (see Table 1) is recorded in Sloane’s data base as A078812.

$\begin{array}[]{r|rrrrrrrrrrrr}\hline\cr n\backslash k&1&2&3&4&5&6&7&8\\\
\hline\cr\hline\cr 1&1&&&&&&&\\\ 2&2&1&&&&&&\\\ 3&3&4&1&&&&&\\\
4&4&10&6&1&&&&\\\ 5&5&20&21&8&1&&\\\ 6&6&35&56&36&10&1&&\\\
7&7&56&126&120&55&12&1&\\\ 8&8&84&252&330&220&78&14&1\\\ \hline\cr\end{array}$
Table 1.  Coefficients $A_{n,k}$ of $Q_{n}(x)$.

###  1.1. Asymptotic normality by singularity analysis

Our first result, using techniques from singularity analysis (see Canfield
[Ca15], Section 3.6), states that the double sequence $A_{n,k}$ is
asymptotically normal with asymptotic mean $a_{n}=\frac{1}{\sqrt{5}}\,\,n$ and
asymptotic variance $b_{n}^{2}=\frac{2}{5}\,\frac{1}{\sqrt{5}}\,\,n$.

######  Theorem 1.1.

There exist real sequences $(a_{n})_{n}$ and $(b_{n})_{n}$ with $b_{n}>0$,
such that

(1.3) $\lim_{n\to\infty}\mathop{\rm
sup}_{x\in\mathbb{R}}\left|{\frac{1}{\sum_{k=0}^{n}A_{n,k}}\sum_{k\leq
a_{n}+xb_{n}}A_{n,k}-\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{t^{2}}{2}}\,dt}\right|=0.$

The sequences given by $a_{n}:=n\,a$ and $b_{n}^{2}:=n\,b^{2}$ satisfy (1.3)
with

$a:=\frac{1}{\sqrt{5}}\text{ and }b^{2}:=\frac{2}{5}\,\frac{1}{\sqrt{5}}.$

Obviously, the sequences are not unique. In singularity analysis they are
frequently called mean and variance (see Bender [Be73]) to indicate the link
to probability theory and the classical Central Limit Theorem (cf. [Fi11],
section 1.1).

In this paper, we provide sequences $(\mu_{n})_{n}$ and
$\left(\sigma_{n}^{2}\right)_{n}$, the expected value and variance of a
suitable sequence of random variable $X_{n}$. We obtain a refined version of
Theorem 1.1. We apply a method of Harper [Ha67] (see also [Ca15] section 3.4:
the method of negative roots). The values $\mu_{n}$ and $\sigma_{n}^{2}$
involve Fibonacci numbers $F_{n}$ and the golden ratio $\varphi$. They
characterize the peaks and modes of the coefficients of the generating series.
These are real-rooted polynomials, which are related to the Jonquière function
(also called polylogarithm $\mathop{\rm Li}_{-1}(t)$) given by
${G}(t):=\sum_{n=1}^{\infty}n\,t^{n}$. The radius of convergence is $R=1$.

###  1.2. Probabilistic approach

Let $X$ be a random variable. Then $X$ is called normally distributed, if the
probability $P(X\leq x)$ is equal to the normal distribution

$\Phi(x):=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{t^{2}}{2}}\,dt,$

for all real $x$. The distribution $\Phi$ is called normal or Gauss
distribution. Let further, $\phi(x)$ be the underlying density function.

######  Definition.

A sequence $(S_{n})_{n}$ of random variables satisfies a Central Limit
Theorem, if there exist sequences $(a_{n})_{n}$ and $(b_{n})_{n}$ with
$a_{n},b_{n}\in\mathbb{R}$ and $b_{n}>0$, such that the normalized random
variables

$S_{n}^{*}:=\frac{S_{n}-a_{n}}{b_{n}},$

converge in distribution against the normal distribution

$P(S_{n}^{*}\leq x)\stackrel{{\scriptstyle{D}}}{{\longrightarrow}}\Phi(x).$

Consider the sequence of random variables $X_{n}\in\\{0,1,\ldots,n\\}$ defined
by:

(1.4) $P(X_{n}=k):=\frac{A_{n,k}}{\sum_{m=0}^{n}A_{n,m}},\text{ where
}Q_{n}(x)=\sum_{k=0}^{n}A_{n,k}\,x^{k}.$

Let $S_{n}:=\sum_{m=1}^{n}X_{m}$. Then

(1.5) $P(S_{n}^{*}\leq x)=\frac{1}{\sum_{m=0}^{n}A_{n,m}}\sum_{k\leq
a_{n}+x\,b_{n}}A_{n,k},$

leads to a probabilistic interpretation of Theorem 1.1.

######  Theorem 1.2.

Let $S_{n}^{*}$ be the normalized random variables as defined in (1.5) with
$a_{n}$ and $b_{n}^{2}$ provided by the expected value $\mathbb{E}(X_{n})$ and
variance $\mathbb{V}(X_{n})$ of the random variable $X_{n}$, associated with
the double sequence $A_{n,k}$. Then $P(S_{n}^{*}\leq
x)\stackrel{{\scriptstyle{D}}}{{\longrightarrow}}\Phi(x)$. Therefore,

(1.6)
$\lim_{n\to\infty}{\left|\\!\left|\frac{1}{\sum_{m=0}^{n}A_{n,m}}\sum_{k\leq\mathbb{E}(S_{n})+\sqrt{\mathbb{V}(S_{n})}\,x}A_{n,k}-\Phi(x)\right|\\!\right|}_{\mathbb{R}}=0.$

Here ${\left|\\!\left|f(x)\right|\\!\right|}_{M}$ denotes the supremum norm of
$f$ on $M\subset\mathbb{R}$. This theorem is obtained by a method introduced
by Harper [Ha67, Ca15] and a result by Pólya for continuous distribution
functions. We determine the explicit values of the mean and variance of
$X_{n}$. These are expressed in terms of Fibonacci numbers $F_{n}$ and are
related to the peaks and plateaus of the unimodal sequence

$A_{n,0}\leq A_{n,1}\leq\ldots\leq A_{n,m}\geq\ldots\geq A_{n,n}$

by Darroch’s Theorem [Da64]. The sequence is log-concave and therefore,
unimodal, since the polynomials $Q_{n}(x)$ are orthogonal polynomials. Indeed,
$Q_{n}(x)=x\,U_{n-1}(x/2+1)$, where $U_{n}(x)$ are the Chebyshev polynomials
of the second kind [HNT20]. We first determine the normalizing factor, the
average over all coefficients.

######  Lemma 1.3.

The normalizing factor is related to Fibonacci numbers. We have

$\sum_{k=0}^{n}A_{n,k}=Q_{n}(1)=F_{2n}.$

######  Proposition 1.4.

The random variable $X_{n}$ defined in (1.4) has the expected value
$\mu_{n}:=\mathbb{E}(X_{n})$ and variance $\sigma_{n}^{2}:=\mathbb{V}(X_{n})$
given by

$\displaystyle\mu_{n}$ $\displaystyle=$
$\displaystyle\frac{2}{5}\left(\frac{F_{2n+1}}{F_{2n}}-\frac{1}{2}+\frac{1}{n}\right)\,n,$
$\displaystyle\sigma_{n}^{2}$ $\displaystyle=$
$\displaystyle\frac{4}{25}\,\left(\frac{F_{2n+1}}{F_{2n}}-\frac{1}{2}-\frac{n}{F_{2n}^{2}}-\frac{1}{2n}\right)\,n.$

######  Remarks.

a) Note that $\mu_{n}$ is never an integer for $n\geq 2$. If $\mu_{n}$ is an
integer then $2\frac{F_{2n+1}}{F_{2n}}n$ is also an integer. Since $F_{2n+1}$
and $F_{2n}$ are coprime, $\frac{2n}{F_{2n}}$ is also an integer. Since
$F_{2n}>2n$ for $n\geq 3$ and the only remaining case is $n=1$.
b) The reciprocal polynomial of $Q_{n}(x)$ has the coefficients
$\binom{N-k-1}{k}$ with $N=2n$. These sequences had been investigated by Tanny
and Zucker [TZ74, TZ78] and Benoumhani [Be03]. It had been shown that the
sequence is strictly log-concave, the smallest mode had been determined and
the indices, at which a double maximum occurs.

Darroch’s theorem [Da64] (see also Benoumhani [Be96]) leads to:

######  Corollary 1.5.

The modes of $Q_{n}(x)$ are located around $\mu_{n}$. Let $m_{n}$ be a mode,
then

$0<\left|\frac{2}{5}\left(\frac{F_{2n+1}}{F_{2n}}-\frac{1}{2}+\frac{1}{n}\right)\,n-m_{n}\right|<1.$

Using elementary calculations and explicit solutions of the Pell–Fermat
equation leads to:

######  Theorem 1.6.

Let $n\in\mathbb{N}$. The polynomials $Q_{n}(x)$ are unimodal and have at most
two modes. The smallest mode $m$ is uniquely determined by

$\frac{\sqrt{5n^{2}+1}-1}{5}\leq m<\frac{\sqrt{5n^{2}+1}+4}{5}.$

The mode is unique if $5m^{2}+2m\neq n^{2}$. We have two modes $m_{k}$ and
$m_{k}+1$, iff

$\left(\begin{array}[]{c}5m_{k}+1\\\
n_{k}\end{array}\right)=\left(\begin{array}[]{cc}161&360\\\
72&161\end{array}\right)^{k}\left(\begin{array}[]{c}1\\\
0\end{array}\right)\qquad\text{ for }k\geq 1.$

######  Remark.

If $5m^{2}+2m=n^{2}$ has positive integer solutions, there are two modes of
$A_{n,m}$ at $m$ and $m+1$. For example, for $n=n_{1}=72$ and $m=m_{1}=32$, we
have $\binom{103}{63}=61218182743304701891431482520=\binom{104}{65}$. The
sequence of these $n_{k}$ starts with

$72,23184,7465176,2403763488,774004377960,249227005939632,\ldots$

and the corresponding $m_{k}$ are

$32,10368,3338528,1074995712,346145280800,111457705421952,\ldots$

The expected values $\mu_{n}$ and variances $\sigma_{n}^{2}$ converge against
the sequences $a_{n}$ and $b_{n}$, obtained in Theorem 1.1 by methods of
singularity analysis.

######  Corollary 1.7.

Let $X_{n}$ be the random variable defined by (1.4) with expected value
$\mu_{n}$ and variance $\sigma_{n}^{2}$. Let $a$ and $b$ be as in Theorem 1.1.
Then

$\displaystyle\lim_{n\to\infty}\frac{\mu_{n}}{n}$ $\displaystyle=$
$\displaystyle\frac{2}{5}\left(\varphi-\frac{1}{2}\right)=\frac{1}{\sqrt{5}}=a,$
$\displaystyle\lim_{n\to\infty}\frac{\sigma_{n}^{2}}{n}$ $\displaystyle=$
$\displaystyle\frac{4}{25}\left(\varphi-\frac{1}{2}\right)=\frac{2}{5}\,\frac{1}{\sqrt{5}}=b^{2},$

where $\varphi$ denotes the golden ratio.

This follows from Proposition 1.4, since we know that

$\lim_{n\to\infty}\frac{F_{2n+1}}{F_{2n}}=\varphi\qquad(\text{Kepler}).$

It is possible to obtain a rate of convergence of (1.6) by utilizing the
Berry–Esseen theorem (see for example [Ca15], section 3.2), which implies
Theorem 1.2, since the variance satisfies

$\lim_{n\to\infty}\sigma_{n}=\infty.$

######  Theorem 1.8.

Let $S_{n}^{*}$ be the normalized random variables as defined in (1.5),
associated with $A_{n,k}$. Then

$|P(S_{n}^{*}\leq x)-\Phi(x)|\,\leq\,C\,\sigma_{n}^{-1}.$

Here $C>0$ can be chosen as $C=0.7975$.

The Central Limit Theorem, to quote Bender [Be73] and Canfield [Ca75],
provides a certain qualitative feel for the numbers $A_{n,k}$. Further
information is provided by local limit theorems (cf. [Ca15], section 3.7). Let
$A_{n,k}^{*}:=A_{n,k}/Q_{n}(1)$. We consider subsets
$K\subset\\{0,1,\ldots,n\\}$ and $k\in K$ with the asymptotic behavior

${A_{n,k}^{*}}\,\,\sim\,\,\left.\frac{\phi(x)}{\sigma_{n}}\right|_{x=(k-\mu_{n})/\sigma_{n}}.$

More generally, we say the doubly indexed sequence $a(n,k)$ satisfies a local
limit theorem on a set $S$ of real numbers provided

$\mathop{\rm sup}_{x\in
S}\Big{|}\frac{\sigma_{n}\,a(n,\left\lfloor\mu_{n}+x\,\sigma_{n}\right\rfloor)}{\sum_{k}a(n,k)}-\phi(x)\Big{|}\longrightarrow
0,$

where $\mu_{n}$ and $\sigma_{n}$ are the expected values and variances
([Ca15], definition 3.7.1). Further, we quote the following result.

######  Theorem 1.9 (Bender, [Ca15], Theorem 3.7.2).

Suppose that $a(n,k)$ are asymptotically normal and
$\sigma_{n}^{2}\rightarrow\infty$. If for each $n$ the sequence $a(n,k)$ is
unimodal in $k$, then $a(n,k)$ satisfies a local limit theorem on the set
$\\{x\,\,:|x|\geq\varepsilon\\}$, for any $\varepsilon>0$. If for each $n$ the
sequence $a(n,k)$ is log-concave in $k$, the $a(n,k)$ satisfies a local limit
theorem on the set $\mathbb{R}$.

We have that $A(n,k)$ are asymptotically normal by Theorem 1.2 and
$\sigma_{n}^{2}\rightarrow\infty$. Further, that $Q_{n}(x)$ has real roots.
This implies by Theorem 1.9:

######  Corollary 1.10.

$\lim_{n\to\infty}{\left|\\!\left|\sigma_{n}\,{A_{n,\left\lfloor\mu_{n}+x\,\sigma_{n}\right\rfloor}^{*}}-\phi(x)\right|\\!\right|}_{\mathbb{R}}=0.$

######  Corollary 1.11.

Further, we have the asymptotic formula

$A_{n,k}\,\,\sim\,\,F_{2n}\frac{\phi(x)}{\sigma_{n}}\text{ for }n\to\infty,$

where $k=\mu_{n}+x\sigma_{n}$ and $x$ is bounded.

##  2\. Singularity analysis

Bender [Be73] provided a criterion for a double sequence $a\left({n,k}\right)$
to be asymptotically normal. Typically the criterion applies, if the
generating function has only one singularity on the circle of convergence. We
recall the approach by Bender (see also [Ca15], section 3.6 Method 4).

######  Theorem 2.1 (Bender).

Let $f(x,t)=\sum_{k,n}\,a(n,k)\,x^{k}\,t^{n}$, with $a(n,k)\geq 0$. Suppose
there exist

  * (i)

a function $A(s)$ continuous and non-zero near $0$,

  * (ii)

a function $r(s)$ with bounded third derivative near $0$,

  * (iii)

a non-negative integer $m$, and

  * (iv)

positive numbers $\varepsilon$ and $\delta$ such that

$\left(1-\frac{t}{r(s)}\right)^{m}\,f\left(e^{s},t\right)-\frac{A(s)}{1-\frac{t}{r(s)}}$

is analytic and bounded for $|s|<\varepsilon$, $|t|<r(0)+\delta$.

Put $a=-r^{\prime}(0)/r(0)$ and $b^{2}=a^{2}-r^{\prime\prime}(0)/r(0)$. If
$b^{2}\neq 0$, then the numbers $a(n,k)$ are asymptotically normal with
$a_{n}=n\,a$ and $b_{n}^{2}=b^{2}\,n$.

###  2.1. Proof of Theorem 1.1

Let $a(n,k)=A_{n,k}$. Let $G(t)=\sum_{n=1}^{\infty}n\,t^{n}$. This power
series has radius of convergence $1$ and has a pole at $t=1$. The generating
series of $A_{n,k}$ (we refer to [HNT20]) is provided by

$f(x,t)=\frac{1}{1-x\,G(t)}=\sum_{n=0}^{\infty}Q_{n}(x)\,t^{n}.$

We have $G(t)=\frac{t}{(t-1)^{2}}$, thus $G(t)=1$ has two real solutions. The
smallest in absolute value is given by $t_{1}=\frac{3-\sqrt{5}}{2}$.

Suppose $x\neq 0$. If $G\left(t\right)=\frac{1}{x}$, then
$xt=\left(1-t\right)^{2}$ i. e. $t^{2}-\left(2+x\right)t+1=0$. Therefore,
$t=1+\frac{x}{2}-\sqrt{\frac{x^{2}}{4}+x}$ with the principle branch of the
square root for $\left|x-1\right|<1$. Let
$r\left(s\right)=1+\frac{e^{s}}{2}-\sqrt{\frac{e^{2s}}{4}+e^{s}}$ for
$\left|s\right|<\ln 2$. Then the pole of $f\left(x,t\right)$ closest to $0$ is
located at $r\left(s\right)$. The other pole is located at
$1/r\left(s\right)=1+\frac{e^{s}}{2}+\sqrt{\frac{e^{2s}}{4}+e^{s}}$.

Since

$\displaystyle f\left(e^{s},t\right)$ $\displaystyle=$
$\displaystyle\frac{1}{1-e^{s}G\left(t\right)}=\frac{\left(1-t\right)^{2}}{1-\left(2+e^{s}\right)t+t^{2}}=1+\frac{e^{s}t}{\left(1-r\left(s\right)t\right)\left(1-t/r\left(s\right)\right)}$
$\displaystyle=$ $\displaystyle
1+\frac{e^{s}}{2\sqrt{\frac{e^{2s}}{4}+e^{s}}}\left(\frac{1}{1-t/r\left(s\right)}-\frac{1}{1-r\left(s\right)t}\right),$

we obtain

(2.1)
$f\left(e^{s},t\right)-\frac{e^{s}}{2\sqrt{\frac{e^{2s}}{4}+e^{s}}}\frac{1}{1-t/r\left(s\right)}=1-\frac{e^{s}}{2\sqrt{\frac{e^{2s}}{4}+e^{s}}}\frac{1}{1-r\left(s\right)t}.$

For $\left|s\right|<\ln 2$ holds $\left|r\left(s\right)\right|<\frac{1}{2}$.
Therefore, $\left|1/r\left(s\right)\right|>2$ and (2.1) is analytic and
bounded for $\left|t\right|<\frac{3}{2}$.

We obtain
$r^{\prime}\left(s\right)=\frac{e^{s}}{2}-\frac{1}{2}\left(\frac{e^{2s}}{4}+e^{s}\right)^{-1/2}\left(\frac{e^{2s}}{2}+e^{s}\right)$
and

$r^{\prime\prime}\left(s\right)=\frac{e^{s}}{2}+\frac{1}{4}\left(\frac{e^{2s}}{4}+e^{s}\right)^{-3/2}\left(\frac{e^{2s}}{2}+e^{s}\right)^{2}-\frac{1}{2}\left(\frac{e^{2s}}{4}+e^{s}\right)^{-1/2}\left(e^{2s}+e^{s}\right).$

Further,

$\displaystyle r^{\prime\prime\prime}\left(s\right)$ $\displaystyle=$
$\displaystyle\frac{e^{s}}{2}-\frac{3}{8}\left(\frac{e^{2s}}{4}+e^{s}\right)^{-5/2}\left(\frac{e^{2s}}{2}+e^{s}\right)^{3}$
$\displaystyle{}+\frac{1}{2}\left(\frac{e^{2s}}{4}+e^{s}\right)^{-3/2}\left(\frac{e^{2s}}{2}+e^{s}\right)\left(e^{2s}+e^{s}\right)$
$\displaystyle{}+\frac{1}{4}\left(\frac{e^{2s}}{4}+e^{s}\right)^{-3/2}\left(e^{2s}+e^{s}\right)^{2}$
$\displaystyle{}-\frac{1}{2}\left(\frac{e^{2s}}{4}+e^{s}\right)^{-1/2}\left(2e^{2s}+e^{s}\right),$

which shows that the third derivative is bounded near $0$,
$r^{\prime}\left(0\right)=\frac{1}{2}-\frac{1}{2}\frac{2}{\sqrt{5}}\frac{3}{2}=\frac{1}{2}-\frac{3}{2\sqrt{5}}$,
and
$r^{\prime\prime}\left(0\right)=\frac{1}{2}+\frac{1}{4}\frac{8}{5\sqrt{5}}\frac{9}{4}-\frac{2}{\sqrt{5}}=\frac{1}{2}-\frac{11}{10\sqrt{5}}$.
With
$r\left(0\right)=\frac{3}{2}-\frac{\sqrt{5}}{2}=\left(\frac{3}{2\sqrt{5}}-\frac{1}{2}\right)\sqrt{5}$.
This yields $a=\frac{1}{\sqrt{5}}$. As
$1/r\left(0\right)=\frac{3}{2}+\frac{\sqrt{5}}{2}$, we obtain

$b=\frac{1}{5}-\left(\frac{3}{2}+\frac{\sqrt{5}}{2}\right)\left(\frac{1}{2}-\frac{11}{10\sqrt{5}}\right)=\frac{1}{5}-\frac{3}{4}+\frac{11}{20}-\frac{\sqrt{5}}{4}+\frac{33}{20\sqrt{5}}=\frac{2}{5\sqrt{5}}.$

##  3\. The Berry–Esseen theorem

Let $X\in\\{0,1,\ldots,n\\}$ be a random variable with expected value
$\mathbb{E}(X)$ and variance $\mathbb{V}(X)$ for $n\in\mathbb{N}$. Let
$P_{X}(x)=\sum_{k=0}^{n}P(X=k)\,x^{k}$ be the probability generating function.
A straightforward calculation leads to

$\displaystyle\mathbb{E}(X)$ $\displaystyle=$ $\displaystyle P^{\prime}(1),$
$\displaystyle\mathbb{V}(X)$ $\displaystyle=$ $\displaystyle
P^{\prime\prime}(1)-\left(P^{\prime}(1)\right)^{2}+P^{\prime}(1).$

There is a Central Limit Theorem for a sequence of independent, but not
necessarily identically distributed random variables. We work this out in the
setting of a triangular array of Bernoulli random variables. To control the
rate of convergence, we utilize a Berry–Esseen theorem (we refer to [Ca15],
section 3).

######  Theorem 3.1.

Let $X_{n,k}$ for $1\leq k\leq n$ be independent random variables with
expected values $\mu_{n,k}$, variances $\sigma_{n,k}^{2}$ and absolute third
central moment

$\rho_{n,k}=\mathbb{E}(|X_{n,k}-\mu_{n,k}|^{3})<\infty.$

Let $\mu_{n}=\sum_{k=1}^{n}\mu_{n,k}$,
$\sigma_{n}^{2}:=\sum_{k=1}^{n}\sigma_{n,k}^{2}$ and
$S_{n}:=\sum_{k=1}^{n}X_{n,k}$. Let $S_{n}^{*}=(S_{n}-\mu_{n})/\sigma_{n}$.
Then

${\left|\\!\left|P(S_{n}^{*}<x)-\Phi(x)\right|\\!\right|}_{\mathbb{R}}\leq
C\,\,\frac{\sum_{k=1}^{n}\rho_{n,k}}{\sigma_{n}^{3}},$

where $C>0$ is a universal constant. This constant can be chosen as $C=0.7975$
[VB72].

###  3.1. Harper’s method

Let $n\in\mathbb{N}$. Let $P_{n}(x)=\sum_{k=0}^{n}a_{n,k}\,x^{k}$ be a monomic
polynomial of degree $n$, $a_{n,k}\geq 0$ and $P_{n}(1)>0$. Suppose the roots
of $P_{n}(x)$ are real and $P_{n}(x)=\prod_{k=1}^{n}\left(x+r_{k}\right)$.
Harper [Ha67] introduced a triangular array of Bernoulli random variables
$X_{n,j}$ with distribution

$P(X_{n,j}=0):=\frac{r_{j}}{1+r_{j}}\text{ and
}P(X_{n,j}=1):=\frac{1}{1+r_{j}}.$

Let $S_{n}:=\sum_{j=1}^{n}X_{n,j}$. Then
$P(S_{n}=k)=\frac{a_{n,k}}{P_{n}(1)}$.

######  Lemma 3.2.

Let $X_{n,j}$ be given. Then

$\mathbb{E}(X_{n,j})=\frac{r_{j}}{1+r_{j}},\,\,\,\mathbb{V}(X_{n,j})=\frac{r_{j}}{(1+r_{j})^{2}},\,\,\,\mathbb{E}(|X_{n,j}-\mathbb{E}(X_{n,j})|^{3})=\frac{r_{j}(1+r_{j}^{2})}{(1+r_{j})^{4}}.$

This implies that
$\mathbb{E}(|X_{n,j}-\mathbb{E}(X_{n,j})|^{3})<\mathbb{V}(X_{n,j})$.

###### Proof of Theorem 1.8.

Due to Lemma, 3.2 we obtain

$\frac{\sum_{k=1}^{n}\rho_{n,k}}{\sigma_{n}^{3}}\leq\frac{\sum_{k=1}^{n}\sigma_{n,k}^{2}}{\sigma_{n}^{3}}=\frac{1}{\sigma_{n}}.$

Finally, Theorem 3.1 gives the result. ∎

###  3.2. Expected values $\mu_{n}$ and variances $\sigma_{n}^{2}$

The family of polynomials $\\{Q_{n}(x)\\}_{n}$, where

$Q_{n}(x)=\sum_{k=0}^{n}A_{n,k}\,x^{k},$

satisfy the three term recurrence relation (1.2). Then
$u_{n}=Q_{n}\left(1\right)$ satisfies the recurrence relation

(3.1) $u_{n}-3u_{n-1}+u_{n-2}=0.$

###### Proof of Lemma 1.3..

The Fibonacci numbers $F_{n}$ satisfy the recurrence relation
$F_{n}-F_{n-1}-F_{n-2}=0$. Therefore,
$F_{n-1}=F_{n}-F_{n-2}=F_{n+1}-2F_{n-1}+F_{n-3}$, which implies
$F_{n+1}-3F_{n-1}+F_{n-3}=0$. Therefore, the sequences constituted by the
$F_{2n}$ or $F_{2n+1}$ satisfy the same recurrence relation as $u_{n}$. Now
$u_{1}=1=F_{2}$ and $u_{2}=3=F_{4}$, so we obtain $u_{n}=F_{2n}$ for $n\geq
1$. ∎

Next we prove the formulas for $\mu_{n}$ and $\sigma_{n}^{2}$.

###### Proof of Proposition 1.4..

To determine $\mu_{n}$ and $\sigma_{n}$, we follow the strategy offered by the
standard method of solving linear non-homogeneous difference equation (we
refer to [El05], Section 2.4). For $v_{n}=Q_{n}^{\prime}\left(1\right)$ holds
$v_{n}-3v_{n-1}+v_{n-2}=u_{n-1}=F_{2n-2}$. The sequences constituted by
$F_{2n}$ or $F_{2n+1}$, resp., are a linearly independent solution of the
homogeneous difference equation (3.1). Therefore, the solution of the
inhomogeneous difference equation

$v_{n}-3v_{n-1}+v_{n-2}=u_{n-1}.$

is a linear combination of $F_{2n}$, $nF_{2n}$, $F_{2n+1}$, $nF_{2n+1}$. From
the initial conditions $v_{1}=1$, $v_{2}=4$, $v_{3}=14$, and $v_{4}=46$, we
can determine the coefficients and obtain
$v_{n}=\frac{2}{5}nF_{2n+1}+\frac{2}{5}F_{2n}-\frac{1}{5}nF_{2n}$ for $n\geq
1$ and observe that this also holds for $n=0$.

Let $w_{n}=Q_{n}^{\prime\prime}\left(1\right)$. Then
$w_{n}-3w_{n-1}+w_{n-2}=v_{n-1}$ and
$N\left(E\right)=\left(p\left(E\right)\right)^{2}$ is an annihilator of the
right hand side. Since the sequences constituted by $F_{2n}$ and $F_{2n+1}$
are a fundamental system of solutions of (3.1), we obtain a fundamental system
of solutions of the inhomogeneous one by $F_{2n}$, $F_{2n+1}$, $nF_{2n}$,
$nF_{2n+1}$, $n^{2}F_{2n}$, and $n^{2}F_{2n+1}$. The coefficients can be
determined from the values $w_{1}=0$, $w_{2}=2$, $w_{3}=14$, $w_{4}=68$,
$w_{5}=282$, and $w_{6}=1068$ and we obtain
$w_{n}=\left(\frac{1}{5}n^{2}-\frac{1}{25}n-\frac{8}{25}\right)F_{2n}+\frac{2}{25}nF_{2n+1}$.
Recall that

$F_{2n}^{2}+F_{2n}F_{2n+1}-F_{2n+1}^{2}=-1\text{ for }n\geq 0.$

This leads to the explicit formula for $\sigma_{n}^{2}$. ∎

##  4\. Location of the modes: proof of Theorem 1.6

###### Proof.

Let $A_{n,m}=\binom{n+m-1}{2m-1}$. Therefore,
$A_{n,m}-A_{n,m+1}=\binom{n+m-1}{2m-1}\left(1-\frac{\left(n+m\right)\left(n-m\right)}{\left(2m+1\right)2m}\right)$.
For the numerator of the expression in brackets, we obtain
$5m^{2}+2m-n^{2}=5\left(m+\frac{1}{5}\right)^{2}-\frac{1}{5}-n^{2}$.
Therefore, $A_{n,m}>A_{n,m+1}$ for $m>\frac{\sqrt{5n^{2}+1}-1}{5}$ and
$A_{n,m-1}<A_{n,m}$ for $m<\frac{\sqrt{5n^{2}+1}-1}{5}+1$.

Obviously, for $n^{2}=5m^{2}+2m$ we have two modes. This equation is
equivalent to $5n^{2}=\left(5m+1\right)^{2}-1$. With $j=5m+1$, we obtain the
Pell–Fermat equation $j^{2}-5n^{2}=1$. All its non-negative solutions are
$\left(\begin{array}[]{c}j\\\
n\end{array}\right)=\left(\begin{array}[]{cc}9&20\\\
4&9\end{array}\right)^{k^{\prime}}\left(\begin{array}[]{c}1\\\
0\end{array}\right)$ for $k^{\prime}\geq 0$. To be a solution to the original
problem, the integer $j$ must satisfy $j\equiv 1\mod 5$. Exactly even powers
$k^{\prime}$ yield such a solution. ∎

###### Proof of Theorem 1.2.

To apply Theorem 1.8, we need to show that the variance $\sigma_{n}^{2}$
proceeds to infinity. But this follows from Corollary 1.7. ∎

##  5\. Local limit theorem: numerical data

We consider Corollary 1.11 for $k=a_{n}+b_{n}$. Here $a_{n}=n\,a$ and
$b_{n}=\sqrt{n}\,\,b$, where we approximated $\mu_{n}$ by $a_{n}$ and
$\sigma_{n}$ by $b_{n}$. This leads to

$\binom{n+\left\lfloor\frac{n}{\sqrt{5}}\right\rfloor-1}{2\,\left\lfloor\frac{n}{\sqrt{5}}\right\rfloor-1}\,\,\sim\,\,\frac{5^{\frac{3}{4}}}{2\sqrt{\pi}\,\,\sqrt{n}}\,F_{2n}.$

A local Berry–Esseen result would suggest a rate of convergence by
$1/\,\sqrt{n}$. Table 2 gives some evidence.

$\begin{array}[]{rrr|rrr|rrr}\hline\cr n&&&n&&&n&&\\\ \hline\cr\hline\cr
2&0&1.4\cdot 10^{0}&20&1.7\cdot 10^{-1}&8.7\cdot 10^{-1}&200&6.6\cdot
10^{-2}&1.0\cdot 10^{-1}\\\ 3&3.7\cdot 10^{-1}&5.3\cdot 10^{-1}&30&1.6\cdot
10^{-1}&2.4\cdot 10^{-1}&300&5.4\cdot 10^{-2}&3.1\cdot 10^{-2}\\\ 4&1.9\cdot
10^{-1}&1.1\cdot 10^{0}&40&1.3\cdot 10^{-1}&5.8\cdot 10^{-1}&400&4.6\cdot
10^{-2}&1.9\cdot 10^{-1}\\\ 5&3.6\cdot 10^{-1}&3.0\cdot 10^{-1}&50&1.3\cdot
10^{-1}&1.6\cdot 10^{-1}&500&4.2\cdot 10^{-2}&9.9\cdot 10^{-2}\\\ 6&2.4\cdot
10^{-1}&9.0\cdot 10^{-1}&60&1.1\cdot 10^{-1}&4.4\cdot 10^{-1}&600&3.8\cdot
10^{-2}&4.2\cdot 10^{-2}\\\ 7&3.3\cdot 10^{-1}&1.6\cdot 10^{-1}&70&1.1\cdot
10^{-1}&1.1\cdot 10^{-1}&700&3.5\cdot 10^{-2}&1.0\cdot 10^{-2}\\\ 8&2.5\cdot
10^{-1}&6.6\cdot 10^{-1}&80&1.0\cdot 10^{-1}&3.4\cdot 10^{-1}&800&3.3\cdot
10^{-2}&1.1\cdot 10^{-1}\\\ 9&3.0\cdot 10^{-1}&7.5\cdot 10^{-2}&90&9.8\cdot
10^{-2}&8.1\cdot 10^{-2}&900&3.1\cdot 10^{-2}&5.6\cdot 10^{-2}\\\ 10&2.5\cdot
10^{-1}&4.7\cdot 10^{-1}&100&9.1\cdot 10^{-2}&2.8\cdot 10^{-1}&1000&2.9\cdot
10^{-2}&2.1\cdot 10^{-2}\\\ \hline\cr\end{array}$ Table 2.  Values of
$\frac{1}{F_{2n}}\binom{n+\left\lfloor\frac{n}{\sqrt{5}}\right\rfloor-1}{2\left\lfloor\frac{n}{\sqrt{5}}\right\rfloor-1}$
(second columns) and
$\left|\frac{2\sqrt{\pi}\sqrt{n}}{5^{\frac{3}{4}}F_{2n}}\binom{n+\left\lfloor\frac{n}{\sqrt{5}}\right\rfloor-1}{2\left\lfloor\frac{n}{\sqrt{5}}\right\rfloor-1}-1\right|\sqrt{n}$
(third columns).

## References

  * [Be73] E. Bender: _Central and local limit theorems applied to asymptotic enumeration._ J. Comb. Theory (A) 15 (1973), 91–111.
  * [Be96] M. Benoumhani: _Sur une propriété des polynômes à racines réelles négatives._ J. Math. Pures Appl. IX. Sér. 75, Number 2 (1996), 85–110.
  * [Be03] M. Benoumhani: _A sequence of binomial coefficients related to Lucas and Fibonacci numbers._ Journal of Integer Seq. 6 (2003), Article 03.2.1.
  * [Ca75] E. R. Canfield: _Asymptotic normality in binomial type enumeration._ Ph. D. Thesis, University of California, San Diego, 1975.
  * [Ca15] E. R. Canfield: _Asymptotic normality in enumeration._ In: Miklós Bóna (ed.) Handbook of Enumeration, CRC Press. Discrete Mathematics and its Applications (2015), 255–280.
  * [Da64] J. N. Darroch: _On the distribution of the number of successes in independent trials._ Ann. Math. Statist. 35, Number 3 (1964), 1317–1321.
  * [El05] S. Elaydi: _Introduction to Difference Equations._ Undergraduate Texts in Mathematics. Springer, New York, third edition, 2005.
  * [Fi11] H. Fischer: _A History of the Central Limit Theorem._ From Classical to Modern Probability Theory. Springer, New York, NY (2011).
  * [Ha67] L. Harper: _Stirling behaviour is asymptotically normal._ Ann. Math. Stat. 38 (1967), 410–414.
  * [HNT20] B. Heim, M. Neuhauser, R. Tröger: _Zeros of recursively defined polynomials._ J. Difference Equ. Appl. 26, Number 4 (2020), 510–531.
  * [Ko01] T. Koshy: _Fibonacci and Lucas Numbers with Applications._ Pure and Applied Mathematics. A Wiley–Interscience Series of Texts, Monographs, and Tracts. 2001.
  * [MV59] A. M. Morgan-Voyce: _Ladder network analysis using Fibonacci numbers._ IRE Trans. on Circuit Theory, CT-6 (Sept. 1959), 321–322.
  * [TZ74] S. Tanny, M. Zucker: _On a unimodal sequence of binomial coefficients._ Discrete Math. 9 (1974), 79–89.
  * [TZ78] S. Tanny, M. Zucker: _Analytic methods applied to a sequence of binomial coefficients._ Discrete Math. 24 (1978), 299–310.
  * [VB72] P. van Beek: _An application of Fourier methods to the problem of sharpening the Berry–Esseen inequality._ Z. Wahrsch. Verw. Gebiete 23 (1972), 187–196.