# The mass radius of the proton

Dmitri E. Kharzeev  Center for Nuclear Theory, Department of Physics and
Astronomy, Stony Brook University, New York 11794-3800, USA  Department of
Physics and RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton,
New York 11973-5000, USA

###### Abstract

The mass radius is a fundamental property of the proton that so far has not
been determined from experiment. Here we show that the mass radius of the
proton can be rigorously defined through the formfactor of the trace of the
energy-momentum tensor (EMT) of QCD in the weak gravitational field
approximation, as appropriate for this problem. We then demonstrate that the
scale anomaly of QCD enables the extraction of the formfactor of the trace of
the EMT from the data on threshold photoproduction of $J/\psi$ and $\Upsilon$
quarkonia, and use the recent GlueX Collaboration data to extract the r.m.s.
mass radius of the proton ${\rm R_{m}=0.55\pm 0.03\ fm}$. The extracted mass
radius is significantly smaller than the r.m.s. charge radius of the proton
${\rm R_{C}=0.8409\pm 0.0004\ fm}$. We attribute this difference to the
interplay of asymptotic freedom and spontaneous breaking of chiral symmetry in
QCD, and outline future measurements needed to determine the mass radius more
precisely.

###### pacs:

12.38.Aw;12.40.Yx;13.60.Le

##  I Introduction

The mass distribution is a fundamental property of a physical object. Yet,
while a lot of information is available about the charge distribution inside
the proton, nothing is known at present about its mass radius. In astrophysics
and cosmology, the study of the mass distribution in galaxies has led to
establishing the presence of Dark Matter that is believed to constitute about
$85\%$ of the total mass of matter in the Universe. Drawing an analogy to the
physics of the proton, the electron scattering experiments reveal the spatial
distribution of quarks (matter visible to photons), but do not directly
constrain the spatial distribution of gluons – “dark matter of QCD” that is
not visible to photons. One may thus fully expect that an experimental
determination of the mass distribution would constitute a big advance in the
understanding of the proton structure.

Because of the extreme weakness of the gravitational field created by a single
proton, its direct measurement at short distances is clearly impossible.
Likewise, a study of graviton–proton scattering is off limits for present
experiments. Does this mean that the mass radius of the proton cannot be
measured? We believe that the possibility to measure the mass distribution
inside the proton is provided by the scale anomaly, reflecting the asymptotic
freedom of QCD Gross and Wilczek (1973); Politzer (1973). In the chiral limit
of massless quarks, the scale anomaly Ellis (1970); Chanowitz and Ellis (1972)
expresses the trace of the energy-momentum tensor (EMT) of QCD in terms of the
scalar gluon operator Collins et al. (1977); Nielsen (1977). It has been
proposed Kharzeev (1996); Kharzeev et al. (1999) that the matrix elements of
this operator (which is largely responsible for the mass of the proton) can be
extracted from the photoproduction of heavy quarkonia near the threshold.
Below, we will show how the formfactor of the scalar gluon operator can be
determined from the recent data on photoproduction of $J/\psi$ near the
threshold recently reported by the GlueX Collaboration Ali et al. (2019). We
will then use this formfactor to extract the mass radius of the proton from
the GlueX data.

##  II The mass distribution and gravitational formfactors

As a first step, let us review how the Newton’s law of gravitation emerges
from the Einstein theory Einstein (1915) in the weak gravitational field, non-
relativistic approximation. The Einstein equation reads

$\displaystyle R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=8\pi G\ T_{\mu\nu},$ (1)

where $g_{\mu\mu}$ is the metric tensor, $R_{\mu\nu}$ is the Ricci curvature
tensor, $R$ is the scalar curvature (Ricci scalar), $G$ is the Newton’s
constant, and $T_{\mu\nu}$ is the EMT. We have omitted the cosmological
constant term that is not relevant for our present discussion, and put the
speed of light $c=1$.

Taking the trace w.r.t. the metric on both sides of (1), we get

$\displaystyle-R=8\pi G\ T,$ (2)

where $T\equiv T_{\mu}^{\mu}$ is the trace of the EMT. This equation relates
the trace of the EMT to the scalar curvature of space-time – in fact, this
relation first appeared in the extension of Nordström’s scalar gravity theory
Nordstrom (1912) proposed by Einstein and Fokker Einstein and Fokker (1914);
see Norton (1984) for a historical overview. In classical Maxwell
electrodynamics, $T_{\rm EM}=0$ without massive charges, so the
electromagnetic field does not curve space-time, and does not gravitate.
Moreover, light does not bend in the presence of massive bodies if they induce
only a scalar curvature – and thus the observation of light bending has ruled
out the scalar gravity in favor of the tensor one (1) proposed by Einstein
Einstein (1915) in 1915.

Nevertheless, in weak gravitational fields the trace of EMT and the temporal
component of the EMT $T_{0}^{0}$ coincide. Therefore the distribution of mass
can be obtained from the formfactor of the trace of the EMT – we will call it
for brevity the “scalar gravitational formfactor”, because this would be the
only formfactor in Nordström’s scalar gravity Nordstrom (1912). To show the
equivalence of the distributions of $T$ and $T_{0}^{0}$ in a weak
gravitational field, let us review how the Newtonian limit emerges from the
Einstein theory, see e.g. Landau and Lifshitz (1951).

In the non-relativistic limit

$\displaystyle g_{00}=1+2\varphi,$ (3)

where $\varphi$ is the gravitational field potential, and the EMT is given by

$\displaystyle T_{\mu}^{\nu}=\mu\ u_{\mu}u^{\nu},$ (4)

where $\mu$ is the mass density, and the 4-velocity of the non-relativistic
body can be chosen as $u_{0}=u^{0}=1$, with all spatial components equal to
zero, $u_{i}=0$. Therefore, in this limit

$\displaystyle T_{0}^{0}=\mu;\hskip 28.45274ptT\equiv
T_{\mu}^{\mu}=T_{0}^{0}=\mu,$ (5)

so the distribution of mass and the distribution of the trace of the EMT
indeed coincide.

The equations (1) and (2) lead to

$\displaystyle R_{\mu\nu}=8\pi G(T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T);$ (6)

the temporal component $\mu=\nu=0$ of this equation is

$\displaystyle R_{0}^{0}=4\pi G\mu,$ (7)

and all other components vanish. For the metric (3), we get

$\displaystyle R_{0}^{0}=\frac{\partial^{2}\varphi}{\partial x^{\mu
2}}\equiv\Delta\varphi,$ (8)

and thus (7) yields the equation describing the gravitational field in
Newtonian mechanics:

$\displaystyle\Delta\varphi=4\pi G\mu.$ (9)

Its solution gives the gravitational field potential created by a distribution
of mass $\mu(R)$:

$\displaystyle\varphi=-G\ \int\frac{\mu(R)\ dV}{R}.$ (10)

For a particle of mass $M=\int\mu\ dV$, the total potential is
$\varphi=-GM/R$, and the force acting on a probe mass $m$ is $F_{g}=-m\
\partial\varphi/\partial R$ which yields the Newton’s law of gravity:

$\displaystyle F_{g}=-G\ \frac{mM}{R^{2}}.$ (11)

The purpose of reviewing this textbook derivation was to show that in the weak
gravitational field, non-relativistic, limit of gravity the distribution of
mass and the distribution of the trace of the EMT are identical, see (5).
Therefore, to measure the mass distribution of a microscopic object with a
weak gravitational field, instead of utilizing graviton scattering, we can
measure its scalar gravitational formfactor (i.e. the formfactor of the trace
of EMT).

We can arrive at the same conclusion by comparing the formfactor of the EMT
with the scalar gravitational formfactor. Let us consider first the formfactor
of the EMT for a proton (spin $1/2$ particle of mass $M$) Pais and Epstein
(1949); Pagels (1966):

$\displaystyle\langle{\bf p}_{1}|T_{\mu\nu}|{\bf p}_{2}\rangle$
$\displaystyle=$ $\displaystyle\left(\frac{M^{2}}{p_{01}\
p_{02}}\right)^{1/2}\frac{1}{4M}\
{\bar{u}}(p_{1},s_{1})\Big{[}G_{1}(q^{2})(p_{\mu}\gamma_{\nu}+p_{\nu}\gamma_{\mu})+G_{2}(q^{2})\frac{p_{\mu}p_{\nu}}{M}+$
$\displaystyle{}+G_{3}(q^{2})\frac{(q^{2}g_{\mu\nu}-q_{\mu}q_{\nu})}{M}\Big{]}u(p_{2},s_{2}),$

where $p_{\mu}=(p_{1}+p_{2})_{\mu}$, $q_{\mu}=(p_{1}-p_{2})_{\mu}$,
$p_{1}^{2}=p_{2}^{2}=M^{2}$, and the four-component spinor $u(p,s)$ satisfies
the free Dirac equation $({\hat{p}}-M)u(p,s)=0$ and is normalized according to
$\sum_{s}{\bar{u}}(p,s)u(p,s)=({\hat{p}}+M)/2M$. The formfactors
$G_{i}(q^{2})$ completely describe the mechanical structure of the spin $1/2$
particle.

The energy-momentum conservation

$\displaystyle\partial^{\mu}T_{\mu\nu}=0$ (12)

implies

$\displaystyle q^{\mu}\langle{\bf p}_{1}|T_{\mu\nu}|{\bf p}_{2}\rangle=0;$
(13)

it is easy to check that (II) satisfies the condition (13) for on-shell
nucleons that obey the free Dirac equation; the symmetry of (II) w.r.t. the
interchange of $p_{1}$ and $p_{2}$ is necessary for (13) to hold.

In the limit of vanishing momentum transfer $q_{\mu}\to 0$, the forward matrix
element of the energy-momentum tensor takes the form

$\displaystyle\langle{\bf p}|T_{\mu\nu}|{\bf
p}\rangle=\left(\frac{M^{2}}{p_{0}^{2}}\right)^{1/2}\ {\bar{u}}(p,s)u(p,s)\
\frac{p_{\mu}p_{\nu}}{M^{2}}\left[G_{1}(0)+G_{2}(0)\right],$ (14)

characterized by the $q^{2}=0$ values of the formfactors $G_{1}(q^{2})$ and
$G_{2}(q^{2})$. The Hamiltonian $H$ is given by the temporal component of the
energy-momentum tensor:

$\displaystyle H=\int d^{3}x\ T_{00}(x);$ (15)

therefore, in the rest frame of the particle, the forward matrix element of
$T_{00}$ should yield the mass of the particle:

$\displaystyle\langle{\bf p}=0|T_{00}|{\bf p}=0\rangle=M.$ (16)

Imposing this constraint on (14), we get the condition

$\displaystyle G_{1}(0)+G_{2}(0)=M.$ (17)

The derivation of this condition is completely analogous to the derivation of
the condition on the electromagnetic formfactor $F(q^{2}=0)=e$, where $e$ is
the electric charge $e$ of the particle.

Let us now consider the matrix element of the trace of the EMT $T\equiv
T_{\mu}^{\mu}$; from (II) we find

$\displaystyle\langle{\bf p}_{1}|T|{\bf
p}_{2}\rangle=\left(\frac{M^{2}}{p_{01}\ p_{02}}\right)^{1/2}\
{\bar{u}}(p_{1},s_{1})u(p_{2},s_{2})\ G(q^{2}),$ (18)

with a new formfactor

$\displaystyle
G(q^{2})=G_{1}(q^{2})+G_{2}(q^{2})\left(1-\frac{q^{2}}{4M^{2}}\right)+G_{3}(q^{2})\frac{3q^{2}}{4M^{2}}.$
(19)

In the rest frame of the particle,

$\displaystyle\langle{\bf p}=0|T|{\bf p}=0\rangle=\langle{\bf p}=0|T_{00}|{\bf
p}=0\rangle=M,$ (20)

therefore

$\displaystyle G(0)=M$ (21)

that is obviously consistent with (19) and (17).

##  III The mass radius

The charge radius of the proton is usually defined Feynman (1973) through the
derivative of its electromagnetic formfactor w.r.t. the momentum transfer
$t=q^{2}\equiv-Q^{2}$ evaluated at $t=0$:

$\displaystyle\langle R_{\rm C}^{2}\rangle=-6\ \frac{dG_{\rm
EM}}{dQ^{2}}\Big{|}_{Q^{2}=0}.$ (22)

To enable a direct comparison to the charge radius, we propose to define the
mass radius analogously through the formfactor of mass density $T_{00}$ given
by (II). Let us compare the derivatives of $T_{00}$ and of the scalar
gravitational formfactor (18) w.r.t. $t$ at $t=0$. Because the formfactor of
$T_{00}$ depends on the reference frame (strictly speaking, we can interpret
it in terms of mass distribution only in the rest frame of the proton), we
have to specify it.

It is natural to choose the Breit frame in which ${\bf p}_{2}=\frac{1}{2}{\bf
q}$, ${\bf p}_{1}=-\frac{1}{2}{\bf q}$. Evaluating the derivatives of the
formfactors of $T_{00}$ and of the trace of the EMT w.r.t. $t$ at $t=0$, we
find that they differ by terms $G_{i}(0)/(4M^{2})$ that have to be compared to
$dG_{i}/dt|_{t=0}\equiv G_{i}(0)/m_{i}^{2}$ that depend on the compositeness
scales $m_{i}^{2}$ of the corresponding formfactors.

This difference results from the frame dependence of the formfactor of
$T_{00}$. Indeed, the relativistic $\gamma$ factor for a nucleon moving with
momentum ${\bf p}=\frac{1}{2}{\bf q}$ is
$\gamma=E/M=\sqrt{M^{2}+(q^{2}/4)}/M=\sqrt{1+q^{2}/(4M^{2})}$, so for
$q\equiv|{\bf q}|\simeq m_{i}$, the nucleon is Lorentz-contracted by
$1/\gamma\simeq(1+m_{i}^{2}/(4M^{2}))^{-1/2}$. In the non-relativistic limit,
when $2M\gg m_{i}$, this is a negligible effect – the size $R\sim 1/m_{i}$ of
a massive non-relativistic body is much bigger than its Compton wavelength
$\lambda\sim 1/M$. Because of this, the mass distribution can be defined
through the scalar gravitational formfactor (19) instead of the formfactor of
$T_{00}$:

$\displaystyle\langle R_{\rm m}^{2}\rangle=\frac{6}{M}\
\frac{dG}{dt}\Big{|}_{t=0},$ (23)

where we took into account the normalization (21). This conclusion agrees with
the arguments given above on the basis of weak gravitational field limit of
the Einstein equation. Later we will verify that the compositeness scale of
the scalar gravitational formfactor $m_{s}^{2}\equiv G(0)/(dG/dt|_{t=0})$
extracted from the experimental data is indeed much smaller than $4M^{2}$.

##  IV Scale anomaly of QCD and the mass distribution

Now that we have stablished that the mass radius of the proton can be
extracted from the formfactor of the trace of the EMT, let us discuss how this
formfactor can be measured. The key to this is the scale anomaly of QCD. In
this theory, quantum effects lead to non-vanishing trace of the EMT even for
massless quarks Collins et al. (1977); Nielsen (1977):

$\displaystyle T\equiv T_{\mu}^{\mu}=\frac{\beta(g)}{2g}\ G^{\mu\nu
a}G_{\mu\nu}^{a}+\sum_{l=u,d,s}m_{l}(1+\gamma_{m_{l}}){\bar{q}}_{l}q_{l}+\sum_{h=c,b,t}m_{h}(1+\gamma_{m_{h}}){\bar{Q}}_{h}Q_{h},$
(24)

where $G^{\mu\nu a}$ is the gluon field strength tensor with color index $a$,
the sum in the second and third terms runs over the light and heavy quark
flavors $q_{l}$ and $Q_{h}$ with masses $m_{l}$ and $m_{h}$ respectively, and
$\gamma_{m}$ are the anomalous mass dimensions. The beta-function of QCD Gross
and Wilczek (1973); Politzer (1973) $\beta(g)=\partial g/\partial(log\mu)$
governs the renormalization group running of the QCD coupling $g$ with scale
$\mu$:

$\displaystyle\beta(g)=-b\frac{g^{3}}{16\pi^{2}}+...,\ \
b=11-\frac{2n_{l}}{3}-\frac{2n_{h}}{3}$ (25)

where the first term in $b$ is due to gluon loops, and the second and third
terms are the contributions from light and heavy quark loops.

To determine the mass radius of the proton, we will be interested in the
matrix element of the operator (24) at small momentum transfer
$|t|<4m_{h}^{2}$; the lightest heavy quark is the charm with mass $m_{c}\simeq
1.25$ GeV, so this inequality implies $|t|\leq 6.25\ {\rm GeV}^{2}$. In this
kinematical region, heavy quarks appear only in virtual $Q{\bar{Q}}$ pairs; as
a result, the heavy quark part of (24) cancels the heavy quark contribution to
the gluon term Shifman et al. (1978):

$\displaystyle\sum_{h=c,b,t}m_{h}(1+\gamma_{m_{h}}){\bar{Q}}_{h}Q_{h}\simeq-\frac{2n_{h}}{3}\
G^{\mu\nu a}G_{\mu\nu}^{a}.$ (26)

As a result, the trace of the energy-momentum tensor that will determine the
mass radius of the proton contains only the contributions from light quarks
and gluons:

$\displaystyle T\equiv T_{\mu}^{\mu}=\frac{{\tilde{\beta}}(g)}{2g}\ G^{\mu\nu
a}G_{\mu\nu}^{a}+\sum_{l=u,d,s}m_{l}(1+\gamma_{m_{l}}){\bar{q}}_{l}q_{l},$
(27)

where ${\tilde{\beta}}$ is the beta-function with $b=11-2n_{l}/3=9$ for three
flavors of light quarks, $u,d$ and $s$.

It is well known that the chiral limit of massless quarks provides an accurate
approximation to the physical world; in this limit, the trace of the EMT (27)
contains only the gluon term. Therefore, since the forward matrix element of
(27) according to (20) yields the mass of the nucleon, we have to conclude
that the mass of the proton in the chiral limit is entirely due to gluons. The
contribution of the second term (“$\sigma$-term”) in (27) for physical values
of light quark masses can be extracted from the experimental data on pion and
kaon scattering amplitudes (for recent work, see Ruiz de Elvira et al. (2018))
or computed in lattice QCD Yang et al. (2016); it contributes about $80\ {\rm
MeV}$, or about $8\%$, to the total proton mass – so the chiral limit is
indeed reasonably accurate.

##  V Quarkonium photoproduction near the threshold

In the chiral limit, the information about the mass radius of the proton is
contained in the matrix element of the scalar gluon operator in (27) at non-
zero momentum transfer. The zero-momentum transfer, forward matrix element of
this operator yields the proton’s mass, and this can be used for evaluating
the scattering length in quarkonium-nucleon interaction Luke et al. (1992);
Kaidalov and Volkovitsky (1992); Kharzeev (1996).

At finite momentum transfer, the matrix element of the scalar gluon operator
in (27) can be measured in photoproduction of vector heavy quarkonium states,
$J/\psi$ and $\Upsilon$ close to the threshold Kharzeev (1996); Kharzeev et
al. (1999). This proposal is based on the following arguments:

  1. 1.

Because $J/\psi$ and $\Upsilon$ are made of a heavy quark and an antiquark,
and the proton at small momentum transfer contains only light quarks, the
corresponding amplitude is dominated by the exchange of gluons.

  2. 2.

Close to the threshold, the characteristic size of the heavy quark-antiquark
pair is $\sim 1/(2m_{h})$; for charm quarks, this is about $0.08\ {\rm fm}$.
Because this size is much smaller than the radius of the proton, the coupling
of gluons to the heavy quark is perturbative, is characterized by a small
coupling constant, and can be described by a local color-neutral gluon
operator of the lowest possible dimension Voloshin (1979); Appelquist and
Fischler (1978); Peskin (1979).

  3. 3.

Because of the vector quantum numbers $J^{PC}=1^{--}$ of $J/\psi$ and
$\Upsilon$, the threshold photoproduction is due to the $t-$channel exchange
of gluons in scalar $0^{++}$ and tensor $2^{++}$ states; the scalar exchange
is described by the operator that is proportional to the first term in (27).
Because of the scale anomaly, its matrix element does not depend on the QCD
coupling constant $g^{2}$, whereas the matrix element of the tensor operator
appears proportional to $g^{2}$, and is sub-leading at weak coupling Luke et
al. (1992); Kaidalov and Volkovitsky (1992); Kharzeev (1996); Kharzeev et al.
(1999); Fujii and Kharzeev (1999).

Let us now examine these arguments in more detail, and use them for the
extraction of the mass radius of the proton. Consider the interaction of the
heavy quark pair with the proton near the threshold, where the velocity of the
quarkonium in the center-of-mass is small: $v_{\psi}\ll c$. The coupling of a
small color-neutral heavy quark-antiquark state to gluons can be described by
the operator $g^{2}{\bf E}^{a2}$, where ${\bf E}^{a}$ is the chromo-electric
field – this is the quadratic QCD Stark effect (the first-order effect is
forbidden by color neutrality). The chromo-magnetic contribution is
proportional to $(v_{\psi}/c)^{2}$ and is suppressed; the operators that
contain covariant derivatives are suppressed by the powers of
$(v_{\psi}/c)^{2}$ as well. The $g^{2}{\bf E}^{a2}$ operator can be
identically represented as a sum of the scalar $0^{++}$ and tensor $2^{++}$
gluon operators Novikov and Shifman (1981):

$\displaystyle g^{2}{\bf E}^{a2}=\frac{g^{2}}{2}({\bf E}^{a2}-{\bf
B}^{a2})+\frac{g^{2}}{2}({\bf E}^{a2}+{\bf B}^{a2})=\frac{8\pi^{2}}{b}\
T+g^{2}T_{00}^{(g)},$ (28)

where $T$ is the trace of the EMT (27) in the chiral limit:

$\displaystyle T=\frac{{\tilde{\beta}}(g)}{2g}\ G^{\mu\nu
a}G_{\mu\nu}^{a}=-\frac{bg^{2}}{32\pi^{2}}G^{\mu\nu a}G_{\mu\nu}^{a},$ (29)

and $T_{00}^{(g)}$ is the temporal component of the gluon part of the EMT of
QCD.

The amplitude of $J/\psi$ photoproduction close to the threshold factorizes
into a short-distance part describing the electric polarizability of the
$c{\bar{c}}$ pair, and the matrix element of the operator (28) over a proton,
see Fig. 1 (left):

$\displaystyle{\cal M}_{\gamma P\to\psi P}(t)=-Qe\ c_{2}\ 2M\ \langle
P^{\prime}|g^{2}{\bf E}^{a2}|P\rangle,$ (30)

where $Qe=2e/3$ describes the coupling of the photon to the electric charge of
the charm quark, $c_{2}$ is the short-distance coefficient describing the
coupling of the chromoelectric fields to the heavy quark pair, and its
transition to the $J/\psi$, $t=(P^{\prime}-P)^{2}$ is the momentum transfer,
and the factor $2M$ is needed to reconcile the relativistic normalization of
states with our normalization of the EMT formfactors (16, 20).

The expression (30) holds only near the threshold, where the scalar gluon
operator dominates over the operators that contain derivatives; they are
suppressed by powers of $J/\psi$ velocity squared, $(v_{\psi}/c)^{2}$. Note
that in this kinematical region the scalar gluon formfactor that enters (30)
cannot be interpreted in terms of the gluon structure functions – indeed, the
gluon structure functions are defined through the matrix elements of traceless
gluon operators in the Operator Product Expansion, and we are interested in
the trace part.

Substituting the relation (28) in (30), we observe that the matrix element of
the first term does not contain the coupling $g^{2}$ as a consequence of scale
anomaly, whereas the second term in (28) is suppressed by $g^{2}$ which is
small at the scale of $Q^{2}=4m_{c}^{2}$; in addition, the first term is
enhanced by a numerical factor. Therefore we can re-write the amplitude (30)
in terms of the scalar gravitational formfactor:

$\displaystyle{\cal M}_{\gamma P\to\psi P}(t)=-Qe\ c_{2}\
\frac{16\pi^{2}M}{b}\ \langle P^{\prime}|T|P\rangle.$ (31)

The differential cross section and the integrated cross section of the
$J/\psi$ photoproduction can now be computed using the standard formulae, see
e.g. Group et al. (2020):

$\displaystyle\frac{d\sigma_{\gamma P\to\psi P}}{dt}=\frac{1}{64\pi s}\
\frac{1}{\left|{\bf p}_{\gamma cm}\right|^{2}}\ \left|{\cal M}_{\gamma
P\to\psi P}(t)\right|^{2},$ (32)

and

$\displaystyle\sigma_{\gamma P\to\psi P}(s)=\int_{t_{min}}^{t_{max}}dt\
\frac{d\sigma_{\gamma P\to\psi P}}{dt},$ (33)

where ${\bf p}_{\gamma cm}$ is the photon momentum in the c.m.s. of the
process, and $s=(p_{\gamma}+P_{p})^{2}$ is the square of the c.m.s. energy. We
expect that the short-distance coefficient $c_{2}$ is on the order of $\pi
r_{c{\bar{c}}}^{2}$, where the size of the $c{\bar{c}}$ pair
$r_{c{\bar{c}}}\simeq 1/2m_{c}\simeq 0.08$ fm. We will fit this parameter to
the GlueX Collaboration data Ali et al. (2019), and then check that it is in
the expected range.

Figure 1: Left: the Feynman diagram of $J/\psi$ photoproduction off a proton.
Right: the differential cross section of $J/\psi$ photoproduction at the
center-of-mass energy $E_{cm}=4.58$ GeV (lab. energy of the photon
$E_{\gamma}=10.72$ GeV); the data is from the GlueX Collaboration Ali et al.
(2019); the theory curve corresponds to the dipole form of the scalar
gravitational formfactor with the parameter $m_{s}=1.24\pm 0.07\ {\rm GeV}$,
corresponding to the mass radius of the proton $R_{\rm m}=0.55\pm 0.03$ fm.

Let us briefly discuss the kinematics of the $\gamma+p\to J/\psi+p$ process.
Because of the large mass of $J/\psi$, close to the threshold the process is
characterized by a sizable minimal momentum transfer $t_{min}$; right at the
threshold, its value is $t_{min}=-M_{\psi}^{2}M/(M_{\psi}+M)\simeq-2.23\ {\rm
GeV}^{2}\simeq-(1.5\ {\rm GeV})^{2}$, where $M_{\psi}\simeq 3.097$ GeV is the
mass of $J/\psi$ and $M\simeq 0.938$ GeV is the mass of the proton. The large
magnitude of $t_{min}$ close to the threshold makes the use of the vector
meson dominance model questionable. On the other hand, $t$ is still much
smaller than $4m_{c}^{2}\simeq 6.25\ {\rm GeV}^{2}$ which justifies the
approach based on Eqs (31) and (32). In this kinematical domain, the factor
$c_{2}$ in (31) can indeed be treated as a constant; when the magnitude of $t$
becomes comparable to $4m_{c}^{2}\simeq 6.25\ {\rm GeV}^{2}$, this factor can
be expected to acquire a significant $t$ dependence. Because $t_{min}$ rapidly
varies with the c.m.s. energy close to the threshold, the energy dependence of
the integrated cross section (33) is sensitive the scalar gravitational
formfactor. However the quantity that is most sensitive to the scalar
gravitational formfactor is the differential cross section (32).

The dominance of the scalar gluon operator over the operators with covariant
derivatives in the QCD multipole expansion is justified by the smallness of
the heavy quark pair velocity in the c.m.s., $v_{\psi}$. The operators with
derivatives are suppressed by even powers of $v_{\psi}$ \- therefore, to limit
their contributions by less than about $10\%$, we have to limit $v_{\psi}\leq
0.3$ – this translates into the limit on the c.m.s. energy $E_{cm}\leq 4.6$
GeV, or the photon lab. frame energy of $E_{\gamma}\leq 11.3$ GeV.

##  VI Extracting the mass radius of the proton

To make a direct comparison with the charge radius of the proton which has
been traditionally extracted by using the dipole form factor Feynman (1973),
we will assume, as a first step, a simple dipole parameterization for the
scalar gravitational formfactor as well:

$\displaystyle G(t)=\frac{M}{\left(1-\frac{t}{m_{s}^{2}}\right)^{2}},$ (34)

where $m_{s}$ is the only adjustable parameter. With the standard definition
(23), this parameter relates to the r.m.s. mass radius of the proton in the
following way:

$\displaystyle\langle R_{\rm m}^{2}\rangle=\frac{12}{m_{s}^{2}}.$ (35)

Using (31) with the formfactor $\langle P^{\prime}|T|P\rangle=G(t)$ given by
(34) to evaluate the differential cross section (32), we can now perform the
fit of the recent data from the GlueX Collaboration Ali et al. (2019) at
Jefferson Lab that is available at $E_{\gamma}\simeq 10.72$ GeV, which is
within our desired kinematical range as described above. The resulting fit of
the data is shown in Fig. 1 (right); the extracted value

$\displaystyle m_{s}=1.24\pm 0.07\ {\rm GeV}$ (36)

corresponding to

$\displaystyle R_{\rm m}\equiv\sqrt{\langle R_{\rm m}^{2}\rangle}=0.55\pm
0.03\ {\rm fm}$ (37)

provides an excellent fit with adjusted ${\bar{r}}^{2}=0.99$
(${\bar{r}}^{2}=1.0$ implies a perfect fit).

The corresponding value of the short-distance coefficient in (31) describing
the coupling of the gluons to the heavy quark pair is $|c_{2}|^{2}=0.043\pm
0.006\ {\rm fm}^{4}$. Each gluon couples to the chromoelectric dipole moment
of the the heavy quark pair, and so $c_{2}$ has dimension of length squared.
As discussed above, we expect that $c_{2}$ is on the order of $\pi
r_{c{\bar{c}}}^{2}$, where the size of the $c{\bar{c}}$ pair
$r_{c{\bar{c}}}\simeq 1/2m_{c}\simeq 0.08$ fm. The extracted value of
$|c_{2}|^{2}=0.043\pm 0.006\ {\rm fm}^{4}$ corresponds to
$r_{c{\bar{c}}}\simeq 0.1$ fm, which is in line with our expectations.

We can compare the value of the mass radius (37) to the proton charge radius
Group et al. (2020)

$\displaystyle R_{\rm C}\equiv\sqrt{\langle R_{\rm C}^{2}\rangle}=0.8409\pm
0.0004\ {\rm fm}$ (38)

that is known with a much better precision, see Bernauer (2020) for a recent
review. It appears that the mass radius of the proton is about $50\%$ smaller
than its charge radius! This observation is statistically significant given
the statistical error bar of our value (37). Of course, there is also a
theoretical systematic uncertainty that we will discuss below; however it does
not appear large enough to explain the observed difference between the charge
and mass radii of the proton.

##  VII The mass radius puzzle

At first glance, this difference may seem surprising – but only if one thinks
of a proton as of a charged ball of a fixed radius, and not as of a quantum
object. First, the charge radius is extracted from the coupling of the photon
to quarks, whereas the mass radius results from the coupling to gluons – and
it appears that the gluon radius of the proton is significantly smaller than
its quark radius. If we write down the spectral representation for the scalar
and charge formfactors (see Fig. 2), the radii of the mass and charge
distributions can be seen to be set by the masses of the lightest physical
states excited from the vacuum by the scalar gluon and vector quark currents,
correspondingly. For the vector quark current, the lightest physical state is
the $\rho$ meson with mass of $m_{\rho}\simeq 770$ MeV, whereas for the scalar
gluon current it is the scalar glueball, with a much larger mass of
$m_{G}\simeq 1600$ MeV, see e.g. Athenodorou and Teper (2020).

While this argument does explain a much smaller mass radius, it is too naïve -
in the QCD vacuum, the scalar gluon current strongly mixes with the scalar
quark one. In fact, there exists a low energy theorem Novikov et al. (1981)
that can be used to relate the off-diagonal correlation function of the scalar
gluon and quark currents to the quark condensate. The analysis Fujii and
Kharzeev (1999); Ellis et al. (1999) based on this low energy theorem and the
experimental hadron spectrum shows a very strong mixing of the scalar glueball
state with the $f_{0}(980)$ meson and a broad $\sigma(500)$ resonance in the
$\pi\pi$ spectrum, which is lighter than the $\rho$ meson. Therefore, the
hadron spectrum in the scalar and vector channels alone cannot explain the
difference between the mass and charge radii.

The reason for the smallness of the mass radius in our opinion is the
interplay of scale anomaly and spontaneously broken chiral symmetry. The
dominant contribution at large distances in the dispersion representation of
the proton formfactor stems from the scalar $\pi\pi$ exchange, see Fig. 2,
left. Because the trace of the EMT is invariant under the Renormalization
Group (RG), its matrix element at small momentum transfer (responsible for the
long-range tail of the mass distribution and thus for the mass radius) can be
calculated Voloshin and Zakharov (1980) using the effective chiral theory.
Matching onto the chiral perturbation theory allows to evaluate the matrix
element of the trace of the EMT in a model-independent way; in the chiral
limit $T=-(\partial_{\mu}\pi)^{2}+...$, so we get Voloshin and Zakharov (1980)

$\displaystyle\langle 0|T|\pi^{+}\pi^{-}\rangle=q^{2},$ (39)

where $q^{2}$ is the invariant mass of the $\pi\pi$ pair. In QCD with $N_{f}$
light quark flavors, the r.h.s. of (39) gets multiplied by $(N_{f}^{2}-1)$
after summing over all possible pion combinations, so the corresponding
spectral density of the correlation function of $T$ is

$\displaystyle\rho^{\pi\pi}(q^{2})=\frac{N_{f}^{2}-1}{32\pi^{2}}\ q^{4}.$ (40)

The matrix element (39) and the spectral density (40) are strongly suppressed
at small invariant masses, corresponding to long-range tail of the proton mass
distribution – as a result, the spectral density of the correlation function
of $T$ peaks at masses about $\sim 1$ GeV Fujii and Kharzeev (1999),
consistent with our finding (36). The underlying reason for the suppression of
the matrix element (39) at low $q^{2}$ is the fact that Goldstone bosons
decouple from the scalar curvature induced by the scale anomaly – this is
analogous to the absence of light bending in scalar gravity discussed in
Section II.

At short distances, where the invariant masses in the spectral representation
are large, the relevant matrix element of the trace of the EMT can be computed
using the QCD perturbation theory; in this case, the two-gluon state dominates
(see Fig. 2, right), with the spectral density Fujii and Kharzeev (1999)

$\displaystyle\rho^{pert}(q^{2})=\left(\frac{bg^{2}}{32\pi^{2}}\right)^{2}\
\frac{N_{c}^{2}-1}{4\pi^{2}}\ q^{4},$ (41)

where $b=(11N_{c}-2N_{f})/3$ is the coefficient of the QCD beta-function. We
thus see that the mass distribution at short distances is governed by the
Renormalization Group flow, and at large distances – by the interplay of scale
anomaly and spontaneously broken chiral symmetry. Because pions are Goldstone
bosons, their couplings involve derivatives of the pion field – otherwise they
would not be invariant w.r.t. the chiral rotations. Because of this, at small
momenta pions decouple from the trace of the EMT, and this confines the mass
distribution of the proton to shorter distances, where it is dominated by the
RG flow of QCD.

Figure 2: Spectral representation of the scalar gravitational formfactor.
Left: large-distance, small invariant mass region, where the dominant
contribution is from the pion pair. Right: short-distance, large invariant
mass region that is dominated by gluon pairs.

It would be interesting to decompose the mass distribution into the quark and
gluon ones. The mass decomposition for the proton Ji (1995); Kharzeev (1996);
Yang et al. (2018); Polyakov and Schweitzer (2018); Lorcé (2018); Metz et al.
(2021); Alexandrou et al. (2020) is a subject of a lively debate at present,
and is subject to the renormalization scale and scheme dependence, as well as
to the frame dependence in the case of $T_{00}$. When quark and gluon
contributions are separated, an additional term appears in the decomposition
(II) that corresponds to the anomalous gravimagnetic moment and gravielectric
dipole moment. These terms however should cancel each other in the total TEM
of the nucleon to obey the Einstein equivalence principle, and this is why we
did not consider them.

The formfactor of the gluon term in the TEM decomposition has been recently
evaluated in lattice QCD Shanahan and Detmold (2019). The authors fit this
formfactor by the dipole form, and extract the effective mass of
$m_{s}=1.13\pm 0.06$ GeV – remarkably, this is consistent with the value
$m_{s}=1.24\pm 0.07$ that we have extracted above from the GlueX data Ali et
al. (2019). The scalar and tensor gravitational formfactors that enter the
photoproduction amplitude have also been evaluated in the approaches based on
holography Hatta and Yang (2018); Hatta et al. (2019); Mamo and Zahed (2020).
In this case, the scale of the formfactor is encoded in the dilaton potential
in the bulk that is constructed to reproduce the hadron spectrum and Regge
trajectories.

A number of papers address the proton mass decomposition basing on the vector
meson dominance Kharzeev (1996); Kharzeev et al. (1999); Gryniuk and
Vanderhaeghen (2016); Strakovsky et al. (2020); Gryniuk et al. (2020); Wang et
al. (2020). Near the threshold, the scattering amplitude of quarkonium
possesses a large real part directly related to the scale anomaly that
strongly affects the cross section, which thus becomes sensitive to the mass
decomposition of the nucleon Kharzeev (1996); Kharzeev et al. (1999). However
a big problem of this approach stems from a large value of $t_{min}$ near the
threshold that necessitates taking account of excited $c{\bar{c}}$ states,
with unknown and interfering scattering amplitudes.

Perturbative approaches to the threshold photoproduction of quarkonium
introduce the two-gluon formfactors and relate them to gluon distributions
Frankfurt and Strikman (2002); Zeng et al. (2020). While these approaches are
close to ours in terms of phenomenology, we stress that the scalar
gravitational formfactor that we have considered cannot be interpreted in
terms of gluon structure functions, that are defined through the matrix
elements of traceless gluon operators; instead, they originate from the trace
operator that dominates near the threshold. In the language of Operator
Product Expansion (OPE), the trace terms correspond to the target mass
corrections to the parton model, see e.g. Kharzeev et al. (1996). Non-
perturbative instanton contributions to scalar formfactors have been recently
evaluated in Shuryak and Zahed (2020).

The contribution from $\Lambda{\bar{D}}^{*}$ $t$-channel exchanges has been
studied in Du et al. (2020); in our approach based on the OPE, this
contribution corresponds to higher dimensional operators that should be
suppressed by the heavy quark mass.

##  VIII Theoretical uncertainties and an outlook

The error bar in the mass radius (37) originates only from the precision of
the current data Ali et al. (2019). But what is the “theoretical systematic
error” involved in its extraction? We can categorize the sources of the
uncertainty in (37) as follows:

  * •

The contribution of gluon operators with derivatives: as mentioned above, the
contamination from these operators is suppressed by powers of the $J/\psi$
velocity in the c.m.s., $v_{\psi}$. The GlueX data Ali et al. (2019) on the
differential cross section that we used are at the energy of $E_{\gamma}\leq
11.3$ GeV, corresponding to $v_{\psi}\simeq 0.3$. This means that the
potential contribution of operators with derivatives is about $10\%$. This
contribution can be further reduced by measuring the differential cross
section at a lower energy – for example, at $E_{\gamma}\leq 9.4$ GeV it should
be less than $5\%$.

  * •

Operators of higher dimension, e.g. the quark-gluon operators: the
contribution from these operators is suppressed by the powers of
$\Lambda^{2}/(4m_{c}^{2})$, where $\Lambda$ is a non-perturbative parameter.
Assuming $\Lambda\sim 1$ GeV, we expect these contributions to be on the order
of $\sim 10-15\ \%$. Note that in threshold photoproduction, the convergence
of the multipole expansion is much better than for the interaction of on-shell
quarkonia where it is governed not by the heavy quark mass but by the
quarkonium binding energy Peskin (1979); Kharzeev (1996). The contributions of
higher-dimensional operators can be further suppressed by extending the
studies of threshold photoproduction to the $\Upsilon$ states, that should
become possible with the advent of the Electron Ion Collider. It could also be
possible to study this process in ultra-peripheral collisions at RHIC and LHC.

  * •

The extrapolation in momentum transfer $t$ and the $t$-dependence of the short
distance coefficient $c_{2}$: we have assumed that the coefficient $c_{2}$ is
$t$-independent. This is justified by the fact that the momentum transfer in
the threshold photoproduction is not large enough to resolve the internal
structure of the produced heavy quark pair with the size $\sim 1/(2m_{c})$,
$t\ll 4m_{c}^{2}$. The range of extrapolation to $t=0$ from $t_{min}$ is also
much smaller than $4m_{c}^{2}$, which justifies the assumption of a constant
$c_{2}$.

A careful evaluation and reduction of the uncertainty in (37) will require a
lot of dedicated theoretical and experimental effort. However the mass
distribution is definitely a fundamental property of the proton. Therefore the
measurements of this distribution, combined with measurements of other
“mechanical” properties of the proton, such as the pressure distribution
Burkert et al. (2018), will definitely advance our understanding of the
quantum origin of mass.

Acknowledgements

This work was first presented at the “Origin of the Proton Mass” workshop
organized at ANL by I. Cloët, X. Ji, Z. Meziani, and J. Qiu; I am grateful to
the organizers and participants of this workshop for stimulating discussions.
I also thank J. Bernauer and  Z. Meziani for useful comments on the
manuscript. This research was supported by the U.S. Department of Energy under
awards DE-FG88ER40388 and DE-SC0012704.

## References

  * Gross and Wilczek (1973) D. J. Gross and F. Wilczek, Phys. Rev. Lett. 30, 1343 (1973).
  * Politzer (1973) H. D. Politzer, Phys. Rev. Lett. 30, 1346 (1973).
  * Ellis (1970) J. R. Ellis, Nucl. Phys. B 22, 478 (1970), [Erratum: Nucl.Phys.B 25, 639–639 (1971)].
  * Chanowitz and Ellis (1972) M. S. Chanowitz and J. R. Ellis, Phys. Lett. B 40, 397 (1972).
  * Collins et al. (1977) J. C. Collins, A. Duncan, and S. D. Joglekar, Phys. Rev. D 16, 438 (1977).
  * Nielsen (1977) N. K. Nielsen, Nucl. Phys. B 120, 212 (1977).
  * Kharzeev (1996) D. Kharzeev, Proc. Int. Sch. Phys. Fermi 130, 105 (1996), eprint nucl-th/9601029.
  * Kharzeev et al. (1999) D. Kharzeev, H. Satz, A. Syamtomov, and G. Zinovjev, Eur. Phys. J. C 9, 459 (1999), eprint hep-ph/9901375.
  * Ali et al. (2019) A. Ali, M. Amaryan, E. Anassontzis, A. Austregesilo, M. Baalouch, F. Barbosa, J. Barlow, A. Barnes, E. Barriga, T. Beattie, et al., Physical review letters 123, 072001 (2019).
  * Einstein (1915) A. Einstein, Sitzungsberichte der Preussischen Akademie der Wissenschaften zu Berlin pp. 844–847 (1915).
  * Nordstrom (1912) G. Nordstrom, Physikalische Zeitschrift 13, 1126 (1912).
  * Einstein and Fokker (1914) A. Einstein and D. Fokker, Adriann, Annalen der Physik 44, 321 (1914).
  * Norton (1984) J. Norton, Historical studies in the physical sciences 14, 253 (1984).
  * Landau and Lifshitz (1951) L. Landau and E. Lifshitz, 2 (1951).
  * Pais and Epstein (1949) A. Pais and S. T. Epstein, Rev. Mod. Phys. 21, 445 (1949).
  * Pagels (1966) H. Pagels, Phys. Rev. 144, 1250 (1966).
  * Feynman (1973) R. P. Feynman, _Photon-hadron interactions_ (W.A. Benjamin, Inc., Reading, Massachusetts, 1973).
  * Shifman et al. (1978) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Phys. Lett. B 78, 443 (1978).
  * Ruiz de Elvira et al. (2018) J. Ruiz de Elvira, M. Hoferichter, B. Kubis, and U.-G. Meißner, J. Phys. G 45, 024001 (2018), eprint 1706.01465.
  * Yang et al. (2016) Y.-B. Yang, A. Alexandru, T. Draper, J. Liang, and K.-F. Liu (xQCD), Phys. Rev. D 94, 054503 (2016), eprint 1511.09089.
  * Luke et al. (1992) M. E. Luke, A. V. Manohar, and M. J. Savage, Phys. Lett. B 288, 355 (1992), eprint hep-ph/9204219.
  * Kaidalov and Volkovitsky (1992) A. B. Kaidalov and P. E. Volkovitsky, Phys. Rev. Lett. 69, 3155 (1992).
  * Voloshin (1979) M. B. Voloshin, Nucl. Phys. B 154, 365 (1979).
  * Appelquist and Fischler (1978) T. Appelquist and W. Fischler, Phys. Lett. B 77, 405 (1978).
  * Peskin (1979) M. E. Peskin, Nucl. Phys. B 156, 365 (1979).
  * Fujii and Kharzeev (1999) H. Fujii and D. Kharzeev, Phys. Rev. D 60, 114039 (1999), eprint hep-ph/9903495.
  * Novikov and Shifman (1981) V. A. Novikov and M. A. Shifman, Z. Phys. C 8, 43 (1981).
  * Group et al. (2020) P. D. Group, P. Zyla, R. Barnett, J. Beringer, O. Dahl, D. Dwyer, D. Groom, C.-J. Lin, K. Lugovsky, E. Pianori, et al., Progress of Theoretical and Experimental Physics 2020, 083C01 (2020).
  * Bernauer (2020) J. C. Bernauer, EPJ Web Conf. 234, 01001 (2020).
  * Athenodorou and Teper (2020) A. Athenodorou and M. Teper, Journal of High Energy Physics 2020, 1 (2020).
  * Novikov et al. (1981) V. A. Novikov, M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 191, 301 (1981).
  * Ellis et al. (1999) J. R. Ellis, H. Fujii, and D. Kharzeev (1999), eprint hep-ph/9909322.
  * Voloshin and Zakharov (1980) M. B. Voloshin and V. I. Zakharov, Phys. Rev. Lett. 45, 688 (1980).
  * Ji (1995) X.-D. Ji, Phys. Rev. Lett. 74, 1071 (1995), eprint hep-ph/9410274.
  * Yang et al. (2018) Y.-B. Yang, J. Liang, Y.-J. Bi, Y. Chen, T. Draper, K.-F. Liu, and Z. Liu, Phys. Rev. Lett. 121, 212001 (2018), eprint 1808.08677.
  * Polyakov and Schweitzer (2018) M. V. Polyakov and P. Schweitzer, Int. J. Mod. Phys. A 33, 1830025 (2018), eprint 1805.06596.
  * Lorcé (2018) C. Lorcé, Eur. Phys. J. C 78, 120 (2018), eprint 1706.05853.
  * Metz et al. (2021) A. Metz, B. Pasquini, and S. Rodini, Phys. Rev. D 102, 114042 (2021), eprint 2006.11171.
  * Alexandrou et al. (2020) C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, H. Panagopoulos, and G. Spanoudes, Phys. Rev. D 101, 094513 (2020), eprint 2003.08486.
  * Shanahan and Detmold (2019) P. Shanahan and W. Detmold, Physical Review D 99, 014511 (2019).
  * Hatta and Yang (2018) Y. Hatta and D.-L. Yang, Physical Review D 98, 074003 (2018).
  * Hatta et al. (2019) Y. Hatta, A. Rajan, and D.-L. Yang, Physical Review D 100, 014032 (2019).
  * Mamo and Zahed (2020) K. A. Mamo and I. Zahed, Physical Review D 101, 086003 (2020).
  * Gryniuk and Vanderhaeghen (2016) O. Gryniuk and M. Vanderhaeghen, Physical Review D 94, 074001 (2016).
  * Strakovsky et al. (2020) I. I. Strakovsky, D. Epifanov, and L. Pentchev, Physical Review C 101, 042201 (2020).
  * Gryniuk et al. (2020) O. Gryniuk, S. Joosten, Z.-E. Meziani, and M. Vanderhaeghen, Physical Review D 102, 014016 (2020).
  * Wang et al. (2020) R. Wang, X. Chen, and J. Evslin, The European Physical Journal C 80, 1 (2020).
  * Frankfurt and Strikman (2002) L. Frankfurt and M. Strikman, Phys. Rev. D 66, 031502 (2002), eprint hep-ph/0205223.
  * Zeng et al. (2020) F. Zeng, X.-Y. Wang, L. Zhang, Y.-P. Xie, R. Wang, and X. Chen, The European Physical Journal C 80, 1 (2020).
  * Kharzeev et al. (1996) D. Kharzeev, H. Satz, A. Syamtomov, and G. Zinovev, Phys. Lett. B 389, 595 (1996), eprint hep-ph/9605448.
  * Shuryak and Zahed (2020) E. Shuryak and I. Zahed (2020), eprint 2008.06169.
  * Du et al. (2020) M.-L. Du, V. Baru, F.-K. Guo, C. Hanhart, U.-G. Meißner, A. Nefediev, and I. Strakovsky, The European Physical Journal C 80, 1 (2020).
  * Burkert et al. (2018) V. D. Burkert, L. Elouadrhiri, and F. X. Girod, Nature 557, 396 (2018).


# Color induction in equiluminant flashed stimuli

Xim Cerda-Company and Xavier Otazu
Computer Vision Center, Computer Science Department,
Universitat Autonoma de Barcelona, Spain

(June 2018)

###### Abstract

Color induction is the influence of the surrounding color (inducer) on the
perceived color of a central region. There are two different types of color
induction: color contrast (the color of the central region shifts away from
that of the inducer) and color assimilation (the color shifts towards the
color of the inducer). Several studies on these effects used uniform and
striped surrounds, reporting color contrast and color assimilation,
respectively. Other authors (Kaneko and Murakami, J Vision, 2012) studied
color induction using flashed uniform surrounds, reporting that the contrast
was higher for shorter flash duration. Extending their work, we present new
psychophysical results using both flashed and static (i.e., non-flashed)
equiluminant stimuli for both striped and uniform surround. Similarly to them,
for uniform surround stimuli we observed color contrast, but we did not obtain
the maximum contrast for the shortest ($10\,ms$) flashed stimuli, but for
$40\,ms$. We only observed this maximum contrast for red, green and lime
inducers, while for a purple inducer we obtained an asymptotic profile along
flash duration. For striped stimuli, we observed color assimilation only for
the static (infinite flash duration) red-green surround inducers (red 1st
inducer, green 2nd inducer). For the other inducers’ configurations, we
observed color contrast or no induction. Since other works showed that non-
equiluminant striped static stimuli induce color assimilation, our results
also suggest that luminance differences could be a key factor to induce it.

##  1 Introduction

The color appearance of objects can be influenced by the color of the
surrounding objects. This effect is called color induction and it can be
classified in two different types: color contrast and color assimilation.
Color contrast occurs when the perceived color of the object shifts away from
the surrounding color (inducer), while color assimilation occurs when the
perceived color shifts toward the inducer. Both color contrast and color
assimilation have been studied since the 19th century [13, 2]. Although
assimilation is more common than contrast in daily life [38], contrast has
been more widely studied, mainly using static stimuli (i.e., the stimuli did
not change along time).

Several authors [26, 27, 30] have shown that, in general, uniform surrounds
induce color contrast and striped surrounds tend to induce color assimilation.
In striped surround stimuli, spatial frequency of stripes is a key factor to
induce color assimilation [14, 35], with higher spatial frequencies
($>2.74\,cpd$) leading to a stronger assimilation [10, 30]. Nevertheless,
Smith et al. [35] found that thick ($<0.7\,cpd$) stripes can induce color
contrast.

Color induction has also been studied using dynamic and flashed stimuli [1,
24, 39, 34, 23]. In dynamic stimuli, the inducer is modulated along time,
being the temporal frequency of the surround modulation an important factor
for color induction, e.g. stronger induction at low temporal frequencies,
falling down beyond $2-3\,Hz$ [39, 34]. In flashed stimuli, the target
stimulus is presented during a brief time (a ’blank’ frame is usually shown
when the target stimulus is not presented) [23]. Some of these studies
measured the color induction of afterimages [1, 24]. They showed that color
contrast can produce afterimages and, conversely, color afterimages can induce
color contrast [1]. Furthermore, Kelly and Martinez-Uriegas [24] concluded
that isoluminant chromatic stimuli create isoluminant chromatic afterimages.

Recently, Kaneko and Murakami [23] published an extensive work in color
induction using equiluminant flashed color stimuli with uniform surrounds.
They measured the color induction at different flash durations (from $10\,ms$
to $640\,ms$) and showed that color contrast significantly depends on the
duration of the flash. They concluded that the shorter the flash duration, the
stronger the contrast. Since they used uniform surrounds, only color contrast
was reported.

In this work, we extended Kaneko and Murakami’s work [23], presenting a new
study of color induction using both uniform and striped surrounds and both
static and flashed stimuli (see several static stimuli examples in Figure 1).
Similarly to other works [26, 27, 30, 10, 14, 35, 23, 1, 24, 39, 34], we
expect to observe color contrast in uniform surround stimuli and color
assimilation in striped surround stimuli. Furthermore, we expect to reproduce
the Kaneko and Murakami’s results [23] for flashed uniform surrounds and to
analyze whether color induction depends on flash duration for both striped and
uniform surrounds. In previous papers [31, 30, 32] one author of this paper
simultaneously reproduced psychophysical results of both color and brightness
induction using a Wavelet model and a neurodynamical model of V1. These models
suggest that color contrast and color assimilation could be the result of the
same mechanism (lateral connections) [45, 44, 32].

Thus, although uniform and striped surrounds could induce opposite effects
(color contrast and color assimilation, respectively), our hypothesis is that,
when the chromatic surrounds are flashed, the temporal behavior of the induced
color would be similar in both cases.

Figure 1: Examples of some stimuli used in this work (see Table 1). In this
figure, we can observe that the achromatic ring (surrounded by 8 black dots)
is perceived differently depending on the surround chromaticity. This effect
is called color induction. For instance, on the top-left figure, the
achromatic ring surrounded by a green inducer is perceived as reddish, while
the achromatic ring surrounded by a purple inducer (bottom-left figure) is
perceived as lime-ish.

##  2 Methods

###  2.1 Apparatus

All experiments were conducted in a dark room, on a calibrated
$21^{\prime\prime}$ SONY GDM-F500R CRT monitor ($1024\times 768\,px$, $100Hz$)
with a viewable image size of $19.8^{\prime\prime}$. The display was viewed
binocularly and freely from an approximate distance of $132\,cm$, subtending
around $17.3\times 13.0\,deg$ of visual angle for the observer. The monitor
was connected to a Wildcat Realizm R500 PCI Express graphics card through a
digital video processor (Cambridge Research Systems ViSaGe MKII Stimulus
Generator) capable of displaying 14-bit color depths. The monitor was
calibrated via a customary software for the stimulus generator (Cambridge
Research Systems, Ltd., Rochester, UK) and a ColorCal (Minolta sensor)
suction-cup colorimeter.

###  2.2 Stimuli

The software was implemented in Matlab (The MathWorks, Inc.), and the video
processor was managed using the Cambridge Research System custom-made toolbox.
We used the same spatial configuration of visual stimuli as Otazu et al. [30],
which was inspired by Monnier and Shevell [26, 27]. In this work, we added a
temporal component, following Kaneko and Murakami [23] flash duration values.

All stimuli were defined in the MacLeod and Boynton color space [5], which is
based in the Smith and Pokorny cone fundamentals [36]. In this opponent space,
the $l$ axis represents the red-green opponency (i.e., ’L vs M’ cone
opponency) and the $s$ axis represents the purple-lime opponency (i.e., ’S vs
(L+M)’ cone opponency), where $s$ is normalized to unity equal-energy white
[5].

####  2.2.1 Spatial Configuration

Several stimuli examples are shown in Figure 1 and a schematic of the
stimuli’s spatial configuration is shown in Figure 2. The test frame was
composed by two circularly symmetric patterns (i.e., the test and the
comparison stimuli) presented side by side and separated by $8.68\,deg$ of
visual angle from the observer’s point of view. Similarly to Otazu et al.
[30], the stimuli’s background was dark and the test ring (in the test
stimulus) was achromatic and surrounded by concentric rings (inducers) of
spatially alternating colors (the 1st and 2nd inducer, according to the
physical distance to the test ring). When these two inducers had the same
chromaticity, the surround was a uniform region. The striped surround was
built with 11 circular stripes (stripes visual frequency was $1.94\,cpd$)
because, as observed in Otazu et al. [30], they produce more color induction
than 5 stripes and are not as thin as 17 stripes. In fact, in the 17 stripes
case, observers reported that for the shortest flashes they could not detect
the test ring. To make the detection of the test ring easier, we drew 8 black
dots of 1 pixel size: 4 dots in the inner radius of the ring and 4 points in
the outer radius (at 0, 90, 180 and 270°). The comparison ring, on the right
side of the frame, was always surrounded by a uniform achromatic disk
approximately metameric to equal-energy white ($l=0.66$, $s=0.98$ and
$Y=20\,cd/m^{2}$) [27].

All chromaticities lay down on the subject’s equiluminant plane. The
calibration of subjects’ equiluminant plane was performed using the Minimally
Distinct Border (MDB) [6, 21, 40, 22, 37, 8, 4]. Thus, for each subject, all
the colors of the stimuli were shown on his/her perceptually equiluminant
surface.

Figure 2: Spatial configuration of the stimuli. Subjects had to adjust the
chromaticity of the comparison ring in order to match perceptually the
chromaticity of the test ring. The colors in this figure are illustrative. All
experimental conditions are described in Table 1.

####  2.2.2 Temporal Configuration

In Figure 3, we show the temporal configuration of the stimuli. Two different
frames were defined: the test frame and the blank frame. During the blank
frame, the test stimulus was an achromatic disk of the same intensity as the
test ring. During the test frame, the test stimulus (either striped or
uniform) was flashed while the achromatic test ring was not modified. That is,
the achromatic test ring remained constant along the experiment and only the
chromatic surrounding rings were flashed. The time duration of our blank frame
was $1\,s$ and the flash duration took values from $10$ to $320\,ms$ in a
dyadic temporal frequency sequence (i.e., $10$, $20$, $40$, $80$, $160$ and
$320\,ms$, $N_{tc}=7$, see temporal conditions in Table 1). The stimulus
sequence (i.e., blank and test frames) was repeated until subjects finished
the task (see Section 2.4).

In addition, we also used static stimuli, which was equivalent to infinite
flash duration of the test frame.

Figure 3: Temporal configuration of the stimuli. The blank frame was shown
during $1\,s$ and the duration of the flash (when the test frame was
presented) depended on the temporal condition. All experimental conditions are
described in Table 1.

All the experimental conditions (both spatial and temporal) are shown in Table
1. The color of the inducers were chosen to have a good representation of the
color space.

Table 1: Summary of the experimental conditions (both spatial and temporal).
In the spatial conditions, we detail the chromaticity sets in MacLeod and
Boynton color space. Only $l$ and $s$ chromatic axes are reported because all
chromaticities lay down on the subject’s equiluminant plane
($Y=20\,cd/m^{2}$). In the temporal conditions, we indicate the flash duration
in $ms$. During the flash, the test frame was presented. In particular, in the
static condition only the test frame was presented (unlimited duration).
Experiment U (Uniform surround)
---
Spatial Conditions
| Test Ring |  | 1st Inducer |  | 2nd Inducer
| l | s |  | l | s |  | l | s
1 | $0.66$ | $0.98$ |  | $0.69$ | $0.98$ |  | $0.69$ | $0.98$
2 | $0.66$ | $0.98$ |  | $0.66$ | $1.38$ |  | $0.66$ | $1.38$
3 | $0.66$ | $0.98$ |  | $0.63$ | $0.98$ |  | $0.63$ | $0.98$
4 | $0.66$ | $0.98$ |  | $0.66$ | $0.58$ |  | $0.66$ | $0.58$
Temporal Conditions (in $ms$)
| 10 | 20 | 40 | 80 | 160 | 320 | static |
|  |  |  |  |  |  |  |
Experiment S (Striped surround)
Spatial Conditions
| Test Ring |  | 1st Inducer |  | 2nd Inducer
| l | s |  | l | s |  | l | s
1 | $0.66$ | $0.98$ |  | $0.69$ | $0.98$ |  | $0.63$ | $0.98$
2 | $0.66$ | $0.98$ |  | $0.66$ | $1.38$ |  | $0.66$ | $0.58$
3 | $0.66$ | $0.98$ |  | $0.63$ | $0.98$ |  | $0.69$ | $0.98$
4 | $0.66$ | $0.98$ |  | $0.66$ | $0.58$ |  | $0.66$ | $1.38$
Temporal Conditions (in $ms$)
| 10 | 20 | 40 | 80 | 160 | 320 | static |

###  2.3 Subjects

The experiments were done by six observers ($N_{sub}=6$), four of them from
our lab (AA, DB, LR and XO), who were familiar with color spaces, and two
others who were not related to the lab (BG and CM). All of them were informed
of the aim of the experiments and consented to participate in the
experimentation. Five observers were completely naïve (AA, BG, CM, DB and LR),
while one of them is one of the authors of the paper (XO). Three of them are
male (AA, DB and XO) and the other three are female (BG, CM and LR). The ages
were comprised between 22 and 45 years old. All of them had normal or
corrected-to-normal vision, tested using the Ishihara [18] and the D-15
Farnsworth Dichotomous Test [15]. To learn the experimental procedure, all
observers did a one day training session before starting the experiments.

###  2.4 Experimental Procedure

The subjects’ task was to adjust the chromaticity of the comparison ring until
they perceived it equal to the test ring. They were instructed to do the
matching according to the test ring color perceived during the test frame,
ignoring the after-effect produced during the blank frame.

We conducted two different experiments (Experiments U and S) in accordance
with the Code of Ethics of the World Medical Association (Declaration of
Helsinki) to study how the surround chromaticity and the flash duration
influence the color perception. Since color induction strongly depends on the
surround type, we divided the experiments according to it. In Experiment U,
the test ring had a uniform surround and, in Experiment S, the test ring had a
striped surround.

To reduce the available color space to one-dimension, similarly to Kaneko and
Murakami [23], we performed a previous experiment where subjects were able to
adjust the chromaticity in a MacLeod and Boynton two-dimensional color space.
We observed that in both experiments (Experiment U and S) the observations
were approximately on the cardinal axis which includes the test and inducers’
chromaticities. Thus, we reduced the available color space to one-dimension,
i.e., the observers only changed the comparison ring chromaticity along $l$ or
$s$ axes, depending on the experimental spatial condition (1-D Experiment).

For each 1-D Experiment (Experiments U and S) we had 28 different experimental
conditions (4 different spatial conditions -$N_{sc}=4$\- and 7 temporal
conditions -$N_{tc}=7$-, see Table 1). Each subject evaluated each
experimental condition 10 different times (i.e., 10 different observations,
$N_{obs}=10$). An experimental condition was determined by the combination of
a spatial condition and a temporal condition. Each run started with 3 minutes
of dark adaptation [23] and subjects performed 48 different observations (8
different spatial conditions, 3 different temporal conditions, and 2
repetitions). Subjects did not have any time restriction, but they were
advised not to take more than 1 minute for each experimental condition. On
average, each run took about 30 minutes. Subjects evaluated the static
conditions of both experiments on two days apart. In that case, they performed
24 observations in each run (8 different spatial conditions, 1 temporal
condition -static- and 3 repetitions, except for the last run, which was 4
repetitions), taking about 15 minutes.

To minimize memory effect of the experiments, we defined a pseudo-random order
of the temporal conditions in each run, while the spatial conditions were
randomized.

##  3 Results

###  3.1 Metric

To represent the induction strength, we define the perceptual color induction
metric

$\Delta C_{i}=\frac{C^{c}_{i}-C^{t}_{i}}{C^{s}_{i}-C^{t}_{i}}\;,$ (1)

where $i$ is the chromatic axis of the MacLeod-Boynton color space along which
the induction was measured ($i=[l,s]$). $C^{c}_{i}$ is the scalar value of the
$i$th axis in the lsY color space of the comparison ring (the observation done
by the subject); $C^{t}_{i}$ is the scalar value of the $i$th axis of the test
ring; and $C^{s}_{i}$ is the scalar value of the $i$th axis of the 1st
inducer. Note that, since it is a 1-D Experiment, all the terms of the
equation are scalar because it only considers the axis where the chromaticity
could be adjusted. This metric is sensitive to both color contrast and color
assimilation: when $\Delta C$ is positive, color assimilation is induced
because $C^{c}_{i}$ is shifted towards $C^{s}_{i}$; and when $\Delta C$ is
negative, color contrast is induced because $C^{t}_{i}$ is shifted away from
$C^{s}_{i}$. Nevertheless, $|\Delta C|$ has to be greater than the just
noticeable difference (JND) to consider that color induction is produced. The
JND region was computed in the CIELab color space ($\Delta E=1$), which is an
approximately perceptually uniform color space, and then transformed to the
MacLeod-Boynton color space. Thus, when $-JND<\Delta C<JND$ we consider that
there is no induction. This metric does not include the 2nd inducer because
only the 1st inducer determines the type of color induction (color contrast or
color assimilation).

For each experimental condition, we averaged all the 10 observations of each
subject and computed the average and the standard error of means (SEM) of all
6 subjects ($N_{sub}=6$).

###  3.2 Experiment U

In this experiment, we studied the color induced by equiluminant uniform
surrounds on an achromatic test ring at different temporal conditions (see
Table 1). The results (Figure 4) show that chromatic contrast was induced in
all experimental conditions, except at the $10\,ms$ flash in the purple-lime
opponent axis. In fact, under these experimental conditions (Spatial
Conditions 2 and 4 -purple and lime inducers- flashed during $10\,ms$),
observers pointed out that the test frame detection was very difficult.

We used a nested ANOVA analysis, which had 6 and 378 degrees of freedom
($N_{tc}-1$ and $N_{tc}N_{sub}(N_{obs}-1)$, respectively) to observe whether
temporal conditions affected to color induction. Once the nested ANOVA
indicated that there were significant differences, a Fisher’s Least
Significant Differences post-hoc analysis (Fisher’s LSD) was performed to
group the temporal conditions according to the color induction they induced.
These analyses showed that for Spatial Conditions 1, 3 and 4, maximum color
contrast induction was produced by $40\,ms$ flash, while in the Spatial
Condition 2 there is no peak at $40\,ms$. Furthermore, for all spatial
conditions there are no significant differences between the perceived colors
at $80$ and $160\,ms$. Moreover, the induction produced by a $320\,ms$ flash
was similar to the induction produced by a static stimulus. All the ANOVA
statistics’ details are shown in Table 2 and the letters below error bars in
Figure 4 show the temporal conditions that induced statistically similar
colors (Fisher’s LSD post-hoc analysis’ results). The temporal conditions
which have the same letter can be considered that induce the same perceptual
color.

Figure 4: Results of Experiment U. The averaged (mean of 6 observers) color
induction metric ($\Delta C$) is plotted against the flash duration (7
different temporal conditions). Separate plots correspond to the different
spatial conditions (see Experiment U in Table 1). The gray area is the region
where subjects could not perceive chromatic differences ($\Delta
C\in[-JND,JND]$). Error bars indicate $\pm 1\,SEM$. An ANOVA analysis of the
data (see Table 2) showed that there were significant differences in color
induction strength for different temporal conditions. Fisher’s post-hoc
analysis (letters below the error bars), which allowed us to measure which
temporal conditions were significantly different, stressed that the peak of
color induction was always perceived at $40\,ms$, except for Spatial Condition
2. Furthermore, static stimulation induced the same color as the longest flash
duration ($320\,ms$) and the perceived color at $80$ and $160\,ms$ did not
vary. In all chromatic conditions, chromatic contrast ($\Delta C<0$) or no
induction was induced.

###  3.3 Experiment S

In this experiment, we studied the color induced by equiluminant striped
surrounds on an achromatic test ring under different temporal conditions (see
Experiment S in Table 1). In Figure 5, we can observe that Spatial Conditions
1 and 4 did not induce any color induction when the stimuli were flashed. Only
one out of 6 subjects observed assimilation in flashed Spatial Condition 1,
and 2 out of 6 perceived assimilation in flashed Spatial Condition 4. By
contrast, static stimulus of Spatial Condition 1 induced chromatic
assimilation ($\Delta C>0$). In Spatial Conditions 2 and 3 only color contrast
($\Delta C<0$) was perceived.

Subjects pointed out that, again, they left a gray color for Spatial
Conditions 2 and 4 (purple-lime axis) flashed during 10 ms because they were
not able to see the test frame.

Similarly to Experiment U, ANOVA analysis showed significant differences
between the induction produced by flashes of different durations in all
spatial conditions (see Experiment S in Table 2). We have not observed any
grouping or behavioral pattern over all spatial conditions in the Fisher’s LSD
post-hoc analysis. Since the profile of the results of Figures 4 and 5 are
completely different and the strongest induction in this experiment is not
observed at the same flash duration as in the previous one, our initial
hypothesis should be rejected.

Figure 5: Results of Experiment S. The averaged (mean of 6 observers) color
induction metric ($\Delta C$) is plotted against the flash duration. Separate
plots correspond to the different spatial conditions (see Experiment S in
Table 1). The gray region indicates where induction could not be perceived
($[-JND,JND]$). Error bars indicate $\pm 1\,SEM$. An ANOVA analysis of the
data (see Experiment S in Table 2) showed that there were significant
differences between color induction for different temporal conditions. Like in
Experiment U, we did a Fisher’s LSD post-hoc analysis (letters below error
bars) to know the temporal conditions that induced significantly different
colors. In this Experiment, we observed that the static stimuli induction had
significant differences with respect to the induction produced by a $320\,ms$
flash. No color assimilation had been induced, except by static stimulus in
Spatial Condition 1. Table 2: Summary of the nested ANOVA results. These
results show that, in all spatial conditions, color induction depended on the
temporal condition of the stimulus (i.e., the flash duration). The nested
ANOVAs had 6 and 378 degrees of freedom, corresponding to $N_{tc}-1$ and
$N_{tc}N_{sub}(N_{obs}-1)$, respectively. Experiment U (Uniform surround)
---
Spatial | $F_{6,378}$ | p
Condition
1 | $18.74$ | $<0.001$
2 | $160.80$ | $<0.001$
3 | $11.95$ | $<0.001$
4 | $41.43$ | $<0.001$
|  |
Experiment S (Striped surround)
Spatial | $F_{6,378}$ | p
Condition
1 | $84.26$ | $<0.001$
2 | $27.19$ | $<0.001$
3 | $20.56$ | $<0.001$
4 | $10.88$ | $<0.001$

##  4 Discussion

We have divided the discussion into two parts, according to the two
experiments (Experiment U and Experiment S).

###  4.1 Uniform surround (Experiment U)

Our results from uniform surround stimuli (see Figure 4) show that static
uniform surrounds induced color contrast in all experimental conditions, in
line with previous works [26, 27, 17, 30, 23]. Concretely, Gordon and Shapley
[17] used uniform surrounds to study how the luminance and brightness of the
test region affect color induction, observing color contrast in all
conditions. In accordance to Kirshcmann’s Third Law [25], they concluded that
brightness, but not luminance, is crucial in the effect on color induction.
They found that color contrast is maximal when the stimulus is equibrightness
(but not equiluminant) and, as brightness contrast is increased, the color
contrast is reduced. In their work, Kaneko and Murakami [23] also used
equiluminant stimuli (using the heterochromatic flicker photometry technique)
and they also found color contrast under all conditions. They used flashed
stimuli and measured the color contrast induced by color surrounds presented
in different temporal conditions. They concluded that the shorter the flash
duration, the stronger the color contrast. By contrast, we observed a clear
peak of color contrast when stimuli were flashed during $40\,ms$. Despite our
work shares several features with Kaneko and Murakami [23] such as temporal
conditions, temporal configuration, equiluminant stimuli and methodology,
there are some differences between the current study and their work. They
measured color induction on a central disk and we measured it on a ring
similar to the one used by Monnier and Shevell [26] and Otazu et al. [30].
Thus, the visual angle of the evaluated feature was different (a disk of
$1\,deg$ and a ring width of $15.5\,min$ of visual angle). In addition, they
introduced a thin black ring around the central disk, i.e., a border of lower
luminance between the central disk and its color surround, which could lead to
different results [43]. Carefully analyzing their raw data, we can observe
that in some spatial conditions there is not a clear peak of color contrast at
the shortest flash ($10\,ms$), but there is around $20$ and $40\,ms$.
Moreover, they showed two subjects (subjects MS and YY) who seemed to obtain
similar results to ours: they had a peak of induction during the short flashes
(around $20\,ms$), but not at the shortest ($10\,ms$). Thus, all these reasons
could explain the dissimilarity between our and their results.

It is assumed that color induction (both color contrast and assimilation) is
the result of neural mechanisms in V1 [45, 44, 41, 10] and stimuli on $l$ (red
and green surrounds) and $s$ (purple and lime surrounds) axes are nearly
independently processed at the first stages of the Human Visual System (HVS),
i.e., in the retina, LGN and V1 [33]. Since stimuli on the $l$ axis are
processed by the parvocellular pathway and stimuli on $s$ axis are processed
by the koniocellular pathway [28], we expected to obtain different results
using stimuli that independently activate these different visual pathways. In
particular, we expected to observe different temporal behaviors for each
pathway because parvocellular and koniocellular pathways have different
processing speeds [11, 7]. From the obtained results (see Figure 4), we can
see that the color contrast when the inducer is purple (Spatial Condition 2 of
Table 1) is completely different to that with other color surrounds. This
spatial condition induced color contrast, i.e., the achromatic test ring was
perceived as lime, when flashed longer than $20\,ms$. Moreover, it is the only
color of the surround with no induction maximum, except at infinite flash
duration (i.e., static temporal condition). This stimulus activates the S-OFF
channel of koniocellular pathway, which directly projects to layer 4A of V1
[12, 9]. By contrast, the S-ON channel of the koniocellular, and both
parvocellular and magnocellular pathways project to layers 2/3, 4C$\beta$, and
4C$\alpha$ of V1, respectively, and all of them converge into layer 2/3 of V1
[33]. The different processing of S-ON and S-OFF channels of koniocellular
pathway could explain the dissimilar psychophysical results (see Figure 4) on
color contrast when inducers had lime (S-ON on test ring) or purple (S-OFF on
test ring) chromaticities. In addition, since all channels of parvocellular
pathway (L-ON, L-OFF, M-ON, and M-OFF) are processed in the same layers of V1
(first in layer 4C$\beta$, 4B and finally in 2/3), and S-ON channel of
koniocellular pathway is mainly processed in layer 2/3, it could explain this
similarity in color induction when inducers had red, green or lime
chromaticities (Spatial Conditions 1, 2 and 4).

###  4.2 Striped surround (Experiment S)

Our results from striped surround stimuli (see Figure 5) show that, similarly
to other authors [14, 35], striped surrounds can induce color contrast, but we
have not observed it under all experimental conditions. Only one out of 28
experimental conditions induced color assimilation, namely when the stimulus
was static, its 1st inducer was red and the 2nd one was green. In contrast,
Monnier and Shevell [26, 27] and Otazu et al. [30] observed that striped
surrounds induce color assimilation in all spatial conditions and never found
color contrast. Although our work has very similar features to these works,
such as spatial configuration and chromaticities (see Section 2.2), they used
non-equiluminant stimuli and we used equiluminant stimuli. Thus, our
hypothesis is that this luminance difference explains the difference between
our results and the ones obtained by these authors. Monnier and Shevell [26,
27], and Otazu et al. [30] introduced a luminance difference between the test
ring and its surround, and they found chromatic assimilation for striped
surrounds. De Weert and Spillmann [42] found no color induction in
equiluminant stimuli, but found color induction when the 1st inducer had lower
luminance than the central region. Extending this work, Cao and Shevell [10]
showed that assimilation in the $l$ axis of MacLeod and Boynton color space
[5] was found when the inducer luminance was lower than the central region
luminance, but not when it was higher, observation that was also reported by
De Weert and Spillmann [42]. In the $s$ axis, they showed that color
assimilation does not depend on the inducing luminance (i.e., induction was
observed when the inducing luminance was either lower or higher than the
central region luminance), but depends on the spatial configuration of the
inducers (i.e., on both spatial frequency and inducer’s spatial separation).
Thus, our hypothesis is supported by De Weert and Spillmann [42] and Cao and
Shevell [10] results. Considering that chromatic assimilation mainly appears
when stimulus is not equiluminant, i.e., only when magnocellular pathway is
activated, it could suggest that magnocellular pathway could act as a switch-
like signaling system activating or deactivating assimilation in both
parvocellular and koniocellular pathways in layer 2/3 (where all the pathways
converge). In recent neurophysiological and psychophysical studies [43, 29],
the authors concluded that brightness and color interact in V1. In particular,
their work supported the hypothesis that the color appearance depends on
brightness contrast [25, 17, 16, 3] because there is a mutual-suppression
(i.e., color assimilation) between color-responsive cells and luminance-
responsive cells. Furthermore, they proposed that these interactions are
driven by double-opponent cells, which respond to both luminance and color
differences (Color-Lum neurons) [19, 20]. The non-opponent neurons, or Lum
neurons, are inactive when the stimulus is equiluminant, and single-opponent
neurons, or Color neurons, respond to large areas of color and do not respond
to luminance differences. These neurophysiological observations, also supports
our hypothesis that luminance difference between the test ring and its
surround, which activates Lum neurons, could be a key factor to induce color
assimilation. In addition, psychophysical studies by Fach and Sharpe [14], and
Smith, Jin and Pokorny [35] concluded that spatial frequency is another
crucial factor to induce color assimilation in equiluminant striped stimuli.
In particular, they observed that very thin stripes ($>9\,cpd$) induce color
assimilation and thick stripes ($<0.7\,cpd$) induce color contrast, with a
transition point from assimilation to contrast around $4\,cpd$ [35].
Considering that the spatial frequency of our stimuli is $1.94\,cpd$, we agree
with them [14, 35, 10]: in that range both color contrast and color
assimilation could be induced.

Comparing the results from flashed and static stimuli, we can observe that in
almost all spatial conditions (except for Spatial Condition 1), the type of
color induction, e.g. assimilation or contrast, did not vary between flashed
and static stimuli, when the flash was longer than $40\,ms$.

In both types of surround configurations, e.g., uniform and striped surrounds
(see Figures 4 and 5), we can see that at the shortest flash duration
($10\,ms$) subjects only perceived color induction when the surrounding colors
were on the $l$ axis of MacLeod and Boynton color space, but did not perceive
any induction when the surround was on the $s$ axis (purple or lime colors).
In fact, when subjects finished the experiment they pointed out that under
these experimental conditions (purple and lime surround colors flashed for
$10\,ms$) they were not able to see the test frame and, therefore, left an
achromatic comparison ring. This consideration goes in line with the
physiological observation that the koniocellular pathway is slower than the
parvocellular pathway [11, 7], i.e., the purple and lime colors are processed
more slowly than the red and green colors.

##  5 Conclusions

Taking into account that we only observed color contrast (except for red-green
inducer in static striped stimuli) for the two different types of surround
(uniform and striped) and that temporal behavior of color induction depends on
visual pathways (see Figures 4 and 5), we can conclude:

  * •

The strongest color contrast is induced by a uniform surround stimulus flashed
for $40\,ms$.

  * •

Purple inducer (which induces a lime chromaticity and, thus, activates the
S-OFF channel in layer 4A) induces a temporal response that is completely
different to the temporal response induced by other inducers such as red,
green and lime (which activates other channels in layer 2/3).

  * •

Striped equiluminant stimuli do not induce color assimilation (except for red-
green inducers).

  * •

The test frame cannot be perceived during flashes shorter than $20\,ms$ when
the colors of the surrounding are on the s axis of the MacLeod and Boynton
color space.

  * •

Our initial hypothesis, i.e., flashed uniform and striped surrounds would
induce opposite colors but with a similar temporal behavior, should be
rejected.

Considering previous studies, we can also conclude that luminance could be a
key factor for color assimilation. In particular, assimilation only appears in
non-equiluminant stimuli, or in equiluminant striped stimuli with a very high
spatial frequency. This could suggest that color contrast and color
assimilation effects are the result of different mechanisms or, at least, the
result of the same mechanism which needs an interaction with luminance to
induce color assimilation.

##  6 Future Work

A more detailed study of the specific contribution of luminance to color
assimilation is needed. There exist some works in this direction [42, 10], but
they did not perform a systematic study of color induction depending on
luminance differences (they used few arbitrary luminance values, and Cao and
Shevell [10] did not use equiluminant stimuli). In particular, it is
interesting to study how the $\Delta C$ values (see Section 3.1) change
depending on test ring, 1st and 2nd inducers’ luminance differences. It would
help to study the precise dependence of chromatic induction on both
equiluminant and non-equiluminant stimuli. It would help to unify the current
work and several previous works [26, 27, 30] under a common framework.

##  7 Funding Information

This work is partially supported by: Spanish Ministry of Economy, Industry and
Competitiveness (DPI2017-89867-C2-1-R); Agencia de Gestio d’Ajuts
Universitaris i de Recerca (AGAUR) (2017-SGR-649); CERCA Programme /
Generalitat de Catalunya.

##  8 Acknowledgements

The authors would like to thank Javier Retana for his useful comments on
statistical analysis procedures, Olivier Penacchio for his comments on the
work and C. Alejandro Parraga for his comments on the laboratory setup. The
authors also want to thank all the subjects for their valuable time.

## References

  * [1] S. Anstis, B. Rogers, and J. Henry.  Interactions between simultaneous contrast and coloured afterimages.  Vision Research, 18:899–911, 1978.
  * [2] W. Von Bezold.  The theory of color and its relation to art and art-industry.  L. Prang and Company, Boston, 1876.
  * [3] D.L. Bimler, G.V. Paramei, and C.A. Izmailov.  Hue and saturation shifts from spatially induced blackness.  Journal of the Optical Society of America A, 26(1):163–172, 2009\.
  * [4] R.M. Boynton.  Implications of the minimally distinct border.  Journal of the Optical Society of America, 63(9):1037–1043, 1973\.
  * [5] R.M. Boynton.  A system of photometry and colorimetry based on cone excitations.  Color Research and Applications, 11:244–252, 1986.
  * [6] R.M. Boynton and P.K. Kaiser.  Vision: The additivity law made to work for heterochromatic photometry with bipartite fields.  Science, 161(3839):366–368, 1968.
  * [7] F. Briggs and W.M. Usrey.  Parallel processing in the corticogeniculate pathway of the macaque monkey.  Neuron, 62:135–146, 2009.
  * [8] M.H. Brill.  Minimally distinct border.  In Ronnier Luo, editor, Encyclopedia of Color Science and Technology, pages 1–3. Springer Berlin Heidelberg, Berlin, Heidelberg, 2014\.
  * [9] E.M. Callaway.  Cells types and local circuits in primary visual cortex of the macaque monkey.  In J.S. Werner and L.M. Chalupa, editors, The New Visual Neurosciences, pages 353–365. The MIT Press, Cambridge, Massachusetts & London, England, 2014.
  * [10] D. Cao and S.K. Shevell.  Chromatic assimilation: spread light or neural mechanism?  Vision Research, 45:1031–1045, 2005.
  * [11] V. Casagrande, F. Yazar, K. Jones, and Y. Ding.  The morphology of the koniocellular axon pathway in the macaque monkey.  Cerebral Cortex, 17(10):2334–2345, 2007.
  * [12] S. Chatterjee and E.M. Callaway.  Parallel colour-opponent pathways to primary visual cortex.  Nature, 426:668–671, 2003.
  * [13] M.E. Chevreul.  De la loi du contraste simultane des couleurs.  Chez Pitois-Levrault et Ce, Paris, 1839.
  * [14] C. Fach and L.T. Sharpe.  Assimilative hue shifts in color gratings depend on bar width.  Perception & Psychophysics, 40(6):412–418, 1986.
  * [15] D. Farnsworth.  The Farnsworth dichotomous test for color blindness: panel D-15, 1947.
  * [16] F. Faul, V. Ekroll, and G. Wendt.  Color appearance: the limited role of chromatic surround variance in the ’gamut expansion effect’.  Journal of Vision, 8(30):1–20, 2008.
  * [17] J. Gordon and R. Shapley.  Brightness contrast inhibits color induction: evidence for a new kind of color theory.  Spatial Vision, 19:133–146, 2006.
  * [18] S. Ishihara.  Tests for colour-blindness.  Kanehara Shippan Co., Ltd., Tokyo, Japan, 1972.
  * [19] E.N. Johnson, M.J. Hawken, and R. Shapley.  The spatial transformation of color in the primary visual cortex of the macaque monkey.  Nature Neuroscience, 4(4):409–416, 2001.
  * [20] E.N. Johnson, M.J. Hawken, and R. Shapley.  The orientation selectivity of color-responsive neurons in macaque v1.  The Journal of Neuroscience, 28(32):8096–8106, 2008.
  * [21] P. Kaiser.  Minimally distinct border as a preferred psychophysical criterion in visual heterochromatic photometry.  Journal of the Optical Society of America A, 61(7):966–971, 1971\.
  * [22] P. Kaiser, B. Lee, P. Martin, and A. Valberg.  The physiological basis of the minimally distinct border demonstrated in the ganglion cells of the macaque retina.  Journal of Physiology, 422:153–183, 1990.
  * [23] S. Kaneko and I. Murakami.  Flashed stimulation produces strong simultaneous brightness and color contrast.  Journal of Vision, 12(1):1–18, 2012.
  * [24] D.H. Kelly and E. Martinez-Uriegas.  Measurements of chromatic and achromatic afterimages.  Journal of Optical Society of America A, 10(1):29–37, 1993.
  * [25] A. Kirschmann.  Ueber die quantitativen verhältnisse des simultanen helligkeits- und farben-contrastes.  Philosophische Studien, 6:417–491, 1891.
  * [26] P. Monnier and S. Shevell.  Large shifts in color appearance from patterned chromatic backgrounds.  Nature Neuroscience, 6:801–802, 2003.
  * [27] P. Monnier and S. Shevell.  Chromatic induction from s-cone patterns.  Vision Research, 44:849–856, 2004.
  * [28] J.J. Nassi and E.M. Callaway.  Parallel processing strategies of the primate visual system.  Nature Reviews Neuroscience, 10:360–372, 2009.
  * [29] V. Nunez, R.M. Shapley, and J. Gordon.  Cortical double-opponent cells in color perception: perceptual scaling and chromatic visual evoked potentials.  i-Perception, 9(1):1–16, 2018.
  * [30] X. Otazu, C.A. Parraga, and M. Vanrell.  Toward a unified chromatic induction model.  Journal of Vision, 10(12):1–24, 2010.
  * [31] X. Otazu, M. Vanrell, and C.A. Parraga.  Multiresolution wavelet framework models brightness induction effects.  Vision Research, 48(5):733–751, 2008.
  * [32] O. Penacchio, X. Otazu, and L. Dempere-Marco.  A neurodynamical model of brightness indcution in v1.  PLOS One, 8(5):1–14, 2013.
  * [33] L.C. Sincich and J.C. Horton.  The circuitry of v1 and v2: integration of color, form, and motion.  Annual Review of Neuroscience, 28:303–326, 2005.
  * [34] B. Singer and M. D’Zmura.  Color contrast induction.  Vision Research, 34(23):3111–3126, 1994.
  * [35] V.C. Smith, P.Q. Jin, and J. Pokorny.  The role of spatial frequency in color induction.  Vision Research, 41(8):1007–1021, 2001.
  * [36] V.C. Smith and J. Pokorny.  Spectral sensitivity of the foveal cone photopigments between 400 and 500 nm.  Vision Research, 15:161–171, 1975.
  * [37] R.L. De Valois and K.K. De Valois.  Sensitivity to Color Variations, chapter 7, pages 212–239.  Oxford Psychology Series. Oxford University Press, Clarendon Press, New York, Oxford, 1988.
  * [38] R.L. De Valois and K.K. De Valois.  Spatial Vision.  Oxford Psychology Series. Oxford University Press, Clarendon Press, New York, Oxford, 1988.
  * [39] R.L. De Valois, M.A. Webster, and K.K. De Valois.  Temporal properties of the brightness and color induction.  Vision Research, 26(6):887–897, 1986.
  * [40] G. Wagner and R. Boynton.  Comparison of four methods of heterochromatic photometry.  Journal of Optical Society of America A, 62(12):1508–1515, 1972\.
  * [41] C.M. De Weert and N.A. Kruysbergen.  Assimilation: central and periferal effects.  Perception, 26(10):1217–1224, 1997.
  * [42] C.M. De Weert and L. Spillmann.  Assimilation: Asymmetry between brightness and darkness?  Vision Research, 35(10):1413–1419, 1995.
  * [43] D. Xing, A. Ouni, S. Chen, H. Sahmoud, J. Gordon, and R. Shapley.  Brightness-color interactions in human early visual cortex.  The Journal of Neuroscience, 35(5):2226–2232, 2015.
  * [44] Q. Zaidi.  Color and brightness induction: from mach bands to three-dimensional configurations.  In K. R. Gegenfurtner and L. T. Sharpe, editors, Color Vision: From genes to perception, pages 317–344. Cambridge University Press, Cambridge, UK, 1999.
  * [45] Q. Zaidi, B. Yoshimi, N. Flanigan, and A. Casanova.  Lateral interactions within color mechanism in simultaneous induced contrast.  Vision Research, 32(9):1695–1707, 1992.