# Electromagnetic properties of neutrinos from scattering on bound electrons
in atom

Junu Jeong(1), Jihn E. Kim(2), and Sungwoo Youn(1) (1)Center for Axion and
Precision Physics Research, IBS, Daejeon 34051, Republic of Korea
(2)Department of Physics, Kyung Hee University, 26 Gyungheedaero, Dongdaemun-
Gu, Seoul 02447, Republic of Korea

###### Abstract

We consider the effects of bound atomic electrons scattered by solar neutrinos
due to the electromagnetic properties of neutrinos. This necessiate
considering the recoil of atomic nucleus, which should be considered in the
momentum conservation, but that effect to the energy conservation is
negligible. This effect changes the kinematic behavior of the scattered
electron compared to that scattered on free electrons. We apply this effect to
the recent XENON1T data, but the bounds obtained from this is not very
restrictive. We obtained the bounds: the (transition) magnetic moment
$|f_{\alpha\beta}|\leq 0.86\times 10^{-7}$ (times the electron Bohr magneton)
and the charge radius $|\tilde{r}|<4.30\times 10^{-17\,}{\rm cm}$. For a non-
vanishing millicharge ($\varepsilon$), the allowed bound is shown in the
$\tilde{r}^{2}-\varepsilon$ plane.

Neutrino Charge Radius, Neutrino Magnetic Moment, Millicharge, XENON1T

###### pacs:

13.15.+g, 12.15.Mm, 13.40.Gp, 13.66.-a

##  I Introduction

Ever since the discovery of neutrinos, it has been a great challenge to find
the electromagnetic(EM) properties of neutrinos Pauli30 ; Reines57 ; LeeTD63 ,
which continued in the gauge theory era Clark73 ; Kim74 ; Kim78 ; Okun86 ;
Barbieri87 . To have the EM properties of neutrinos at the observable level,
particles in the beyond-the-standard-model(BSM) had to be introduced Kim76 ;
Yanagida88 ; Giunti15 ; KimMoscow19 . Two relevant EM properties of neutrinos
are the magnetic moment(MM) and the charge radius(CR). The MM, coupling to the
EM field strength, is gauge invariant. On the other hand, the CR, coupling to
the EM field itself, is not gauge invariant. Gauge invariance is needed to
maintain the renormalizability. Therefore, if the probing wave length in the
effective theory is much larger than the cutoff length-scale in the
renormalization, as in the Fermi weak interaction much above the electroweak
length-scale $10^{-16\,}$cm, CR can be also considered as a useful physical
parameter. The cutoff energy scale for the CR is the mass of a heavy BSM
particle whose exchange produces the relevant charge radius. In this paper, we
take this viewpoint considering the electromagnetic properties of neutrinos.

The published experimental limit on the MM of $\nu_{e}$ in units of the
electron Bohr magneton is BdNuMagMom17 , $|f|\lesssim 2.8\times 10^{-11}$, and
the limit on the squared CR, $\tilde{r}^{2}$ is BdNuChargeR10 ,
$\tilde{r}^{2}=[-2.1,+3.3]\times 10^{-32}{\rm cm^{2}}$. The standard model(SM)
prediction Fujikawa80  on the neutrino MM is much smaller than the upper limit
presented above, by a factor of O($10^{-8}$). Note that the millicharge limit
of neutron is O($10^{-21\,}$cm) ChargeNeutron88 , which however cannot be
directly used for a limit on the millicharges of neutrinos.

Recently, the XENON group considered the possibility of MMs of solar neutrinos
for the excess events in their data, for an exposure of 0.65 tonne-year with
an extremely low background rate of $76\pm 2(\rm
stat.)\,$(Events)/${\rm(t\cdot y\cdot keV)}$ Xenon1T20 . The plot on these
excess events starts around the electron recoil energy near 2–3 keV, and ends
around 30 keV as shown in Fig. 1 (a). Out of the fitted 42,179.4 events, the
estimated SM background by the solar neutrinos is 220.8 events Xenon1T20 ,
only 0.52%. For estimating non-vanishing electromagnetic properties of solar
neutrinos, therefore, one may neglect the SM background and associate the bulk
of 42,179.4 events with the electromagnetic properties of neutrinos in the
scattering process. In this way, Ref. Xenon1T20  obtained the bound on the
neutrino MM, $[1.4,2.9]\times 10^{-11}$ times the electron Bohr magneton. But,
the assumed cross section by the XENON group has the reciprocal dependence on
the recoil electron energy. This reciprocal dependence is obtained by assuming
the two-body scattering Vogel89 , $\nu_{e}+e$, which may not be a correct
method. The reason is that in kicking out an electron from a heavy atom, the
atomic nucleus can carry away some momentum and the momentum conservation used
in the $\nu_{e}+e_{\rm in\,atom}\to\nu_{e}^{\prime}+e^{\prime}$ is not exact.

In Sec. II, we express the XENON1T unit “${\rm Events/t\cdot y\cdot keV}$” in
terms of the cross section unit “$\,\mathrm{MeV}^{3}$” to compare with the
theoretical prediction. In Sec. III, we discuss the possibility of including
electromagnetic properties of neutrinos in the scattering process. Here, the
three-body phase space with a non-relativistic atom is discussed. In Subsec.
III.4, we briefly comment on the possibility of neutrino milli-charges. In
Sec. IV, from XENON1T data we obtain the CR and millicharge bounds of
electron-type neutrino. Section V is a brief summary on the EM properties of
neutrinos. In Appendix A, we estimate the parameters used in the calculation.

##  II Event Rate and Cross Section in the XENON1T Data

(a) (b)

Figure 1: (a) Cross section vs. the data curve for $\mu_{\nu}=7\times
10^{-11}\mu_{\rm B}$ reproducing the RHS panel of Fig. 1 of Xenon1T20 , and
(b) a conversion of the event rate to the cross section.

The XENON group presents data in units of ${\rm Events/t\cdot y\cdot keV}$. To
convert it to cross section, we use the right-hand side (RHS) panel of Fig. 1
of Xenon1T20 , and followed their analysis using the neutrino magnetic moment
$\mu_{\nu}=7\times 10^{-11}\mu_{\rm B}$ and the solar $pp$ neutrino flux of
Bahcall04 . They used the inverse-$T_{e}$ rule presented in Vogel89 . In our
paper here, the $pp$ neutrino energy flux of Bahcall04  in the region,
$E_{\nu}\geq 4\,\,\mathrm{MeV}$, is fitted by $F^{\rm solar}_{\rm pp}(E)$. The
analytic form for the $pp$-flux function is presented in Vitagliano17 ,

$F^{\rm solar}_{pp}(E)=P_{0}E^{2}(Q+m_{e}-E)\sqrt{(Q+m_{e}-E)^{2}-m_{e}^{2}},$
(1)

where $Q\approx 423.41\,\textrm{keV}$ is the end point energy of neutrino, and
$m_{e}\approx 0.511\,\textrm{MeV}$ is the electron mass.
$P_{0}=188.143\,\textrm{MeV}^{-5}$ is the normalization constant. From this
form, some averages of energy functions are

$\begin{split}\langle E^{-2}\rangle=&~{}30.0256\,\textrm{MeV}^{-2},\\\ \langle
E^{-1}\rangle=&~{}4.46821\,\textrm{MeV}^{-1},\\\ \langle
E^{1}\rangle=&~{}0.267946\,\textrm{MeV}^{1},\\\ \langle
E^{2}\rangle=&~{}0.0794926\,\textrm{MeV}^{2},\\\ \langle
E^{3}\rangle=&~{}0.0251769\,\textrm{MeV}^{3},\\\ \langle\ln{E/Q_{\rm
min}}&\rangle=-0.528546.\\\ \end{split}$ (2)

Figure 1 (a) is the overlap of the cross section and the RHS panel data of
Fig. 1 of Ref. Xenon1T20 . Here, we used $Q_{\rm min}$ as the minimum of
$|{\bf p}_{e}^{\prime}|^{2}/(2m_{e})$. Figure 1 (b) is the cross section
$d\sigma/dT_{e}$ in the MeV-3 units, obtained from the XENON1T data using the
inverse-$T_{e}$ rule. Then, our calculation of $d\sigma/dT_{e}$ resulting from
some electromagnetic properties of neutrinos, not using the inverse-$T_{e}$
rule but using our formula, will be compared with Fig. 1 (b).

##  III Cross Sections with Electromagnetic Form Factors of Neutrino

The photon vertex does not change the chirality. For the mass and MM,
therefore, a chirality change must take place. So, the MM can arise by
attaching a photon vertex in the mass generating diagrams. For the CR, it is
even simpler because the CR has the same chiral property as that given by the
EM coupling. For the neutrino, the gauge coupling is expressed as the form
factor $F^{\nu}_{1}(q^{2})$ with $F^{\nu}_{1}(0)=0$. The CR is defined by the
first term in the expansion of $F^{\nu}_{1}(q^{2})$ in terms of momentum
transfer $q^{2}$. As noted in the beginning, we use the charge radius for the
scale $q^{2}\ll$(cutoff scale)2. In addition, if a (almost) massless extra
photon (ex-photon) in the BSM sector is present, there is a possibility that
neutrinos can carry millicharges Holdom86 ; ParkJC07 . The chiral property due
to millicharges is the same as that due to the charge radius. The vector
couplings of the charge radii and millicharges of neutrinos can mix with the
$G_{F}$-order SM couplings, but the interference between these vector and
axial-vector interactions and the SM amplitude has not played an important
role because of the tiny contribution in the data of XENON1T if it were from
solar neutrinos.

###  III.1 SM cross section and beyond

The effective interaction in the SM for the electron-type neutrino scattering
on an electron is given by the interaction tHooft71 ; KimRMP81 ,

$\displaystyle\frac{G_{F}}{\sqrt{2}}\bar{\nu}_{e}\gamma^{\mu}(1+\gamma_{5})\nu_{e}\,\bar{e}\gamma_{\mu}(g_{V}+g_{A}\gamma_{5})e.$
(3)

When we consider the neutrino scattering, kicking out an electron from an atom
with atomic number $Z$, the following six-fermion interaction with parameter
$M_{\rm eff}$ can be considered,

$\begin{split}\frac{1}{M_{\rm
eff}^{3}}\cdot\frac{G_{F}}{\sqrt{2}}\bar{\nu}_{e}\gamma^{\mu}(1+\gamma_{5})\nu_{e}\,\bar{e}\gamma_{\mu}(g_{V}+g_{A}\gamma_{5})e\,\overline{A}A,\end{split}$
(4)

where $A$ is treated as a quantum field with $Z-1$ atomic electrons.111$A$ may
be considered as a boson, which gives the same result. A strategy to estimate
$M_{\rm eff}$ for four-fermion interaction is given in Appendix A. We
interpret $M_{\rm eff}$ as representing the electromagnetic process occurring
in the atomic orbits. For four-fermion interaction, it is a kind of definition
for the atomic process. For six-fermion interaction of (4), we use the same
$M_{\rm eff}$ determined in Appendix A.

In the SM, we have

$\begin{split}g_{V}=\frac{1}{2}+2\sin^{2}\theta_{W},~{}g_{A}=\frac{1}{2},\end{split}$
(5)

where $\sin^{2}\theta_{W}\simeq 0.231$ sin2LHC18 ; KimRMP81  is the weak
mixing angle, and $e^{2}=g^{2}\sin^{2}\theta_{W}$. The EM couplings of
neutrinos are

$\begin{split}\overline{\nu}_{\beta}\gamma^{\mu}eF_{1\,\alpha\beta}^{\nu}\nu_{\alpha}A^{\rm
em}_{\mu}+\overline{\nu}_{\beta}\sigma^{\mu\nu}\frac{ef_{\alpha\beta}}{2M}\nu_{\alpha}(q_{\mu}A^{\rm
em}_{\nu}-q_{\nu}A^{\rm em}_{\mu})\end{split}$ (6)

where $F_{1\,\alpha\beta}$ and $\frac{f_{\alpha\beta}}{2M}$ are matrices
between weak eigenstates of neutrinos, $\nu_{\beta}(k_{\nu}^{\prime})$ and
$\nu_{\alpha}(k_{\nu})$. $f_{\alpha\beta}$ is the MM of neutrino in units of
the Bohr magneton of the mass $M$ particle. For $\alpha=\beta$, Eqs. (3) and
(6) combine, and in particular for the electron-type neutrino $\alpha=\beta$
we have

$\begin{split}\frac{G_{F}}{\sqrt{2}}\bar{\nu}_{e}\gamma^{\mu}(1+\gamma_{5})\nu_{e}\,\bar{e}\left(\left[g_{V}+\frac{F_{1}^{\nu_{e}}}{\sqrt{2}G_{F}q^{2}}\right]\gamma_{\mu}+g_{A}\gamma_{\mu}\gamma_{5}\right)e.\end{split}$
(7)

The SM neutrino $\nu_{\alpha}$ is a two-component spinor. Therefore, in the
cross section calculation the spin sums for a Majorana neutrino is
$\frac{1}{2}$ of that of a Dirac neutrino, which will be taken into account.

For calculating scattering cross section on electrons bound in the orbits of a
rest nucleus, the wave function of electrons in the electron cloud around the
nucleus must be considered. Let the uncertainty of the outgoing electron
position be $\Delta r$. Then, the uncertainty of the momentum is $\Delta
p\approx 1/\Delta r$. We consider the recoil energy greater than 2.5 keV
(glimpsing Fig. 1 (a) ), corresponding to to $\Delta r<3.9\times 10^{-9}$cm
with $\Delta p=\sqrt{2m_{e}T_{e}}$ for the electron recoil energy $T_{e}$.
But, in the neutrino scattering with the incoming solar neutrino energy less
than $0.45\,\,\mathrm{MeV}$ Bahcall04 , the uncertainty of neutrino position
is less than $4.38\times 10^{-11}$cm which is much shorter than the orbit
radius of the $n$-th level of Xenon, $10^{-10}n$ cm. So, in the neutrino
scattering we can safely use the two-body scattering formula as described in
Appendix A.

The $3\to n$ transition rate is

$\begin{split}Rate=(2\pi)^{4-3n}\frac{d^{3}{\bf p}_{1}^{\prime}\cdots
d^{3}{\bf
p}_{n}^{\prime}}{2E_{1}2E_{2}2E_{3}2E_{1}^{\prime}2E_{2}^{\prime}\cdots
2E_{n}^{\prime}}\frac{\delta^{4}(p_{1}+p_{2}+p_{3}-p_{1}^{\prime}\cdots-
p_{n}^{\prime})}{V^{2}}|\langle p_{1}^{\prime}\cdots
p_{n}^{\prime}|T|p_{1}p_{2}p_{3}\rangle|^{2}\end{split}$ (8)

from which we calculate the three body scattering cross section. A flux of
incoming particles 1 (neutrino in our case) scatter on two particles 2
(electron in our case) and 3 (atom in our case). Thus, the cross section
consists of two parts, one with a flux factor $1/|{\bf v}_{1}-{\bf v}_{2}|$
and the other $1/|{\bf v}_{1}-{\bf v}_{3}|$. Taking ${\bf v}_{3}=0$, the
neutrino-electron scattering is

$\begin{split}d\sigma=\frac{(2\pi)^{4-9}}{|{\bf v}_{1}-{\bf
v}_{2}|}\frac{d^{3}{\bf p}_{1}^{\prime}d^{3}{\bf p}_{2}^{\prime}d^{3}{\bf
p}_{3}^{\prime}}{2E_{1}2E_{2}2E_{3}2E_{1}^{\prime}2E_{2}^{\prime}2E_{3}^{\prime}}\frac{\delta^{4}(\sum
p_{i}-\sum p_{f}^{\prime})}{V}|\langle
p_{1}^{\prime}p_{2}^{\prime}p_{3}^{\prime}|T|p_{1}p_{2}p_{3}\rangle|^{2}\end{split}$
(9)

where we can take $V$ as the atomic volume to which the flux of neutrinos
sweep. For $1/V$, the sum $B$ in Eq. (42) is used. The $T$-matrix squared in
Eq. (9) is summed over the spins

$\begin{split}\frac{1}{(2s_{e}+1)(2s_{\nu}+1)}\sum_{s_{\nu},s_{e}}\sum_{s^{\prime}_{\nu},s^{\prime}_{e}}|{\rm
Amp}|^{2}=(|T|^{2}_{\rm Q}).\end{split}$ (10)

If we consider the MM interaction also, then we will have the factor
$(|T|^{2}_{\rm MM}+|T|^{2}_{\rm Q})$ instead of $(|T|^{2}_{\rm Q})$. For the
volume $V$, we use the Xenon volume for each principal quantum number $n$
separately. The $n$-th shell volume is $V=(4\pi/3)(na_{Z})^{3}$ where $a_{Z}$
is the K-shell radius of Xenon atom. Then, for the
$\nu_{\alpha}(k_{\nu})+e(p_{e})+A(P_{A})\to\nu^{\prime}_{\beta}(k_{\nu}^{\prime})+e^{\prime}(p_{e}^{\prime})+A^{\prime}(P_{A}^{\prime})$
scattering, Eq. (9) becomes

$\begin{split}d\sigma=\frac{3Z^{3}\alpha^{3}_{\rm
em}}{2^{12}\pi^{4}}\,\frac{m_{e}^{2}|{\bf p}_{e}^{\prime}|^{2}d|{\bf
p}_{e}^{\prime}|E_{\nu}^{\prime}dE_{\nu}^{\prime}}{(m_{e}+T_{e})E_{\nu}M_{A}^{2}}\delta(E_{\nu}+m_{e}-\delta_{B}+M_{A}-E_{\nu}^{\prime}-(m_{e}+T_{e})-M_{A})|\langle
p_{1}^{\prime}p_{2}^{\prime}p_{3}^{\prime}|T|p_{1}p_{2}p_{3}\rangle|^{2}\,x\end{split}$
(11)

where $\delta_{B}$ is the binding energy and $T_{e}$ is the kinetic energy of
the final electron $e^{\prime}$, and $x$ is given in Appendix A for the
process in consideration. The following scalar products will be useful,

$\begin{split}&q^{2}=-2E_{\nu}E_{\nu}^{\prime}(1-\cos\theta),~{}q=k_{\nu}-k_{\nu}^{\prime},\\\
&k_{\nu}\cdot k_{\nu}^{\prime}=-k_{\nu}\cdot q=k_{\nu}^{\prime}\cdot
q=-\frac{q^{2}}{2},\\\ &p_{e}\cdot
p_{e}^{\prime}=m_{e}^{2}-\frac{q^{2}}{2}~{}\textrm{assuming atom at rest},\\\
&k_{\nu}^{\prime}\cdot p_{e}=m_{e}E_{\nu}^{\prime}=k_{\nu}\cdot
p_{e}^{\prime}=E_{\nu}E_{e}^{\prime}-E_{\nu}|{\bf
p}_{e}^{\prime}|\cos\theta_{e},\end{split}$ (12)

where we used massless neutrinos, $\theta_{e}$ is the polar angle of the
outgoing electron direction relative to the incoming neutrino direction,222The
polar angle made by the outgoing neutrino direction relative to the incoming
neutrino direction will be denoted as $\theta$. and the non-relativistic
approximation $E_{e}^{\prime}\simeq m_{e}+\frac{|{\bf
p}_{e}^{\prime}|^{2}}{2m_{e}}=m_{e}+T_{e}$ is used.

Firstly, let us show the inverse-$T_{e}$ rule, used in Sec. II. Note the
kinematic variables for
$\nu_{e}(k_{\nu})+e(p_{e})\to\nu_{e}^{\prime}(k^{\prime}_{\nu})+e^{\prime}(p^{\prime}_{e})$,
for (almost) massless neutrinos and in the limit $M_{A}=\infty$,

$\begin{split}&E_{\nu}^{\prime}=|{\bf k}_{\nu}^{\prime}|=\sqrt{E_{e}^{\prime
2}+E_{\nu}^{2}-m_{e}^{2}},\\\ &2k_{\nu}\cdot
k_{\nu}^{\prime}=2E_{\nu}E_{\nu}^{\prime}(1-\cos\theta),\end{split}$ (13)

where $\theta$ is the angle of the outgoing netrino relative to the incoming
neutrino. It is possible to express $\theta$ in terms of $E_{e}^{\prime}$
Vogel89 ,
$\cos\theta=(E_{\nu}+m_{e})E_{\nu}^{-1}T_{e}^{1/2}(T_{e}+2m_{e})^{-1/2}$ such
that $|{\bf p}^{\prime}_{\nu}|\simeq E_{\nu}^{\prime}$ becomes

$\begin{split}E_{\nu}^{\prime}\simeq\sqrt{E_{\nu}^{2}+2m_{e}T_{e}-2E_{\nu}\sqrt{2m_{e}T_{e}}(E_{\nu}+m_{e})E_{\nu}^{-1}T_{e}^{1/2}(T_{e}+2m_{e})^{-1/2}}.\end{split}$
(14)

Now, the $|{\bf k}_{\nu}^{\prime}|d|{\bf
k}_{\nu}^{\prime}|=E_{\nu}^{\prime}dE_{\nu}^{\prime}=-E_{\nu}^{\prime}dT_{e}$
integration with the energy conservation delta function $\delta(f(T_{e}))$,
with the approximation $T_{e}\ll m_{e}$ and $E_{\nu}^{\prime}\simeq E_{\nu}$
in the final formula, we obtain

$\begin{split}\frac{d\sigma}{dT_{e}}\simeq\frac{1}{2^{5}\pi
m_{e}^{2}T_{e}}.\end{split}$ (15)

The reason for the $1/T_{e}$ dependence is the specific $T_{e}$ powers in Eq.
(14), and if it is slightly violated then the $1/T_{e}$ dependence would not
result. Thus, we anticipate the inverse-$T_{e}$ rule is not applicable if the
momentum carried by the scattered Xenon atom is considered.

Let us now proceed the three-body scattering given in Eq. (47). The four
momenta delta function

$\begin{split}\delta^{(3)}({\bf k}_{\nu}-{\bf k}^{\prime}_{\nu}-{\bf
p}^{\prime}_{e}-{\bf
P}^{\prime}_{A})\delta(E_{\nu}-E_{\nu}^{\prime}-\Delta-\frac{|{\bf
p}_{e}^{\prime}|^{2}}{2m_{e}}-\frac{|{\bf
P}_{A}^{\prime}|^{2}}{2M_{A}})\end{split}$ (16)

where ${\bf p}_{e}^{\prime}$ and ${\bf P}_{A}^{\prime}$ are the momenta of the
outing $e^{\prime}$ and outgoing atom $A^{\prime}$ with $m_{e}$ and $M_{A}$,
respectively, and $E_{\nu}-E_{\nu}^{\prime}=\Delta=T_{e}+\delta_{B}$ with the
binding energy $\delta_{B}$ of the electron in the atom. The energy transfer
to the atom is neglected in the energy-conservation $\delta$ function, since
it is so small. The averaged spin sums of $|T|^{2}$ from Eq. (4) is

$\begin{split}&\frac{G_{F}^{2}}{2M_{\rm eff}^{6}}{\rm
Tr\,}\gamma^{0}(1+\gamma_{5})\gamma^{\rho}\gamma^{0}k\hskip-4.26773pt/^{\prime}\gamma^{\alpha}(1+\gamma_{5})k\hskip-4.26773pt/\cdot{\rm
Tr\,}\gamma^{0}(g_{V}^{*}+g_{A}^{*}\gamma_{5})\gamma_{\rho}\gamma_{0}(p\hskip-4.26773pt/^{\prime}+m_{e})\gamma^{0}\gamma_{\alpha}(g_{V}+g_{A}\gamma_{5})(p\hskip-4.26773pt/+m_{e})\cdot{\rm
Tr\,}(P\hskip-6.25963pt/^{\prime}+M_{A})(P\hskip-6.25963pt/+M_{A})\\\
&~{}=\frac{G_{F}^{2}}{2M_{\rm eff}^{6}}{\rm
Tr\,}2(1-\gamma_{5})\gamma^{0}\gamma^{\rho}\gamma^{0}k\hskip-4.26773pt/^{\prime}\gamma^{\alpha}k\hskip-4.26773pt/\cdot{\rm
Tr\,}\gamma^{0}\gamma_{\rho}\gamma^{0}(p\hskip-4.26773pt/^{\prime}+m_{e})\gamma_{\alpha}\\\
&\qquad\qquad\qquad\big{[}(|g_{V}|^{2}+|g_{A}|^{2})p\hskip-4.26773pt/-(g_{V}g_{A}^{*}+g_{V}^{*}g_{A})p\hskip-4.26773pt/\gamma_{5}+(|g_{V}|^{2}-|g_{A}|^{2})m_{e}-(g_{V}g_{A}^{*}-g_{V}^{*}g_{A})m_{e}\gamma_{5}\big{]}\cdot(8M_{A}^{2})\\\
&~{}=\frac{G_{F}^{2}}{2M_{\rm eff}^{6}}2\cdot
4(k^{\prime\rho}k^{\alpha}-k\cdot
k^{\prime}g^{\rho\alpha}+k^{\prime\alpha}k^{\rho}-i\varepsilon^{\rho\mu\alpha\nu}k^{\prime}_{\mu}k_{\nu})\times
4\big{(}(|g_{V}|^{2}+|g_{A}|^{2})(p^{\prime}_{\rho}p_{\alpha}-p\cdot
p^{\prime}g_{\rho\alpha}+p^{\prime}_{\alpha}p_{\rho})\\\
&\qquad~{}+(|g_{V}|^{2}-|g_{A}|^{2})m_{e}^{2}g_{\rho\alpha}-i(g_{V}g_{A}^{*}+g_{V}^{*}g_{A})\varepsilon_{\rho\kappa\alpha\eta}p^{\prime\kappa}p^{\eta}\big{)}\cdot(8M_{A}^{2})\\\
&~{}=\frac{G_{F}^{2}}{2M_{\rm eff}^{6}}2\cdot
4\Big{(}4(|g_{V}|^{2}+|g_{A}|^{2})(2k\cdot pk^{\prime}\cdot p^{\prime}+2k\cdot
p^{\prime}k^{\prime}\cdot
p)+4\big{[}(|g_{V}|^{2}-|g_{A}|^{2})m_{e}^{2}(-2k\cdot k^{\prime})\\\
&\qquad-2{\rm Re\,}(g_{V}g_{A}^{*})(k\cdot pk^{\prime}\cdot p^{\prime}-k\cdot
p^{\prime}k^{\prime}\cdot p)\big{]}\Big{)}\cdot(8M_{A}^{2})\\\
&~{}=\frac{2^{8}G_{F}^{2}M_{A}^{2}m_{e}^{2}E_{\nu}^{2}}{M_{\rm
eff}^{6}}\Big{(}(|g_{V}|^{2}+|g_{A}|^{2})\big{[}(1-\frac{|{\bf
p}_{e}^{\prime}|^{2}}{2m_{e}E_{\nu}})+(\frac{E_{e}^{\prime}}{m_{e}}-\frac{|{\bf
p}_{e}^{\prime}|}{m_{e}}\cos\theta_{e})^{2}\big{]}+(|g_{V}|^{2}-|g_{A}|^{2})m_{e}^{2}(-q^{2}/2E_{\nu}^{2})\\\
&\qquad-{\rm Re\,}(g_{V}g_{A}^{*})\big{[}(1-\frac{|{\bf
p}_{e}^{\prime}|^{2}}{2m_{e}E_{\nu}})+(\frac{E_{e}^{\prime}}{m_{e}}-\frac{|{\bf
p}_{e}^{\prime}|}{m_{e}}\cos\theta_{e})^{2}\big{]}\Big{)}\end{split}$ (17)

where we will use $M_{\rm eff}^{6}$ given in Eq. (44) of Appendix A and
$\sum_{s}u_{e_{\rm atom}}(p)\bar{u}_{e_{\rm atom}}(p)\simeq
p\hskip-4.26773pt/+m_{e}$. Here, however, the three momentum conservation is
treated accurately. On the other hand, most people use the two body scattering
and the uncertaintly here is the probability in finding a prospective bound
electron in the atom and the subsequent integration over the probability
function of the electron cloud in the atom, which will be very complicated.
This calculation is not found in the literature.

###  III.2 Magnetic moment contribution $|T|^{2}_{\rm MM}$

If neutrinos have magnetic moment, we consider the second term in Eq. (6). In
this case, summing over the spins, viz. Eq. (7), we have

$\begin{split}|T|^{2}_{\rm MM}&=\frac{1}{4}\left(\frac{e^{2}f_{M}}{2MM_{\rm
eff}^{3}}\right)^{2}\frac{q_{\mu}q_{\alpha}}{q^{4}}\,{\rm
Tr}\,\sigma^{\mu\nu}(\frac{k\hskip-4.26773pt/_{i}}{2})\sigma^{\alpha\beta}(\frac{k\hskip-4.26773pt/_{f}}{2})\cdot{\rm
Tr}\,\gamma_{\nu}(p\hskip-4.26773pt/_{i}+m_{e})\gamma_{\beta}(p\hskip-4.26773pt/_{f}+m_{e})\cdot{\rm
Tr\,}(P\hskip-6.25963pt/^{\prime}+M_{A})(P\hskip-6.25963pt/+M_{A})\\\ \\\
&=\frac{1}{64}\left(\frac{e^{2}f_{M}}{2MM_{\rm
eff}^{3}}\right)^{2}\frac{1}{q^{4}}{\rm
Tr}\,(q\hskip-4.26773pt/\gamma^{\nu}-\gamma^{\nu}q\hskip-4.26773pt/)k\hskip-4.26773pt/_{i}(q\hskip-4.26773pt/\gamma^{x}-\gamma^{x}q\hskip-4.26773pt/)k\hskip-4.26773pt/_{f}\cdot{\rm
Tr}\gamma_{\nu}(p\hskip-4.26773pt/_{i}+m_{e})\gamma_{\beta}(p\hskip-4.26773pt/_{f}+m_{e})\cdot(8M_{A}^{2})\\\
&=\frac{1}{64}\left(\frac{e^{2}f_{M}}{2MM_{\rm
eff}^{3}}\right)^{2}\frac{1}{q^{4}}\,Q^{\nu x}L_{\nu
x}\cdot(8M_{A}^{2})\end{split}$ (18)

where we neglected the neutrino mass. We obtain333The trace of six gamma
matrices, $(1/4){\rm
Tr}\,\gamma^{\rho}a\hskip-4.26773pt/b\hskip-4.26773pt/\gamma^{\sigma}c\hskip-4.26773pt/d\hskip-4.26773pt/=g^{\rho\sigma}(a\cdot
b\,c\cdot d-a\cdot c\,b\cdot d+a\cdot d\,b\cdot c)+a\cdot
c(b^{\rho}d^{\sigma}+b^{\sigma}d^{\rho})+b\cdot
d(a^{\rho}c^{\sigma}+a^{\sigma}c^{\rho})\\\ -a\cdot
d(b^{\rho}c^{\sigma}+b^{\sigma}c^{\rho})-b\cdot
c(a^{\rho}d^{\sigma}+a^{\sigma}d^{\rho})+a\cdot
b(-c^{\rho}d^{\sigma}+c^{\sigma}d^{\rho})+c\cdot
d(a^{\rho}b^{\sigma}-a^{\sigma}b^{\rho})$, is used.

$\begin{split}Q^{\nu x}&\equiv{\rm
Tr}(q\hskip-4.26773pt/\gamma^{\nu}-\gamma^{\nu}q\hskip-4.26773pt/)(k\hskip-4.26773pt/_{i}+m_{\nu})(q\hskip-4.26773pt/\gamma^{x}-\gamma^{x}q\hskip-4.26773pt/)(k\hskip-4.26773pt/_{f}+m_{\delta})\\\
&=4g^{\nu x}(2q\cdot k_{i}\,q\cdot k_{f})+4k_{i}\cdot
k_{f}(q^{\nu}q^{x}-q^{2}\,g^{\nu
x})+4q^{2}(k_{i}^{\nu}k_{f}^{x}+k_{i}^{x}k_{f}^{\nu})\\\
&\quad-4q^{\nu}(k_{f}^{x}\,q\cdot k_{i}+k_{i}^{x}\,q\cdot
k_{f})-4q^{x}k_{f}^{\nu}\,q\cdot k_{i}+4q^{x}k_{i}^{\nu}\,q\cdot k_{f},\\\
L_{\nu x}&\equiv 4(p^{i}_{\nu}p^{f}_{x}+p^{f}_{\nu}p^{i}_{x}-p^{i}\cdot
p^{f}g_{\nu x}+m_{e}^{2}g_{\nu x}),\end{split}$ (19)

from which we obtain, neglecting O($m_{e}^{2}/E_{\nu}^{2}$), using Eq. (12),

$\begin{split}\frac{1}{q^{4}}Q^{\nu x}L_{\nu
x}=&\frac{1}{q^{4}}\Big{(}32k_{i}\cdot k_{f}\,q\cdot p_{i}\,q\cdot
p_{f}+16q^{2}(-k_{i}\cdot k_{f}\,p_{i}\cdot p_{f}+2k_{i}\cdot
p_{i}\,k_{f}\cdot p_{f}+2k_{i}\cdot p_{f}\,k_{f}\cdot p_{i})\\\ &-32q\cdot
k_{i}(q\cdot p_{i}\,k_{f}\cdot p_{f}+q\cdot p_{f}\,k_{f}\cdot p_{i})-32q\cdot
k_{f}(q\cdot p_{i}\,k_{i}\cdot p_{f}+q\cdot p_{f}\,k_{i}\cdot p_{i})\\\
&+16m_{e}^{2}(4q\cdot k_{i}\,q\cdot k_{f}-q^{2}k_{i}\cdot k_{f})\Big{)}\\\
=&\frac{1}{q^{4}}\Big{(}-4q^{6}+16q^{2}(k_{i}\cdot p_{i}+k_{f}\cdot
p_{i})(k_{i}\cdot p_{f}+k_{f}\cdot p_{f})\Big{)}\\\
&=4Q^{2}-\frac{16}{Q^{2}}(k_{i}\cdot p_{i}+k_{f}\cdot p_{f})(k_{i}\cdot
p_{f}+k_{f}\cdot
p_{i})=8E_{\nu}E_{\nu}^{\prime}(1-\cos\theta)-\frac{8m_{e}^{2}}{(1-\cos\theta)}\end{split}$
(20)

where $q^{2}=-Q^{2}$, and $|T|^{2}_{\rm MM}$ is

$\begin{split}|T|^{2}_{\rm
MM}=\frac{\pi^{2}}{2}\left(\frac{f_{\alpha\beta}}{m_{e}M_{\rm
eff}^{3}}\right)^{2}\alpha^{2}_{\rm
em}\Big{(}E_{\nu}E_{\nu}^{\prime}(1-\cos\theta)-\frac{m_{e}^{2}}{(1-\cos\theta)}\Big{)}(8M_{A}^{2}).\end{split}$
(21)

Using Eq. (11), we obtain the cross section for the incident neutrino
$\nu_{\alpha}$ going to $\nu_{\beta}$ with a MM (including transition MM
also), $f_{\alpha\beta}$ (in units of the electron Bohr magneton),444We use
the incoherent cross sector for $\Delta p\gg 2.5\,\,\mathrm{eV}$.

$\begin{split}\frac{d\sigma}{dT_{e}d\cos\theta}&=\frac{3Z^{3}\alpha_{\rm
em}^{3}m_{e}^{3}E_{\nu}^{\prime}|{\bf
p}_{e}^{\prime}|}{2^{10}\pi^{4}E_{\nu}(m_{e}+T_{e})}4\pi^{2}\left(\frac{f_{\alpha\beta}}{m_{e}M_{\rm
eff}^{3}}\right)^{2}\alpha^{2}_{\rm
em}\Big{(}E_{\nu}E_{\nu}^{\prime}(1-\cos\theta)-\frac{m_{e}^{2}}{(1-\cos\theta)}\Big{)}\\\
&=\frac{3Z^{3}\alpha_{\rm em}^{5}m_{e}^{3}}{2^{8}\pi^{2}M_{\rm
eff}^{6}}\frac{\sqrt{2T_{e}/m_{e}}}{(1+T_{e}/m_{e})^{5/2}}(1-\frac{\delta_{B}}{m_{e}}-\frac{T_{e}}{m_{e}})^{2}\Big{(}(1-\cos\theta)-\frac{m_{e}^{2}}{E_{\nu}E_{\nu}^{\prime}(1-\cos\theta)}\Big{)}\left(f_{\alpha\beta}\right)^{2}\\\
&=0.883\times
10^{-5}\,\mathrm{MeV}^{-3}\frac{\sqrt{2T_{e}/m_{e}}}{(1+T_{e}/m_{e})^{5/2}}(1-\frac{\delta_{B}}{m_{e}}-\frac{T_{e}}{m_{e}})^{2}\Big{(}(1-\cos\theta)-\frac{m_{e}^{2}}{E_{\nu}E_{\nu}^{\prime}(1-\cos\theta)}\Big{)}\left(f_{\alpha\beta}\right)^{2}\end{split}$
(22)

leading to

$\begin{split}\frac{d\sigma}{dT_{e}}=0.883\times
10^{-5}\,\mathrm{MeV}^{-3}\frac{\sqrt{2T_{e}/m_{e}}}{(1+T_{e}/m_{e})^{5/2}}(1-\frac{\delta_{B}}{m_{e}}-\frac{T_{e}}{m_{e}})^{2}\Big{(}2+\frac{2m_{e}^{2}}{E_{\nu}E_{\nu}^{\prime}}\ln\frac{1}{\delta_{\rm
inf}}\Big{)}\left(f_{\alpha\beta}\right)^{2}\end{split}$ (23)

where $\theta_{\rm min}$ is near 0 is used, i.e. $\cos\theta|_{\rm
max}=\frac{1}{2}-\frac{\sqrt{2m_{e}T_{e}}}{<E_{\nu}^{\prime}>}$, which is
integrated over $T_{e}$ and the energy flux of solar $pp$ neutrinos Bahcall04
to give

$\begin{split}\textrm{For MM}:~{}\frac{d\sigma}{dT_{e}}&\simeq 1.071\times
10^{-5}\,\mathrm{MeV}^{-3}\left(2-8.288\right)f_{\alpha\beta}^{2}.\end{split}$
(24)

The infrared divergence in the last term of Eq. (22) is cured by the soft
photon emission processes Bloch37 , and Eq. (24) is used for $T_{e}\geq
2.5\,\mathrm{keV}$. Equation (24) gives a number $-6.73\times
10^{-5}\,\mathrm{MeV}^{-3}f_{\alpha\beta}^{2}$. Comparing with the almost flat
central region ($T_{e}=15-20\,\,\mathrm{keV}$) of Fig. 1 (b), we obtain
$|f_{\alpha\beta}|\leq 0.86\times 10^{-7}$ which is much larger than that
$2.8\times 10^{-11}$ obtained by the Borexino collaboration BdNuMagMom17 . The
huge difference between our estimation and previous ones seems to be from the
over estimation of cross sections in the previous analyses, using the inverse
$T_{e}$ rule and the coherent scattering.

###  III.3 Contributions from vector and axial-vector charges $|T|^{2}_{Q}$

The electromagnetic $F_{1}^{\nu}$ form factor of a neutrino is identified from
Giunti15

$\begin{split}&\Lambda_{\mu}(q)=f_{Q}(q^{2})\gamma_{\mu}-f_{M}(q^{2})\,i\sigma_{\mu\nu}q^{\nu}+f_{E}(q^{2})\sigma_{\mu\nu}q^{\nu}\gamma_{5}+f_{A}(q^{2})(q^{2}\gamma_{\mu}-q_{\mu}q\hskip-4.26773pt/)\\\
&F_{1}^{\nu}(q^{2})=f_{Q}(q^{2})=F_{1}^{\nu}(0)+q^{2}\frac{dF_{1}^{\nu}(q^{2})}{dq^{2}}\Big{|}_{q^{2}=0}+\cdots.\end{split}$
(25)

The charge radius is the $q^{2}$ dependent part,

$\begin{split}\tilde{r}^{2}\equiv\langle
r^{2}\rangle=6\frac{dF_{1}^{\nu}(q^{2})}{dq^{2}}\Big{|}_{q^{2}=0}.\end{split}$
(26)

The elastic cross section includes the interference of the SM weak interaction
and the vector and axial vector couplings of neutrinos with

$\begin{split}G_{V}=\frac{1}{2}+2\sin^{2}\theta_{W}+e^{2}\left(\frac{{\tilde{r}}^{2}/6}{\sqrt{2}G_{F}}-\frac{\varepsilon}{\sqrt{2}G_{F}Q^{2}}\right),\textrm{
and}~{}g_{A}=\frac{1}{2},\end{split}$ (27)

where $Q^{2}=-q^{2}>0$ for the $t$-channel scattering. However, we consider
the 3-body scattering to take into account the momentum of the recoiling atom.
The parameter describing the atomic electromagnetic effects is described in
Appendix A in terms of $B$ and $C$. Inserting (17) into (11), we obtain

$\begin{split}\frac{d\sigma}{dT_{e}}&=\frac{3Z^{3}\alpha_{\rm
em}^{3}G_{F}^{2}m_{e}^{4}}{16\pi^{4}}\cdot\frac{\sqrt{2m_{e}T_{e}}E_{\nu}(E_{\nu}-\delta_{B}-T_{e})}{M_{\rm
eff}^{6}(1+T_{e}/m_{e})}\approx 1.22\times
10^{-23}\cdot\frac{\sqrt{2m_{e}T_{e}}E_{\nu}(E_{\nu}-\delta_{B}-T_{e})}{(1+T_{e}/m_{e})}\cdot\frac{1}{M_{\rm
eff}^{6}}\\\ &=6.33\times
10^{-6}\,\mathrm{MeV}^{-3}\cdot\frac{\big{[}\frac{E_{\nu}}{m_{e}}(\frac{E_{\nu}}{m_{e}}-\frac{\delta_{B}+T_{e}}{m_{e}})\big{]}}{\big{(}1+\frac{T_{e}}{m_{e}}\big{)}^{4}}\sqrt{\frac{2T_{e}}{m_{e}}}\end{split}$
(28)

with $B\simeq 1.144\times 10^{4}$, where the following $|T|^{2}_{Q}$ for the
SM value tHooft71  is

$\begin{split}|T|^{2}_{Q}=\frac{8G_{F}^{2}m_{e}^{2}}{M_{\rm
eff}^{2}}\,\left\\{(G_{V}+g_{A})^{2}E_{\nu}^{2}+(G_{V}-g_{A})^{2}E_{\nu}E_{\nu}^{\prime}+m_{e}E_{e}^{\prime}(g_{A}^{2}-G_{V}^{2})\right\\},\end{split}$
(29)

and we used $1/M_{\rm eff}^{6}$ given in Eq. (44).

###  III.4 Millicharge

Figure 2: The neutrino coupling to photon through the kinetic mixing with ex-
photon $\gamma_{\rm ex}$.

There can be another tiny EM property of neutrinos if there exist another
(almost) massless U(1) gauge boson. If neutrinos couple to this ex-photon,
$\gamma_{\rm ex}$, neutrinos can carry millicharges via the kinetic mixing
Holdom86 , for which the Feynman diagram is shown in Fig. 2. Since the SM
neutrino is in a SM doublet, L-handed electron $e_{L}$ must carry the same
millicharge under the ex-photon gauge group. Then, the magnitude of the
$\nu_{\mu}$ coupling to the ex-photon is limited from the $\nu_{\mu}+e$
scattering data (through neutral current coupling by $Z_{\mu}$),

$\begin{split}(\delta_{\rm ex}e)_{\rm
from\,\nu\,vertex}\frac{1}{q^{2}}(\delta_{\rm ex}e)_{\rm
from\,e\,vertex}\leq\frac{1.167\times
10^{-5}}{\sqrt{2}\,\textrm{GeV}^{2}}\varepsilon(e_{L})=0.413\times
10^{-5}/\,\textrm{GeV}^{2}\end{split}$ (30)

where $\varepsilon(e_{L})(=\varepsilon(\nu_{L}))$ is the neutral current
coupling of $e_{L}$ to the $Z$ boson in the SM KimRMP81 , i.e. $0.5$. Since
there is no way to pinpoint $\delta_{\rm ex}$, we set it as 1, transfering its
uncertainty to the kinetic mixing value $\varepsilon=F_{1}^{\nu}(0)$.

##  IV Fitting

The contributions from the MM and CR of neutrinos add incoherently. For the MM
contribution, we obtained already Eq. (23). For the CR contribution, we should
consider the interference with the weak amplitude in the SM, which is given in
Eq. (28). From Fig. 1 (b), the cross section is less than $5\times
10^{-19}\,\,\mathrm{MeV}^{-3}$ for $T_{e}>17$ keV. For
$\sin^{2}\theta_{W}\simeq 0.231$ sin2LHC18 ; KimRMP81 , Eq. (28) gives, with
$\langle E\rangle=0.2668$ and $\langle E^{2}\rangle=0.0788\,\mathrm{MeV}^{2}$
of solar $pp$ neutrinos, at ${\tilde{r}}^{2}=\varepsilon=0$,

$\begin{split}\textrm{For CR}:~{}\simeq&6.33\times
10^{-6}\,\mathrm{MeV}^{-3}\cdot\Big{\\{}\big{(}0.3018-1.02\,T_{e,\rm
MeV}\big{)}\Big{[}1.462+\Gamma\Big{]}^{2}\\\ &+\big{(}0.3018-1.02\,T_{e,\rm
MeV}\big{)}\Big{[}0.462+\Gamma\Big{]}^{2}\Big{\\}}\frac{2.768\,T_{e,\rm
MeV}^{1/2}}{\big{(}1+1.957\,T_{e,\rm MeV}\big{)}^{4}},\end{split}$ (31)

where $T_{e,\rm MeV}=T_{e}/\,\mathrm{MeV}$, which becomes at $T_{e}=17$ keV

$\begin{split}&0.5703\times
10^{-6}\,\mathrm{MeV}^{-3}\cdot\Big{\\{}\Big{[}1.462+\Gamma\Big{]}^{2}+\Big{[}0.462+\Gamma\Big{]}^{2}\Big{\\}}\to
1.341\times 10^{-6}\,\mathrm{MeV}^{-3}{\rm~{}for~{}}\Gamma=0.\end{split}$ (32)

The ratio of $5\times 10^{-19}\,\,\mathrm{MeV}^{-3}$ and (32) is $\approx
0.790\times 10^{-13}$. For solar neutrino energies of order 400 keV, $Q^{2}$
given in (12) is of order $(0.3-0.4\,\,\mathrm{MeV})^{2}$. Taking $\langle
Q^{2}\rangle\simeq 0.1\,\mathrm{MeV}^{2}$, the charge radius limit from Eq.
(26) is $2.18\times 10^{-6}\,\mathrm{MeV}^{-1}$, or $4.30\times 10^{-17}$cm.
Note that $\Gamma$ in Eq. (31) is

$\begin{split}\Gamma&=4\pi\alpha_{\rm
em}\left(\frac{{\tilde{r}}^{2}/6}{\sqrt{2}G_{F}}-\frac{\varepsilon}{\sqrt{2}G_{F}Q^{2}}\right)=0.9265\times
10^{9}\left(\tilde{r}^{2}\,\mathrm{MeV}^{2}-6\varepsilon_{\rm MeV}\right)\\\
&=2.387\times 10^{30}\left(\tilde{r}^{2}{\rm cm}^{-2}-2.329\times
10^{-21}\varepsilon_{\rm MeV}\right),\end{split}$ (33)

with $\varepsilon_{\rm MeV}=\varepsilon/Q_{\rm MeV}^{2}$ (where $Q_{\rm MeV}$
is $Q$ in units of MeV).

Figure 3: The green bars are the $1\sigma_{\rm stat}$-level backgrounds. The
systematic error is a factor 3 smaller than the statistical error Aprile20 .

If we allow millicharges, the $1\sigma$ error (taken as the $\pm 2$ events out
of 76 in Fig. 1 (a)) is shown as the yellow band in the
$\tilde{r}^{2}-\varepsilon$ plane of Fig. 4. There already exist bounds on the
millicharges of light dark matter from the stellar evolution data, viz. Fig.
12 of CMBeps1  and Fig. 1 of CMBeps2 , presenting $\varepsilon<10^{-12}$ in
the dark matter mass range less than $10^{4}$ eV, which however cannot be used
for our SM neutrinos. The difference is that we studied in this paper the
(almost) massless SM neutrinos which cannot be dark matter in the Universe and
we considered the EM properties of neutrinos only from solar neutrinos.

Figure 4: The bounds on the charge radius (vertical bar) at $\varepsilon=0$
and the $1\sigma$ yellow band in the charge radius vs. millicharge plane,
where $\varepsilon_{\rm MeV}=\varepsilon/Q_{\rm MeV}^{2}$ (where $Q_{\rm MeV}$
is $Q$ in units of MeV) depends on the effective momentum transfer in the
process.

##  V Conclusion

We obtained the bounds on the electromagnetic properties of neutrinos implied
by the XENON1T data: a bound on the magnetic moment $|f_{\alpha\beta}|\leq
0.86\times 10^{-7}$(times the electron Bohr magneton), the charge radius
$|\tilde{r}|<4.30\times 10^{-17\,}{\rm cm}$, and the millicharge bound (Fig.
4) if there exists a massless extra photon.

##  Appendix A Estimate of probability amplitude

If an electron is kicked out from an atomic orbit, an electron in the outer
orbit will fill that vacancy immediately. The transition is occurring
electromagnetically via the $E1$ or $M1$ transition. The the $E1$ transition
rate between the eigenstates $|k\rangle$ and $|n\rangle$ is presented a long
time ago Schiff55 ; Gottfried03

$\begin{split}\frac{(2\pi)^{2}\alpha_{\rm em}}{3}\omega_{kn}^{3}|\langle
k|{\bf r}|n\rangle|^{2},\end{split}$ (34)

where ${\bf r}$ is the coordinate of the EM field. We can estimate $\langle
k|{\bf r}|n\rangle$ for the ground state ($k=n=1$) as $a_{B}/Z$, with
$a_{B}=0.529\times 10^{-8\,}$cm=268.1 MeV-1, and for the excited state $n$ as
$na_{B}/Z$. We use $a_{Z}=14.48\,\mathrm{keV}^{-1}$ for $Z=54$. Then, the rest
of the matrix element takes into account only the spherical harmonics and
gives selection rules between angular momentum eigenfunctions. Now we consider
kicking out the electron outside the atom and the selection rule always
permits it. So, we will consider
$\omega_{nn^{\prime}}=\omega_{n\infty}=|E_{n}/\hbar|$. Thus, the transition
rate to kick out an electron to outside is Gottfried03 ,

$\begin{split}\frac{4\alpha_{\rm em}k^{3}}{3}\frac{1}{2j_{i}+1}|\langle
n_{f}j_{f}||D||n_{i}j_{i}\rangle|^{2}\end{split}$ (35)

Considering the Hydrogenic radial wave function $C_{n}r^{n-1}e^{-r/na}$ with
$C_{n}=(2/n)^{2n+1}/\sqrt{4\pi a_{Z}^{3}}$ with $a_{Z}=a_{B}/Z$, we can
estimate the last factor, between the plane wave and the bound state
$|n\rangle$, $|2^{2n}(n+1)!a_{Z}^{4}/n^{n-1}\sqrt{\pi a_{Z}^{3}V}|^{2}$. Thus,
the transition rate to kick out an electron to outside by electromagnetic
interactions in the atom is

$\begin{split}\frac{4\alpha_{\rm em}k^{3}}{3}\frac{a_{Z}^{5}}{\pi
V}\left(\frac{2^{2n+1}[(n+1)!]}{n^{n-1}}\right)^{2}.\end{split}$ (36)

Let us introduce a parameter $M_{\rm eff}$ reproducing Eq. (36), using the
following effective interaction

$\begin{split}\frac{1}{M^{2}_{\rm
eff}}\bar{e}\gamma^{\mu}e\bar{A}\gamma_{\mu}A,\end{split}$ (37)

and we use particle physicists’ $2\to n$ transition rate,

$\begin{split}Rate=(2\pi)^{4-3n}\frac{d^{3}{\bf p}_{1}^{\prime}\cdots
d^{3}{\bf p}_{n}^{\prime}}{2E_{1}2E_{2}2E_{1}^{\prime}2E_{2}^{\prime}\cdots
2E_{n}^{\prime}}\frac{\delta^{4}(p_{1}+p_{2}-p_{1}^{\prime}\cdots-
p_{n}^{\prime})}{V}|\langle p_{1}^{\prime}\cdots
p_{n}^{\prime}|T|p_{1}p_{2}\rangle|^{2}.\end{split}$ (38)

The two-body scattering rate from Eq. (38), using (37) and averaging over a
flux of incoming particles on a rest target, is

$\begin{split}&(2\pi)^{-2}\frac{d^{3}{\bf
p}_{e}^{\prime}}{2^{4}m_{e}E_{e}^{\prime}M_{A}^{2}}\frac{\delta(\sum_{i}E_{i}-\sum_{f}E_{f})}{V}\frac{2^{6}}{M_{\rm
eff}^{4}}m_{e}E_{e}^{\prime}M_{A}^{2}\\\ &=\frac{4|{\bf
p}_{e}^{\prime}|^{2}d|{\bf
p}_{e}^{\prime}|}{}\frac{\delta(E_{\nu}^{\prime}+m_{e}+\frac{|{\bf
p}_{e}^{\prime}|^{2}}{2m_{e}}+M_{A}-E_{\nu}-m_{e}+\delta_{B}-M_{A})}{\pi
V}\frac{1}{M_{\rm eff}^{4}}=4\pi^{-1}\frac{m_{e}|{\bf p}_{e}^{\prime}|}{M_{\rm
eff}^{4}V}\end{split}$ (39)

where we used

$\begin{split}16(p_{e}^{\mu}{p^{\prime}_{e}}^{\nu}-p^{\prime}_{e}\cdot
p_{e}g^{\mu\nu}+p_{e}^{\nu}{p^{\prime}_{e}}^{\mu}+m_{e}^{2}g^{\mu\nu})(P_{A}^{\mu}{P^{\prime}_{A}}_{\nu}-P^{\prime}_{A}\cdot
P_{A}g_{\mu\nu}+P_{A}^{\nu}{P^{\prime}_{A}}_{\mu}+M_{A}^{2}g_{\mu\nu})=\\\
16(2p_{e}\cdot P_{A}p_{e}^{\prime}\cdot P_{A}^{\prime}-2p^{\prime}_{e}\cdot
p_{e}P^{\prime}_{A}\cdot P_{A}+2p_{e}\cdot P_{A}^{\prime}p_{e}^{\prime}\cdot
P_{A}+2m_{e}^{2}P_{A}\cdot P_{A}^{\prime}\\\ -(P_{A}\cdot
P_{A}^{\prime})(-2p_{e}\cdot p_{e}^{\prime}+4m_{e}^{2})+M_{A}^{2}(-2p_{e}\cdot
p_{e}^{\prime}+4m_{e}^{2})\\\ \to\frac{32}{M_{\rm eff}^{4}}(p_{e}\cdot
P_{A}p_{e}^{\prime}\cdot P_{A}^{\prime}+p_{e}\cdot
P_{A}^{\prime}p_{e}^{\prime}\cdot P_{A}-M_{A}^{2}p_{e}\cdot
p_{e}^{\prime}+m_{e}^{2}M_{A}^{2})\\\ \simeq\frac{32}{M_{\rm
eff}^{2}}\left(3m_{e}^{2}M_{A}^{2}-m_{e}E_{e}^{\prime}M_{A}^{2}\right)\simeq\frac{2^{6}}{M_{\rm
eff}^{4}}m_{e}E_{e}^{\prime}M_{A}^{2}.\end{split}$ (40)

This can be compared to (36), and we obtain

$\begin{split}\frac{1}{M_{\rm eff}^{2}}=\sqrt{\frac{\alpha_{\rm
em}}{3}}\frac{a_{Z}^{5/2}}{\sqrt{m_{e}(m_{e}+T_{e})}}k_{n}^{3/2}\left(\frac{2^{2n+1}[(n+1)!]}{n^{n-1}}\right),\end{split}$
(41)

where we use $k_{n}$ for the threshold value, $k_{n}\simeq|E_{n}(Xe)|/\hbar$.
In practice, let us use Eq. (36) summing over the electrons from each orbit of
the Xenon atom with $k_{n}^{3}$,

$\begin{split}&A\equiv\sum_{n=1}^{5}\,\frac{1}{n^{3}}f(n)\left(\frac{2^{2n+1}[(n+1)!]}{n^{n-1}}\right)^{2}=7.028\times
10^{5},\\\
&B\equiv\sum_{n=1}^{5}\,\frac{1}{n^{6}}f(n)\left(\frac{2^{2n+1}[(n+1)!]}{n^{n-1}}\right)^{2}\simeq
1.144\times 10^{4},\end{split}$ (42)

where $f(n)=2,8,18,18,8$ for $n=1,2,\cdots,5$, respectively. This takes into
account the number of electrons in the Xenon atom $Z$, corresponding to the
incoherent process. $B$ appears when we takes into account $V$ in the
denominator, which is the case in the 3-body scattering. Thus, we estimate,
for the Xenon atom of $Z=54$ and $M_{A}\simeq 7.028\times
10^{5}\,\mathrm{MeV}$,

$\begin{split}\frac{1}{M_{\rm eff}^{2}}=\sqrt{\frac{\alpha_{\rm
em}}{3}}\frac{\sqrt{A}}{Z^{5/2}m_{e}^{2}\sqrt{1+T_{e}/m_{e}}}=0.150(1+T_{e}/m_{e})^{-1/2}\,\,\mathrm{MeV}^{-2}.\end{split}$
(43)

If we use (43), then

$\begin{split}\frac{1}{M_{\rm eff}^{6}}=5.69\times
10^{4}\frac{\,\mathrm{keV}^{-15/2}}{[m_{e}(m_{e}+T_{e})]^{3/2}}\sum_{n}k_{n}^{9/2}\left(\frac{2^{2n+1}[(n+1)!]}{n^{n-1}}\right)^{3}\equiv
2.27\times
10^{5}\frac{\,\mathrm{MeV}^{-3}}{[m_{e}(m_{e}+T_{e})]^{3/2}}\,C\end{split}$
(44)

where $C$ becomes555If we use $Z_{\rm eff}(n)$, including electron
interactions, we obtain $C=5.136\times 10^{12}$, and ${1}/{M_{\rm
eff}^{6}}=1.17\times 10^{18}[m_{e}(m_{e}+T_{e})]^{-3/2}\,\,\mathrm{MeV}^{-3}$.

$\begin{split}C\simeq 2.287\times 10^{12}.\end{split}$ (45)

The $3\to n$ transition rate is

$\begin{split}Rate=(2\pi)^{4-3n}\frac{d^{3}{\bf p}_{1}^{\prime}\cdots
d^{3}{\bf
p}_{n}^{\prime}}{2E_{1}2E_{2}2E_{3}2E_{1}^{\prime}2E_{2}^{\prime}\cdots
2E_{n}^{\prime}}\frac{\delta^{4}(p_{1}+p_{2}+p_{3}-p_{1}^{\prime}\cdots-
p_{n}^{\prime})}{V^{2}}|\langle p_{1}^{\prime}\cdots
p_{n}^{\prime}|T|p_{1}p_{2}p_{3}\rangle|^{2}\end{split}$ (46)

from which we calculate the three body scattering cross section. A flux of
incoming particles 1 (neutrino in our case) scatter on two particles 2
(electron in our case) and 3 (atom in our case). Thus, the cross section
consists of two parts, one with a flux factor $1/|{\bf v}_{1}-{\bf v}_{2}|$
and the other $1/|{\bf v}_{1}-{\bf v}_{3}|$. Taking ${\bf v}_{3}=0$, the
neutrino-electron scattering is

$\begin{split}d\sigma=\frac{(2\pi)^{4-9}}{|{\bf v}_{1}-{\bf
v}_{2}|}\frac{d^{3}{\bf p}_{1}^{\prime}d^{3}{\bf p}_{2}^{\prime}d^{3}{\bf
p}_{3}^{\prime}}{2E_{1}2E_{2}2E_{3}2E_{1}^{\prime}2E_{2}^{\prime}2E_{3}^{\prime}}\frac{\delta^{4}(\sum
p_{i}-\sum p_{f}^{\prime})}{V}|\langle
p_{1}^{\prime}p_{2}^{\prime}p_{3}^{\prime}|T|p_{1}p_{2}p_{3}\rangle|^{2}\end{split}$
(47)

where we can take $V$ as the atomic volume to which the flux of neutrinos
sweep. For $1/V$, the sum $B$ in Eq. (42) will be used below. The $T$-matrix
squared in Eq. (47) is summed over the spins

$\begin{split}\frac{1}{(2s_{e}+1)(2s_{\nu}+1)}\sum_{s_{\nu},s_{e}}\sum_{s^{\prime}_{\nu},s^{\prime}_{e}}|{\rm
Amp}|^{2}=(|T|^{2}_{\rm Q}).\end{split}$ (48)

If we consider the MM interaction also, then we will have a factor
$(|T|^{2}_{\rm MM}+|T|^{2}_{\rm Q})$ instead of $(|T|^{2}_{\rm Q})$. For the
volume $V$, we use the Xenon volume for each $n$ separately. The $n$-th shell
volume we used $V=(4\pi/3)(na_{Z})^{3}$.

###### Acknowledgements.

We thank Kyungwhan Ahn, Elena Aprile, Ki-Young Choi, Paul Frampton, Wonho Jhe,
Junu Jeong, Sin-Kyu Kang, Sejin Kim, and Myungbo Shim for useful discussions.
J.E.K. thanks the APCTP, for the hospitality extended to him during his visit,
where this work has been started. J.E.K. is supported in part by the National
Research Foundation (NRF) grant NRF-2018R1A2A3074631, and S.Y. is supported in
part by the Institute for Basic Science
(IBS-R017-D1-2020-a00/IBS-R017-Y1-2020- 345 a00).

## References

  * (1)
  * (2) W. Pauli, Letter to the T$\ddot{\rm u}$bingen meeting of the German Physical Society, December 4, 1930.
  * (3) F. Reines and C. L. Cowan, _Neutrino magnetic moment upper limit_ , Phys. Rev. 107, 528 (1957) [doi:10.1103/PhysRev.107.528].
  * (4) J. Bernstein and T. D. Lee, _Electromagnetic form factor of the neutrinos_ , Phys. Rev. Lett. 11, 512 (1963) [doi:10.1103/PhysRevLett.11.512].
  * (5) R. B. Clark and R. D. Pedigo, _Forward-peaked $\nu_{e}\,e$ scattering and the solar-neutrino problem_, Phys. Rev. D 8, 2261 (1973) [doi: 10.1103/PhysRevD.8.2261].
  * (6) J. E. Kim, V. S. Mathur, and S. Okubo, _Electromagnetic properties of the neutrino from neutral-current experiments_ , Phys. Rev. D 9, 3050 (1974) [doi: 10.1103/PhysRevD.9.3050].
  * (7) J. E. Kim, _Effects of the transition magnetic moment of the neutrino_ , Phys. Rev. Lett. 41, 360 (1978) [doi: 10.1103/PhysRevLett.41.360].
  * (8) L. B. Okun, M. B. Voloshin, and M. Y. Vysotsky, Sov. J. Nucl. Phys. 44 (1986) 440, Yad. Fiz. 44 (1986) 677–680 [Sov. Phys. JETP 64 (1986) 446–452], and Zh. Eksp. Teor. Fiz. 91 (1986) 754–765.
  * (9) R. Barbieri, R. N. Mohapatra and T. Yanagida, _The magnetic moment of the neutrino and its implications for neutrino signal from SN1987a_ , Phys. Lett. B 213, 69 (1988) [doi: 10.1016/0370-2693(88)91049-0].
  * (10) J. E. Kim, _Neutrino magnetic moment_ , Phys. Rev. D 14, 3000 (1976) [doi: 10.1103/PhysRevD.14.3000].
  * (11) M. Fukugita and T. Yanagida, _Neutrino oscillation in matter induced by charged scalar particles_ , Phys. Lett. B 206, 93 (1988) [doi: 10.1016/0370-2693(88)91268-3].
  * (12) See, for example, C. Giunti and A. Studenikin, _Neutrino electromagnetic interactions: A window to new physics_ , Rev. Mod. Phys. 87, 531 (2015) [e-Print: 1403.6344 [hep-ph]].
  * (13) For a review, see, J. E. Kim, _History of neutrino magnetic moment_ , e-Print: 1911.06883 [hep-ph].
  * (14) Borexino Collaboration (M. Agostini et al. ), _Limiting neutrino magnetic moments with Borexino Phase-II Solar neutrino data_ , Phys. Rev. D 96, 091103 (2017) [e-Print: 1707.09355 [hep-ex]].
  * (15) TEXONO Collaboration (M. Deniz et al. ), _Measurement of $\bar{\nu}_{e}$-electron scattering cross-section with a CsI(Tl) scintillating crystal array at the Kuo-Sheng nuclear power reactor_, Phys. Rev. D 81, 072001 (2010) [e-Print: 0911.1597 [hep-ex]].
  * (16) K. Fujikawa and R. Shrock, _The magnetic moment of a massive neutrino and neutrino spin rotation_ , Phys. Rev. Lett. 45, 963 (1980) [doi: 10.1103/PhysRevLett.45.963].
  * (17) J. Baumann, R. Gaehler, J. Kalus, and W. Mampe, _Experimental limit for the charge of the free neutron_ , Phys. Rev. D 37, 3107 (1988) [doi: 10.1103/PhysRevD.37.3107].
  * (18) E. Aprile et al. (XENON Collaboration), _Excess electronic recoil events in XENON1T_ , Phys. Rev. D 102, 072004 (2020) [e-print:2006.097121 [hep-ex]].
  * (19) J. N. Bahcall and C. Peña-Garay, _Solar models and solar neutrino oscillations_ , New J. Phys. 6, 63 (2004).
  * (20) See, for example, P. Vogel and J. Engel, _Neutrino electromagnetic form factors_ , Phys. Rev. D 39, 3378 (1989) [doi: 10.1103/PhysRevD.39.3378].
  * (21) B. Holdom, _Two U(1)s and $\varepsilon$ charge shift_, Phys. Lett. B 166, 196 (1986).
  * (22) For more references, see, J.-H. Huh, J. E. Kim, J.-C. Park, and S. C. Park, _Galactic 511 keV line from MeV milli-charged dark matter_ , Phys. Rev. D 77, 123503 (2007) [eprint: 0711.3528 [astro-ph]].
  * (23) J. E. Kim, P. Langacker, M. Levine, and H. H. Williams, _A theoretical and experimental review of the weak neutral current: A determination of its structure and limits on deviations from the minimal SU(2) ${}_{L}\times$U(1) electroweak theory_, Rev. Mod. Phys. 53, 211 (1981) [doi:10.1103/RevModPhys.53.211].
  * (24) G. ’t Hooft, _Prediction for neutrino-electron cross-sections in Weinberg’s model of weak interactions_ , Phys. Lett. B 37, 195 (1971) [doi:10.1016/0370-2693(71)90050-5].
  * (25) A. Bodek, _Measurement of the effective weak mixing angle $\sin^{2}\theta_{eff}^{lept}$ from the forward-backward asymmetry of Drell-Yan events at CMS_, [eprint:1808.03170 [hep-ex]]. It has not changed much from the value of 1980’s KimRMP81 .
  * (26) W. S. Cho, J.-H. Huh, I.-W. Kim, J. E. Kim, and B. Kyae, _Constraining WIMP magnetic moment from CDMS II experiment_ , Phys. Lett. B 687, 6 (2010) and 694, 496 (E) (2011) [e-Print: 1001.0579 [hep-ph]].
  * (27) J. D. Vergados, _Searching for light WIMPs via their interaction with electrons_ , J. Phys. G 47, 095007 (2020) [eprint:1911.12956 [hep-ph]].
  * (28) F. Bloch and A. Nordsieck, _Note on the Radiation Field of the Electron_ , Phys. Rev. 52, 54 (1937) [doi:10.1103/PhysRev.52.54].
  * (29) E. Aprile, a private communication.
  * (30) W. L. Xu, C. Dvorkin, and A. Chael, _Probing sub-GeV dark matter-baryon scattering with cosmological observables_ , Phys. Rev. D 97, 103530 (2018) [eprint: 1802.06788 [astro-ph.CO]].
  * (31) E. D. Kovetz, V. Poulin, V. Gluscevic, K. K. Boddy, R. Barkana, and M. Kamionkowski, _Tighter limits on dark matter explanations of the anomalous EDGES 21cm signal_ , Phys. Rev. D 98, 103529 (2018) [eprint: 1807.11482 [astro-ph.CO]].
  * (32) L. I. Schiff, Quantum Mechanics (McGraw-Hill Book Co., New York, 1955) p.403, Eq. (45.22).
  * (33) K. Gottfried and T-M. Yan, _Quantum Mechanics: Fundamentals_ (Springer, Heidelberg, Germany, 2003), p.466.
  * (34) E. Vitagliano, J. Redondo, and G. Raffelt, _Solar neutrino flux at keV energies_ , JCAP 12, 010 (2017) [e-print: 1708.02248 [hep-ph]].


1 1 0 Received: date; Accepted: date; Published: date Stability and evolution
of synonyms and homonyms in signaling game Dorota Lipowska 1,†,∗ and Adam
Lipowski 2,† Adam Lipowski, Dorota Lipowska Correspondence:
<EMAIL_ADDRESS>These authors contributed equally to this work.

###### Abstract

Synonyms and homonyms appear in all natural languages. We analyse their
evolution within the framework of the signaling game. Agents in our model use
reinforcement learning, where probabilities of selection of a communicated
word or of its interpretation depend on weights equal to the number of
accumulated successful communications. When the probabilities increase
linearly with weights, synonyms appear to be very stable and homonyms decline
relatively fast. Such behaviour seems to be at odds with linguistic
observations. A better agreement is obtained when probabilities increase
faster than linearly with weights. Our results may suggest that a certain
positive feedback, the so-called Metcalfe’s Law, possibly drives some
linguistic processes. Evolution of synonyms and homonyms in our model can be
approximately described using a certain nonlinear urn model.

###### keywords:

multi-agent modeling, signaling game, language formation, synonyms, homonyms,
urn model

##  1 Introduction

Evolution and structure of language is often analysed using computational
modelling Cangelosi and Parisi (2002); Nolfi and Mirolli (2010); D’Ulizia et
al. (2020). A particularly appealing research paradigm is inspired by the idea
that language might have spontaneously appeared in a population of
communicating individuals, possibly with some adaptive features Pinker and
Bloom (1990). This standpoint prompted numerous analysis of multi-agent
models, which mimic such communication and try to infer the properties of the
emerging language and its possible further evolution Steels (2012); Gong et
al. (2014); Kirby et al. (2014).

In certain models of this kind, language emergence and evolution is studied
using the signaling game Lewis (2002), where communicating agents must decide
which signal (i.e., a word) to send or how to interpret the signal they have
received. To cope with this, agents very often use some form of the
reinforcement learning Skyrms (2010); Lenaerts et al. (2005); Barrett (2006);
Franke (2016); Mühlenbernd and Franke (2012); Lipowska and Lipowski (2018);
Van Eecke and Beuls (2020). Language that emerges in such models may provide a
unique form-meaning mapping (in a signaling game terminology, it is a
signaling system), but there are also other possibilities. In some cases,
synonyms or homonyms can emerge, destroying thus the unambiguity of the form-
meaning mapping. Neglecting some linguistic nuances Ravin and Leacock (2000),
synonymy means that a single concept can be expressed by different words while
in the case of homonymy, one word carries different meanings.

Of course, synonyms and homonyms should not be considered undesirable or
unrealistic as they exist in virtually every natural language, and are
actually quite common. For example, out of approximately 60,000 entries in
Webster’s Seventh New Collegiate Dictionary, about 35% are homonyms with two
or more meanings Byrd et al. (1987). The thesaurus based on Corpus of
Contemporary American English (containing 385 million words) points to 30,000
words with a total of 370,000 synonyms Davies (2009). Let us notice, however,
that in most cases synonymous words correspond to similar, not the same,
meanings Jackson and Amvela (2007).

It seems quite plausible that synonyms or homonyms change in time. For
example, the frequency of their usage may change and gradually one form will
be preferred over the other, and the latter eventually can even disappear.
Although linguistic data are difficult to interpret, there are some
indications that in natural languages synonyms are quite rare in contrast to
homonyms, which appear to be more common Hurford (2003); Clark (1990). Some
linguists even insist that true synonyms do not exist or at best are very rare
compared to homonyms Lyons (1981); Goldberg (1995). There are some arguments
that the difference in the frequency of synonyms and homonyms may be due to
evolutionary pressures favouring speakers rather than hearers Hurford (2003),
or to language acquisition in childhood Markman (1989).

Let us notice that synonyms actually compete in a quite different way from
homonyms, which can be demonstrated already within the framework of the
signaling game. While synonymous words compete for being selected by a
speaker, for a homonymous word, it is the hearer’s role to assign an
appropriate interpretation. It is thus possible that such a difference can
affect an overall dynamics of synonyms and homonyms and eventually result in
different degrees of their prevalence.

In the present paper, we examine evolution of synonyms and homonyms within the
framework of the signaling game. We show that this evolution depends on the
selection mechanism that is implemented in signaling game. Our results suggest
that to be consistent with linguistic observations, the dynamics of our model
should implement a certain positive feedback, known in some marketing or
economic contexts as Metcalfe’s Law Shapiro and Varian (1998). Qualitatively,
some of our numerical results can be better understood by referring to certain
urn models.

##  2 Model

First, we briefly describe a certain urn model that will help us to understand
some aspects of our mult-agent signaling game.

###  2.1 Nonlinear urn model

Urn models were introduced by Pólya Eggenberger and Pólya (1923) and
intensively studied since then Pemantle (2007). In such models, one considers
an urn with white and black balls. At each step, one of the balls is drawn
randomly from the urn and its colour observed. It is then placed back in the
urn together with an additional ball of the same colour, and the process is
repeated. In the simplest version, the probability to select a ball of a given
colour is proportional to the number of such balls in the urn. Particularly
interesting for us, however, is a nonlinear version of such a model, where the
selection probability is proportional to the number of balls raised to a
certain power $\alpha$ Drinea et al. (2002). The case $\alpha>1$ can be
interpreted as a positive feedback, commonly referred to as Metcalfe’s law
Shapiro and Varian (1998). In this case, the urn becomes dominated by balls of
one colour. For $\alpha<1$, one might say that there is a negative feedback
and the urn tends to the state with an equal number of balls of each colour.
The orignal Pólya urn model is equivalent to the $\alpha=1$ case and it
separates these two different regimes.

To develop a heuristic understanding of the behaviour of such a nonlinear urn
model Drinea et al. (2002), first we denote the number of white and black
balls by $x(t)$ and $y(t)$, where $t$ is the total number of balls in the
system. The probability that the ball selected at time $t+1$ is white equals
$\frac{x(t)^{\alpha}}{x(t)^{\alpha}+y(t)^{\alpha}}$. Thus, the expected change
of $x(t)$ might be written as

$\Delta
x(t)=E[x(t+1)-x(t)]=\frac{x(t)^{\alpha}}{x(t)^{\alpha}+y(t)^{\alpha}},$ (1)

and similarly for $y(t)$. Using the heuristic approximation $\frac{\Delta
y(t)}{\Delta x(t)}=\frac{dy}{dx}$ we obtain

$\frac{dy}{dx}=\frac{y^{\alpha}}{x^{\alpha}}.$ (2)

When $\alpha=1$, the solution of Eq. (2) is $y=cx$, where $c$ is a certain
constant. Integrating Eq. (2) for $\alpha<1$, one obtains that $y/x$ goes to 1
in the long run. For $\alpha>1$, one obtains that $y/x$ goes to 0 or infinity
Drinea et al. (2002). As we will suggest in the following, such a nonlinear
urn model can help us understand the stability and evolution of synonyms and
homonyms, at least within the framework of the signaling game.

###  2.2 Multi-agent signaling game with reinforcement learning

In our model, we have a population of $N$ agents, which play a variant of the
signaling game trying to establish names for $n_{o}$ objects. Agents are
placed at the sites of a network, and in the present paper we assume that it
is a complete graph, where each site is connected, and thus can communicate,
with all the remaining sites. Each agent ($A$) for each object ($o$) has an
inventory, where it stores $n_{w}$ words ($i$) with their corresponding
weights:

$w(A)_{i,o}\ \ \ A=1,\ldots,N\ \ i=1,\ldots,n_{w}\ \ o=1,\ldots,n_{o}$ (3)

In an elementary step, a randomly selected agent (the speaker) communicates
with one of its neighbours (the hearer). The speaker ($S$) chooses an object
($o$) and from a corresponding inventory selects a word ($i$), with a
probability of selection ($p(S)_{i,o}$) proportional to the weight of this
word ($w(S)_{i,o}$) raised to a certain power $\alpha$. More precisely,

$p(S)_{i,o}=\frac{w(S)^{\alpha}_{i,o}}{\sum^{n_{w}}_{j=1}w(S)^{\alpha}_{j,o}}$
(4)

To interpret the communicated word ($i$), i.e., to select an appropriate
object ($o$), the hearer ($H$) takes into account the weights of the
communicated word in all inventories and a similar form of the probability of
selection is used:

$p(H)_{i,o}=\frac{w(H)^{\alpha}_{i,o}}{\sum^{n_{o}}_{q=1}w(H)^{\alpha}_{i,q}}$
(5)

When hearer’s interpretation (i.e., the object it selects) agrees with
speaker’s choice, the reward for their communicative success is that they both
increase by unity the weights of the communicated word in their respective
inventories, which improves the chances of choosing successful words in future
communication attempts (reinforcement learning). To avoid an excessive
increase of weights, which could result in stable configurations, we apply a
population renewal, which means that the selected agent with a certain (small)
probability $p$ does not become a speaker but instead is replaced with a newly
created agent with all weights $w(A)_{i,o}=1$. The population renewal is not
the only way to avoid a standstill as one can use some other techniques such
as memory loss or lateral inhibition Spike et al. (2017).

To examine a time dependance in our model, a unit of time $t=1$ is defined as
$N\cdot n_{o}$ of elementary steps, i.e., in a unit of time for each agent
each object is on average once selected by a speaker.

The model described above has already been used to examine the role of the
network structure in the emergence of linguistic coherence Lipowska and
Lipowski (2018). It was observed that $\alpha>1$ leads to a faster convergence
than $\alpha=1$ but not necessarily to a linguistic coherence state (i.e., a
signaling system). A moderate intensity of a population renewal helps to reach
such a coherence, but too large impedes the process. However, some other
factors, such as the structure of the network of interactions or parameters
$n_{o}$ and $n_{w}$, also affect the evolution of the model Lipowska and
Lipowski (2018). We do not examine our model for $\alpha<1$ because in this
case the probabilities of selection increase sublinearly with accumulated
weights, which inhibits reaching the linguistic coherence. As we show in the
following sections, parameter $\alpha$ qualitatively affects the evolution of
synonyms and homonyms in our model.

The nonlinear urn model described in the previous subsection plays only an
auxiliary role. The signaling game model that we defined in this section can
be considered as a multi-urn generalization of such a single-urn model with
weights $w(A)_{i,o}$ (Eq. (3)) corresponding to the number of balls. Let us
also note that the selection probabilites Eq. (4) and Eq. (5) are determined
by the same set of weights. To select a word for an object, the speaker takes
into account the weights of all words in the corresponding repository (the
summation in Eq. (4) is over words). On the other hand, to interpret the
communicated word, the hearer takes into account the weights of this word in
all repositories (the summation in Eq. (5) is over objects). In some
implementations of the signaling game, separate inventories are used for the
selection of communicated words and for their interpretation Mühlenbernd
(2014).

##  3 Synonyms

To examine the evolution of synonyms, we prepared an initial configuration
with synonyms, where for each object ($i$), agents have a set pair of words
($2i$ and $2i-1$) of weights larger than those of other words in the
respective inventory (for simplicity, we chose two consecutive words for
subsequent objects). An example of such configuration for $n_{o}=2$ and
$n_{w}=5$ is shown in Fig. 1. Initially, thus, agents are inclined to use one
of those two words for an object. To monitor to what extent the synonymous
structure persists in the system, we measured the parameter $m$ defined as

$m=\frac{1}{Nn_{o}}\sum_{A=1}^{N}\sum_{i=1}^{n_{o}}\frac{|w(A)_{2i,i}-w(A)_{2i-1,i}|}{w(A)_{2i,i}+w(A)_{2i-1,i}}$
(6)

This parameter measures the relative difference of weights of synonymous
words. When synonyms persist in the system and have approximately equal
weights, then $m$ is close to 0. When synonyms are eliminated, typically one
of the synonymous words loses importance (i.e., its weight) and $m$ increases
toward unity.

Figure 1: Two examplary agent’s inventories of five words (left column) with
weights (right column; $n_{o}=2$ and $n_{w}=5$). To examine the evolution of
synonyms, we used an initial configuration where in each agent’s inventory,
there is a pair of words with larger weights. When talking about object
$O_{1}$ (or $O_{2}$), agents are initially inclined to communicate using
synonymous words $r_{1}$ and $r_{2}$ (or $r_{3}$ and $r_{4}$).

Numerical simulations show that persistence of synonyms strongly depends on
$\alpha$ (Fig. 2). As simulations with $n_{o}=20$ and $n_{w}=60$ demonstrate,
the characteristic life-time of synonyms increases for decreasing $\alpha$ and
it seems to diverge for $\alpha\rightarrow 1$. For $\alpha=1$ synonyms seem to
be very stable (even up to $t=10^{6}$) and such a behaviour was observed for
several sets of parameters $n_{o}$ and $n_{w}$.

Figure 2:  Time evolution of synonyms as measured with the parameter $m$ (Eq.
(6)). Simulations were done for $N=10^{3}$ agents and the presented results
are averages over 100 independent runs. Note the logarithmic scale on the
horizontal axis.

Such a behaviour of synonyms can be understood referring to the nonlinear urn
model described in section 2.1. Assuming that for a given object the two
competing synonyms have dominant weights, we can neglect other words in this
inventory. Moreover, we can assume that these two words have rather small
weights in inventories corresponding to other objects. The weights of the
first and second synonyms at time $t$ are denoted by $x(t)$ and $y(t)$,
respectively. The average change of $x(t)$ in one elementary step is equal to
the probability that this synonym will be selected by the speaker (which is
$\frac{x(t)^{\alpha}}{x(t)^{\alpha}+y(t)^{\alpha}}$, because the two
synonymous words dominate in the inventory) times the probability that this
word will be correctly interpreted by the hearer (which equals 1, as we
assumed that this word does not have large weights in any other inventory).
Thus the average increase $\Delta x(t)$ has exactly the same form as in Eq.
(1). Since the analogous equation can be written for $\Delta y(t)$, we obtain
that the competition of synonyms in our model is equivalent to the nonlinear
urn model (subject to some simplifying assumptions that we made). From the
analysis of the latter, it is thus clear that for $\alpha>1$, one of the
synonyms gets eliminated while for $\alpha=1$, both of them are likely to
persist.

To illustrate the evolution of the model with synonymous initial conditions,
we plotted the weights of given objects and words (Fig. 3, Fig. 4). The
weights are normalized with the total weight (the sum of weights over all
objects and words). Initally, for each object two words have their weights
larger (than other words in the inventory) and mainly such synonyms are in
use. For $\alpha=1.5$ (Fig. 3), after approximately $25\cdot 10^{3}$ steps, a
certain asymmetry is noticeable and gradually some of these synonymous words
start to dominate. Further evolution leads to the formation of a unique
object-word mapping (i.e, signaling system), where a single word is almost
always used by a speaker to communicate a given object and this word is almost
always correctly interpreted by a hearer. Since $\alpha=1.5$ is not much
greater than unity, the process of elimination of synonyms is relatively slow.

Figure 3:  Evolution of the weights (averaged over $N=10^{3}$ agents) for
given objects and words. Simulations were done for $n_{o}=10$, $n_{w}=30$,
$\alpha=1.5$, and $p=10^{-5}$. The initial configuration of weights includes
10 pairs of synonymous words. After $25\cdot 10^{3}$ steps, some asymmetry in
weights can be noticed and at $t=5\cdot 10^{4}$, an almost perfect unique
object-word mapping can be seen.

Similar calculations for $\alpha=1$ (Fig. 4) show that in this case the pairs
of synonyms exist basically unchanged for at least $t=10^{6}$ steps, in
agreement with the above analysis based on the nonlinear urn model (Fig. 4).

Figure 4:  Evolution of the weights (averaged over $N=10^{3}$ agents) for
given objects and words. Simulations were done for $n_{o}=10$, $n_{w}=30$,
$\alpha=1$, and $p=10^{-5}$. The initial configuration of weights includes 10
pairs of synonymous words. Evolution of the model shows that in this case
synonyms are very stable.

##  4 Homonyms

To examine the evolution of homonyms, we prepared the initial configuration
with homonyms, i.e., this time one word of a large weight is associated with a
chosen pair of objects. An example of such a configuration is shown in Fig. 5.

Figure 5:  Two examplary agent’s inventories of two words (left column) with
weights (right column; $n_{o}=2$ and $n_{w}=2$). To examine the evolution of
homonyms, we used an initial configuration where agent’s inventories are
grouped in pairs containing the same word with a weight larger than those of
other words in these two inventories. When talking about any of the objects
$O_{1}$ or $O_{2}$, agents are initially inclined to communicate using the
word $r_{1}$. However, as a homonym, this word in both inventories has equal
weights and the probability of its correct interpretation is only 0.5.

In the case of homonyms, the communicated word can be incorrectly interpreted
by the hearer, which reduces the success rate $s$ of communicating agents,
where $s$ is defined as a fraction of successfull communication attempts.
Thus, such a parameter is suitable for monitoring the stability of homonyms.

We made numerical simulations for the simplest possible case of $n_{o}=2$ and
$n_{w}=2$ (a single homonym) but we believe that our conclusions concerning
the stability of homonyms may also hold in other cases. We examined several
values of $\alpha$ and the population renewal probability $p=10^{-5}$ or
$10^{-6}$ (Fig. 6). As expected, initially $s$ is close to 0.5 (see Fig. 5),
however, it increases over simulation time, which indicates that homonymy is
gradually eliminated and a unique object-word mapping is formed.

Figure 6:  Time evolution of the success rate $s$ for the initial
configuration with homonyms. Simulations were done for $N=10^{3}$ agents, and
$n_{o}=2$, $n_{w}=2$. The presented results are averages over 100 independent
runs.

Let us notice the nonmonotonic behaviour with respect to the parameter
$\alpha$. For $\alpha=1$, the success rate $s$ rapidly increases, but the
increase is much smaller for $\alpha=1.2$ and 1.5. For $\alpha\geq 2$, the
increase of $s$ is again relatively fast. Since an increase in the success
rate $s$ is related to elimination of homonyms, our simulations show that
homonyms are the most stable for $\alpha=1.2$ and $1.5$.

However, it is in our opinion somewhat surprising that for $\alpha=1$, we
observe a fast increase of $s$, which is much in contrast with the stability
of synonyms in this case. Some understanding of such behaviour can be inferred
from the analysis of a certain urn model, which should approximately reflect
the beginning of homonym elimination. The process starts when in a certain
inventory a non-homonymous word acquires a sufficiently large weight and
further evolution makes it a dominant word, which suppresses a homonym. For
simplicity, let us consider only the inventory, where such a process takes
place (first inventory in Fig. 7). Denoting the weights of the competing words
by $x(t)$ (homonymous word) and $y(t)$ (non-homonymous word), we can write the
following equations that describe their average change in a single step:

$\Delta x(t)=p_{h}\frac{x(t)^{\alpha}}{x(t)^{\alpha}+y(t)^{\alpha}},\ \ \Delta
y(t)=p_{nh}\frac{y(t)^{\alpha}}{x(t)^{\alpha}+y(t)^{\alpha}}$ (7)

The factors $p_{h}$ and $p_{nh}$ are the probabilities that homonymous and
non-homonymous words, respectively, will be correctly interpreted. The terms
$\frac{x(t)^{\alpha}}{x(t)^{\alpha}+y(t)^{\alpha}}$ and
$\frac{y(t)^{\alpha}}{x(t)^{\alpha}+y(t)^{\alpha}}$ in Eq. (7) are the
probabilites to select homonymous and non-homonymous words, respectively. From
the above equations and with the same heuristic replacement as in Eq. (2), we
obtain:

$\frac{dx}{dy}=\frac{p_{h}}{p_{nh}}\frac{x^{\alpha}}{y^{\alpha}}$ (8)

Figure 7:  The non-homonymous word $r_{2}$ acquired a weight larger in the
first inventory (for the object $O_{1}$) than in the second one (for the
object $O_{2}$). When it happen to be communicated by a speaker who has chosen
the object $O_{1}$, it is likely to be correctly interpreted by a hearer,
which will result in a further increase in its weight in the first inventory.

Let us notice that when a non-homonymous word that acquired a large weight in
a certain inventory is communicated, it is very likely that it will be
correctly interpreted (as illustrated in Fig. 7). Thus we expect that $p_{nh}$
is close to 1 and $p_{h}$ is close to 0.5 because at the beginning of the
homonym elimination, both homonyms are symmetric (i.e., of equal weights).
Integrating Eq. (8), for $\alpha=1$ we obtain $x=Cy^{p_{h}/p_{nh}}$, where $C$
is a certain constant. Thus, for $t\rightarrow\infty$ and $p_{h}/p_{nh}<1$
(which is the case, as justified above), the ratio $x/y$ equals 0. In other
words, the weight of a non-homonymous word will dominate that of a homonymous
word. It means that the non-homonymous word, even with a much smaller weight
initially, will overtake a homonym in the long run.

In our opinion, this result is by no means obvious and because the
differential equation approach (Eq. (8)) is based on some heuristic
assumptions, we made simulations of the urn model with the rules corresponding
to Eq. (7) with $\alpha=1$. Namely, there are $x(t)$ white balls and $y(t)$
black balls in the urn. We select randomly one of them. If it is a white one,
then with probability 1/2 we add a white ball to the urn (and with probability
1/2 we do nothing). If the selected ball is black, we add one black ball to
the urn. Numerical results for this urn model support the above analysis (Fig.
8). We can see that even when initial conditions strongly favour $x(t)$
(homonymous word), it is $y(t)$ (non-homonymous word) that dominates in the
long run.

Figure 8:  Time evolution of the average weights $x(t)$ and $y(t)$ for the urn
model with a stochastic addition (see the text). The presented results are
averages over 20 independent runs.

When $\alpha>1$, the prefactor $\frac{p_{h}}{p_{nh}}$ in Eq. (8) is
unimportant and similarly as for the nonlinear urn model Eqs. (1-2), we obtain
that $y/x$ is equal to 0 or $\infty$. However, in the context of the evolution
of homonyms, we have the initial condition that favours $x(t)$ and thus for
$\alpha>1$ one gets $y/x\rightarrow 0$ and homonyms appear to be stable
against words that due to fluctuations acquired a relatively large weight. We
would like to emphasize that this conclusion is based on the heuristic
analysis of the nonlinear urn model. The signaling game with its stochastic
fluctuations and additional factors such as population renewal only to some
extent agrees with this behaviour. Indeed, as demonstrated in Fig. 6, for
$\alpha>1$ homonyms seem to be more stable than for $\alpha=1$, but in the
long run, they are also gradually eliminated.

Finally, we present the evolution of weights in simulations with a larger
number of objects and words. The initial weights detrmine 5 homonyms (each
representing two objects). For $\alpha=1.5$ (Fig. 9) at $t=15\cdot 10^{3}$, we
can notice that for some objects non-homonymous words acquired substantial
weights and they gradually suppress the homonym. Around $t=5\cdot 10^{4}$ all
homonyms are eliminated and a unique object-word mapping is formed. However,
in this case homonyms are relatively long-lived (see Fig. 6). We made similar
calculations for $\alpha=1$ (Fig. 10), in which case homonyms seem to be
eliminated a bit faster. The resulting language contains also some synonyms,
and because $\alpha=1$, they are likely to persist.

Figure 9:  Evolution of the weights (averaged over $N=10^{3}$ agents) for
given objects and words. Simulations were done for $n_{o}=10$, $n_{w}=30$,
$\alpha=1.5$, and $p=10^{-5}$ and the initial configuration of weights
corresponds to 5 homonyms, each with 2 meanings (objects). After $15\cdot
10^{3}$ steps, non-homonymous words for some objects acquire increasing
weights and gradually dominate homonymous ones. Around $5\cdot 10^{4}$, we can
see an almost perfectly unique object-word mapping.  Figure 10:  Evolution of
the weights (averaged over $N=10^{3}$ agents) for given objects and words.
Simulations were done for $n_{o}=10$, $n_{w}=30$, $\alpha=1$ and $p=10^{-5}$
and the initial configuration of weights corresponds to 5 homonyms, each with
2 meanings (objects). Around $3\cdot 10^{5}$, homonyms are basically
eliminated but the emerging language contains synonyms (even multiple).

Fig. 4 and Fig. 10 show that the signaling game with $\alpha=1$ is
characterized by very stable synonyms and a relatively fast elimination of
homonyms, which seems to be different from linguisitc observations of natural
languages. With moderate nonlinearity (e.g., $\alpha=1.5$), both synonyms and
homonyms decay on comparable time scales (Fig. 3, Fig. 9), which seems to be
more adequate.

Finally, let us notice that the initial configurations that we used have
maximal numbers of synonyms or homonyms and differ, of course, from existing
languages, where the frequency of such structures is much smaller. Such a
choice, however, simplifies computational analysis and enables us to monitor
the evolution of syno- and homonyms. We expect that a signaling game with a
more realistic initial configuration should exhibit an analogous behaviour but
its computational analysis would be much more demanding, e.g., due to a much
larger number of objects ($n_{o}$) that would need to be considered.

##  5 Conclusions

In the present paper, we suggest that the presence of synonoyms and homonyms
in natural languages may give us some valuable clues as to the nature of the
mechanisms that drive linguistic processes. Within the framework of the
signaling game, we argue that the reinforcement learning should operate in the
super-linear regime ($\alpha>1$) with probabilities of selections increasing
faster than linearly with the accumulated weights. The linear regime
($\alpha=1$) would, instead, lead to languages with very stable synonyms and
relatively fast decaying homonyms, which is probably at odds with some
linguistic observations Hurford (2003).

We have argued that certain aspects of the dynamics of synonyms and homonyms
in the signaling game can be understood due to the similarity to the nonlinear
urn model Drinea et al. (2002). The relationship is not entirely
straightforward: the monopolistic regime that appears in the nonlinear urn
model translates as a removal of synonyms but for homonyms, it actually
protects their stability. The signaling game is of course more complex than
the urn model and an analysis of the latter provides only a limited insight
into the evolution of synonyms and homonyms in the signaling game. Let us also
note that our model contains several parameters such as $N$, $n_{o}$, $n_{w}$,
$\alpha$, or $p$. Together with the choice of the network of interactions (it
is a complete graph in our study) or the initial values of weights, these are
all the factors influencing the evolution of the model. We did only very
limited numerical analyses for the values of parameters that, in our opinion,
correspond to some generic (within some bounds) behaviour of the model. Of
course, more complete analyses of the model would be desirable.

As we have already mentioned, synonyms, contrary to homonyms, do not reduce
the communicative efficiency in the signaling game. However, such a property
should not be considered as implying that synonyms are stable and homonyms are
not. Although for $\alpha=1$, we indeed observe stable synonyms and unstable
homonyms, the reverse situation takes place for $\alpha>1$ (in which case the
impact of synonyms and homonyms on the communicative efficiency is the same).
As a possible further extension, one could consider the possibility that
synonyms are less stable if there is a cost related to the size of the
inventory used, i.e., when there is a pressure against using more words than
necessary for perfect communication. Moreover, homonymy can be evolutionarily
stable if the signaling game is modified, for example, when there are
additional sources of information that can be used for disambiguation by the
hearer O’Connor (2015); Santana (2014); Mühlenbernd (2021).

Perhaps an interesting question is why nonlinear ($\alpha>1$) rather than
(maybe naively expected) linear ($\alpha=1$) feedback drives linguistic
processes. In marketing or economic contexts, a competition between, e.g.,
video recording formats, operating systems, and even keyboard types is often
discussed, in which cases the value of the system grows probably faster than
linearly with the number of users (Metcalfe’s Law) Shapiro and Varian (1998);
Arthur (1994). In the signaling game, it would mean that a benefit of using a
certain word (and thus a probability of its future selection) increases faster
than linearly with the number of successful communications. Considering the
complexities of language evolution, with its various social, biological, and
cognitive aspects, it seems quite likely.

conceptualization, D.L. ; methodology, A.L.; software, A.L. and D.L.;
validation, A.L. and D.L.; investigation, A.L. and D.L..; writing–review and
editing, A.L. and D.L.; visualization, A.L. All authors have read and agreed
to the published version of the manuscript.”, please turn to the CRediT
taxonomy for the term explanation. Authorship must be limited to those who
have contributed substantially to the work reported.

The authors declare no conflict of interest. References

## References

  * Cangelosi and Parisi (2002) Cangelosi, A.; Parisi, D., Eds.  Simulating the Evolution of Language; Springer-Verlag: London, 2002\.  doi:black10.1007/978-1-4471-0663-0.
  * Nolfi and Mirolli (2010) Nolfi, S.; Mirolli, M., Eds.  Evolution of Communication and Language in Embodied Agents; Springer-Verlag: Berlin, Heidelberg, 2010.  doi:black10.1007/978-3-642-01250-1.
  * D’Ulizia et al. (2020) D’Ulizia, A.; Ferri, F.; Grifoni, P.  A survey on modeling language evolution in the new millennium.  New Gener. Comput. 2020, 38, 97–124.
  * Pinker and Bloom (1990) Pinker, S.; Bloom, P.  Natural language and natural selection.  Behav. Brain Sci. 1990, 13, 707–727.  doi:black10.1017/S0140525X00081061.
  * Steels (2012) Steels, L., Ed.  Experiments in Cultural Language Evolution; Vol. 3, Advances in Interaction Studies, John Benjamins: Amsterdam, 2012.  doi:black10.1075/ais.3.
  * Gong et al. (2014) Gong, T.; Shuai, L.; Zhang, M.  Modelling language evolution: Examples and predictions.  Phys. Life Rev. 2014, 11, 280–302.  doi:black/10.1016/j.plrev.2013.11.009.
  * Kirby et al. (2014) Kirby, S.; Griffiths, T.; Smith, K.  Iterated learning and the evolution of language.  Curr. Opin. Neurobiol. 2014, 28, 108–114.  doi:black10.1016/j.conb.2014.07.014.
  * Lewis (2002) Lewis, D.K.  Convention: A Philosophical Study; Blackwell: Oxford, UK, 2002\.  doi:black10.1002/9780470693711.
  * Skyrms (2010) Skyrms, B.  Signals: Evolution, Learning, and Information; Oxford University Press: Oxford, UK, 2010.  doi:black10.1093/acprof:oso/9780199580828.001.0001.
  * Lenaerts et al. (2005) Lenaerts, T.; Jansen, B.; Tuyls, K.; De Vylder, B.  The evolutionary language game: An orthogonal approach.  J. Theor. Biol. 2005, 235, 566–582.  doi:black10.1016/j.jtbi.2005.02.009.
  * Barrett (2006) Barrett, J.A.  Numerical Simulations of the Lewis Signaling Game: Learning Strategies, Pooling Equilibria, and the Evolution of Grammar.  UC Irvine: Institute for Mathematical Behavioral Sciences, 2006.
  * Franke (2016) Franke, M.  The evolution of compositionality in signaling games.  J. Log. Lang. Inf. 2016, 25, 355–377.  doi:black10.1007/s10849-015-9232-5.
  * Mühlenbernd and Franke (2012) Mühlenbernd, R.; Franke, M.  Simulating the emergence of conventions in small-world networks.  Proceedings of the 21st International Conference on Behavior Representation in Modeling & Simulation (BRIMS 2012), 2012, pp. 37–42.
  * Lipowska and Lipowski (2018) Lipowska, D.; Lipowski, A.  Emergence of linguistic conventions in multi-agent reinforcement learning.  PLoS One 2018, 13, e0208095.
  * Van Eecke and Beuls (2020) Van Eecke, P.; Beuls, K.  Re-conceptualising the Language Game Paradigm in the Framework of Multi-Agent Reinforcement Learning, 2020, [arXiv:cs.AI/2004.04722].
  * Ravin and Leacock (2000) Ravin, Y.; Leacock, C.  Polysemy: Theoretical and Computational Approaches; Oxford University Press: Oxford, UK, 2000.
  * Byrd et al. (1987) Byrd, R.J.; Calzolari, N.; Chodorow, M.; Klavans, J.L.; Neff, M.S.; Rizk, O.A.  Tools and methods for computational linguistics.  Comput. Linguist. 1987, 13, 219–240.
  * Davies (2009) Davies, M.  The 385+ million word Corpus of Contemporary American English (1990–2008+): Design, architecture, and linguistic insights.  Int. J. Corpus Linguist. 2009, 14, 159–190.
  * Jackson and Amvela (2007) Jackson, H.; Amvela, E.Z.  Words, Meaning and Vocabulary: An Introduction to Modern English Lexicology; Bloomsbury Publishing, 2007.
  * Hurford (2003) Hurford, J.R.  Why synonymy is rare: Fitness is in the speaker.  Advances in Artificial Life. Proceedings of the Seventh European Conference on Artificial Life ECAL 2003; Banzhaf, W.; Christaller, T.; Dittrich, P.; Kim, J.T.; Ziegler, J., Eds.; Springer-Verlag: Berlin, Heidelberg, 2003; Vol. 2801, Lecture Notes in Artificial Intelligence, pp. 442–451.  doi:black10.1007/978-3-540-39432-7_47.
  * Clark (1990) Clark, E.V.  On the pragmatics of contrast.  J. Child Lang. 1990, 17, 417–431.
  * Lyons (1981) Lyons, J.  Language, Meaning and Context; Fontana: London, 1981.
  * Goldberg (1995) Goldberg, A.E.  Constructions: A construction Grammar Approach to Argument Structure; University of Chicago Press: Chicago, 1995.
  * Markman (1989) Markman, E.M.  Categorization and Naming in Children: Problems of Induction; MIT Press, 1989.
  * Shapiro and Varian (1998) Shapiro, C.; Varian, H.R.  Information Rules: A Strategic Guide to the Network Economy; Harvard Business Press, 1998.
  * Eggenberger and Pólya (1923) Eggenberger, F.; Pólya, G.  Über die Statistik vorketter vorgänge.  Zeit. Angew. Math. Mech. 1923, 3, 279–289.
  * Pemantle (2007) Pemantle, R.  A survey of random processes with reinforcement.  Probab. Surv. 2007, 4, 1–79.  doi:black10.1214/07-PS094.
  * Drinea et al. (2002) Drinea, E.; Frieze, A.; Mitzenmacher, M.  Balls and bins models with feedback.  Proceedings of the thirteenth annual ACM-SIAM symposium on Discrete algorithms (SODA’02). Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2002, pp. 308–315.
  * Spike et al. (2017) Spike, M.; Stadler, K.; Kirby, S.; Smith, K.  Minimal requirements for the emergence of learned signaling.  Cogn. Sci. 2017, 41, 623–658.
  * Mühlenbernd (2014) Mühlenbernd, R.  Structure and Variation of Signaling Conventions in Scale-free Networks.  Proceedings of the European Conference on Social Intelligence (ECSI-2014); Herzig, A.; Lorini, E., Eds., 2014, pp. 347–358.
  * O’Connor (2015) O’Connor, C.  Ambiguity is kinda good sometimes.  Philos. Sci. 2015, 82, 110–121.
  * Santana (2014) Santana, C.  Ambiguity in cooperative signaling.  Philos. Sci. 2014, 81, 398–422.
  * Mühlenbernd (2021) Mühlenbernd, R.  Evolutionary stability of ambiguity in context signaling games.  Synthese 2021, 198, 11725–11753.
  * Arthur (1994) Arthur, W.B.  Increasing Returns and Path Dependence in the Economy; University of Michigan Press, 1994.