# Repulsive and attractive Casimir interactions in liquids

Anh D. Phan  University of South Florida, Tampa, USA<EMAIL_ADDRESS>N.
A. Viet  Institute of Physics, Hanoi, Vietnam

###### Abstract

The Casimir interactions in the solid-liquid-solid systems as a function of
separation distance have been studied by the Lifshitz theory. The dielectric
permittivity functions for a wide range of materials are described by Drude,
Drude-Lorentz and oscillator models. We find that the Casimir forces between
gold and silica or MgO materials are both the repulsive and attractive. We
also find the stable forms for the systems. Our studies would provide a good
guidance for the future experimental studies on the dispersion interactions.

###### pacs:

Valid PACS appear here

##  I Introduction

The dispersion interactions, the Casimir force, between neutral objects have
brought attraction for many years. There are a lot of factors affecting on the
value of force, such as geometry and material properties. Each of them gives
rise to hot subjects of ongoing investigation. Some experiments have examined
the influence of the dielectric properties of objects on the Casimir force 1 ;
2 ; 3 ; 6 . A number of settings used to study the interaction in terms of
theory are ideal metals, real metals and semiconductors 3 ; 4 ; 6 ,
metamaterials, and two objects placed in liquids 1 ; 2 ; 5 . These studies
have significantly advanced our understanding of the subtle effect of geometry
and material on the Casimir-Lifshitz interactions, especially for designing
nanodevices and nanotechnologies.

In the Lifshitz theory, the dispersion interactions primarily depend on
dielectric permittivity functions of materials. Changing dielectric function
alters the Casimir interactions. There are some ways to modify dielectric
functions, including illuminating a light on the silicon 7 ; 8 , which make
drifting carriers on semiconductor materials. In principle, there are some
models to describe dielectric response functions of real materials, for
example, plasma and Drude models for metals 6 ; 9 ; 11 , Drude-Lorentz and
oscillator models for liquids 2 ; 11 , oxides and others 10 ; 11 ; 12 . Based
on these models, the Casimir forces were obtained by numerical integrations
and series expansion methods 19 .

It has been theoretically shown that the attractive Casimir interaction always
occurs between two (non-magnetic) dielectric bodies related by reflection.
Therefore, the repulsive force is a striking feature creating inspiration for
scientists to make accurate measurements of nano electromechanical machines
where the repulsive force plays an important role and might resolve the
stiction problems. The repulsive Casimir forces can be observed in systems
which have the presence of liquids 2 , metamaterials and metallic geometries
16 . Recent experiments have pointed out that there a repulsive force exists
between a gold sphere and a silica plate, separated by bromobenzen 2 . As a
matter of fact, the repulsive Casimir forces between solids arise when the
dielectric of material surfaces 1 and 2 and an intervening liquid obey the
relation
$\varepsilon_{1}(i\xi)>\varepsilon_{liquid}(i\xi)>\varepsilon_{2}(i\xi)$ over
a wide imaginary frequency range $\xi$.

A previous theoretical 14  study has noticed that it is difficult to establish
an equilibrium configuration of sytems in a vacuum medium. In the reference 15
, the authors showed that they were able to form some stable configurations of
Teflon-Si and Silica-Si immersed in ethanol. The equilibrium is explicitly
explained by dispersion properties. In the present work, our theoretical
studies have shown that the equilibra can be obtained by placing Au-MgO,
Silica-MgO and Au-Silica sytems in bromobenzen.

In this paper, the Casimir-Lifshitz forces between material plate systems made
in oxides and metals immersed in bromobenzen are calculated. The combination
between these results and the proximity force approximation (PFA) method allow
us to compute the Casimir interactions in different configurations. We find
that the magnitude of the Casimir force between two dielectric bodies depends
on the configuration and distance between two bodies. The shape usually used
in experiments is a combination of a sphere and a plate because one can avoid
the problem of alignment and easily control the distance between them. The
energy interactions between a plate-plate system per unit area can be obtained
by using the relationship between the Casimir energy of two plannar objects
and the dispersion force of a sphere-plate system.

The rest of the paper is organized as follows: In Sec. II the theoretical
formulations of Casimir-Lifshitz force interaction are introduced. In section
III, the numerical results for the Casimir force between two bodies are
presented. Important conclusions and discussions are finally given in section
IV.

##  II Lifshitz theory for force calculations

For the force calculations, we used Lifshitz theory without considering effect
of temperature. The separations used here were less than 1 $\mu$m, therefore
thermal corrections at $T=300K$ are not significant. As previously noted in 2
; 3 ; 20 ; 21 , the Lifshitz formula at zero temperature for the Casimir force
acting between between two parallel flat bodies per unit area, separated by a
distance $d$ are given by

$\displaystyle F(d)$ $\displaystyle=$
$\displaystyle-\dfrac{\hbar}{2\pi^{2}}\int_{0}^{\infty}qk_{\perp}dk_{\perp}\int_{0}^{\infty}d\xi$
(1)
$\displaystyle\times\left(\dfrac{r_{TM}^{(1)}r_{TM}^{(2)}}{e^{2qd}-r_{TM}^{(1)}r_{TM}^{(2)}}+\dfrac{r_{TE}^{(1)}r_{TE}^{(2)}}{e^{2qd}-r_{TE}^{(1)}r_{TE}^{(2)}}\right).$

Here the reflection coefficients $r_{TM,TE}^{(1)}$ and $r_{TM,TE}^{(2)}$ for
two independent polarizations of the electromagnetic field (transverse
magnetic and transverse electric fields) are

$\displaystyle
r_{TM}^{(p)}=r_{TM}^{(p)}\left({\xi,k_{\bot}}\right)=\frac{{\varepsilon^{(p)}\left({i\xi}\right)q-\varepsilon^{(2)}\left({i\xi}\right)k^{(p)}}}{{\varepsilon^{(p)}\left({i\xi}\right)q+\varepsilon^{(2)}\left({i\xi}\right)k^{(p)}}},$
(2) $\displaystyle
r_{TE}^{(p)}=r_{TE}^{(p)}\left({\xi,k_{\bot}}\right)=\frac{{\mu^{(2)}(i\xi)k^{(p)}-\mu^{(p)}(i\xi)q}}{{\mu^{(2)}(i\xi)k^{(p)}+\mu^{(p)}(i\xi)q}},$
(3)

where

$\displaystyle
q=\sqrt{k_{\bot}^{2}+\varepsilon^{(2)}\left({i\xi}\right)\mu^{(2)}\left({i\xi}\right)\frac{{\xi^{2}}}{{c^{2}}}},$
(4) $\displaystyle
k^{(p)}=\sqrt{k_{\bot}^{2}+\varepsilon^{(p)}\left({i\xi}\right)\mu^{(p)}\left({i\xi}\right)\frac{{\xi^{2}}}{{c^{2}}}}.$
(5)

in which $\varepsilon^{(p)}\left({i\xi}\right)$ and
$\mu^{(p)}\left({i\xi}\right)$ are the dielectric permittivity and the
magnetic permeability of the first body (p = 1) and the second body (p=3),
respectively. $\varepsilon^{(2)}\left(\omega\right)$ and
$\mu^{(2)}\left({i\xi}\right)$ are the dielectric function and the
permeability of a liquid filled between two bodies. Here, medium ‘2’ selected
is a bromobenzen so $\mu^{(2)}\left({i\xi}\right)=1$. Moreover, in this paper,
the non-magnetic materials used such as germanium, gold and oxides have also
$\mu^{(p)}\left({i\xi}\right)=1$. $k_{\bot}$ magnitude of the wave vector
component perpendicular on the plate, is frequency variable along the
imaginary axis ($\omega=i\xi$).

We recall that Lifshitz formula, routinely used to interpret current
experiments, express the Casimir force between two parallel plates as an
integral over imaginary frequencies $i\xi$ of a quantity involving the
dielectric permittivities of the plates $\omega=i\xi$. It is important to note
that, in principle, recourse to imaginary frequencies is not mandatory because
it is possible to rewrite Lifshitz formula in a mathematically equivalent
form, involving an integral over the real frequency axis. In this case,
however, the integrand becomes a rapidly oscillating function of the
frequency, which hampers any possibility of numerical evaluation. Another
remarkable point is that occurrence of imaginary frequencies in the expression
of the Casimir force is a general feature of all recent formalisms hence
extending Lifshitz theory to non-planar geometries 17 ; 18 . The problem is
that the electric permittivity $\varepsilon(i\xi)$ at imaginary frequencies
cannot be measured directly by any experiment. The only way to determine it by
means of dispersion relations, which allow the expression of
$\varepsilon(i\xi)$ in terms of the observable real-frequency electric
permittivity $\varepsilon(i\xi)$. In the standard works on the Casimir effect,
$\varepsilon(i\xi)$ is expressed with the Kramers-Kronig relation in terms of
an integral of a quantity involving the imaginary part of the electric
permittivity 13

$\displaystyle\varepsilon(i\xi)=1+\frac{2}{\pi}\int\limits_{0}^{\infty}{d\omega\frac{{\omega{\mathop{\rm
Im}\nolimits}\varepsilon(\omega)}}{{\omega^{2}+\xi^{2}}}},$ (6)

where ${\mathop{\rm Im}\nolimits}\varepsilon(\omega)$ is calculated using the
tabulated optical data for the complex index of refraction.

The well-known dielectric function described for gold is the Drude model 13

$\displaystyle\varepsilon(i\xi)=1+\frac{{\omega_{p}^{2}}}{{\xi(\xi+\gamma)}},$
(7)

where $\omega_{p}=9.0$ eV, $\gamma=0.035$ eV are the plasma frequency and the
relaxation parameter of Au, respectively.

The imaginary part of the resulting dielectric function at 6 and 295 K of pure
MgO are shown in 12 . The optical features have been fitted to a classical
oscillator model using the complex dielectric function

$\displaystyle\varepsilon(\omega)=\varepsilon_{\infty}+\sum\limits_{j}{\frac{{\omega_{p,j}^{2}}}{{\omega_{TO,j}^{2}-\omega^{2}-i2\omega\gamma_{i}}}},$
(8)

where $\varepsilon_{\infty}$ is a high-frequency contribution, and
$\omega_{TO,j},{2\gamma_{i}}$ and $\omega_{p,j}$ are the frequency, full width
and effective plasma frequency of the $j$th vibration. The values of these
parameters can be found in 12 . Of course with such simple model for the
permittivity of MgO, there is no need to use dispersion relations to obtain
the expression of $\varepsilon(i\xi)$, for this can be simply done by the
substitution $\omega\to i\xi$ in the r.h.s of Eq.(8) 24

$\displaystyle\varepsilon(i\xi)=\varepsilon_{\infty}+\sum\limits_{j}{\frac{{\omega_{p,j}^{2}}}{{\omega_{TO,j}^{2}+\xi^{2}+2\xi\gamma_{j}}}}.$
(9)

In the case of bromobenzen and silica, it has recently been used for
measurement of repulsive forces between gold and silica surfaces. Extremely
weak repulsion was measured, indicating that the dielectric functions of
bromobenzen and silica are very similar in magnitude. In fact, oscillator
models are constructed to represent the dielectric function at imaginary
frequencies. The form of the oscillator model is given by

$\displaystyle\varepsilon(i\xi)=1+\sum_{i}{\frac{{C_{i}}}{{1+\xi^{2}/\omega_{i}^{2}}}},$
(10)

where the coefficients $C_{i}$ are the oscillator’s strengths corresponding to
(resonance) frequencies $\omega_{i}$ 1 ; 2 ; 25 . The dielectric data was
fitted in a wide frequency range 2 . They are more accurate in comparison with
other simple oscillator models. Moreover, many older references used limited
dielectric data, so the oscillator models with second or third order may lead
to the difference in Casimir force calculations. The parameters we used in the
present paper for bromobenzen and also silica come from 2 .

##  III Numerical results and discussions

The Casimir attractive force usually occurs in experiments and theoretical
calculations. When bromobenzen is filled in the gap between two bodies, the
Casimir force is attractive if the dielectric functions do not satisfy one
condition
$\varepsilon_{1}(i\xi)>\varepsilon_{liquid}(i\xi)>\varepsilon_{2}(i\xi)$ for
all frequencies $\xi$. Therefore, by describing Fig. 1 as the dielectric
response function as a function of the frequency gives us some predictions of
repulsive and attractive forces.

Figure 1: (Color online) The dielectric function of various materials plotted
at imaginary frequencies $\xi$.

This graph shows that
$\varepsilon_{Au}(i\xi)>\varepsilon_{MgO}(i\xi)>\varepsilon_{liquid}(i\xi)$
and
$\varepsilon_{MgO}(i\xi)>\varepsilon_{Au}(i\xi)>\varepsilon_{liquid}(i\xi)$ at
$\xi<6.5$ eV, thus the interactions between Au and MgO body immersed in
bromobenzen liquid and in vacuum are attractive in this range. In the range of
$\xi>6.5$ eV,
$\varepsilon_{MgO}(i\xi)>\varepsilon_{liquid}(i\xi)>\varepsilon_{Au}(i\xi)$,
it causes the repulsive interaction. Similarly, in the gold-bromobenzen-silica
system, at extremely small frequencies , the forces are attractive. In the
larger frequency region, the Casimir forces are repulsive. Besides, the
similar explainations are applied to understand the interaction in the MgO-
bromobenzen-Au system. The numerical calculations of the normalized Casimir
force are provided in Fig. 2

Figure 2: (Color online) Relative Casimir force between two semi-infinite
plates normalized by the perfect metal force $F_{o}(d)=-\pi^{2}\hbar
c/240d^{4}$. The liquid used in this calculation is bromobenzen.

In the MgO-bromobenzen-Silica system, it can be clearly seen that there are
two positions in each curve where the Casimir force is equal to zero. The
first points are corresponding to unstable equilibria $d_{us}^{(1)}\approx 13$
nm because the interaction force changes from the attactive force to the
repulsive force, the second point $d_{s}^{(1)}\approx 110$ nm is a stable
position. There is only one position in the Au-bromobenzen-Silica system and
the Au-bromobenzen-MgO sytem, The interaction forces disappear at
$d_{s}^{(2)}\approx 275$ nm and $d_{s}^{(3)}\approx 5.5$ nm, stable position
of each sytem, respectively.

Figure 3: (Color online) Schematic picture of the setting considered in our
calculations. A sphere is located in bromobenzen at a distance $d$ away from a
material plate.

In order to consider the Casimir interactions between a spherical body and a
plate at a distance of close approach $d$ at a temperature T = 300 K, it is
very useful to utilize the PFA method to calculate. Experimental results for
the Casimir force in the plane-sphere geometry are usually compared with PFA-
based theoretical models. The spherical surface is assumed to be nearly flat
over the scale of $d$. Although the Casimir force is not additive, PFA is
often expected to provide an accurate description when $R\gg d$. Here, the
radius of Au sphere that is used in configurations is $R=40$ $\mu$m in order
to calculate Casimir interaction by the proximity force approximation (PFA)
method because the ratio of $d$ to $R$ is small enough to PFA results becoming
enormously accurate. It can be described by Fig. 3. In this approach, the
surfaces of the bodies are treated as a superposition of infinitesimal
parallel plates 22 .

$\displaystyle
F_{sp}^{PFA}(d)=\int\limits_{0}^{R}{F_{pp}(d+R-\sqrt{R^{2}-r^{2}})}2\pi rdr.$
(11)

here ${F_{pp}}$ is the Casimir force for two parallel plates of unit area.

When using the PFA method, one important point is that the interactions
between a gold sphere or a magnesium oxide sphere and a silica plate are equal
to the interactions between a magnesium oxide plate, which has the same
radius, and a gold plate or a silica plate. There is no difference in
calculations and results as well because the PFA method does not consider a
structure of bodies when their shape is modified or is spherical or
cylindrical shape. The equivalent situations occur in other materials. The
resulting Casimir forces are shown in Fig. 4.

Figure 4: (Color online) The Casimir forces of various sphere-bromobenzen-
plate systems are estimated as a function of separation, described in the text
with the spherical redius $R=40$ $\mu$m.

In the reference 2 , the authors experimentally measured and theoretically
calculated the Casimir interaction between a gold sphere and a silica plate
immersed in bromobenzen in the range from $20$ nm to $60$ nm. Our results in
this range are the same for this range. But when we extend the considered
range of distance, the attractive-repulsive transition occurs at approximately
$190$ nm. This position makes this system stable. Another consequence of Fig.
4 demonstrates that stable position of the Au-bromobenzen-MgO sytem moves to
$3.5$ nm to balance between the attractive and repulsive forces. It can be
explained that increasing the separation distance of infinitesimal parallel
plates causes the fast reduction of the dispersion interaction. At the same
minimal separation distance $d$, the attractive force acting on a sphere is
less than that of a plate in the same effective area. Finally, in the system
of a MgO sphere and a silica plate embeded in bromobenzen, there is only one
presence of non-interaction posion at nearly $10$ nm. It is an unstable
position.

Figure 5: (Color online) The Casimir energy is calculated as a function of
separation for different materials .

In addition, PFA formula and Eq.(11) can allow us to estimate the Casimir
energy per unit area between two plate bodies illustrated in Fig. 5. The
Casimir energy is approximated by 22

$\displaystyle F_{sp}^{PFA}(d)=2\pi RE(d),$ (12)

where $E(d)$ is the Casimir energy per a unit area for planar bodies.

We have also applied the PFA method to calculate the Casimir force in sphere-
sphere systems, we continue to calculate by the PFA method. The formula for
this calculation is given

$\displaystyle F_{ss}^{PFA}(d)=2\pi\int\limits_{0}^{R_{2}}rdr\times$
$\displaystyle
F_{pp}(d+R_{1}-\sqrt{R_{1}^{2}-r^{2}}+R_{2}-\sqrt{R_{2}^{2}-r^{2}}),$ (13)

where the radii of two spherical objects are ${R_{1}}$ and ${R_{2}}$,
respectively. It is assumed that $R_{2}<R_{1}$. In this study, we consider
$R_{1}=40$ $\mu$m and the case of $R_{2}=R_{1}$, $R_{1}=2R_{2}$ and
$R_{1}=2R_{2}$.

Here, having calculated $F_{ss}^{PFA}(d)$ in a sphere-sphere system using
Eq.(13) and $F_{sp}^{PFA}(d)$ in a sphere-plate system using Eq.(11). These
results obtained show that, when increasing $d$, the ratio
$F_{ss}^{PFA}(d)/F_{sp}^{PFA}(d)$ does not depend on the distance $d$. It is a
constant with its magnitude as a function of the radius of two spheres,
$F_{ss}^{PFA}(d)/F_{sp}^{PFA}(d)=1/2$ when $R_{1}=R_{2}$ and
$F_{ss}^{PFA}(d)/F_{sp}^{PFA}(d)=1/3$ when $R_{1}=2R_{2}$. Generalizing this
ratio, if $R_{1}=nR_{2}$, we have $F_{ss}^{PFA}(d)/F_{sp}^{PFA}(d)=1/(n+1)$.
This character is likely to be explained by the results in 22 . When the
second sphere is extremely small in comparison with the first one, the
interaction force goes to zero. In this case, the Lifshitz formula used to
calculate the Casimir force should be transformed to the Casimir-Polder
formula describing the interaction between an atom and a microscopic object.
Moreover, the PFA method is not accurate because this approach is useful if
the size of objects is much larger than the separation between them. On the
other hand, we have

$\displaystyle F_{ss}^{PFA}(d)=\frac{1}{n+1}F_{sp}^{PFA},$ (14)

where $R_{1}=nR_{2}$. If $F_{sp}^{PFA}=0$, $F_{ss}^{PFA}$ must be zero.
Therefore, the unstable and stable positions are constant and unchanged when a
radius of a second sphere varies.

One demonstrated that the Casimir force between two objects embedded in
liquids can be derived from the well-known Lifshitz formula at least if the
object is not made of nonabsorbing materials 5 . That explains why the
Lifshitz expressions is used in order to calculate the Casimir force and
compare with experimental data. Nothing changes in the dielectric functions of
bodies immersed in liquids. On the other hand, several experiments verified
that when metals are placed in liquids, there is a variation of Drude
parameters in the metal 23 . The discrepacy of the interaction between “dry”
and “wet” can reach to 15 $\%$ in this case. But they measured the Casimir
force between two metal plates and got the error. Besides, maybe the
dielectric of liquids and low index materials play much more important role in
Casimir force. In the reference 2 , the change of Drude parameters are not
taken into account but the theoretical calculations are close to the
experimental data curves when we have liquids and the low index materials.

##  IV Conclusions

In this work, we have extended the Lifshitz theory to calculate the Casimir
force. Liquid, silica and magnesium oxide are represented by oscillator
models. Although further studies are required to determine the repulsive
Casimir force accurately, our results show that MgO and silica is a good
candidate for the demonstration of quantum levitation. The contribution of
bromobenzen is important because it is an important factor making the purely
repulsive force or the repulsive-attractive transition. After calculating the
Casimir force between two bodies per unit area and associating proximity force
approximation method, it is easy to compute the interaction of different
material plates with a material sphere. Based on the formula of the Casimir
force between a sphere and a plate, it is convenient to estimate the free
energy interaction of bodies. The result is a prediction for further
experimental studies.

###### Acknowledgements.

This work was supported by the Nafosted Grant No. 103.02.57.09. We thank Prof.
Lilia M. Woods and Prof. P. J. van Zwol for helpful discussions and comments.

## References

  * (1) J. N. Munday, F. Capasso, and V. A. Parsegian, Nature 457, 170 (2009).
  * (2) P. J. van Zwol, and G. Palasantzas, Phys. Rev. A 81, 062502 (2010).
  * (3) F. Chen, U. Mohideen, G. L. Klimchitskaya, and V. M. Mostepanenko, Phys. Rev. A 74, 022103 (2006).
  * (4) I. Pirozhenko, and A. Lambrecht, Phys. Rev. A 77, 013811 (2008).
  * (5) Christian Raabe, and Drik-Gunnar Welsch, Phys. Rev. A 71, 013811 (2005).
  * (6) V.B. Svetovoy, P.J. van Zwol, G. Palasantzas, and J. Th. M. De Hosson, Phys. Rev. B 77, 035439 (2008).
  * (7) Norio INUI, J. Phys. Soc. Jpn. 75, 024004 (2006).
  * (8) F. Chen, G.L. Klimchitskaya, V.M. Mostepanenko, and U. Mohideen, Phys. Rev. B 76, 035338 (2007).
  * (9) R. Esquivel-Sirvent, Phys. Rev. A 77, 042107 (2008).
  * (10) R. Castillo-Garza,C.-C. Chang, D. Jimenez, G. L. Klimchitskaya, V. M. Mostepanenko, and U. Mohideen, Phys. Rev. A 75, 062114 (2007).
  * (11) I. Pirozhenko, and A. Lambrecht, Phys. Rev. A 77, 013811 (2008).
  * (12) Tao Sun, Philip B. Allen, David G. Stahnke, Steven D. Jacobsen, and Christopher C. Homes, Phys. Rev. B 77, 134303 (2008).
  * (13) Giuseppe Bimonte, Phys. Rev. A 83, 042109 (2011).
  * (14) S. J. Rahi, M. Kardar, and T. Emig, Phys. Rev. Lett. 105, 070404 (2010).
  * (15) Alejandro W. Rodriguez, Alexander P. McCauley, David Woolf, F. Capasso, J. D. Joannopoulos, and Steven G. Johnson, Phys. Rev. Lett. 104, 160402 (2010).
  * (16) R. Zhao, J. Zhou1, Th. Koschny, E. N. Economou, and C. M. Soukoulis, Phys. Rev. Lett. 103, 103602 (2009).
  * (17) T. Emig, N. Graham, R. L. Jaffe, and M. Kardar, Phys. Rev. Lett. 99, 170403 (2007); Phys. Rev. D 77, 025005, (2008).
  * (18) O. Kenneth, and I. Klich, Phys. Rev. B 78, 014103 (2008).
  * (19) R. S. Decca, E. Fischbach, G. L. Klimchitskaya, D. E. Krause, D. L ̵́ pez, and V. M. Mostepanenko, Phys. Rev. A 84, 042502 (2011).
  * (20) J. N. Munday, F. Capasso, V. A. Parsegian, and S. M. Bezrukov, Phys. Rev. A 78, 032109 (2008).
  * (21) P. J. van Zwol, G. Palasantzas, and J. Th. M. De Hosson, Phys. Rev. B 79, 195428 (2009).
  * (22) Bo. E. Sernelius, and C. E. Romn-Velzquez, Phys. Rev. A 78, 032111 (2008).
  * (23) R. Esquivel-Sirvent, J. Chem. Phys 132, 194707 (2010).
  * (24) Giuseppe Bimonte, Phys. Rev. A 81, 062501 (2010).
  * (25) I. Pirozhenko, and A. Lambrecht, Phys. Rev. A 80, 042510 (2009).


# Evolution of cooperation in the public goods game with Q-learning

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  Shengfeng Deng  School of
Physics and Information Technology, Shaanxi Normal University, Xi’an 710061,
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

Recent paradigm shifts from imitation learning to reinforcement learning (RL)
is shown to be productive in understanding human behaviors. In the RL
paradigm, individuals search for optimal strategies through interaction with
the environment to make decisions. This implies that gathering, processing,
and utilizing information from their surroundings are crucial. However,
existing studies typically study pairwise games such as the prisoners’ dilemma
and employ a self-regarding setup, where individuals play against one opponent
based solely on their own strategies, neglecting the environmental
information. In this work, we investigate the evolution of cooperation with
the multiplayer game – the public goods game using the Q-learning algorithm by
leveraging the environmental information. Specifically, the decision-making of
players is based upon the cooperation information in their neighborhood. Our
results show that cooperation is more likely to emerge compared to the case of
imitation learning by using Fermi rule. Of particular interest is the
observation of an anomalous non-monotonic dependence which is revealed when
voluntary participation is further introduced. The analysis of the Q-table
explains the mechanisms behind the cooperation evolution. Our findings
indicate the fundamental role of environment information in the RL paradigm to
understand the evolution of cooperation, and human behaviors in general.

##  I 1\. Introduction

Cooperation is crucial for the sustainable development of humans from tackling
climate warming Milinski  _et al._ (2006), fighting global infectious
diseases, and dealing with the melting of glaciers causing the accelerated
extinction of species, where the web of lives that sustains us is being worn
and torn apart. The solution to these issues depends on human cooperation
Kollock (1998). Deciphering the cooperation mechanism is thus a question of
paramount importance, and is crucial for addressing the challenges of these
public tragedies Hardin (1968).

To achieve this aim, researchers investigate the emergence of cooperation by
adopting the framework of evolutionary game theory Nowak and Sigmund (2004);
Nowak (2006), where the prisoners’ dilemma (PD) game Axelrod and Hamilton
(1981) serves as a typical model to study the evolution of cooperation
behavior in pairwise interactions Doebeli and Hauert (2005). In the PD, there
is a potential conflict for the two engaged players between the individual and
their collective interests Rapoport and Chammah (1965); Axelrod and Hamilton
(1981). Importantly, the public goods game (PGG) Hauert and Szabó (2003); Perc
_et al._ (2013); Tanimoto (2013) extends the scenario to multilateral
interactions, where the number of participants could be arbitrary. In PGG,
cooperation contributes to the common goods, but the free-riding behaviors
ruin the benefit of the collective and ultimately lead to the tragedy of
commons if no countermeasure is engaged Dawes (1980); Fischbacher  _et al._
(2001); Wang  _et al._ (2009). Nowadays, the PGG is widely adopted to discuss
many issues we are facing, such as resource management, environmental
protection, and the provision of public services Gintis  _et al._ (2003).

After endeavors for several decades, various mechanisms for the emergence of
cooperation in the PGG have been proposed, including voluntary participation
Szabó and Hauert (2002); Hauert  _et al._ (2002a, b); Semmann  _et al._
(2003), punishment Fehr and Gächter (2002); Boyd  _et al._ (2003); Perc  _et
al._ (2017); Yang  _et al._ (2017) and reward Sigmund  _et al._ (2001);
Szolnoki and Perc ; Fang  _et al._ (2019), social diversity Santos  _et al._
(2008), noise Guan  _et al._ (2007), network heterogeneity Rong and Wu (2009),
reputation and reciprocity Xia  _et al._ (2023), etc. These endeavors provide
important insights into the mechanism of cooperation emergence. However, most
of these findings are mainly based on imitation learning (IL) Nowak and May
(1992); Roca  _et al._ (2009), where players imitate strategies of neighbors
who are better off in terms of payoffs. In essence, IL can be taken as a
simple version of social learning Bandura and Walters (1977), where
individuals learn from others in their socioeconomic activities, through
observation or instruction, which may or may not involve direct experiences.
As a different paradigm, reinforcement learning (RL) Sutton and Barto (2018)
provides a completely distinct framework for understanding human behaviors,
which has garnered increasing attention in recent years, including deciphering
the emergence mystery of cooperation Tanabe and Masuda (2012); Ezaki  _et al._
(2016); Horita  _et al._ (2017); Ding  _et al._ (2019); Fan  _et al._ (2022);
Zhang  _et al._ (2020); Wang  _et al._ (2022, 2023); He  _et al._ (2022); Ding
_et al._ (2023); Geng  _et al._ (2022); Yang  _et al._ (2024); Shi and Rong
(2022); Zhang  _et al._ (2024a); Zhao  _et al._ (2024), trust Zheng  _et al._
(2024), resource allocation Andrecut and Ali (2001); Zhang  _et al._ (2019),
and other collective behaviors for humans Tomov  _et al._ (2021); Shi and Rong
(2022). There players acquire information through interactions with the
environment and adjust behaviors based on rewards provided by the environment,
aiming to maximize their cumulative rewards.

Figure 1: Three different setups for PGG. (Left) PGG with Fermi rule: players
engage in the PGG and update their strategies using the Fermi updating rule.
(Middle) PGG with Q-learning algorithm: players engage in the PGG and update
their strategies using the Q-learning algorithm. (Right) VPGG with Q-learning
algorithm: players engage in the VPGG and update their strategies using the
Q-learning algorithm. The grey circle represents the focal individual, and
red, blue, and pink represent the cooperators, defectors, and loners in the
neighborhood. At the bottom of each panel, the shown are the Fermi rule, and
the corresponding Q-tables for the other two setups using Q-learning,
respectively.

Within the framework of RL, some new insights into the mechanism of
cooperation emergence for the pairwise game have been provided. For example,
by adopting the Bush-Mosteller model as a simple RL setup, Masuda _et al._
reproduce the failure of network reciprocity Ezaki and Masuda (2017) and
discuss the condition for cooperation emergence Masuda and Nakamura (2011).
Jia _et al._ use the Monte Carlo method to study pattern formation and phase
transitions towards cooperation in social dilemmas that are driven by RL, they
show that global players play a decisive role in ensuring cooperation Jia  _et
al._ (2021). Song _et al._ combine the iterated prisoner’s dilemma game with
the Bush-Mosteller reinforcement learning model and show that there exists a
moderate switching dynamics of the interaction intensity that is optimal for
the evolution of cooperation Song  _et al._ (2022). In Ref. Wang  _et al._
(2022), when reinforcement learners playing PD game are subject to Lévy noise,
a particularly positive role of Lévy noise is unveiled within this social
dilemma. Ding _et al._ reveal that a strong memory and a long-sighted
expectation yield the emergence of coordinated optimal policies, where both
players act like “win stay, lose shift” to maintain a high level of
cooperation Ding  _et al._ (2023).

Till now, these studies are primarily based on the PD where only two players
are engaged in each single game. Only recently, few works have started to
discuss the evolution of multiplayer game within the RL paradigm by adopting
PGG. In Ref. Jia  _et al._ (2022), the researchers investigate how cooperative
behavior evolves in a PGG on a network with empty nodes based on the Bush-
Mosteller model; the results show that only the existence of appropriate empty
nodes in the network can stimulate individual cooperative behavior. Wang _et
al._ propose a framework that combines the PGG and adaptive reward mechanisms
to better capture decision-making behaviors in multi-player interactions,
where they are mainly concentrating on the synergistic effects of them in RL
Wang  _et al._ (2023). Ref. Zhang  _et al._ (2024b) introduces the Q-learning
algorithm into the voluntary PGG to explore the evolution of cooperation, and
uncovers when the synergy factor is large and the adjust loner payoff’s
multiply factor is smaller, the number of cooperators becomes gradually
consistent, and the proportion of defectors evolves non-linearly for different
synergy factors.

However, these works based on PGG consider external factors, such as empty
nodes or rewards, it’s unclear what’s evolution of cooperation within the
original PGG game. More importantly, these works only consider the
individual’s own action information or are termed as a self-regarding setup,
where the players neglect the information of surroundings. This setup
simplifies the model but obviously contradicts our daily experiences, where
our decision-making is based upon not only our own actions but also the
environmental information (e.g. neighbors’ actions) we perceived. It is
generally believed that once more information about their surroundings is
incorporated, we could make better decisions to maximize the payoffs Sutton
and Barto (2018).

In this work, we investigate the evolution of cooperation of the original PGG
within the paradigm of RL where the environmental information is incorporated.
Specifically, we empower players with the Q-learning algorithm to revise their
strategies, where each player is guided by a Q-table. The environmental
information is incorporated into the state of the Q-table in a way that keeps
the Q-table simple enough. Our simulation shows that the incorporation of
environmental information is able to promote the cooperation prevalence, and
cooperation is more likely to emerge than the case of IL. Besides, for the PGG
with voluntary participation, the cyclic dominance in the previous studies
Szabó and Hauert (2002) within IL disappears; instead, a non-monotonic
cooperation dependence of the game parameter is observed, which is counter-
intuitive.

Figure 2: (Color online) Cooperation phase transition versus the gain factor
$r$ for three different models. (a) Players update their strategies using the
Fermi rule in PGG with $K$=0.1, and the threshold for cooperation emergence is
around $r_{c}$=4.59. (b) Players update their strategies using the Q-learning
algorithm in PGG, and the threshold for cooperation emergence is around
$r_{c}$=3.29. (c) Players update their strategies using the Q-learning
algorithm in VPGG, and the threshold for cooperation emergence is around
$r_{c}$=2.02. Other parameters: $\epsilon$=0.01, $\alpha$=0.1, $\gamma$=0.9,
$\sigma$=1.0, and the population size $N=100\times 100$.

##  II 2\. Model

Let’s consider a square lattice with periodic boundary condition, where each
node represents an individual engaging in the public goods game. Specifically,
we investigate two models with the Q-learning algorithm: the public goods game
(PGG) and the voluntary public goods game (VPGG). In the first model, the
focal player initiates a public goods game with four of its nearest neighbors,
where each player can either choose cooperation (C) or defection (D). The
choice of C means that an individual invests a certain amount (set to 1) in
the public pool, while D represents a free-riding behavior without any
contribution. Different from PGG, participation in VPGG is voluntary, players
can choose to be “loner” (L) as the third strategy, if they do not want to
play the game at all, and receive a small fixed income denoted as $\sigma$
(fixed at 1 if not stated otherwise). It is noteworthy that Ref. Zhang  _et
al._ (2024b) mainly focus on the impact of the loner strategy with self-
regarding setup on the evolution of cooperation, without considering the
environmental information of individuals and using adjustable payoffs for
loners. Here, we maintain consistency with previous studies on IL Szabó and
Hauert (2002) by fixing the loner’s payoff at 1, while incorporating
individuals’ environmental information into the Q-table setup to study the
evolution of cooperation within the RL framework. This allows for a better
comparison with the results of IL.

Suppose now a player $i$ initiates a PGG, five participants ($n_{i}=5$) are
involved, while the number of participants in VPGG needs to subtract the
number of loners $n_{L}$ in their neighborhood. After the investment, the
amount of money in the public pool is multiplied by a gain/synergy factor $r$,
and then evenly divided among all participating players. The expected payoffs
for different types of players $i$ are

$\pi_{i}=\begin{cases}\frac{rn_{C}}{n_{C}+n_{D}}-1,&\text{ if }s_{i}=C\\\
\frac{rn_{C}}{n_{C}+n_{D}},&\text{ if }s_{i}=D\\\ \sigma,&\text{ if
}s_{i}=L\end{cases}$ (1)

where $s_{i}$ represents the strategy of the player $i$, $n_{C,D}$ represents
the number of cooperators and defectors, respectively. By convention, in the
case of $r\in(1,n_{i})$ or $r\in(1,n_{i}-n_{L})$, social dilemmas arise where
the choice D brings a higher payoff than C.

In our models, individuals update their strategies according to the Q-learning
algorithm Watkins (1989, 1992), where each has a Q-table to guide their
decision-making. The Q-table is a two-dimensional table expanded by states and
actions. The action set consists of all available strategies, which is
$\mathbb{A}=\\{$C,D$\\}$ in PGG, and $\mathbb{A}=\\{$C,D,L$\\}$ in VPGG. The
state set serves to characterize the environment perceived by individuals. In
the previous self-regarding setup, the state only records the action adopted
by the player itself in the last round. Here, the environment information is
captured in a simple manner by comparing the number of cooperators $n_{C}$ and
defectors $n_{D}$ among the player’s four nearest neighbors. Therefore, there
are three possible states $s_{1,2,3}$ defined by the condition of
$n_{C}\textgreater n_{D}$, $n_{C}=n_{D}$, $n_{C}\textless n_{D}$,
respectively. We use the same state set in both models
$\mathbb{S}=\\{s_{1},s_{2},s_{3}\\}$. The two Q-tables are shown at the bottom
of Fig. 1. The idea of Q-table is score to the value of each action $a$ within
each state $s$. A larger value of $Q_{s,a}$ within a given row is supposed to
be more valuable within that state, which means the action $a$ is more likely
to be chosen when within the state $s$.

Input: $\alpha,\gamma,\epsilon$

1 Initialization;

2 for  _each agent_ do

3 Pick an action randomly from $\mathbb{A}$;

4 Create a Q-table with each item in the matrix near zero;

5

6for  _each agent_ do

7 Generate state: $s$;

8

9repeat

10 for  _each agent_ do

11 if  _rand() $<\epsilon$_ then

12 Pick an action randomly from $\mathbb{A}$;

13 else

14 Takes action according to the Q-table;

15 $a=\arg\max_{a_{i}}Q_{s,a_{i}}$;

16

17

18 for  _each agent_ do

19 Compute rewards according to Eq.(1);

20 Update state and Q-table according to Eq.(2);

21

22

23until  _the termination condition is met_ ;

Algorithm 1 Pseudocode of PGG/VPGG with Q-learning

Without loss of generality, each player is randomly assigned with an action
from the action set $\mathbb{A}$, and the elements $Q_{s,a}$ in the Q-tables
randomly assigned a value between [0, 1) independently. In our Monte Carlo
simulations, the evolution follows a synchronous updating procedure. In each
round $t$, with a probability $\epsilon$, they opt for a random action
$a\in\mathbb{A}$; with probability $1-\epsilon$, they strictly follow the
guidance of their Q-tables by choosing the action with the largest Q-value
given their current states. After all players make their moves, they receive
the reward $\pi_{i}$ following Eq. (1), and are able to compute the new state
$s^{\prime}$. As the next step, each player $i$ tries to draw lessons by
revising the corresponding $Q^{i}_{s,a}$ they just adopted in this round,
according to the Bellman equation Sutton and Barto (2018):

$\displaystyle Q^{i}_{s,a}(t+1)$
$\displaystyle=(1-\alpha)Q^{i}_{s,a}(t)+\alpha\left(\pi_{i}(t)+\gamma\max_{a^{\prime}}Q^{i}_{s^{\prime},a^{\prime}}(t)\right),$
(2)

where $s$ and $a$ represent the current state and action of the focal
individual $i$, and $s^{\prime}$ is the new state at $t+1$. $\alpha\in(0,1]$
is the learning rate determining how fast of the Q-value is revised,
$\gamma\in[0,1]$ is the discount factor, defining the significance of future
rewards. Notice that, $\pi_{i}(t)$ is the payoff for player $i$ when she/he
acts as the focal player initiating a PGG game, not all payoffs received.
Because the item $Q^{i}_{s,a}$ reflects the value of action the player $i$
actively made Zhang  _et al._ (2020). When all players update their Q-tables,
a single MC round is completed. The evolution consists of many MC steps and is
terminated until the system reaches equilibrium or the desired duration is
reached, the pseudocode is provided in Algorithm 1.

Besides the two models using Q-learning proposed above, we also adopt the
Fermi updating rule Szabó and Tőke (1998) in the traditional framework of IL
as the baseline model for comparison. Here, players make decisions by
comparing utility with their neighbors and imitate strategies of those who are
better off in the game. Since this setup has been extensively studied in
previous game-theoretic works, we refer to Ref. Szabó and Hauert (2002);
Flores  _et al._ (2022) for the simulation details. The three involved model
setups are intuitively illustrated in Fig. 1.

Figure 3: (Color online) Typical cooperation snapshots for the three setups.
In (a1-a4), PGG with Fermi rule, and cooperation clusters are formed to resist
the exploitation by defectors. The gain factor $r=4.60$. In (b1-b4), PGG with
the Q-learning algorithm, no apparent clusters of cooperators are seen
instead, and the evolution time takes longer. In (c1-c4), VPGG with the
Q-learning algorithm, there are also no obvious cooperation clusters, but the
resulting cooperation prevalence $f_{c}$ in panel (c4) is much higher than the
panel (b4). $\epsilon=0.01$, $\alpha$=0.1, $\gamma$=0.9, $\sigma$=1.0,
$r=3.32$ are used in both (b1-b4) and (c1-c4). Other parameters: the
population size $N=100\times 100$.  Figure 4: (Color online) Typical time
series in PGG with Q-learning. (a) Time evolution of two strategies, which
differs from the typical “first down and later up” trend often observed as a
sign of network reciprocity in the framework of imitation learning. (b) The
temporal evolution of the fractions of the three states $s_{1,2,3}$. After the
transient period, $s_{3}$ dominates. (c) The evolution of average reward for
two different strategies, where there is a crossover observed. (d-f) Time
evolution of action preference in three different states. The emergence of
cooperation primarily occurs in state $s_{1}$ ($n_{C}<n_{D}$). Other
parameters: $\epsilon=0.01$, $\alpha$=0.1, $\gamma$=0.9, $\sigma$=1.0,
$r$=3.32, and the population size $N=100\times 100$.

##  III 3\. Results

We first report the fraction of different strategies as a function of the gain
factor $r$ in Fig. 2 for the three models. Fig. 2(a) shows the results with
the Fermi updating rule, where we observe that cooperation emerges as the gain
factor is larger than a critical threshold $r>r_{c}\approx 4.59$, in line with
previous studies Flores  _et al._ (2022). When the gain factor $r$ is high
enough, a cooperation phase transition is seen as shown here. Notice that, the
model with a strictly synchronous updating exhibits similar phase transition,
with the threshold being $r_{c}\approx 4.80$ (data now shown), showing a bit
more difficult for the cooperation to emerge. Fig. 2(b) shows the results for
PGG with the Q-learning algorithm. Compared to Fig. 2(a), the cooperation
phase transition is qualitatively the same, but with a reduced threshold
$r_{c}\approx 3.29$. This means that cooperation is more likely to emerge
within the reinforcement learning paradigm than in the case of imitation
learning.

When voluntary participation is introduced, Fig. 2(c) shows that the threshold
is further reduced to $r_{c}\approx 2.02$ where cooperation is most likely to
emerge among the three models we studied here. Detailed examination indicates
that the introduction of “loner” is an efficient manner to considerably
inhibit the growth of defectors. For a lower gain factor $r<r_{c}$, the
fraction of loners dominate. Defectors appear only when cooperators survive as
$r>r_{c}$, but as the gain factor further increases cooperators dominates. The
fraction of defectors remains $f_{D}\lesssim 0.2$ for the whole range. This
overall picture aligns with findings from previous work Szabó and Hauert
(2002) by IL. However, there is a non-monotonic dependence region for
$r\in(2.13,2.50)$ where the increase in $r$ leads to the reduction in $f_{C}$
accompanied by the increase in $f_{D}$. Notice that, the non-vanishing $f_{L}$
at the end of $r\rightarrow 5$ effectively inhibits the invasion of defectors,
guaranteeing a higher level of cooperation compared to the case without loner
[see Fig. 2(b)]. A systematic investigation of the impact of the two learning
parameters is given in Appendix A. It shows that when individuals both focus
on historical experience and have a long-term vision, cooperation is likely to
emerge. Otherwise, either a forgetful property (a large $\alpha$) and/or a
short-term vision (a small $\gamma$) lead to the failure of its emergence.

To develop some intuition for the evolution process for the three scenarios,
we present some typical spatio-temporal snapshots of the system. Starting from
random initial conditions, Figs. 3(a1-a4) show that some cooperation clusters
are formed. This is a typical illustration of network reciprocity Nowak and
May (1992), where the structured population enables cooperators to persist by
forming clusters, thereby avoiding the exploitation by defectors. For the
given parameter $r=4.60$, being slightly larger than the threshold, the
resulting prevalence of cooperation is stabilized around 0.13, while the
cooperation clusters are not frozen but fluctuate all the time.

The PGG with the Q-learning algorithm is illustrated in Fig. 3(b1-b4), where
$r$ is also chosen slightly larger than the threshold. As can be seen, it’s
evident that the patterns are quantitatively distinct from Figs. 3(a1-a4),
where there are no apparent cooperator clusters observed. Over time, the
number of cooperators gradually declines, and upon reaching a steady state,
these cooperators are evenly scattered across the whole domain. Additionally,
the time-scale reaching a steady state seems longer in this case than in Figs.
3(a1-a4).

Finally, we showcase the pattern evolution for the VPGG with Q-learning, with
identical parameter $r=3.32$ as adopted in Fig. 3(b1-b4) for comparison.
Similar to the PGG, no obvious cooperation clusters are seen in Fig. 3(c1-c4).
As time goes by, however, both the numbers of loners and defectors decrease,
and cooperators dominate in the end, with a notably higher cooperation
prevalence compared to Fig. 3(b1-b4). Furthermore, the patterns depicted do
not display the waves with cyclic dominance as observed in IL Szabó and Hauert
(2002). This suggests a different mechanism working here.

Figure 5: (Color online) Two scenarios within state $s_{3}$ where
$n_{C}<n_{D}$. The gray circle represents the focal individual, who is either
to cooperate or to defect. Red and blue circles represent the cooperators and
defectors in neighbors, respectively. Cooperation is almost impossible in
scenario (a), while it may occur in scenario (b).

##  IV 4\. Mechanism analysis

To understand how cooperation evolves when people play PGG with Q-learning
algorithm. In this section, we move to the mechanism analysis of the latter
two cases playing with Q-learning.

###  IV.1 4.1 Mechanism analysis of PGG with Q-learning

We commence our analysis by delving into PGG with Q-learning, where the
cooperation phase transition and pattern evolution have been shown above. To
further observe the temporal evolution, we present time series in Fig. 4(a),
with identical parameters adopted in Fig. 3(b1- b4). Two strategy fractions
remain at the initial level around 0.5 and begin to decline later on and
eventually stabilize. This observation is fundamentally different from the
“first down and later up” trend often regarded as a sign of network
reciprocity in IL Perc  _et al._ (2008); Szolnoki and Perc (2009).

To understand this trend, we need to take into account the surrounding
information and the resulting payoffs. As depicted in Fig. 4(b), with random
initial conditions, the proportions of the three states $s_{1,2,3}$ are
approximately identical. Fig. 4(c) shows that initially free-riders (D)
consistently yield higher payoffs than cooperators (C) since approximately
half of the individuals choose to cooperate. As a result, players tend to
defect by learning, leading to an increasing proportion of the state $s_{3}$
($n_{C}<n_{D}$), and accordingly the other two state fractions of $s_{1,2}$
decrease. Specifically, the state $s_{3}$ includes two different scenarios
(see Fig. 5): (a) The focal individual is surrounded by four defectors, in
which case choosing C yields a negative payoff, which is disadvantageous. (b)
There is one cooperator present in the neighborhood around the focal
individual, choosing C will result in a positive payoff, leading to an
increase in the Q-value in this state. While opting for strategy D yields a
higher payoff but this also regresses the environment back to scenario (a).

Figure 6: (Color online) Time evolution of $\Delta$Q after transient within
state $s_{3}$. Typical time evolution of $\Delta Q=Q_{s_{3},C}-Q_{s_{3},D}$ in
state $s_{3}$ ($n_{C}\textless n_{D}$). The inset shows the frequencies of
occurrences when $\Delta Q>0$ and $\Delta Q<0$ at the same time intervals. Two
parameters $r=3.32$ (black) and $r=3.4$ (orange) are chosen here for
comparison. Other parameters: $\epsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$, and
the population size $N=100\times 100$.

In the process of decision-making using the Q-learning algorithm, individuals
focus on both historical experience and long-term rewards, therefore they try
to avoid falling into scenario (a). Consequently, by factoring in the low
payoff in scenario (a) and the potential for higher payoffs in the new state
$s_{2}$ ($n_{C}=n_{D}$), the preference of choosing C in state $s_{3}$ arises,
leading to the emergence of cooperation. In fact, the choice of C indeed
brings higher payoffs on average as shown by the crossover in Fig. 4(c). In
Fig. 4(d-e), we show the evolution of action preferences in three different
states. It can be observed that the emergence of cooperation primarily occurs
mostly in the state $s_{3}$.

Figure 7: (Color online) Typical time series in VPGG with Q-learning. (a) Time
evolution of three strategies. As the defectors begin to decline, the fraction
of cooperators shifts from a decreasing to an increasing trend, while the
corresponding loners shift from an increasing to a decreasing trend,
indicating the emergence of cooperators within clusters of loners. (b) Time
evolution of the fractions of the three states. After a sufficiently large
number of Monte Carlo steps, it can be observed that the system is most likely
in the state of $s_{1}$, while there exists a certain proportion of states
$s_{2}$ and $s_{1}$ in dynamic equilibrium. (c) The evolution of average
reward for three different strategies. In the early stages, loners
consistently generate higher average payoffs compared to the other strategies.
(d-f) Time evolution of action preference in three states $s_{1,2,3}$. The
emergence of cooperation primarily occurs in state $s_{1}$. Other parameters:
$\epsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$, $r=2.11$ and the population size
$N=100\times 100$.

To further confirm, we monitor the preference evolution captured by the
Q-values of a typical individual within state $s_{3}$, as shown in Fig. 6.
It’s convenient to define $\Delta Q=Q_{C}-Q_{D}$ as the cooperation
preference, a positive value corresponds to the case that cooperation is
preferred by learning and vice versa. With this in mind, we observe that the
value of $\Delta Q$ fluctuates up and down around 0. Most of the time, $\Delta
Q<0$ with intermittent occurrences of $\Delta Q>0$. This means that D is the
preferred mostly but C is also selected from time to time. The two actions
compete and reach a balance as seen in Fig. 4(f). As the gain factor $r$
increases, the payoff advantage of defectors over cooperators diminishes,
indicating an increased likelihood of cooperation emergence. This is validated
in the inset of Fig. 6, which shows an increased frequency of $\Delta Q>0$
when $r$ increases to 3.4.

The above analysis indicates that self-reflective decision-making within
Q-learning does not rely on clusters of cooperators to resist defectors.
Instead, it possesses the potential to escape from the most adverse
environments of full defection, thereby leading to the emergence of
cooperation.

###  IV.2 4.2 Mechanism analysis of VPGG with Q-learning

We next turn to the mechanism in the VPGG with the Q-learning, especially the
impact of voluntary participation, shown in Fig. 7. Starting from random
initial conditions, the fraction of loners (L) starts to increase, and both
defectors and cooperators decrease. But once the fraction of D stabilizes near
zero, cooperators start to increase and loners decrease, and eventually, the
two fractions become stabilized with approximately equal size for the given
parameters, shown in Fig. 7(a). This evolution is also captured by the
fraction of the three states $s_{1,2,3}$, where the fraction of state $s_{1}$
dominates and the fraction of $s_{3}$ diminishes.

Fig. 7(c) shows that the payoff of the loner is much larger than the other two
fractions at the initial stage, which explains the number increases in loners.
Figs. 7(d-f), respectively, show the evolution of action fraction within the
three states $s_{1,2,3}$, where the action preference differs significantly.
To specifically understand the non-trivial evolution trend in Fig. 7(a), we
need to examine the action prevalence in the three states $s_{1,2,3}$ by
incorporating the state fractions shown in Fig. 7(b). At the first stage
($t\lesssim 10^{4}$), the state fraction of $s_{3}$ is slightly higher than
the other two, where the action fraction of loner is highest ($10^{2}\lesssim
t\lesssim 10^{4}$), the three action fractions in $s_{1,2}$ are approximately
identical; these observations then explain the increasing trend in loners and
decreasing trend of both C and D.

At around $t=10^{4}$, the action fraction of L becomes largest in all three
$s_{1,2,3}$, and loner becomes the most preferred choice in all cases by
learning at this moment. This change then leads to the peak of fraction in L
in Fig. 7(a). As $f_{D}\rightarrow 0$, this change provides an ideal
surrounding for players to choose cooperation, which brings higher rewards
than the case of being loners (see Fig. 7(c)). Therefore, the state $s_{1}$
dominates afterward, and the action of C becomes a more preferred action than
L, which is supported by the payoff comparison provided in Table I in Appendix
B. That’s the reason for the crossover for loners and cooperators as seen in
later evolution in Fig. 7(d). Though, within state $s_{2}$, being a loner is a
better choice than being a cooperator (see Table II in Appendix B). Within
state $s_{3}$, being a loner is still the best choice, and being a cooperator
is worst, though the fraction of $s_{3}$ is too small.

Taken together Fig. 7(a-c), players opt to be loners within $s_{2,3}$, and
cooperator is more preferred within $s_{1}$. This is the reason why only
cooperators and loners survive and coexist. The overall picture is that
individuals in this setup choose strategy L from time to time to inhibit the
growth of defectors, laying the groundwork for cooperators to thrive in a sea
of loners.

Figure 8: (Color online) Time evolution of Q value in a stable state. In VPGG,
time series of Q-value for three different strategies in a stable state. Other
parameters: $\epsilon=0.01$, $\alpha=0.1$, $\gamma=0.9$, and the population
size $N=100\times 100$.

By this point, one might question are the coexistence stability of L and C,
and whether the strategy L can entirely be transformed into C or vice versa.
The answer is no. An all-C system is unstable and susceptible to the invasion
of D. Similarly, an all-L system is also unstable and can be invaded by C, and
an all-D system is susceptible to the invasion of L. This cyclic dominance
determines that none could dominate, they have to coexist in the form of a
subtle balance, even though the fraction of D is very small. An example is
shown in Fig. 8, where we present the time series of Q values for all three
strategies in a stable state, randomly selecting an individual. We see that
$Q_{C}$ fluctuates up and down around $Q_{L}$, while $Q_{D}$ shows a
dependence on $Q_{C}$, their orders are reshuffled consistently. This
observation demonstrates that the three strategies are in dynamic equilibrium.

###  IV.3 4.3 Analysis of non-monotonic variation

Figure 9: (Color online) The analysis of non-monotonic interval in VPGG. (a)
The left side of the y-axis shows the fraction of strategies as a function of
the gain factor $r$. The right side shows the tendency of choosing strategy L
based on theoretical analysis, characterized by the proportion of players
where both conditions of $\pi_{L}>\pi_{C}$ and $\pi_{L}>\pi_{D}$ are
satisfied. The inset shows the corresponding values of $\Delta Q$ and
$Q_{L,D}$. (b) The fraction of cooperators as the discount factor $\gamma$
varies for three different $r$ values. Other parameters: $\epsilon=0.01$,
$\alpha=0.1$, $\gamma=0.9$, and the population size $N=100\times 100$.

Till now, the non-monotonic dependence of cooperation level for
$r\in(2.13,2.50)$ in Fig. 2(c) remains unexplained. In this interval, the
fraction of the strategy L remains nearly constant around 0.5, and the
cooperation prevalence decreases as the gain factor $r$ is increased.

To proceed, we estimate the tendency to choose strategy L in theory by
computing the proportion of players where the condition
$\pi_{C}\textless\pi_{L}$ and $\pi_{D}\textless\pi_{L}$ are both satisfied
within all possible states. Here, $\pi_{C}$, $\pi_{D}$, and $\pi_{L}$
represent the payoffs of choosing C, D, and L, respectively. This reveals the
tendency of individuals to choose the strategy L as a function of $r$, as
depicted in Fig. 9(a). As shown in Fig. 9(a), the tendency to choose strategy
L exhibits a stage-wise plateau in the interval $r\in(2.15,2.49)$. Within this
range, the fraction of loners remains fixed, while the fraction of cooperators
shows a decreasing trend. As $r$ further increases to $r\geq 2.5$, the
tendency of loners continues to decrease, approaching a new plateau, from the
results of numerical simulations, we notice a reduction in the fraction of
loners and a simultaneous increase in the fraction of cooperators. Notice
that, these stage-wise plateaus correspond to the present rewards that are
immediately received at present, where the impact of expected payoffs in the
future is not included.

The reason why the plateau well matches the non-monotonic interval is exactly
because the impact of expected payoffs in the future is weakened. As depicted
in Fig. 9(b), three typical cases within the non-monotonic region
$r\in(2.13,2.50)$ are shown, where we observe that the dependence of the
cooperation prevalence on the discount factor $\gamma$ is non-monotonic, there
is an optimal value of $\gamma$ in each case yielding the highest cooperation
prevalence. This is in sharp contrast with the results obtained from previous
observations Ding  _et al._ (2023); Zheng  _et al._ (2024), where a higher
expectation (i.e. a larger $\gamma$) always leads to a better outcome of
evolution; two such examples outside of the non-monotonic region ($r=2.05$ and
2.66) are shown in Fig. 12 in Appendix C C. This difference implies that,
within this non-monotonic interval, a high emphasis on future rewards is less
favorable for the emergence of cooperation, suggesting potential challenges in
maintaining long-term cooperation stability, the underlying cause of this
result is the uncertainty in the environment due to changes in the strategies
of other individuals.

Specifically, as discussed in the above section, L and C are more preferred
than D within the three states $s_{1,2,3}$ in typical scenarios. But this is
not always the case, the inset in Fig. 9(a) shows that the Q-value for L is no
longer larger than the value for D, instead their difference $\Delta
Q=Q_{D}-Q_{L}$ becomes positive and is enlarged as $r$ increases within this
region. This means that players hesitate in the choice of action D or L,
leading to increasing uncertainty of the environment. This explains the non-
monotonic trend observed in Fig. 9(b), where the power of long-term vision
doesn’t work anymore and players prioritize their immediate rewards.

Outside of this non-monotonic region, players no longer get vacillated between
L and D as the environmental information becomes certain, and the expectation
factor for the future again plays its role in the emergence of cooperation. By
further increasing in $r$, the advantage of D over L becomes more obvious, and
the same is also true for $Q_{C}>Q_{D}$, leading to the decreasing fraction of
L and increasing trend in C.

##  V 5\. Conclusion and discussion

In summary, we investigate the original PGG and PGG with voluntary
participation (VPGG) within the paradigm of reinforcement learning, where each
player acts following a Q-learning algorithm. Different from most previous
works where the self-regarding Q-learning algorithm is employed, we
incorporate environmental information into the state by comparing the numbers
of cooperators and defectors in their neighborhoods. We observe that the
thresholds for cooperation emergence in both cases occurs are much reduced
compared to the case in the IL framework. The cooperation emergence is found
to be most likely to occur in VPGG. Unexpectedly, a non-monotonic dependence
of cooperation prevalence on the gain factor is revealed in the VPGG. By the
analysis of the evolution of state, action preference, rewards, and the
evolution of the Q-table, we clarify the mechanism of cooperation in both
cases and the reason behind the non-monotonic dependence in VPGG.

Our research reveals significant differences in the mechanisms of cooperation
emergence between the traditional IL and our RL. Within RL, the emergence of
cooperation does not rely on the cooperator clustering, the so-called network
reciprocity Nowak and May (1992), instead by incorporating the historical
experiences, present rewards, and the expected payoffs in the future, players
avoid complete defection and cooperators emerge in the form of scattered
distribution. Furthermore, in the VPGG, the choice of the loner (non-
participating players) for the individuals significantly inhibits the survival
of defectors, and players turn to cooperators when the loners are dominating.

Unexpectedly, there is a non-monotonic interval where the number of
cooperators decreases with the gain factors $r$. Our analysis reveals that
this is due to the enhanced environmental uncertainty, reducing the
expectation of future rewards, and cooperation is thus degraded as players put
more weight on the present rewards. This uncertainty induced risk aversion
disappears as the gain factor $r$ is further increased where the environment
becomes less uncertain, and the cooperation prevalence continues to increase
with increasing $r$.

Till now, reinforcement learning has been employed to explain various socio-
economic activities beyond cooperation, including trust Zheng  _et al._
(2024), resource allocation Andrecut and Ali (2001); Zhang  _et al._ (2019),
and other human behaviors Zhang  _et al._ (2020); Tomov  _et al._ (2021); Shi
and Rong (2022). In fact, there is abundant experimental evidence in
neuroscience Lee  _et al._ (2012); Rangel  _et al._ (2008) indicating that
many decision-making processes of humans fall into RL paradigm. Notice that,
we do not devalue social learning Olsson  _et al._ (2020). We conjecture that,
in our daily lives both paradigms are relevant, and each works in some
specific scenarios and together they underpin most of our decision-makings. Of
course, this speculation requires further validation through behavioral
experiments. We believe that the research conducted in this work is crucial
for elucidating some basic questions along this research line.

##  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). ZGZ is supported by Excellent Graduate
Training Program of Shaanxi Normal University [Grants No. LHRCTS23064].

##  Appendix A Appendix A: Phase diagrams of the learning parameters

To systematically examine the impact of the learning parameters on the
cooperation evolution, we provide phase diagrams for both PGG and VPGG with
Q-learning algorithm.

In PGG with Q-learning, we find a prominent emergence of cooperation is
observed in a physically meaningful region, see the phase diagram Fig. 10, it
reports the final fractions of cooperators in the learning parameter domain
($\gamma$, $\alpha$). As shown, there is a white region where the cooperation
emerges, where the learning rate $\alpha$ is small and the discount factor
$\gamma$ is large. This observation means that when individuals focus on both
historical experience and long-term vision, cooperation is likely to emerge.
Otherwise, either a forgetful property (a large $\alpha$) and/or a short-term
vision (a small $\gamma$) lead to the failure of its emergence.

Similar results can also be observed in the VPGG, as shown in Fig. 11. The
emergence of cooperators fails with the decrease in $\gamma$ and/or the
increase in $\alpha$. However, the fraction of defectors remains small across
the entire parameter space. This result is easily understandable because, with
the introduction of the loner strategy, individuals tend to choose strategy L
with a higher payoff rather than the all-defect strategy with zero payoffs as
we discuss above.

These observations are in line with our previous findings revealed in the
emergence of cooperation in the prisoners’ dilemma Ding  _et al._ (2023) or
the donation game Zhang  _et al._ (2024a), and the emergence of trust Zheng
_et al._ (2024).

Figure 10: (Color online) The color-coded stationary fraction of cooperators
in the domain ($\gamma$, $\alpha$) for PGG with the Q-learning algorithm.
Cooperation emerges at the corner of large $\gamma$ and small $\alpha$,
indicating that when individuals focus on both historical experience and the
long-term vision, cooperation emerges. Each data is averaged 100 times after a
transient of $5\times 10^{5}$. Other parameters: $\epsilon=0.01$, $r=3.32$,
and the population size $N=100\times 100$.

Figure 11: (Color online) The color-coded stationary fraction of the three
strategies in the domain ($\gamma$, $\alpha$) for VPGG with the Q-learning
algorithm. Cooperation emerges at the corner of large $\gamma$ and small
$\alpha$, showing that when individuals focus on both historical experience
and the long-term vision, cooperation emerges. Each data is averaged 100 times
after a transient of $5\times 10^{5}$. Other parameters: $\epsilon=0.01$,
$r=2.11$, and the population size $N=100\times 100$.

##  Appendix B Appendix B: Payoff analysis in VPGG

For a more intuitive understanding of the transition in individual action
preferences from the perspective of payoffs, we provide payoffs for different
strategies in three states $s_{1,2,3}$, as shown in table 1 to table 3.

Table 1 shows payoffs for individuals in state $s_{1}$ where $n_{C}>n_{D}$.
Only in scenario $(nC,0D)$ does choosing cooperation yields a higher payoff
than choosing the loner. This implies that in an environment where defectors
exist, cooperators are prone to shifting toward being loners.

Table 2 displays payoffs for individuals in state $s_{2}$ where $n_{C}=n_{D}$.
Choosing L consistently yields the highest payoff, except in the scenario
$(0C,0D)$, where all individuals around are loners, facilitating the emergence
of cooperation.

Table 3 presents payoffs for individuals in state $s_{3}$ where $n_{C}<n_{D}$.
Choosing L consistently yields the highest payoff, suggesting individuals
typically opt for higher payoff strategies when in this state.

$s_{1}$ Action | C ($a_{1}$) | D ($a_{2}$) | L ($a_{3}$)
---|---|---|---
1C, 0D | $1.11$ | $1.055$ | $1$
2C, 0D | $1.11$ | $1.406$ | $1$
2C, 1D | $0.5825$ | $1.055$ | $1$
3C, 0D | $1.11$ | $1.5825$ | $1$
3C, 1D | $0.688$ | $1.266$ | $1$
4C, 0D | $1.11$ | $1.688$ | $1$
Table 1: In the VPGG, the average payoffs when choosing different actions in state $s_{1}$, where $n_{C}>n_{D}$. Parameter: $r=2.11$. $s_{2}$ Action | C ($a_{1}$) | D ($a_{2}$) | L ($a_{3}$)
---|---|---|---
0C, 0D | $1.11$ | $0$ | $1$
1C, 1D | $0.406$ | $0.703$ | $1$
2C, 2D | $0.266$ | $0.844$ | $1$
Table 2: In the VPGG, the average payoffs when choosing different actions in state $s_{2}$, where $n_{C}=n_{D}$. Parameter: $r=2.11$. $s_{3}$ Action | C ($a_{1}$) | D ($a_{2}$) | L ($a_{3}$)
---|---|---|---
0C, 4D | $-0.578$ | $0$ | $1$
0C, 3D | $-0.4725$ | $0$ | $1$
0C, 2D | $-0.296$ | $0$ | $1$
0C, 1D | $0.055$ | $0$ | $1$
1C, 2D | $0.055$ | $0.5275$ | $1$
1C, 3D | $-0.156$ | $0.422$ | $1$
Table 3: In the VPGG, the average payoffs when choosing different actions in
state $s_{3}$, where $n_{C}<n_{D}$. Parameter: $r=2.11$.

##  Appendix C Appendix C: Dependence on $\gamma$ outside of the non-monotonic
interval

As a comparison, here we present the fraction of cooperators outside of the
non-monotonic interval (at $r=2.05$ and $r=2.66$) as the discount factor
$\gamma$ varies, as shown in Fig. 12. We can observe that the fraction of
cooperators continues to increase with the increasing $\gamma$, quantitatively
different from the observations made in Fig. 9(b) in the non-monotonic
interval. This indicates that individuals valuing future rewards are conducive
to the emergence of cooperation, which aligns with our previous findings Ding
_et al._ (2023); Zheng  _et al._ (2024), where a higher expectation (i.e. a
larger $\gamma$) always leads to a better outcome of evolution.

Figure 12: (Color online) The fraction of cooperators outside of the non-
monotonic interval. We can see that the fraction of cooperation increases with
the increase of $\gamma$. Other parameters: $\epsilon=0.01$, $\alpha=0.1$, and
the population size $N=100\times 100$.

## References

  * Milinski  _et al._ (2006) M. Milinski, D. Semmann, H.-J. Krambeck, and J. Marotzke, Proceedings of the National Academy of Sciences 103, 3994 (2006).
  * Kollock (1998) P. Kollock, Annual Review of Sociology 24, 183 (1998).
  * Hardin (1968) G. Hardin, Science 162, 1243 (1968).
  * Nowak and Sigmund (2004) M. A. Nowak and K. Sigmund, Science 303, 793 (2004).
  * Nowak (2006) M. A. Nowak, Science 314, 1560 (2006).
  * Axelrod and Hamilton (1981) R. Axelrod and W. D. Hamilton, Science 211, 1390 (1981).
  * Doebeli and Hauert (2005) M. Doebeli and C. Hauert, Ecology Letters 8, 748 (2005).
  * 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).
  * Hauert and Szabó (2003) C. Hauert and G. Szabó, Complexity 8, 31 (2003).
  * Perc  _et al._ (2013) M. Perc, J. Gómez-Gardeñes, A. Szolnoki, L. M. Florĺa, and Y. Moreno, Journal of The Royal Society Interface 10, 20120997 (2013).
  * Tanimoto (2013) J. Tanimoto, Physical Review E 87, 062136 (2013).
  * Dawes (1980) R. M. Dawes, Annual Review of Psychology 31, 169 (1980).
  * Fischbacher  _et al._ (2001) U. Fischbacher, S. Gächter, and E. Fehr, Economics Letters 71, 397 (2001).
  * Wang  _et al._ (2009) J. Wang, F. Fu, T. Wu, and L. Wang, Physical Review E 80, 016101 (2009).
  * Gintis  _et al._ (2003) H. Gintis, S. Bowles, R. Boyd, and E. Fehr, Evolution and Human Behavior 24, 153 (2003).
  * Szabó and Hauert (2002) G. Szabó and C. Hauert, Physical Review Letters 89, 118101 (2002).
  * Hauert  _et al._ (2002a) C. Hauert, S. D. Monte, J. Hofbauer, and K. Sigmund, Science 296, 1129 (2002a).
  * Hauert  _et al._ (2002b) C. Hauert, S. D. Monte, J. Hofbauer, and K. Sigmund, Journal of Theoretical Biology 218, 187 (2002b).
  * Semmann  _et al._ (2003) D. Semmann, H.-J. Krambeck, and M. Milinski, Nature 425, 390 (2003).
  * Fehr and Gächter (2002) E. Fehr and S. Gächter, Nature 415, 137 (2002).
  * Boyd  _et al._ (2003) R. Boyd, H. Gintis, S. Bowles, and P. J. Richerson, Proceedings of the National Academy of Sciences 100, 3531 (2003).
  * Perc  _et al._ (2017) M. Perc, J. J. Jordan, D. G. Rand, Z. Wang, S. Boccaletti, and A. Szolnoki, Physics Reports 687, 1 (2017).
  * Yang  _et al._ (2017) B. Yang, M. Yang, W. Liu, and X. Chen, Acta Physica Sinica 66, 196401 (2017).
  * Sigmund  _et al._ (2001) K. Sigmund, C. Hauert, and M. A.Nowak, Proceedings of the National Academy of Sciences 98, 10757 (2001).
  * (25) A. Szolnoki and M. Perc, Europhysics Letters .
  * Fang  _et al._ (2019) Y. Fang, T. P. Benko, M. Perc, H. Xu, and Q. Tan, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 475, 20190349 (2019).
  * Santos  _et al._ (2008) F. C. Santos, M. D. Santos, and J. M. Pacheco, Nature 454, 213 (2008).
  * Guan  _et al._ (2007) J. Guan, Z. Wu, and Y. Wang, Physical Review E 76, 056101 (2007).
  * Rong and Wu (2009) Z. Rong and Z. Wu, Europhysics Letters 87, 30001 (2009).
  * Xia  _et al._ (2023) C. Xia, J. Wang, M. Perc, and Z. Wang, Physics of Life Reviews 46, 8 (2023).
  * Nowak and May (1992) M. A. Nowak and R. M. May, Nature 359, 826 (1992).
  * Roca  _et al._ (2009) C. P. Roca, J. A. Cuesta, and A. Sánchez, Physics of Life Reviews 6, 208 (2009).
  * Bandura and Walters (1977) A. Bandura and R. H. Walters,  _Social Learning Theory_ , Vol. 1 (Englewood cliffs Prentice Hall, 1977).
  * Sutton and Barto (2018) R. S. Sutton and A. G. Barto,  _Reinforcement Learning: An Introduction_ (MIT press, 2018).
  * Tanabe and Masuda (2012) S. Tanabe and N. Masuda, Journal of Theoretical Biology 293, 151 (2012).
  * Ezaki  _et al._ (2016) T. Ezaki, Y. Horita, M. Takezawa, and N. Masuda, PLoS Computational Biology 12, e1005034 (2016).
  * Horita  _et al._ (2017) Y. Horita, M. Takezawa, K. Inukai, T. Kita, and N. Masuda, Scientific Reports 7, 39275 (2017).
  * Ding  _et al._ (2019) H. Ding, G. Zhang, S. Wang, J. Li, and Z. Wang, Physica A: statistical mechanics and its applications 536, 122551 (2019).
  * Fan  _et al._ (2022) L. Fan, Z. Song, L. Wang, Y. Liu, and Z. Wang, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 123140 (2022).
  * 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).
  * Yang  _et al._ (2024) Z. Yang, L. Zheng, M. Perc, and Y. Li, Applied Mathematics and Computation 463, 128364 (2024).
  * Shi and Rong (2022) Y. Shi and Z. Rong, IEEE Transactions on Circuits and Systems II: Express Briefs 69, 2463 (2022).
  * Zhang  _et al._ (2024a) J. Zhang, Z. Rong, G. Zheng, J. Zhang, and L. Chen, Journal of Physics: Complexity 5, 025006 (2024a).
  * Zhao  _et al._ (2024) C. Zhao, G. Zheng, C. Zhang, J. Zhang, and L. Chen, Chaos: An Interdisciplinary Journal of Nonlinear Science 34, 073123 (2024).
  * 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.-C. Lai, Physical Review E 99, 032302 (2019).
  * Tomov  _et al._ (2021) M. S. Tomov, E. Schulz, and S. J. Gershman, Nature Human Behaviour 5, 764 (2021).
  * Ezaki and Masuda (2017) T. Ezaki and N. Masuda, PloS one 12, e0189220 (2017).
  * Masuda and Nakamura (2011) N. Masuda and M. Nakamura, Journal of Theoretical Biology 278, 55 (2011).
  * Jia  _et al._ (2021) D. Jia, H. Guo, Z. Song, L. Shi, X. Deng, M. Perc, and Z. Wang, New Journal of Physics 23, 083020 (2021).
  * Song  _et al._ (2022) Z. Song, H. Guo, D. Jia, M. Perc, X. Li, and Z. Wang, Neurocomputing 513, 104 (2022).
  * Jia  _et al._ (2022) D. Jia, T. Li, Y. Zhao, X. Zhang, and Z. Wang, Applied Mathematics and Computation 413, 126658 (2022).
  * Zhang  _et al._ (2024b) H. Zhang, T. An, P. Yan, K. Hu, J. An, L. Shi, J. Zhao, and J. Wang, Chaos, Solitons & Fractals 178, 114358 (2024b).
  * Watkins (1989) C. J. C. H. Watkins,  _Learning from delayed rewards (Ph.D. thesis)_ , Ph.D. thesis (1989).
  * Watkins (1992) P. Watkins, Christopher J. C. H.and Dayan, Machine Learning 8, 279 (1992).
  * Szabó and Tőke (1998) G. Szabó and C. Tőke, Phys. Rev. E 58, 69 (1998).
  * Flores  _et al._ (2022) L. S. Flores, M. A. Amaral, M. H. Vainstein, and H. C. Fernandes, Chaos, Solitons & Fractals 164, 112744 (2022).
  * Perc  _et al._ (2008) M. Perc, A. Szolnoki, and G. Szabó, Physical Review E 78, 066101 (2008).
  * Szolnoki and Perc (2009) A. Szolnoki and M. Perc, The European Physical Journal B 67, 337 (2009).
  * 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. Montague, Nature Reviews. Neuroscience 9, 545 (2008).
  * Olsson  _et al._ (2020) A. Olsson, E. Knapska, and B. Lindström, Nature Reviews Neuroscience 21, 197 (2020).