# Surface roughness in finite element meshes

Fabian Loth<EMAIL_ADDRESS>Thomas Kiel  Kurt Busch  Philip Trøst
Kristensen  Humboldt-Universität zu Berlin, Institut für Physik, AG
Theoretische Optik & Photonik, Newtonstraße 15, 12489 Berlin, Germany  Max-
Born-Institut, Max-Born-Straße 2A, 12489 Berlin, Germany

###### Abstract

We present a practical approach for constructing meshes of general rough
surfaces with given autocorrelation functions based on the unstructured meshes
of nominally smooth surfaces. The approach builds on a well-known method to
construct correlated random numbers from white noise using a decomposition of
the autocorrelation matrix. We discuss important details arising in practical
applications to the physicalmodeling of surface roughness and provide a
software implementation to enable use of the approach with a broad range of
numerical methods in various fields of science and engineering.

###### keywords:

surface roughness , finite element method , autocorrelation matrix , Cholesky
decomposition

††journal: Journal of Computational Physics

##  1 Introduction

Surface roughness plays an important role in physical and chemical phenomena
and is, therefore, of great importance in science and engineering [1]. Indeed,
in the relatively broad area of nanophotonics alone with which we have first-
hand modeling experience, it is known that the optical response of any
structure can, in general, be significantly changed by surface roughness. For
instance, this plays a role in light scattering by small particles [2],
influences the optical properties of plasmonic nanostructures [3, 4, 5], leads
to detrimental losses in photonic crystal waveguides [6], modifies the
performance of hyperbolic metamaterials [7] and affects the Casimir force [8].
Moreover, it is by now well recognized, that irregularities and protrusions on
general metallic surfaces can lead to the formation of hot-spots with enormous
optical field enhancements which, in turn, enables surface-enhanced Raman
spectroscopy [9] with applications ranging from biophysics [10] to the
detection of trace materials such as explosives [11]. Furthermore, surface
roughness influences the photon extraction of light emitting diodes [12] and
the quantum efficiency of solar cells [13]. Apart from applications in
nanophotonics, it is known that scattering of phonons at rough surfaces of
materials can reduce the thermal conductivity [14], just as surface roughness
has a decisive influence on qualitative as well as quantitative details of
fluid dynamics, in particular in the vast research field of micro- and
nanofluidics [15]. Lately, it was demonstrated that surface roughness can
increase the catalytic activity of materials [16] and theoretical studies show
that surface roughness affects the electron confinement in semiconductor
quantum dots [17]. As a last example of our non-exhaustive list, we note that
surface roughness plays a key role in the tribology of micro- and
nanoelectromechanical systems [18].

Practical characterizations and measurements of surface roughness at the
micro- and nanoscale can be performed with techniques such as atomic force
microscopy, confocal laser scanning microscopy, scanning interferometry and
scatterometry [1]. This data can then be used in the physical modeling of
experiments, for example, to asses if the measured or expected levels of
surface roughness can explain the deviations from the theoretical predictions
which are usually based on idealized smooth geometries. Similarly, the
characteristics can be used in the design of devices to model if the expected
or measured surface roughness is likely to degrade the performance of the
device. In practice, such modeling is almost exclusively done by numerical
means, and therefore it is of some interest to explore and discuss the
practical generation of finite element meshes for general rough surfaces with
given autocorrelation functions based on unstructured meshes of nominally
smooth surfaces.

Surface roughness is the deviation in normal direction of a real surface from
its nominal form and is typically characterized by the root mean squared (rms)
roughness $R_{\mathrm{q}}$ in the normal direction of the nominal form and the
correlation length $l$ in the lateral (tangential) direction of the nominal
form. Actually, a mesh with given rms roughness can be easily achieved by
simply moving points on the mesh of the nominal surface in the normal
direction by an amount set by uncorrelated random numbers. However, such a
naive approach produces roughness with a vanishing correlation length, and
this is far from a good model for the roughness of typical real surfaces. In
order to introduce a finite correlation length, one option is to perform a
convolution of the uncorrelated random numbers with a filter function of the
desired form [19]. This method, which is known as the spectral method [20], is
suitable for regular grids of rectangular surfaces and produces periodic
roughness. For instance, it has been used for the computation of the
absorption of light by rough metal surfaces [21]. As an alternative to the
spectral method, we discuss in this article how one can instead work directly
with the unstructured meshes of nominal structures. Such unstructured meshes
are typical for finite element techniques.

Specifically, we create correlated numbers by using a well-known method based
on decomposition of the autocorrelation matrix [22, 23]. This method is useful
for a relatively broad range of three-dimensional geometries, including cases
which cannot be constructed from a rectangle without distorting the distances
between the points on the surface. In particular, we apply it to generate
rough spheres, as illustrated in fig. 1, which shows the nominal mesh of a
sphere generated with Gmsh [24] along with a single realization of a
corresponding rough mesh.

Figure 1: Unstructured mesh of a sphere with radius $r$ generated with Gmsh
[24] and a rough sphere with rms roughness $R_{\mathrm{q}}/r=0.05$,
correlation length $l/r=0.15$ and element size $h/r=0.05$, which is called
_characteristic length_ in Gmsh.

This article is organized as follows. In section 2, we present the method for
constructing meshes of general rough surfaces, and section 3 contains details
of the method in order to use it in practice. Finally, section 4 holds the
conclusions.

##  2 Method

Given a nominally smooth surface, we describe the local deviation in the
normal direction as a stationary Gaussian stochastic process $D=D(t)$. Each
realization of the stochastic process yields an ordinary function $d(t)$,
where $t$ denotes the tangential coordinate along the surface. In the discrete
case considered in this work, each $t_{i}$ denotes the position of a vertex in
the nominal mesh and $d(t_{i})=d_{i}$ denotes the corresponding shift in
normal direction. The mean of the Gaussian process $\langle D\rangle=0$
vanishes, as it corresponds to the nominal surface. For simplicity we can set
the variance to one and multiply the deviation with the desired root mean
squared roughness $R_{\mathrm{q}}$ at the end of the calculation. As a
consequence, the autocovariance, the autocorrelation and the second moment are
equivalent for this stochastic process.

We start from the simplest stationary Gaussian stochastic process, which is
known as Gaussian white noise $W$ and which is determined by its
autocorrelation matrix

$\langle WW^{\mathrm{T}}\rangle=I\,,\quad I_{ij}=\delta_{ij}\,.$ (1)

Each realization $w$ is a vector of normally distributed uncorrelated random
numbers, which, in practice, can be conveniently generated with a pseudo-
random number generator, such as the Mersenne Twister [25] used in this work
[26]. One example realization in one dimension, which uses $w$ directly as
height deviation, is depicted in the top panel of fig. 2. It clearly shows
that Gaussian white noise is a rather bad model for most numerical models of
roughness. Not only is the lack of a finite correlation length in lateral
direction in principle inconsistent with reality, but the discrete realization
in practice enforces a correlation length equal to the mesh size. This latter
property fundamentally changes the convergence characteristics of any
numerical scheme relying on refinement of the calculation mesh to increase
accuracy.

Even if Gaussian white noise is a bad model for most, if not all, realistic
cases of roughness, it is a wonderful fact, that we can use it to construct
any other Gaussian random vector [23]. For the examples in this work, we
assume that the height deviation $D$ has a Gaussian autocorrelation matrix of
the form

$\langle DD^{\mathrm{T}}\rangle=R\,,\quad
R_{ij}=\exp\left(-\frac{|t_{i}-t_{j}|^{2}}{2l^{2}}\right)\,,$ (2)

but we note that other forms of correlation, e.g. exponential correlation
[27], can be treated immediately by simply inserting the desired functional
form for $R_{ij}$.

Figure 2: Discrete Gaussian white noise $w$ with $n=100$ points generated with
a random number generator (top) and the height deviation $d$ of a one-
dimensional rough surface created using a Gaussian autocorrelation with
correlation length $l/c=0.03$ (bottom). The tangential coordinate $t$ and the
amplitudes $w$ and $d$ are normalized to the total length $c$ and the root
mean squared roughness $R_{\mathrm{q}}$, respectively.

To construct correlated random numbers from Gaussian white noise, we seek a
matrix $L$ such that $D=LW$. Substitution into the autocorrelation yields:

$R=\langle DD^{\mathrm{T}}\rangle=\langle(LW)(LW)^{\mathrm{T}}\rangle=\langle
LWW^{\mathrm{T}}L^{\mathrm{T}}\rangle=L\langle WW^{\mathrm{T}}\rangle
L^{\mathrm{T}}=LL^{\mathrm{T}}\,,$ (3)

from which it follows, that any decomposition of the desired autocorrelation
matrix into a matrix $L$ and its transpose $L^{\mathrm{T}}$ will suffice. As
stated in the introduction, this is a well-known approach for generating
correlated numbers in general [22, 23], and we apply it here to generate
surface roughness in finite element meshes. The matrix $L$ can be obtained by
the Cholesky decomposition $R=LL^{\mathrm{T}}$, where $R$ and $L$ are,
respectively, symmetric positive-definite and lower triangular matrices [28].
We note that $R$ is symmetric by construction and for now we assume that it is
also positive-definite; this assumption is further discussed in section 3.2.
The method can be directly applied to unstructured meshes of general surfaces
by defining $|t_{i}-t_{j}|$ to be the distance between vertices $i$ and $j$.
It proceeds in eight distinct steps as follows:

  1. 1.

Generate a mesh of the nominal surface using a mesh generator.

  2. 2.

Compute the pairwise distances between all vertices.

  3. 3.

Construct the desired autocorrelation matrix $R$.

  4. 4.

Compute the Cholesky decomposition $R=LL^{\mathrm{T}}$.

  5. 5.

Generate Gaussian white noise $w$ using a random number generator.

  6. 6.

Calculate the height deviation $d=Lw$ with the desired autocorrelation.

  7. 7.

Multiply $d$ by the desired rms roughness to obtain
$\tilde{d}=R_{\mathrm{q}}d$.

  8. 8.

Shift all vertices in normal direction according to the elements of
$\tilde{d}$.

This approach, which is implemented in a collection of Python scripts and
published together with this article [29], was used to generate the rough
sphere in fig. 1. The effect of different roughness parameters is illustrated
in fig. 3. It shows, that the root mean squared roughness $R_{\mathrm{q}}$
controls the amplitude of the roughness while the correlation length $l$ is
proportional to the typical distance between bumps on the surface. Moreover,
it illustrates that one typically needs a finer mesh for shorter correlation
lengths.

Figure 3: Meshes of spheres with nominal radius $r$ and different roughness.
The rms roughness is doubled from $R_{\mathrm{q}}/r=0.04$ (left) to
$R_{\mathrm{q}}/r=0.08$ (right) and the correlation length is doubled from
$l/r=0.1$ (top) to $l/r=0.2$ (bottom). The element size is $h=l/2.5$.

##  3 Details

Even if the principles of the suggested approach are relatively simple, its
practical implementation calls for the discussion of a number of details.

###  3.1 Sampling

Describing surface roughness as a discrete stochastic process is an
approximation to the actual continuous rough surface. We examine the quality
of this approximation by means of the sampling theorem [30], which states that
a function which is bandwidth-limited to frequencies smaller than
$f_{\mathrm{max}}$ can be reconstructed exactly from its equidistantly sampled
values if the sampling frequency is larger than $2f_{\mathrm{max}}$.

By means of the Wiener-Khinchin theorem, the spectral density $S(f)$ of a
stationary stochastic process is given by the Fourier transform of the
autocorrelation function [31, 32]. Thus, for a Gaussian autocorrelation the
spectral density $S(f)$ is Gaussian and not strictly bandwidth-limited.
Nevertheless, we may define a frequency $f_{\mathrm{max}}$, such that $S(f)$
is negligibly small for $f>f_{\mathrm{max}}$. Then, for a sampling frequency
of $2f_{\mathrm{max}}$, the surface in principle can be reconstructed with
negligible aliasing using the Whittaker-Shannon interpolation formula [30]. In
this work, however, we are interested in the construction of meshes for use
with numerical methods such as finite element methods, wherefore we are
typically restricted to linear interpolations of the calculation domain into
triangular or tetrahedral elements. In this case the discretization in general
needs to be finer to represent high frequency fluctuations correctly, cf. fig.
3. Since a finer mesh increases the size of the numerical problem, one needs
to find a compromise between accuracy and computation time depending on the
problem at hand. Even if this is a general characteristic of the numerical
calculations, the introduction of roughness will typically lead to an inherent
uncertainty in the final results and thereby add a natural bound on the
calculation accuracy beyond which it does not make sense to refine the
calculation mesh further. Indeed, if one finds that the introduction of
roughness leads to a $10$ percent variation in the figure of merit, one can
safely perform calculations with an estimated error of one percent. In the
examples in figs. 1 and 3, the element size is chosen as $h=l/2.5$, where $l$
is the correlation length.

###  3.2 Eigenvalues of the autocorrelation matrix

The Cholesky decomposition of a symmetric matrix exists if and only if the
matrix is positive-definite meaning that all eigenvalues are strictly
positive. Thus, if the numerical Cholesky decomposition of the autocorrelation
matrix succeeds, then the matrix is positive-definite. Furthermore, Schoenberg
proved that the function $\exp(-|x|^{p})$ with $x\in\mathbb{R}$ is positive
definite if $0<p\leq 2$ and not positive definite if $p>2$ [33]. We have
found, however, that under certain conditions the numerical Cholesky
decomposition fails even for $p=2$, which corresponds to the Gaussian
autocorrelation given by eq. 2. We examine this phenomenon by computing the
eigenvalues of the Gaussian autocorrelation matrix with different correlation
lengths $l$ for a one-dimensional surface of length $c$ sampled with fixed
element size $h/c=0.01$, see fig. 4, left.

Figure 4: Logarithm of the computed eigenvalues $\lambda_{i}$ for the Gaussian
autocorrelation matrix with different correlation lengths $l$ in descending
order for an open (left) and a periodically closed (right) one-dimensional
surface of length $c$ and sampled with fixed element size $h/c=0.01$. For
correlation lengths $l/h\gtrsim 3$, points on the right side of the minimum
correspond to negative eigenvalues.

First, we notice that the symmetric eigenvalue problem is well-conditioned
[34] even if the matrix itself is ill-conditioned as it is the case for the
autocorrelation matrix. Because of Szegö’s theorem [35], we expect the
eigenvalues of the autocorrelation matrix to have the same distribution as the
spectral density. Indeed, the eigenvalue spectrum has a normal distribution
which shows a parabola in the logarithmic plot in fig. 4. As the eigenvalues
reach the machine precision $\epsilon\approx 10^{-16}$, however, all smaller
eigenvalues are distributed around zero with an absolute value on the order of
$\epsilon$. This happens for $l/h\gtrsim 3$, and in these cases, the Cholesky
decomposition fails. Instead, we compute the eigendecomposition

$R=Q\,\Lambda\,Q^{\mathrm{T}}\,,$ (4)

where $Q$ is orthogonal and $\Lambda$ is diagonal and contains the
eigenvalues. We take the positive part $\Lambda^{+}$ by setting all negative
eigenvalues to zero, and use $Q\sqrt{\Lambda^{+}}$ instead of the Cholesky
factor $L$.

An additional complication may arise in the case of periodic boundary
conditions, as seen in the right panel of fig. 4. In this case, the theorem by
Schoenberg does not apply, thus, it is not guaranteed that the autocorrelation
matrix is positive-definite. The relevant quantity for the following
discussion is the ratio of the correlation length $l$ and the circumference
$c$. For the cases $l/c\gtrsim 0.07$ (which corresponds in fig. 4 to
$l/h\gtrsim 7$), the Gaussian eigenvalue spectrum is cut off at a value
$\lambda_{\mathrm{cut}}\gg\epsilon$, after which the eigenvalues are
distributed around zero with absolute values of the order of
$\lambda_{\mathrm{cut}}$. This effect is visible in the spectra as a number of
distinct floors, at which the magnitude of the eigenvalues are approximately
constant. It turns out that this value is of the same order as the minimal
correlation, i.e. the correlation between most distant points in the mesh. To
avoid this problem, we need to chose the correlation length such that the
minimal correlation is smaller than the desired precision. Fortunately, this
is inherently unproblematic, as a strong correlation of all points would
essentially constitute a different geometry than that of a rough surface.
Clearly, these considerations hold also for two-dimensional surfaces, e.g. a
sphere (fig. 5). As a side note, we see that the eigenvalue spectrum has steps
of almost the same values. The width of these steps increases such that the
spectrum is linear in the logarithmic plot. This corresponds to the fact that
for each vertex, the distance to surrounding vertices is roughly an integer
multiple of the element size and the number of vertices with distance $d$
increases as $d$ increases itself.

Figure 5: Logarithm of the computed eigenvalues $\lambda_{i}$ of the Gaussian
autocorrelation matrix with different correlation length $l$ in descending
order for a spherical mesh with circumference $c$ and element size $h/c=0.01$.
For correlation lengths $l/h\gtrsim 3$, points on the right side of the
minimum correspond to negative eigenvalues.

###  3.3 Distance matrix

In order to construct the autocorrelation matrix one needs to compute the
pairwise distances between all vertices of the nominal surface mesh. On curved
surfaces, it is arguably reasonable to interpret this as the geodesic
distance, i.e. the shortest path along the surface. Analytic expressions for
the geodesic distance, however, are only available for simple geometries, such
as spheres. For other geometries, we use the Euclidean distance as an
approximation. This approximation is expected to be best in the limit of short
correlation lengths relative to the curvature of the surface. In particular,
sharp edges and corners should be avoided. To handle such geometries, one may
have to use an additional step to first round the corners, before applying the
roughness. We note, that the geodesic distance between points on an arbitrary
triangular mesh may be computed using an algorithm by Mitchell, Mount and
Papadimitrou [36] and variants of it [37]. However, it is much more costly
than the simple Euclidean distance and thus it was not used in our present
work.

##  4 Conclusion

We have presented a practical approach for constructing finite element meshes
of general rough surfaces with a desired autocorrelation based on a
decomposition of the autocorrelation matrix. In addition, we have discussed
important details of the method and especially explained and solved the
problem of negative eigenvalues of the autocorrelation matrix. The approach
can be directly used with the code provided together with this article [29].

## Acknowledgment

This work was supported by the German Federal Ministry of Education and
Research through the funding program Photonics Research Germany (Project
13N14149). F.L. and K.B. further acknowledge the support by Deutsche
Forschungsgemeinschaft (DFG) within the priority program SPP 1839 ”Tailored
Disorder” (project BU 1107/10-2).

## References

  * [1] Y. Gong, J. Xu, and R.C. Buchanan.  Surface roughness: A review of its measurement at micro-/nano-scale.  Physical Sciences Reviews, 3(1), 2018.
  * [2] C. Li, G.W. Kattawar, and P. Yang.  Effects of surface roughness on light scattering by small particles.  J. Quant. Spectrosc. Radiat. Transfer, 89(1-4):123–131, 2004.
  * [3] A. Trügler, J.-C. Tinguely, J.R. Krenn, A. Hohenau, and U. Hohenester.  Influence of surface roughness on the optical properties of plasmonic nanoparticles.  Phys. Rev. B, 83(8):081412, 2011.
  * [4] A. Trügler, J.-C. Tinguely, G. Jakopic, U. Hohenester, J.R. Krenn, and A. Hohenau.  Near-field and sers enhancement from rough plasmonic nanoparticles.  Phys. Rev. B, 89(16):165409, 2014.
  * [5] Y.-W. Lu, L.-Y. Li, and J.-F. Liu.  Influence of surface roughness on strong light-matter interaction of a quantum emitter-metallic nanoparticle system.  Sci. Rep., 8(1), 2018.
  * [6] S.G. Johnson, M.L. Povinelli, M. Soljacic, A. Karalis, S. Jacobs, and J.D. Joannopoulos.  Roughness losses and volume-current methods in photonic-crystal waveguides.  Applied Physics B, 81(2-3):283–293, 2005.
  * [7] S. Kozik, M.A. Binhussain, A. Smirnov, N. Khilo, and V. Agabekov.  Investigation of surface roughness influence on hyperbolic metamaterial performance.  Advanced Electromagnetics, 3(2), 2014.
  * [8] P. J. van Zwol, V. B. Svetovoy, and G. Palasantzas.  Characterization of optical properties and surface roughness profiles: The casimir force between real materials.  In Casimir Physics, pages 311–343. Springer Berlin Heidelberg, 2011\.
  * [9] M. Fleischmann, P.J. Hendra, and A.J. McQuillan.  Raman spectra of pyridine adsorbed at a silver electrode.  Chemical Physics Letters, 26(2):163 – 166, 1974.
  * [10] Katrin Kneipp, Harald Kneipp, Irving Itzkan, Ramachandra R Dasari, and Michael S Feld.  Surface-enhanced raman scattering and biophysics.  Journal of Physics: Condensed Matter, 14(18):R597–R624, 2002.
  * [11] Aron Hakonen, Per Ola Andersson, Michael Stenbæk Schmidt, Tomas Rindzevicius, and Mikael Käll.  Explosive and chemical threat detection by surface-enhanced raman scattering: A review.  Analytica Chimica Acta, 893:1 – 13, 2015.
  * [12] T. Fujii, Y. Gao, R. Sharma, E. L. Hu, S. P. DenBaars, and S. Nakamura.  Increase in the extraction efficiency of GaN-based light-emitting diodes via surface roughening.  Applied Physics Letters, 84(6):855–857, February 2004.
  * [13] J. Krč, M. Zeman, O. Kluth, F. Smole, and M. Topič.  Effect of surface roughness of ZnO:al films on light scattering in hydrogenated amorphous silicon solar cells.  Thin Solid Films, 426(1-2):296–304, February 2003.
  * [14] D. H. Santamore and M. C. Cross.  Effect of surface roughness on the universal thermal conductance.  Physical Review B, 63(18), April 2001.
  * [15] James B. Taylor, Andres L. Carrano, and Satish G. Kandlikar.  Characterization of the effect of surface roughness and texture on fluid flow: Past, present, and future (keynote).  In ASME 3rd International Conference on Microchannels and Minichannels, Parts A and B. ASME, 2005.
  * [16] Ikhyun Kim, Gisu Park, and Jae Jeong Na.  Experimental study of surface roughness effect on oxygen catalytic recombination.  International Journal of Heat and Mass Transfer, 138:916–922, August 2019.
  * [17] R. Macêdo, M. S. Sena, J. Costa e Silva, A. Chaves, and J. A. P. da Costa.  The role of surface roughness on the electron confinement in semiconductor quantum dots.  In Latin America Optics and Photonics Conference. OSA, 2012.
  * [18] Bharat Bhushan.  Micro/nanotribology of MEMS/NEMS materials and devices.  In Nanotribology and Nanomechanics, pages 1031–1089. Springer-Verlag, 2005.
  * [19] N. Garcia and E. Stoll.  Monte carlo calculation for electromagnetic-wave scattering from random rough surfaces.  Phys. Rev. Lett., 52(20):1798–1801, 1984.
  * [20] Karl F Warnick and Weng Cho Chew.  Numerical simulation methods for rough surface scattering.  Waves in Random Media, 11(1):R1–R30, January 2001.
  * [21] D. Bergström, J. Powell, and A. F.H. Kaplan.  The absorption of light by rough metal surfaces—a three-dimensional ray-tracing analysis.  Journal of Applied Physics, 103(10):103515, May 2008.
  * [22] H.F. Kaiser and K. Dickman.  Sample and population score matrices and sample correlation matrices from an arbitrary population correlation matrix.  Psychometrika, 27(2):179–182, 1962.
  * [23] Robert G. Gallager.  Stochastic Processes.  Cambridge University Press, December 2013.
  * [24] C. Geuzaine and J.-F. Remacle.  Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities.  International Journal for Numerical Methods in Engineering, 79(11):1309–1331, 2009.
  * [25] Makoto Matsumoto and Takuji Nishimura.  Mersenne twister: A 623-dimensionally equidistributed uniform pseudo-random number generator.  ACM Transactions on Modeling and Computer Simulation, 8(1):3–30, jan 1998.
  * [26] T. E. Oliphant.  Python for scientific computing.  Computing in Science Engineering, 9(3):10–20, May 2007.
  * [27] J.A. Ogilvy and J.R. Foster.  Rough surfaces: gaussian or exponential statistics?  Journal of Physics D: Applied Physics, 22(9):1243–1251, 1989.
  * [28] N. J. Higham.  Cholesky factorization.  WIREs Comp. Stat., 1:251–254, 2009.
  * [29] F. Loth.  roughmesh.  https://github.com/fabian-loth/roughmesh, 2020.
  * [30] C.E. Shannon.  Communication in the presence of noise.  Proc. IRE, 37(1):10–21, 1949.
  * [31] Norbert Wiener.  Generalized harmonic analysis.  Acta Mathematica, 55(0):117–258, 1930.
  * [32] A. Khintchine.  Korrelationstheorie der stationären stochastischen prozesse.  Mathematische Annalen, 109(1):604–615, December 1934.
  * [33] I.J. Schoenberg.  Metric spaces and positive definite functions.  Transactions of the American Mathematical Society, 44(3):522–522, March 1938.
  * [34] F. L. Bauer and C. T. Fike.  Norms and exclusion theorems.  Numerische Mathematik, 2(1):137–141, December 1960.
  * [35] T.K. Moon and W.C. Sterling.  Mathematical Methods and Algorithms for Signal Processing.  Prentice Hall, 2000.
  * [36] J.S.B. Mitchell, D.M. Mount, and C.H. Papadimitriou.  The discrete geodesic problem.  SIAM J. of Computing, 16(4):647–668, 1987.
  * [37] V. Surazhsky, T. Surazhsky, D. Kirsanov, S.J. Gortler, and H. Hoppe.  Fast exact and approximate geodesics on meshes.  ACM Trans. Graph., 24(3):553–560, 2005.


# Decoding trust: A reinforcement learning perspective

Guozhong Zheng  School of Physics and Information Technology, Shaanxi Normal
University, Xi’an 710061, P. R. China  Jiqiang Zhang  School of Physics and
Electronic-Electrical Engineering, Ningxia University, Yinchuan 750021, P. R.
China  Jing Zhang  School of Physics and Information Technology, Shaanxi
Normal University, Xi’an 710061, P. R. China  College of Information Science
and Technology, Donghua University, Shanghai 201620, 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

Behavioral experiments on the trust game have shown that trust and
trustworthiness are universal among human beings, contradicting the prediction
by assuming _Homo economicus_ in orthodox Economics. This means some mechanism
must be at work that favors their emergence. Most previous explanations
however need to resort to some factors based upon imitative learning, a simple
version of social learning. Here, we turn to the paradigm of reinforcement
learning, where individuals update their strategies by evaluating the long-
term return through accumulated experience. Specifically, we investigate the
trust game with the Q-learning algorithm, where each participant is associated
with two evolving Q-tables that guide one’s decision making as trustor and
trustee respectively. In the pairwise scenario, we reveal that high levels of
trust and trustworthiness emerge when individuals appreciate both their
historical experience and returns in the future. Mechanistically, the
evolution of the Q-tables shows a crossover that resembles human’s
psychological changes. We also provide the phase diagram for the game
parameters, where the boundary analysis is conducted. These findings are
robust when the scenario is extended to a latticed population. Our results
thus provide a natural explanation for the emergence of trust and
trustworthiness without external factors involved. More importantly, the
proposed paradigm shows the potential in deciphering many puzzles in human
behaviors.

##  I 1\. Introduction

Trust and trustworthiness are central components of our human civilization
Hardin (2002), especially required in dealing with many pressing threats such
as climate change, pandemic threats, international conflicts, energy crises
etc. As “a lubricant for social system” Arrow (1974), trust can facilitate
cooperation and contribute to the economic growth Zak and Knack (2001); Algan
and Cahuc (2013). A large-scale survey shows many of our human beings give a
positive answer when facing the statement “most people can be trusted”, though
a recent survey (see “Using Data to Understand Our World”) Ortiz-Ospina and
Roser (2016) shows that over the past four decades fewer people say they
“trust each other”, which is an alarming signal. As fundamental questions,
what is the mechanism for the emergence of trust and the associated
trustworthiness, and under what conditions they are likely to sustain, have
attracted considerable attention from different areas in the past decades.

Trust is the willingness of an agent (the trustor) to act in a way that
benefits the other (the trustee) with the expectation that the trustee will
return part of the profit afterwards. But, the trustor has no control over the
trustee’s action, the action of trust thus puts herself/himself in a
vulnerable position. Without any mechanism to enforce the returns, the trustee
tends to maximize her/his interest by simply walking away without
reciprocation; with the backward induction, the trustor should not trust the
trustee due to the expected loss of investment. That’s the prediction from the
assumption of  _Homo economicus_ in orthodox Economics Solow and Simon (1957);
Fukuyama (1995); Samuelson and Nordhaus (2005), where individuals are supposed
to be self-interested, their behaviors are guided by the maximization of
payoff.

A large number of laboratory experiments on the canonical model — the trust
game Berg  _et al._ (1995), however, reveals strikingly different observations
Berg  _et al._ (1995); Johnson and Mislin (2011). Across different countries
with varied experimental protocols, the trustors are willing to send an
average of $50\%$ of their initial endowments to trustees, and the trustees
return an average of $37\%$ of the earnings back to trustors. Despite the
variations in details, it turns out that we human beings remarkably prefer
trusting others and are trustworthy to a large extent.

Figure 1: Emergence of trust and trustworthiness. The color-coded stationary
fractions of the four strategies in the domain $(\gamma,\alpha)$ values with
both $\gamma$, $\alpha$ from 0 to 1. A large fraction of TR combination is
seen in the corner of large $\gamma$ and small $\alpha$. Each data is averaged
over 100 realizations. Other parameters: $\epsilon=0.01$, $g=3$, $x=1.0$, and
$\omega=0.5$.

Many efforts subsequently have been made to understand the above
discrepancies. Zucker Zucker (1986) systematically discussed three trust-
producing modes: trust is tied to the past or expected exchange, to social
characteristics, and to formal societal structures. Ref. Buchan  _et al._
(2004) explored the impact of culture and its relationship with indirect
reciprocity Nowak (2005). There are other explanations where the emergence of
trust and trustworthiness is attributed to further factors, such as social
awareness Mcnamara  _et al._ (2008), reputation Xia  _et al._ (2022),
information Manapat  _et al._ (2013), delayed Manapat and Rand (2012) or
partial information Rauwolf and Bryson (2018) etc. While network reciprocity
is found to potentially promote many altruistic behaviors, such as cooperation
Nowak (2006), fairness Page and Sigmund (2000), and honesty Capraro  _et al._
(2020), Ref. Kumar  _et al._ (2020) however reveals that its impact on the
evolution of trust and trustworthiness is marginal for whatever types of
underlying networks of the population. Note that, most of aforementioned
studies follow the imitation learning rule Roca  _et al._ (2009), such as the
Moran or Fermi rule Szabó and Tőke (1998), where individuals imitate the
strategies of their neighbors who may have higher payoffs. In essence,
imitative learning can be taken as a simple version of social learning Bandura
and Walters (1977), where individuals learn from others in their socio-
economic activities, through observation or instruction, which may or may not
involve direct experiences.

While social learning is ubiquitous in nature and our society, there is
however a different paradigm of learning that has been largely ignored — the
reinforcement learning (RL) Sutton and Barto (2018). As one of the three main
classes in machine learning algorithms, RL specializes in decision-making
based upon experience and has achieved tremendous success in many aspects of
science and technology, especially after the marriage to deep learning LeCun
_et al._ (2015). Yet, only recently, RL starts to be applied to the
evolutionary game theory to help understand the emergence of cooperation Zhang
_et al._ (2020a); Song  _et al._ (2022); Ding  _et al._ (2023), the minority
game Andrecut and Ali (2001); Zhang  _et al._ (2019), and other collective
behaviors in complex systems Zhang  _et al._ (2020b); Tomov  _et al._ (2021);
He  _et al._ (2022); Shi and Rong (2022).

In fact, the RL has a solid foundation in neuroscience and the existing
evidences have shown that the working area and the physiological processes for
RL is different from social learning, manifesting itself as a fundamentally
distinct learning paradigm Lee  _et al._ (2012); Rangel  _et al._ (2008);
Olsson  _et al._ (2020). In the social learning paradigm, the observations of
utility are first made, and then their decisions are made by utility
comparison with others, e.g. imitating the strategies of those peers who have
higher utilities. As a kind of value-based learning, individuals score
different actions in RL and action is chosen probabilistically based on these
scores. Here, the most important distinction is that social learning is based
upon the utility-comparison rule within its neighborhood, while each
individual with RL develops a unique policy by self-reflection. The rule
remains the same in social learning, but the policies in RL are coevolving
with their surroundings and could be unique for each player. This makes the
two learning paradigms fundamentally different. Within the paradigm of RL, we
are interested in the following questions: _Is such endogenous self-reflective
learning way capable of providing an new paradigm to decipher the emergence of
trust and trustworthiness?_ It could provide fundamentally different insights
from the current social learning paradigm by exogenous comparison with peers.

In this work, we investigate the evolution of trust and trustworthiness within
the paradigm of reinforcement learning. Specifically, we adopt a Q-learning
algorithm Watkins (1989, 1992) to study the trust game, where each person is
guided by two Q-tables, respectively for the role of trustor and trustee. In
the two-person scenario, we find that a high level of trust and
trustworthiness emerges when individuals both respect the historical
experience and have a long-term vision. Similar phenomena are observed when
the two-person scenario is extended to a population level. The analysis of the
Q-tables reveals underlying mechanism behind where the individuals’
preferences change at different stages. We also compute the phase diagram
within the associated game parameters, and several boundaries are identified
for the onset of the emergence and for regions of different prevalences.

The paper is organized as follows: we introduce our Q-learning setup for the
trust game in Sec. 2II. In Sec. 3III, we show results for the two-person
scenario, and provide a mechanistic analysis. The impact of gain factor on
trust evolution is studied, and the boundary analysis is conducted in Sec.
4IV. In Sec. 5E, we extend our study of the evolution of trust within a
population on a 1d lattice. Finally, we conclude our work together with
discussions in Sec. 6VI.

##  II 2\. Methods

We start by introducing the trust game for two-player scenario, where two
agents act as a trustor and a trustee in turn. Initially, the trustor is
endowed with one monetary unit, who can choose either to stay the status quo
or to invest the trustee with trust by transferring an investment fraction
$x\in(0,1]$. With the status quo, denoted as not trust (N), the game is over.
With trust (T), the investment is multiplied by a gain factor $g$, the game
enters into the second stage. The trustee has to make a choice between
reciprocity (R) and betrayal (B). In the former, the trustee transfers a
return fraction of its earnings $\omega\in(0,1]$ to the trustor as return,
i.e. $\omega gx$. The choice of betrayal means no return, the trustee walks
away without reciprocity. At the end of the second session, the final payoffs
for the trustor and trustee are, respectively, $1-x+\omega gx$ and
$(1-\omega)gx$ if R is chosen, whereas the payoffs correspond to $1-x$ and
$gx$ if B is selected. Here, the investment of the trustor is a manifestation
of trust, and the transfer back to the trustor can be interpreted as
trustworthiness. Table 1 shows the payoff matrix for the two agents.

$i$ (trustor) $j$ (trustee) | R |  B
---|---|---
T | ${1+(g\omega-1)x,(1-\omega)gx}$ |  ${1-x,gx}$
N | ${1,0}$ |  ${1,0}$

Table 1: The payoff matrix in the trust game. Within each item, the first
payoff is for the trustor (row player) and the second for the trustee (column
player).

In the one-shot anonymous scenario, it’s obvious that the trustee prefers to
walk away, and betrayal is the reasonable choice to maximize its payoff.
Similarly, the trustor is supposed not to trust as no return is expected.
Therefore, the rational solution for the trust game ends with the solution (N,
B). Notice that, from the payoff structure in Table 1, the trust game bears a
strong resemblance to the Prisoner’s Dilemma, but the two are fundamentally
different because the two agents in the latter are symmetrical, while they
play different roles and act in sequence in the trust game.

For simplicity, we consider the temporal evolution of this two-player scenario
where only two individuals are engaged, and they play the role of trustor and
trustee in turn. Specifically, they adopt a Q-learning algorithm Watkins
(1989, 1992), a model-free reinforcement learning. In this algorithm, each
player has two Q-tables denoted as $Q^{A}_{s_{n}}$ and $Q^{B}_{s_{n}}$ in
hands that guide their decision-making for the role of trustor and trustee,
respectively, see Table 2. One takes an action from the action set $\\{T,N\\}$
as a trustor or from $\\{R,B\\}$ as a trustee. The system is in one of four
states denoted as $\mathbb{S}=\\{s_{1},...,s_{4}\\}$ with
$s_{1}=TR,s_{2}=TB,s_{3}=NR,s_{4}=NB$ as shown in $Q^{A}_{s_{n}}$, whereas the
states in $Q^{B}_{s_{n}}$ are the same but just reshuffled,
$s_{1}=RT,s_{2}=RN,s_{3}=BT,s_{4}=BN$. The elements $Q_{s,a}$ in the Q-tables
is a value function used to measure the value of the corresponding action in a
given state, which is constantly updated over time. Q-learning aims to finds
the optimal policy in the sense of maximizing the expected accumulated reward
starting from a given state.

Without loss of generality, the two individuals are initially assigned with a
random strategy from $\\{$TR, TB, NR, NB$\\}$ with equal probabilities.
Following the common practice, the two roles are switched in turn for the two
agents from round to round, noting that the way of randomly assigned role does
not change our findings. In each round $t$, with a probability $\epsilon$, the
two agents randomly choose an action from the corresponding action set to
conduct a trial-and-error exploration; With $1-\epsilon$, they choose the
action with the larger Q value in the associated Q-table to maximize their
rewards. This completes the stage of decision-making for their actions. The
corresponding element in the Q-table is updated as follows

$\displaystyle Q_{s,a}(t+1)$
$\displaystyle=Q_{s,a}(t)+\alpha\left(R+\gamma\max_{a^{\prime}}Q_{s^{\prime},a^{\prime}}(t)-Q_{s,a}(t)\right)$
(1)
$\displaystyle=(1-\alpha)Q_{s,a}(t)+\alpha\left(R+\gamma\max_{a^{\prime}}Q_{s^{\prime},a^{\prime}}(t)\right),$

where $s,a$ represents the current state and action of the focal individual,
and $s^{\prime}$ is the new state at $t+1$. $\alpha\in(0,1]$ is the learning
rate, which captures the contribution of current step. A larger value of
$\alpha$ means that the agent is more forgetful as old Q values tend to be
more rapidly modified. $\gamma\in[0,1]$ is the discount factor, measuring the
weight of future rewards, as $\max_{a^{\prime}}Q_{s^{\prime},a^{\prime}}(t)$
is the maximal value expected in the new state. This completes the stage of
Q-table updating, and a single round is done. The evolution protocol is
summarized in Fig. 8 in Appendix AA for clarity.

State Action | T ($a_{1}$) | N ($a_{2}$) |  State Action | R ($a_{1}$) | B ($a_{2}$)
---|---|---|---|---|---
(T,R) ($s_{1}$) | $Q_{s_{1}a_{1}}$ | $Q_{s_{1}a_{2}}$ |  (R,T) ($s_{1}$) | $Q_{s_{1}a_{1}}$ | $Q_{s_{1}a_{2}}$
(T,B) ($s_{2}$) | $Q_{s_{2}a_{1}}$ | $Q_{s_{2}a_{2}}$ |  (R,N) ($s_{2}$) | $Q_{s_{2}a_{1}}$ | $Q_{s_{2}a_{2}}$
(N,R) ($s_{3}$) | $Q_{s_{3}a_{1}}$ | $Q_{s_{3}a_{2}}$ |  (B,T) ($s_{3}$) | $Q_{s_{3}a_{1}}$ | $Q_{s_{3}a_{2}}$
(N,B) ($s_{4}$) | $Q_{s_{4}a_{1}}$ | $Q_{s_{4}a_{2}}$ |  (B,N) ($s_{4}$) | $Q_{s_{4}a_{1}}$ | $Q_{s_{4}a_{2}}$

Table 2: The two Q-tables for each individual in the 2-player scenario, the
left $Q^{A}_{s_{n}}$ and right $Q^{B}_{s_{n}}$ are Q-tables respectively for
the role of trustor and trustee.

Figure 2:  Crossover in time series. Typical time series in the two-agent
scenario for the learning parameter combination $\alpha=0.1$ and $\gamma=0.9$.
(a) Four time series at the individual level; (b) Four time series at the
pairwise level; (c) The evolution of net payoffs to the trustor; (d) The
evolution of payoffs corresponding to different actions chosen by the trustee.
The vertical dashed lines mark the same transition moment where the level of
trust and trustworthiness starts to rise. Each data is averaged over 500
realizations in (a-d), besides a sliding window average of 60 steps is
conducted in (c-d). Other parameters: $\epsilon=0.01$, $g=3$, $x=1.0$, and
$w=0.5$.

In our practice, we focus on the case with fixed the game parameters $g=3$,
$x=1.0$, and $\omega=0.5$ Kumar  _et al._ (2020), this corresponds to such a
scenario: the trustor transfers all the money to trustee, and after the money
is tripled, half of the profit is returned to the trustor. Unless stated
otherwise, the exploration probability $\epsilon$ is chosen to be 0.01, to
balance the random exploration and the exploitation of the Q-table. Note that,
we adopt a time-varying form for both learning rate and exploration rate at
the beginning Auer  _et al._ (2002), i.e. $\epsilon=\frac{4.0}{\sqrt{t}}$,
$\alpha=\frac{6.0}{\sqrt{t}}$. This means that the two values starts with a
large value, but as time goes by they get fixed once approach the desired
value (will be set in the following studies), see Fig. 3(c). Though this setup
does change the evolutionary outcome in the long term compare to the case for
all-fixed parameters. Details see Appendix AB. For each case, the simulation
is run over $10^{8}$ time steps to guarantee the evolution reaches a stable
state.

##  III 3\. Results

###  III.1 3.1. Emergence of trust and trustworthiness

We find that a prominent emergence of both trust and trustworthiness is
observed in a physically meaningful region, see the phase diagram Fig. 1. It
reports the the final fractions of the four strategies $\\{$TR, TB, NR,
NB$\\}$ in the learning parameters domain ($\gamma$, $\alpha$). As shown,
there is a red region where the fraction of TR dominates ($f_{TR}>0.6$), which
happens when the learning rate $\alpha$ is small and the discount factor
$\gamma$ is large. This observation indicates that when individuals focus on
both historical experience and the long-term vision, both levels of trust and
trustworthiness rise. Otherwise, either the strongly forgetful property (a
large $\alpha$) and/or the short-term vision (a small $\gamma$) leads to the
failure of their emergence. Notice that, the fraction of the strategy TB
remains small across the whole domain, meaning that once the trust is adopted,
the agent also shows trustworthiness when acting in the role of trustee —
trustors never betray. But once non-trust is chosen, the reciprocity and
betrayal are equally likely to be chosen.

Figure 3:  Evolution of Q-tables. (a) The evolution of $\Delta Q^{B}_{s_{n}}$
for the trustee and (b) $\Delta Q^{A}_{s_{n}}$ for the trustor in all four
possible states. If $\Delta Q^{A}_{s_{n}}>0$, the action of T is preferred
within the state $s_{n}$ when in the role of trustor, otherwise N is flavored.
Similarly, if $\Delta Q^{B}_{s_{n}}>0$ , the action R is flavored when acting
as a trustee, otherwise B is preferred. The vertical dashed lines in (a) and
(b) seperate the two stages, NB is preferred in the first stage, but is then
followed by a turnover, where a positive feedback is formed to promote both
the trust and trustworthiness prevalences. (c) The two curves $\Delta
Q^{A}_{TR}$ in (a) and $\Delta Q^{B}_{RT}$ in (b) are put together for
clarity, along with the exploration rate $\epsilon$ and the learning rate
$\alpha$. The vertical line in (c) marks the approximate transitions at which
the individuals continue to choose to trust and show trustworthiness. Each
data is averaged over 500 realizations. Parameters: $\epsilon=0.01$, $g=3$,
$\alpha=0.1$, $\gamma=0.9$. $x=1.0$, $w=0.5$.

To understand the emergence of trust and trustworthiness, we show the time
evolution of the four fractions for the case of $\alpha=0.1$ and $\gamma=0.9$,
a typical combination of parameters located in the red region in Fig. 1, see
Fig. 2(a). We can see that their fractions start from around 0.25 due to
random initialization, and then the two fractions of distrust (NR and NB)
rise, but as time goes by these two fractions turn down, all three fractions
except the strategy TR (the red line) decline, the fraction of the strategy NB
(the black line) almost vanishes, and the system becomes stable in the end.

The reason for the initial decline in trust is straightforward since the level
of reciprocity is low at the beginning (even below $50\%$), putting the
trustor at a disadvantageous position, see Fig. 2(b). The disadvantage of
trustors is more clearly seen in Fig. 2(c), where their net earnings are
negative, being in a loss state. At the same time, the trustee is inclined to
betrayal because as long as the trustor is willing to invest, betrayal tends
to yield higher short-term payoffs compare to reciprocity. This is in line
with the prediction from _Homo economicus_ perspective. This leads to a slight
increase in the preference of betrayal, and consequently enhances the trend in
the preference decline of trust.

Unexpectedly, the fraction declines in trust ceases after around a hundred
rounds and starts to turn up. This crossover can be understood by monitoring
simultaneously the fractions of the two players’ strategies [Fig. 2(b)]. As
can be seen, after the transient, the trustor tends not to invest in a trustee
who has betrayed in the last round, but learns to invest to who showed the
reciprocity, which explains the fraction of TB declines but TR starts to rise.
This observation is also explained in Fig. 2(d), where the advantage in the
payoff by choosing B over R is lost at the moment as indicated by the dashed
line. Intuitively, as the trustee learns that the action of betrayal incurs no
investment, one starts to reciprocate the trustor instead of walking away.
Once this trend starts, a higher level of reciprocity is preferred for the
trustee, which in turn enhances the preference in action T for the trustor.
This then forms a positive feedback that finally yields a high preference in
both trust and trustworthiness.

###  III.2 3.2. Evolution of Q-tables

For a deeper understanding of the mechanism, we direct our attention to the
evolution of the two Q-tables. It shows the preference in NB like a _Homo
economicus_ is only present at the very beginning, the evolution in the later
stage aiming for maximizing payoffs in the long term forces them to turn to
TR, which resembles human’s psychological changes.

Specifically, we focus on the relative magnitude of each row of the two
Q-tables, represented as $\Delta
Q^{A,B}_{s_{n}}=Q^{A,B}_{s_{n}a_{1}}-Q^{A,B}_{s_{n}a_{2}}$, which determines
the preferred action within the given state $s_{n}$ according to the idea of
Q-learning. For example, if $\Delta Q^{A}_{s_{n}}\textgreater 0$, this means
that the action T is considered to be the optimal choice within the state
$s_{n}$ when acting as a trustor; otherwise N is supposed to be a better
choice. Likewise, $\Delta Q^{B}_{s_{n}}\textgreater 0$ means that action R is
the optimal choice when playing as a trustee, otherwise B is preferred. In our
study, the evolution of the $\Delta Q^{A,B}_{s_{n}}$ values for both
individuals are found the same, as their learning parameters are identical,
implicating they have quite similar cognitive processes. Therefore, we only
focus on the evolution of the $\Delta Q^{A,B}_{s_{n}}$ values for one of two
individuals, as illustrated in Fig. 3.

Fig. 3(a) and 3(b) show respectively the time evolution of $\Delta
Q^{B}_{s_{n}}$ and $\Delta Q^{A}_{s_{n}}$ for all four states. At the initial
stage ($t\lesssim 100$, before the marked dashed line), $\Delta
Q^{A,B}_{s_{n}}$ in both Q-tables mostly become more negative by learning, and
the strategy for the individual converges to NB, where as a trustor one is
unwilling to invest, and as a trustee one opts to betray. This is reasonable
because betrayal for the trustee is better off than reciprocity in the short-
term, since no investment can avoid potential money loss. As a result, the
dominating state in the system is NB for a trustor, and BN for a trustee, both
players act indeed as a _Homo economicus_.

As time comes to $t\approx 100$, the advantage of betrayal over reciprocity
diminishes, because no investment comes from the trustor. This then causes a
reversal of $\Delta Q^{B}_{BN}$, and the action of R is then preferred. This
critical transition, however, does not immediately leads to the boom of trust
or trustworthiness, because all four $\Delta Q^{A}_{s_{n}}$ at the moment are
averagely all negative, and distrust is still dominating. Actually, the action
of reciprocity is unstable, since for $\Delta Q^{B}_{BN}>0$ and $\Delta
Q^{B}_{RN}<0$, the state for the trustee oscillates between RN and BN.
Therefore, at this stage, still no driving force for either trust or
trustworthiness is seen, as shown in Fig. 3(c), where TR is unstable for
either trustor or trustee since the associated $\Delta Q$ are both negative at
the left to the vertical dashed line.

An important change unfolds afterwards (i.e., $t\gtrsim 100$), as can be seen
in Fig. 3(c). There, as the learning rate $\alpha$ and the exploration
$\epsilon$ decreases, some optimal or nearly optimal policies have been learnt
by individuals, and they rely more on their historical experience and conduct
less random explorations. As shown, the trustee becomes gradually inclined to
choose reciprocity rather than betrayal, as $\Delta Q$ for RT stops decreasing
and starts to grow. When this turnover is detected, the trustor also starts to
trust as a response, where a turnover is also present for $\Delta Q$ for TR,
but with a time delay.

Once both players turn to TR, this forms a positive feedback loop, both are
well paid, leading to an increasing $\Delta Q$ that are both positive in the
end. When the system becomes stable in the long run ($t\gtrsim 10^{6}$), the
trustee always chooses R since all its $\Delta Q>0$ [Fig. 3(a)]. On the
trustor side, $\Delta Q>0$ for the action of reciprocity (i.e. for $s_{n}=$ TR
and NR), but the trustor chooses not to trust ($\Delta Q<0$ for $s_{n}=$ TB
and NB) when the trustee betrayed [Fig. 3(b)]. This punishment-like policy
further forces the occasional betrayals back to the reciprocity. This then
results in stable emergence of trust, where a decent level of trust and
trustworthiness is seen.

Note that the final rise of $\Delta Q$ for both Q-tables also benefits from
the diminishing exploration rate [Fig. 3(c)], where both trustor and trustee
choose TR with a large probability. In a noisy scenario (a large $\epsilon$),
which can be interpreted as many misunderstandings or “trembling hands”
Fundenberg and Maskin (1990), such misunderstanding can considerably suppress
the level of the trust and trustworthiness, see Appendix C.

The above analysis shows that the emergence of trust and trustworthiness is
caused by the preference transition from NB to TR. To further confirm the
analysis, we compute the joint probability for two consecutive states
$P({s_{t},s_{t+1}})$, all state transitions at different stages are shown in
Fig. 4. In the early stage of evolution ($0\\!<\\!t\\!<\\!1000$), almost all
mode transitions are detected, however, some modes are more likely to happen,
such RN-BN, BN-RN, BN-BN, BN-RT, RT-BN, as shown in Fig. 4(a). This means that
the state BN is the main state at this stage, in line with our above argument.
As analyzed above, once the trust and trustee start to form a positive
feedback loop, the strategy of TR starts to dominate, as can be seen in Fig.
4(b). Apart from the dominating TR-TR mode, the other two bars BN-RT and RT-BN
are also present, meaning that the strategy pair RT flips to BN from time to
time, but the flipping back to RT is equally likely, as the two bars are of
nearly the same height. In Fig. 4(c-d), the mode of RT-RT becomes even higher,
though the other two bars (i.e. BN-RT and RT-BN) are still present due to the
non-vanishing $\epsilon$. Till then, the evolution of trust and
trustworthiness becomes stable.

Figure 4:  Joint probabilities for two consecutive states. The joint
probabilities $P({s_{t},s_{t+1}})$ of all sixteen state transitions in four
typical stages: (a) $0\sim 1000$, (b) $1.0\times 10^{3}\sim 1.0\times 10^{5}$,
(c) $1.0\times 10^{5}\sim 1.0\times 10^{6}$, (d) $2.0\times 10^{6}\sim
3.0\times 10^{6}$ time steps. In the end, the system mainly stays in TR state.
Each data is averaged over 500 realizations. Parameters: $\epsilon=0.01$,
$g=3$, $\alpha=0.1$, $\gamma=0.9$. $x=1.0$, and $w=0.5$.

Based on the above analysis, the emergence of trust can be roughly divided
into three stages:

  * 1)

Initially, the trustee finds that betrayal is more profitable than
reciprocating the trustor, and is thus inclined to choose B. Betrayal
gradually becomes prevalence. As a consequence, the trustor chooses to be non
trusting as the net earning of investing is negative. The level of trust and
trustworthiness both decrease.

  * 2)

As less investment is detected, the tendency towards betrayal was reversed,
the trustee starts to be inclined to reciprocating the investment. Once this
turnover is on site, the trustor also starts to invest.

  * 3)

The preference change in TR for the two players then forms a positive feedback
that strengthens the advantage of T and R in their Q-tables. In the end, the
trustee prefers R in almost all scenarios, and the trustor trusts those
reciprocating trustees but invests no money to the betrayed trustee as a
punishment. This guarantees a decent level of trust and trustworthiness.

However, as the two learning parameters ($\gamma$, $\alpha$) deviate from the
red region in Fig. 1(a), the mechanism behind stages 2) and 3) could be
ruined. Analyses based on the evolution of Q-tables indicate that for a large
$\alpha$, the past experience is rapidly washed out, so that no lesson can be
drawn from the history. In this case, the evolution of the game degrades to
the classic iterated scenario in the absence of Q-learning, where the trustor
is not willing to invest, and trust and trustworthiness fail to emerge Kumar
_et al._ (2020). Meanwhile, for the case of small $\gamma$, the positive
feedback between T and R fails to establish without confidence of the future
reward. For a detailed analysis see Appendix DD.

##  IV 4\. Impact of gain factor and boundary analysis

The revealed mechanism is robust against the game parameter $x$, but shows
intricate dependence on the return fraction $\omega$ and the gain factor $g$.
In the original model Berg  _et al._ (1995), the investment is tripled on the
trustee side, while in some other work the gain factor $g$ is set to be 2 or
other values Gereke  _et al._ (2018); Woolley and Fishbach (2017). To
systematically investigate the impact of the gain factor $g$ on the emergence
of trust and trustworthiness, here we first show the dependence of the TR
fraction on the investment fraction $x$ and the return fraction $\omega$ for
$g=2,3,4,5$, shown in Fig. 5. As can be seen, the fractions of TR show no any
dependence on the investment fraction $x$ in all four cases, which is
reasonable since the investment fraction only affects the absolute payoff, but
not the relative values in the Q-tables and thus the level of trust and
trustworthiness shows no dependence on $x$. However, a larger gain factor
widens the region where the trust and the trustworthiness emerge. We find that
there is a simple relationship determining the left boundary. For a trustor,
the bottom line to invest is not to lose money, conditioned by $gx\omega\geq
x$, leading to the following inequality

$\displaystyle g\omega\geq 1.$ (2)

This immediately gives the left boundaries in Fig. 5, i.e.
$\omega_{c}^{(1)}=1/g$, which is well confirmed and explains the independence
on $x$.

This theoretical argument is better validated in Fig. 6 within the parameter
domain $\omega-g$, with a fixed $x$ due to its irrelevance. We see that the
hyperbolic boundary fits perfectly the simulation results within a wide range
of the gain factor $g$. A closer lookup shows that the boundary slightly
shifts to the right as $g$ increases compared with the theoretic prediction.
This implies that when the gain becomes so high, the trustor would take some
risk of losing money to invest the trustee. Actually, the relationship
revealed in Eq. (2) is supported by previous experiments Evans and Krueger
(2011); Lenton and Mosley (2011); Johnson and Mislin (2011), where they found
that as the gain factor increases, the promotion in trust is seen.

Interestingly, the right boundary seen in Fig. 5 seems independent on the gain
factor $g$, where $\omega_{c}^{(2)}\approx 0.8$. This means that for a
trustee, one would only show trustworthiness only if she can keep at least
$20\%$ of the pie, otherwise one just walks away even if the gain factor is
large enough that a positive earning is expected. This observation implies
that the emergence of trustworthiness is built upon fairness, individuals in
the game seemingly also desire a fair division. In fact, as we can further see
in Fig. 6, the ideal division $\omega^{(3)}=0.5$ seemingly marks as another
threshold, below which a considerably high level of TR is then possible,
especially for a large gain factor. This means that a high level of
reciprocity further requires that the trustee can keep at least half of the
pie.

Figure 5:  Impact of gain factor. The color-coded fraction of the strategy TR
in the parameter domain of $x-\omega$, where the gain factor $g=2,3,4$ and 5,
respectively in (a-d). Each data is averaged $2\times 10^{4}$ times after a
transient $2\times 10^{7}$. Other parameters: $\epsilon=0.01$, $\alpha=0.1$,
and $\gamma=0.9$.  Figure 6:  Phase diagram and boundaries. The color-coded
fraction of the strategy TR in the parameter domain of $\omega-g$, where
$x=0.5$. Three suggested thresholds (dashed line) are respectively
corresponding to $\omega^{(1,3,2)}=1/g,0.5$, and 0.8 from left to right. Each
data is averaged over 100 realizations. Other parameters: $\epsilon=0.01$,
$\alpha=0.1$, and $\gamma=0.9$.

##  V 5\. 1-dimensional Lattice

In fact, the findings are not restricted to the above 2-player scenario, they
are robust and can also be seen at the population case. An example of
1-dimensional lattice with the size $N=50$ and $k=2$ is shown in Fig. 7. Since
there are two nearest neighbors for each individual in this scenario, the
Q-table has to be expanded accordingly, detailed settings can be seen in
Appendix EE. We find a qualitatively similar phenomenon that trust and
trustworthiness tend to arise when both historical experiences and the long-
term vision are emphasized, and the fraction of TR stabilizes at around 0.7 in
the long run. We can see that at the early stage, the nontrusting strategy is
dominating against their two nearest neighbors when they act as a trustor ,
but as time goes by, they gradually tend to invest. Also, the level of
reciprocity is low at the beginning but rises to be high later on. These
observations are similar to the two-player scenario. Notice that, due to the
presence of more neighbors, a smaller exploration rate $\epsilon$ and a higher
expectation of future reward $\gamma$ are preferred to maintain a stable
surrounding of a high level of TR, compared with the two-player scenario.

Figure 7:  Spatiotemporal evolution of 1d latticed population. (a) and (b)
correspond to scenarios where the individual $i$ acts as a trustor or a
trustee, respectively. In each scenario, for actions are possible, e.g., the
strategy (N,T) in (a) means as a trustor, the player chooses not to trust the
left nearest neighbor, but trust the right nearest neighbor. Parameters:
$N=50$, $k=2$, $\epsilon=0.001$, $\alpha=0.1$, and $\gamma=0.98$.

Actually for a pair of players, they act according to the their two Q-tables
regarding the other person, relying only on the agreement reached between
them. Therefore, even when an individual has multiple neighbors, trust in one
neighbor does not imply trust in others, i.e. the occurrence of trust and
trustworthiness is by nature pairwise between the two involved individuals,
independent of the rest. As a result, when extended to other complex networked
populations, we expect that the mechanism and phenomena are similar to the
2-individual scenario, which remain to be validated in the future.

##  VI 6\. Discussion

In summary, we have investigated the trust game within the paradigm of
reinforcement learning, each player acts following a Q-learning algorithm, and
we focus on the evolution of trust and trustworthiness. Surprisingly, high
levels the trust and trustworthiness emerges in the two-player scenario when
players both care about the historical experience and have long-term vision.
The evolution of the associated two Q-tables reveals a crossover in the action
preference, which resembles the psychological transition when we human beings
playing the game. Our boundary analysis shows that a high level of trust and
trustworthiness requires that the net earnings for the investment for the
trustor is positive and the trustee can keep half of the earnings in hand.
Furthermore, if the action choice deviates much from learnt Q-tables, this
“trembling hand” effect undermines the evolution of trust and trustworthiness,
where the desired relationships are broken down. Finally, the emergence of
trust and trustworthiness can also be seen when the scenario is changed into a
latticed population.

Interestingly, part of these observations were also seen in a series of
experiments by Engle-Warnick and Slonim Engle-Warnick and Slonim (2004, 2006,
2001). They explored the evolution of strategies in repeated trust game
experiments, and revealed that experience, attitude toward the future,
institutions, and other factors have an important influence on the strategy
selection in the repeated trust game. The results confirm that concerns for
the future of repeated interactions and past experiences in game history are
important for the persistence of trust.

Most importantly, we do not resort to any external factor as assumed in most
previous work with social learning. Our reinforcement learning provides a
natural explanation for the emergence of trust and trustworthiness. In fact,
existing efforts within this paradigm show that it can also provide
explanations for understanding cooperation Zhang  _et al._ (2020a); Song  _et
al._ (2022); Ding  _et al._ (2023), resource coordination Andrecut and Ali
(2001); Zhang  _et al._ (2019) etc. These work suggest that reinforcement
learning demonstrates its power in explaining human behaviors, where the
evolution of Q-table provides a uniform mechanism framework. Given the
consistent experimental deviation from the predictions of _Homo economicus_
regarding different altruistic behaviors Camerer and Thaler (1995); Camerer
_et al._ (2004); Henrich  _et al._ (2005); Johnson and Mislin (2011), such as
cooperation, fairness, trust and trustworthiness, the reinforcement learning
may provide a uniform paradigm to decipher complexities of human psychology,
shedding new insights into the understanding of moral behaviors Capraro and
Perc (2018, 2021).

Although we adopt reinforcement learning as our framework, we do not deny the
value of social learning paradigm that has been widely used in previous
studies. In fact, the two paradigms are not contradictory, but complementary
to each other. Till now, there are plenty of experimental evidences in
neuroscience that the decision making for both learning ways have solid neural
bases Lee  _et al._ (2012); Rangel  _et al._ (2008); Olsson  _et al._ (2020),
and indicate that they are suitable for different scenarios. The learning
processes in the real world are probably a mixture of the two. An important
question in this direction is to infer the type of learning from the
behavioral experiments. Only when the learning paradigms in realities are
clarified, we are on the right track to understand many important issues the
human society such as cooperation, fairness, honesty, and so on.

##  VII Acknowledgments

This work was supported by the National Natural Science Foundation of China
[Grants Nos. 12075144,12165014]. ZGZ is supported by Excellent Graduate
Training Program of Shaanxi Normal University [Grants No. LHRCTS23064].

##  Appendix A Appendix A: Our protocol of Q-learning

For clarity, our Q-learning setup of trust game is summarized in Fig. 8, which
is a synchronous updating protocol. (1) Initialize all Q-tables with small
random number within $(0,1)$ to mimic the unawareness of agents to the
environment, and also a random strategy within $\\{TR,TB,NR,NB\\}$ for each
agent. (2) In the game process, each agent chooses an action by either
exploitation or exploration, afterwards their rewards are obtained by
collecting payoffs. (3) The learning process is to update their Q-tables, and
their states also need to be updated. Steps (2) and (3) iterate until the
system reaches statistically stable or the desired time duration.

Figure 8:  The flowchart of protocol of the trust game.

##  Appendix B Appendix B: Setup of exploration rate and learning rate

To speed up the evolution, we adopt a time-varying form for the exploration
rate $\epsilon$ and learning rate $\alpha$ at the initial stage, which is a
common practice Auer  _et al._ (2002). To be specific,
$\epsilon=\frac{4.0}{\sqrt{t}}$, $\alpha=\frac{6.0}{\sqrt{t}}$, and remained
unchanged after dropping to the desired values (e.g. 0.01 and 0.1 in Fig 2,
respectively). The logic behind $\alpha=\frac{6.0}{\sqrt{t}}$ can be
interpreted as that the naive individuals at the beginning put more weight on
the future reward since they know that they have almost no experience
initially; but as time goes by, as more experiences have been learnt, they put
more weight on the history with a decreasing $\alpha$ till to the desired
learning rate. The same logic applies to the exploration rate $\epsilon$,
where more explorations are conducted at beginning when no policy is formed,
but the explorations decrease to a minimal level when optimal policies are
approached. This practice generally is able to accelerate the convergence of
the evolution, but does not change the system evolution in the long run.

As an comparison, Fig. 11 adopt the all-fixed the exploration rate $\epsilon$
and learning rate $\alpha$ at the very beginning. It shows that the four
fractions are consistent with the results in Fig. 2(a) using the time-varying
form, but indeed the transient takes a much longer time than the evolution in
Fig. 2(a).

##  Appendix C Appendix C: The “trembling hands” effect

The search of optimal policies benefits from the presence of the exploration
rate $\epsilon$. But, once the evolution becomes stable, the players can be
considered accumulated sufficient experiences that they have found their
optimal policies. In that case, the exploration rate $\epsilon$ may be taken
as the “trembling hands” Fundenberg and Maskin (1990); Nowak (2006), players
erroneously take actions, deviated from the guidance from their Q-tables. Fig.
9 shows that as the increase of $\epsilon$, the fraction of TR monotonically
decreases, the fraction of NB rises accordingly. This suggests that a decent
level trust and trustworthiness requires a weak “trembling hands” effect, too
many misunderstandings would undermine and break down the foundation of the
two.

Figure 9:  The “trembling hands” effect in the two-person scenario. The
fractions of TR and NB versus the exploration rate $\epsilon$. Each data is
averaged over 50 realizations. Other parameters: $\alpha=0.1$, $\gamma=0.9$.

##  Appendix D Appendix D: Failed cases of trust emergence

As the two key learning parameters $(\gamma,\alpha)$ deviate the ideal
combination, the emergence of trust and trustworthiness fails. Three cases
with typical parameter combinations are investigated, both the time series of
the four fractions and their corresponding $\Delta Q_{s_{n}}^{A,B}$ are shown
in Fig. 10.

The failures can be attributed to two aspects. When the learning rate $\alpha$
becomes large [e.g. $\alpha=0.9$ in Fig. 10(d-f) and (g-i)], this means that
the historical experiences of players is removed immediately, almost no lesson
is kept in the Q-table. In this scenario, the Q-learning algorithm loses its
strength and the evolution degenerates to the traditional iterated trust game,
where the trustor is not willing to invest in the trustee. Note that, once no
investment is made, the reciprocating behaviors from the trustee is pointless,
therefore, no decent level of trust and trustworthiness is seen, as shown in
Fig. 10(d) and (g). When the discount factor $\gamma$ becomes smaller, this
causes another problem. Without confidence of the future reward, it’s hard for
the trustor to foresee potential reciprocity from the trustee to select to
trust, and vice versa. As a result, the positive feedback between trust and
reciprocity fails to form. A decent level of trust and trustworthiness is
still hard to see [e.g. $\gamma=0.1$ in Fig. 10(a-c) and (d-f)]. With these
observations, it’s reasonable to see why the emergence of trust and
trustworthiness is observed for the combination of small $\alpha$ and large
$\gamma$ in Fig. 1.

Figure 10:  Three cases where the emergence of trust and trustworthiness is
failed to see. The parameter combinations are: (a-c) $\alpha=0.1$,
$\gamma=0.1$; (d-f) $\alpha=0.9$, $\gamma=0.1$; (g-i) $\alpha=0.9$,
$\gamma=0.9$. Other parameters are the same as Fig. 3

##  Appendix E Appendix E: Setup of Q-table in 1D latticed population

When we extend the two-individual scenario to the one-dimensional latticed
population, where each player connects 2 nearest-neighbors, i.e. the degree
$k=2$. The action set for the trustor is then extended to be
$\mathbb{A}_{r}=\\{(T,T),(T,N),(N,T),(N,N)\\}$, where the first and the second
are respectively the strategy playing against to its left and right nearest-
neighbor. Similarly, the action set for the trustee is
$\mathbb{A}_{e}=\\{(R,R),(R,B),(B,R),(B,B)\\}$. The state for the trustor is
expanded as $\mathbb{S}={s_{1},...,s_{16}}$, with $s_{1}=(TT,RR)$,
$s_{2}=(TT,RB)$,…, $s_{16}=(NN,BB)$. For example, $s_{2}=(TT,RB)$ means that
the player has invested her two nearest-neighbors in the last round, and the
left one showed reciprocity but the right one betrayed. The two Q-tables
$Q^{A,B}_{s_{n}}$ are illustrated in Table 3 and Table 4, respectively for the
trustor and trustee.

State Action | (T, T) ($a_{1}$) | (T, N) ($a_{2}$) | (N, T) ($a_{3}$) | (N, N) ($a_{4}$)
---|---|---|---|---
(TT, RR) ($s_{1}$) | $Q_{s_{1}a_{1}}$ | $Q_{s_{1}a_{2}}$ | $Q_{s_{1}a_{3}}$ | $Q_{s_{1}a_{4}}$
(TT, RB) ($s_{2}$) | $Q_{s_{2}a_{1}}$ | $Q_{s_{2}a_{2}}$ | $Q_{s_{2}a_{3}}$ | $Q_{s_{2}a_{4}}$
(TT, BR) ($s_{3}$) | $Q_{s_{3}a_{1}}$ | $Q_{s_{3}a_{2}}$ | $Q_{s_{3}a_{3}}$ | $Q_{s_{3}a_{4}}$
(TT, BB) ($s_{4}$) | $Q_{s_{4}a_{1}}$ | $Q_{s_{4}a_{2}}$ | $Q_{s_{4}a_{3}}$ | $Q_{s_{4}a_{4}}$
(TN, RR) ($s_{5}$) | $Q_{s_{5}a_{1}}$ | $Q_{s_{5}a_{2}}$ | $Q_{s_{5}a_{3}}$ | $Q_{s_{5}a_{4}}$
… | … | … | … | …
(NT, RR) ($s_{9}$) | $Q_{s_{9}a_{1}}$ | $Q_{s_{9}a_{2}}$ | $Q_{s_{9}a_{3}}$ | $Q_{s_{9}a_{4}}$
… | … | … | … | …
(NN, RR) ($s_{13}$) | $Q_{s_{13}a_{1}}$ | $Q_{s_{13}a_{2}}$ | $Q_{s_{13}a_{3}}$ | $Q_{s_{13}a_{4}}$
… | … | … | … | …
(NN, BB) ($s_{16}$) | $Q_{s_{16}a_{1}}$ | $Q_{s_{16}a_{2}}$ | $Q_{s_{16}a_{3}}$ | $Q_{s_{16}a_{4}}$

Table 3: The Q-table for the trustor in the one-dimensional latticed
population, where each player connects 2 nearest-neighbors.

State Action | (R, R) ($a_{1}$) | (R, B) ($a_{2}$) | (B, R) ($a_{3}$) | (B, B) ($a_{4}$)
---|---|---|---|---
(RR, TT) ($s_{1}$) | $Q_{s_{1}a_{1}}$ | $Q_{s_{1}a_{2}}$ | $Q_{s_{1}a_{3}}$ | $Q_{s_{1}a_{4}}$
(TT, RB) ($s_{2}$) | $Q_{s_{2}a_{1}}$ | $Q_{s_{2}a_{2}}$ | $Q_{s_{2}a_{3}}$ | $Q_{s_{2}a_{4}}$
(RR, TN) ($s_{3}$) | $Q_{s_{3}a_{1}}$ | $Q_{s_{3}a_{2}}$ | $Q_{s_{3}a_{3}}$ | $Q_{s_{3}a_{4}}$
(TT, BB) ($s_{4}$) | $Q_{s_{4}a_{1}}$ | $Q_{s_{4}a_{2}}$ | $Q_{s_{4}a_{3}}$ | $Q_{s_{4}a_{4}}$
(RR, NT) ($s_{5}$) | $Q_{s_{5}a_{1}}$ | $Q_{s_{5}a_{2}}$ | $Q_{s_{5}a_{3}}$ | $Q_{s_{5}a_{4}}$
… | … | … | … | …
(RB, TT) ($s_{9}$) | $Q_{s_{9}a_{1}}$ | $Q_{s_{9}a_{2}}$ | $Q_{s_{9}a_{3}}$ | $Q_{s_{9}a_{4}}$
… | … | … | … | …
(BR, TT) ($s_{13}$) | $Q_{s_{13}a_{1}}$ | $Q_{s_{13}a_{2}}$ | $Q_{s_{13}a_{3}}$ | $Q_{s_{13}a_{4}}$
… | … | … | … | …
(BB, NN) ($s_{16}$) | $Q_{s_{16}a_{1}}$ | $Q_{s_{16}a_{2}}$ | $Q_{s_{16}a_{3}}$ | $Q_{s_{16}a_{4}}$

Table 4: The Q-table for the trustee in the one-dimensional latticed
population, where each player connects 2 nearest-neighbors. Figure 11:  The
evolution of four strategies fractions at an all-fixed exploration rate
$\epsilon=0.01$. Each data is averaged over 500 realizations. Other
parameters: $\alpha=0.1$, $\gamma=0.9$.

## References

  * Hardin (2002) R. Hardin,  _Trust and Trustworthiness_ (Russell Sage Foundation, 2002).
  * Arrow (1974) K. J. Arrow,  _The Limits of Organization_ (Norton & Company, 1974).
  * Zak and Knack (2001) P. J. Zak and S. Knack, The Economic Journal 111, 295 (2001).
  * Algan and Cahuc (2013) Y. Algan and P. Cahuc, Annual Review of Economics 5, 521 (2013).
  * Ortiz-Ospina and Roser (2016) E. Ortiz-Ospina and M. Roser, Our World in Data (2016).
  * Solow and Simon (1957) R. M. Solow and H. A. Simon,  _Models of Man: Social and Rational_ (Willey, 1957).
  * Fukuyama (1995) F. Fukuyama,  _Trust: The Social Virtues and the Creation of Prosperity_ (New York: Free Press, 1995).
  * Samuelson and Nordhaus (2005) P. Samuelson and W. Nordhaus,  _Economics (18th edition)_ (McGraw-Hil Education, 2005).
  * Berg  _et al._ (1995) J. Berg, J. Dickhaut, and K. McCabe, Games and Economic Behavior 10, 122 (1995).
  * Johnson and Mislin (2011) N. D. Johnson and A. A. Mislin, Journal of Economic Psychology 32, 865 (2011).
  * Zucker (1986) L. G. Zucker, Research in Organizational Behavior 8, 53 (1986).
  * Buchan  _et al._ (2004) N. R. Buchan, R. T. A. Croson, and R. M. Dawes, American Journal of Sociology 108, 168 (2004).
  * Nowak (2005) M. A. Nowak, Nature 437, 1291 (2005).
  * Mcnamara  _et al._ (2008) J. Mcnamara, P. Stephens, S. Dall, and A. Houston, Proceedings of the Royal Society B: Biological Sciences 276, 605 (2008).
  * Xia  _et al._ (2022) C. Xia, Z. Hu, and D. Zhao, New Journal of Physics 24, 083041 (2022).
  * Manapat  _et al._ (2013) M. L. Manapat, M. A. Nowak, and D. G. Rand, Journal of Economic Behavior $\&$ Organization 90, S57 (2013).
  * Manapat and Rand (2012) M. Manapat and D. Rand, Dynamic Games and Applications 2, 401 (2012).
  * Rauwolf and Bryson (2018) P. Rauwolf and J. J. Bryson, Dynamic Games and Applications 8, 891 (2018).
  * Nowak (2006) M. A. Nowak, Science 314, 1560 (2006).
  * Page and Sigmund (2000) K. M. Page and N. K. Sigmund, Proceedings Biological Sciences 267, 2177 (2000).
  * Capraro  _et al._ (2020) V. Capraro, M. Perc, and D. Vilone, Physical Review E 101 (2020).
  * Kumar  _et al._ (2020) A. Kumar, V. Capraro, and M. Perc, Journal of The Royal Society Interface 17 (2020).
  * Roca  _et al._ (2009) C. P. Roca, J. A. Cuesta, and A. Sánchez, Physics of Life Reviews 6, 208 (2009).
  * Szabó and Tőke (1998) G. Szabó and C. Tőke, Phys. Rev. E 58, 69 (1998).
  * 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).
  * LeCun  _et al._ (2015) Y. LeCun, Y. Bengio, and G. Hinton, Nature 521, 436 (2015).
  * Zhang  _et al._ (2020a) J.-Q. Zhang, S.-P. Zhang, L. Chen, and X.-D. Liu, Physical Review E 101, 042402 (2020a).
  * Song  _et al._ (2022) Z. Song, H. Guo, D. Jia, M. Perc, X. Li, and Z. Wang, Neurocomputing 513, 104 (2022).
  * Ding  _et al._ (2023) Z.-W. Ding, G.-Z. Zheng, C.-R. Cai, W.-R. Cai, L. Chen, J.-Q. Zhang, and X.-M. Wang, arXiv:2307.04612 (2023).
  * Andrecut and Ali (2001) M. Andrecut and M. Ali, Physical Review E 64, 067103 (2001).
  * Zhang  _et al._ (2019) S.-P. Zhang, J.-Q. Dong, L. Liu, Z.-G. Huang, L. Huang, and Y.-C. Lai, Physical Review E 99, 032302 (2019).
  * Zhang  _et al._ (2020b) S.-P. Zhang, J.-Q. Zhang, L. Chen, and X.-D. Liu, Nonlinear Dynamics 99, 3301 (2020b).
  * Tomov  _et al._ (2021) M. S. Tomov, E. Schulz, and S. J. Gershman, Nature Human Behaviour 5, 764 (2021).
  * He  _et al._ (2022) Z. He, Y. Geng, C. Du, L. Shi, and Z. Wang, New Journal of Physics 24, 123038 (2022).
  * Shi and Rong (2022) Y. Shi and Z. Rong, IEEE Transactions on Circuits and Systems II: Express Briefs 69, 2463 (2022).
  * 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).
  * Olsson  _et al._ (2020) A. Olsson, E. Knapska, and B. Lindström, Nature Reviews Neuroscience 21, 197 (2020).
  * 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).
  * Auer  _et al._ (2002) P. Auer, N. Cesa-Bianchi, and P. Fischer, Machine Learning 47, 235 (2002).
  * Fundenberg and Maskin (1990) D. Fundenberg and E. Maskin, The American Economic Review 80, 274 (1990).
  * Gereke  _et al._ (2018) J. Gereke, M. Schaub, and D. Baldassarri, PloS one 13, e0199834 (2018).
  * Woolley and Fishbach (2017) K. Woolley and A. Fishbach, Journal of Consumer Psychology 27, 1 (2017).
  * Evans and Krueger (2011) A. M. Evans and J. I. Krueger, Journal of Experimental Social Psychology 47, 171 (2011).
  * Lenton and Mosley (2011) P. Lenton and P. Mosley, Journal of Economic Psychology 32, 890 (2011).
  * Engle-Warnick and Slonim (2004) J. Engle-Warnick and R. L. Slonim, Journal of Economic Behavior $\&$ Organization 55, 553 (2004).
  * Engle-Warnick and Slonim (2006) J. Engle-Warnick and R. L. Slonim, Games and Economic Behavior 54, 95 (2006).
  * Engle-Warnick and Slonim (2001) J. Engle-Warnick and R. L. Slonim, Economics Papers (2001).
  * Camerer and Thaler (1995) C. Camerer and R. H. Thaler, Journal of Economic perspectives 9, 209 (1995).
  * Camerer  _et al._ (2004) C. F. Camerer, G. Loewenstein, and M. Rabin,  _Advances in Behavioral Economics_ (Princeton university press, 2004).
  * Henrich  _et al._ (2005) J. Henrich, R. Boyd, S. Bowles, C. Camerer, E. Fehr, H. Gintis, R. McElreath, M. Alvard, A. Barr, J. Ensminger, _et al._ , Behavioral and Brain Sciences 28, 795 (2005).
  * Capraro and Perc (2018) V. Capraro and M. Perc, Frontiers in Physics 6, 107 (2018).
  * Capraro and Perc (2021) V. Capraro and M. Perc, Journal of the Royal Society interface 18, 20200880 (2021).