# Quark anomalous magnetic moment in an extreme magnetic background from pQCD

Eduardo S. Fraga<EMAIL_ADDRESS>Instituto de Física, Universidade Federal
do Rio de Janeiro, Caixa Postal 68528, 21941-972, Rio de Janeiro, RJ, Brazil
Leticia F. Palhares<EMAIL_ADDRESS>Universidade do Estado do Rio de
Janeiro, Instituto de Física-Departamento de Física Teórica-Rua São Francisco
Xavier 524, 20550-013 Maracanã, Rio de Janeiro, Brasil  Cristian Villavicencio
<EMAIL_ADDRESS>Centro de Ciencias Exactas and Departamento de
Ciencias Básicas, Facultad de Ciencias, Universidad del Bío-Bío, Casilla 447,
Chillán, Chile

###### Abstract

We compute the one-loop QCD correction to the photon-quark-antiquark vertex in
an extremely strong magnetic background, i.e., one in which $\sqrt{eB}$ is
much larger than all other mass scales. We resort to the lowest-Landau level
approximation, and consider on-shell fermions. We find that the the total
magnetic moment is such that that the anomalous magnetic moment contributes to
the electric part $\sim iE_{3}\sigma_{3}$, the other contributions to
$\bm{B}\cdot\bm{\sigma}$ and $i\bm{E}\cdot\bm{\sigma}_{\perp}$ being
completely suppressed. We show results for the anomalous magnetic moment of
quarks up and down, scaled by their values in the vacuum, as a function of the
external magnetic field, considering dressed gluons and a renormalization
scale $\mu_{\mathrm{\tiny\overline{MS}}}=\sqrt{|eB|}$.

##  I Introduction

The original one-loop calculation of the anomalous magnetic moment [1],
together with the spectacular agreement with experimental data that followed
refined higher-order perturbative calculations (see Ref. [2] for a review),
represents a major landmark of theoretical physics.

The one-loop correction to the anomalous magnetic moment (AMM) of a fermion
that moves in the presence of an external magnetic background will, however,
deviate from the original prediction by Schwinger of $\alpha/2\pi$ [1, 3],
where $\alpha\approx 1/137$ is the QED fine-structure constant. This was
already known by Schwinger, since what is currently known as the Schwinger
phase appears in each fermion propagator [4]. The functional dependence of
this correction on the magnitude of the external magnetic field has first been
addressed almost half a century ago [5, 6].

Motivated by the ultra-strong magnetic fields that can be achieved in current
experiments at the Large Hadron Collider (LHC), we compute the one-loop
perturbative QCD correction to the photon-quark-antiquark vertex in an
extremely strong magnetic background, i.e., one in which $\sqrt{eB}$ is much
larger than all other mass scales. Here $e$ is the fundamental electric charge
and $B$ is the magnetic field strength. For such high fields, we can resort to
the lowest-Landau level (LLL) approximation, and consider on-shell fermions.

The AMM of quarks in the presence of an external magnetic field can be related
to the phenomenon of magnetic catalysis [7, 8, 9] and, in the presence of a
hot plasma, might play a role in the so-called chiral magnetic effect [10].
The ultra-high magnetic fields that can be produced in high-energy heavy ion
collisions [11, 12, 13] bring the question of whether the AMM might prove to
be experimentally relevant in this physical setting. Therefore, effects from
the AMM have been under investigation in quark and hadronic matter for over
two decades [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32]

Usually, estimates of magnetic corrections to the fermion AMM are based on the
Schwinger ansatz [4], which relates the AMM with the tensor spinor structure
of the selfenergy. On the other hand, not much has been done in obtaining the
AMM from radiative corrections to the fermion-photon coupling. The AMM has
been calculated from the fermion-photon vertex, including magnetic field
corrections, using the Schwinger proper-time method for low magnetic fields
[33].

In the paper at hand, conversely, we explore the magnetic effects on the quark
AMM in the presence of a very intense magnetic background. We find that the
total magnetic moment is such that that the AMM in the presence of a very
strong magnetic background contributes to the electric part $\sim
iE_{3}\sigma_{3}$, which should not be confused with an anomalous electric
dipole moment [34]. The other AMM contributions to $\bm{B}\cdot\bm{\sigma}$
and $i\bm{E}\cdot\bm{\sigma}_{\perp}$ are completely suppressed.

We also show results for the anomalous magnetic moment of quarks up and down,
scaled by their values in the vacuum, as a function of the external magnetic
field, considering dressed gluons and a renormalization scale
$\mu_{\mathrm{\tiny\overline{MS}}}=\sqrt{|eB|}$. We consider two two
scenarios: one with a constant gluon mass, $m_{g}=0.3$ GeV, and one with the
gluon mass extracted from the one-loop correction to its polarization tensor
in presence of a large external magnetic field.

The paper is organized as follows. In Section II we present the general
formalism. First, we outline the computation of the anomalous correction in
QED. Then, we compute the AMM of QCD in a strong magnetic background, with
special emphasis on the QCD correction to the vertex triangular diagram that
contributes to the structure function $F_{2}$. In Section III we present and
discuss our results for the QCD correction to the AMM as a function of the
magnetic field in a few phenomenological scenarios. Finally, Section IV
contains our summary and outlook.

##  II Anomalous magnetic moment: formalism

###  II.1 Outline of AMM in one-loop QED

In order to set the notation and define the main physical quantities, we start
with a brief discussion of well-known textbook results for the anomalous
magnetic moment in QED (see, e.g., Ref. [34]).

The Dirac equation

$[i\not{D}-m]\psi=0$ (1)

provides, at tree level, the fermion-antifermion-photon vertex
$e_{q}\overline{u}(p^{\prime})\gamma^{\mu}u(p)$, where $m$ is the spinor mass,
$e_{q}$ its electric charge ($e_{q}=-e$ for the electron), and
$D_{\mu}=\partial_{\mu}+ie_{q}A_{\mu}$ the covariant derivative. Here, we use
the Weyl representation for the Dirac gamma matrices in Minkowski space. The
equation of motion naturally exhibits the AMM when written in quadratic form,
so that charged spinors in an electromagnetic field $F_{\mu\nu}$ satisfy

$\left[D^{2}+m^{2}-\frac{e_{q}}{2}F_{\mu\nu}\sigma^{\mu\nu}\right]\psi=0\,,$
(2)

where $[D_{\mu},D_{\nu}]=-ie_{q}F_{\mu\nu}$, and
$\sigma^{\mu\nu}=\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$. One-loop corrections
bring about other structures, with terms proportional to $(p^{\prime}+p)$ and
$(p^{\prime}-p)$. In this case, one can resort to the Gordon identity

$\overline{u}(p^{\prime})(p^{\prime\mu}+p^{\mu})u(p)=\overline{u}(p^{\prime})\left[2m\,\gamma^{\mu}-i\sigma^{\mu\nu}(p^{\prime}_{\nu}-p_{\nu})\right]u(p)\,,$
(3)

where the fermions are on-shell. The photon momentum $q=p^{\prime}-p$ can be
recognized in the previous equation, so the contribution to the fermion-photon
vertex $\gamma_{\mu}$ is interpreted as a correction to the fermion charge,
while the term $\sigma_{\mu\nu}q^{\nu}$ is identified with the AMM.

At any loop order, the vertex is constrained by Lorentz symmetry, the Dirac
equation for spinors, $\not{k}_{i}u(k_{i})=m\,u(k_{i})$, the Ward identity,
and the Gordon identity. So, one can write it in the following standard form:

${\cal
M}^{\mu}=e_{q}\overline{u}(p^{\prime})\left[F_{1}\left(\frac{q^{2}}{m^{2}}\right)\gamma^{\mu}+F_{2}\left(\frac{q^{2}}{m^{2}}\right)\frac{i\sigma^{\mu\nu}}{2m}q_{\nu}\right]u(p)\,.$
(4)

On one hand, one can see that the function $F_{1}$ encodes charge
renormalization, which brings a scale dependence to $e_{q}$. On the other
hand, $F_{2}$ has the structure of a magnetic moment, and indeed

$g=2\left[1+F_{2}\left(\frac{q^{2}}{m^{2}}\right)\right]\,.$ (5)

To one loop, only the triangular diagram of Fig. 1 contributes to $F_{2}$ and,
consequently, to $(g-2)$. Other diagrams correspond to particle-reducible
corrections to the external legs and contribute only to $F_{1}$. The one-loop
contribution to $F_{2}$ is thus given by

$\begin{split}i\delta{\cal
M}^{\mu}=(ie_{q})^{3}\int_{k}\,D_{\mu\nu}(k)\,\overline{u}(p^{\prime})\gamma^{\nu}S(k+p^{\prime})\\\
\gamma^{\mu}S(k+p)\gamma^{\rho}\,u(p)\,,\end{split}$ (6)

where the free fermion propagator in momentum space is
$S(k)=i/(\not{k}-m+i\epsilon)$ and the gluon free propagator in momentum space
in the Feynman gauge is $D_{\mu\nu}(k)=-ig_{\mu\nu}D(k)$, with
$D(k)=i/(k^{2}+i\epsilon)$. Here, we use the notation
$\int_{k}\equiv\int\frac{d^{4}k}{(2\pi)^{4}}$.

Using the standard Feynman trick of writing

$\frac{1}{ABC}=\int_{xyz}\frac{2}{\left[xA+yB+zC\right]^{3}}\,,$ (7)

where we use the following notation for the integrals in Feynman parameters,

$\int_{xyz}=\int_{0}^{1}dx~{}dy~{}dz~{}\delta(x+y+z-1),$ (8)

and performing usual quantum field theory manipulations, yields

$F_{2}(q^{2})=\frac{\alpha}{\pi}m^{2}\int_{xyz}\frac{z(1-z)}{(1-z)^{2}m^{2}-xyq^{2}}\,.$
(9)

For an on-shell photon ($q^{2}=0$), this leads to the well-known result
$F_{2}(0)=\alpha/2\pi$, so that $g=2+\alpha/2\pi$.

Figure 1: Vertex triangular diagram from QED that contributes to $F_{2}$.

###  II.2 One-loop QCD correction to the photon-quark-antiquark vertex in a
strong magnetic background

Let us now consider the triangle diagram depicted in Fig. 2, encoding the
first QCD correction to the quark-photon vertex. In the vacuum, its result is
analogous to the one discussed in the last section, except for the inclusion
of the color structure of the $SU(N_{c})$ group and the swapping of the
electromagnetic coupling for the strong coupling $\alpha_{s}=g_{s}^{2}/4\pi$.
This gives rise to a AMM correction of the following form:

$a_{q}^{\mathrm{vac}}\equiv F_{2}^{\text{\tiny
vac}}(0)=\frac{\alpha_{s}}{2\pi}\frac{N_{c}^{2}-1}{2N_{c}}\,.$ (10)

This kind of triangle diagrams describing corrections to some fermion-
antifermion-boson interaction has been calculated in the presence of external
magnetic field with different methods [35, 33, 36, 37, 38, 39].

In the presence of a strong magnetic field $\bm{B}=B\hat{\bm{z}}$, the quark
fields are quantized in Landau levels. Here we adopt the Landau gauge
representation [40, 41], with $A^{\mu}=(0,0,Bx,0)$. So, the matrix elements
for the quark field operator have the form

$\langle
0|\psi(\xi)|p_{\parallel},p_{y},n\rangle=f_{n}(p_{y},\xi_{\perp})u_{n}(p_{\parallel})\,e^{-ip\cdot\xi_{\parallel}}\,,$
(11)

where $\xi$ is a spacetime coordinate coordinate four-vector and
$p^{\mu}=(p^{0},p^{1},p^{2},p^{3})$ is the four-momentum. The parallel and
perpendicular vectors are defined with respect to the magnetic field
direction, $\bm{B}=B\hat{\bm{z}}$: $v_{\parallel}=(v^{0},0,0,v^{3})$ and
$v_{\perp}=(0,v^{1},v^{2},0)$. Also, $n$ denotes the index of Landau levels,
$u$ is the spinor, and $f_{n}$ is a diagonal matrix in spinorial space related
to the solution of the Dirac equation in the presence of a magnetic background
[41]. For the external photon leg, the associated photon field matrix element
reads:

$\langle 0|{\cal
A}_{\mu}(\zeta)|q\rangle=\varepsilon_{\mu}(q)\,e^{-iq\cdot\zeta}\,,$ (12)

where $\zeta_{\mu}$ is a spacetime coordinate, $q_{\mu}$ is the four-momentum
and $\varepsilon_{\mu}$ is the polarization vector.

We use the Minkowski metric with an explicit separation of parallel and
perpendicular components defined as
$g_{\mu\nu}=g_{\mu\nu}^{\parallel}+g_{\mu\nu}^{\perp}$, where
$g_{\mu\nu}^{\parallel}\equiv{\rm diag}(1,0,0,-1)$ and
$g_{\mu\nu}^{\perp}\equiv{\rm diag}(0,-1,-1,0)$. This means that for any
vector $v_{\perp}^{2}=-\bm{v}_{\perp}^{2}$. Also, we use a general notation
for four-vector coordinates
$(v^{0},v^{1},v^{2},v^{3})=(v_{t},v_{x},v_{y},v_{z})$. The temporal components
of the relevant external momenta are given by

$\displaystyle p_{0}$ $\displaystyle=\sqrt{p_{z}^{2}+2n\ell_{q}^{-2}+m^{2}},$
(13) $\displaystyle p_{0}^{\prime}$ $\displaystyle=\sqrt{p_{z}^{\prime
2}+2n^{\prime}\ell_{q}^{-2}+m^{2}},$ (14) $\displaystyle q_{0}$
$\displaystyle=|\bm{q}|,$ (15)

where $\ell_{q}=|e_{q}B|^{-1/2}$ is the magnetic length and $m$ is the quark
mass.

The one-loop vertex transfer matrix associated with the QCD correction to the
quark-photon vertex is then

$\begin{split}{\cal
T}=-e_{q}g_{s}^{2}\sum_{l,l^{\prime}}\int_{\xi\xi^{\prime}\zeta}e^{i\left(p^{\prime}\cdot\xi_{\parallel}^{\prime}-p\cdot\xi_{\parallel}-q\cdot\zeta\right)}\,\bar{D}_{\nu\rho}(\xi-\xi^{\prime})\\\
\bar{u}_{n^{\prime}}(p_{\parallel}^{\prime})f^{\dagger}_{n^{\prime}}(p_{y}^{\prime},\xi_{\perp}^{\prime})\,\gamma^{\nu}t^{a}\,\bar{G}_{l}(\zeta,\xi^{\prime})\,\\\
\not{\varepsilon}\,\bar{G}_{l^{\prime}}(\xi,\zeta)\,\gamma^{\rho}t^{a}\,f_{n}(p_{y},\xi_{\perp})u_{n}(p_{\parallel}),\end{split}$
(16)

where $g_{s}$ is the gauge coupling constant, $t^{a}$ are Gell-mann matrices,
and $\bar{G}_{n}$ is the fermion Green’s function in coordinate space in the
presence of the external magnetic field for a given Landau level $n$ [8].
Here, we adopt the Feynman gauge for the gluon and its propagator
$\bar{D}_{\nu\rho}(x)$ in coordinate space can be expressed in terms of the
Fourier transformation as

$\bar{D}_{\nu\rho}(\xi-\xi^{\prime})=g_{\nu\rho}\int_{k}e^{-ik\cdot(\xi-\xi^{\prime})}D(k)\,.$
(17)

The quark propagator can be expressed in terms of the Fourier transformation
of the parallel coordinates

$\bar{G}_{n}(\xi,\zeta)=\int_{k_{\parallel}}e^{-ik\cdot(\xi-\zeta)_{\parallel}}G_{n}(k_{\parallel};\xi_{\perp},\zeta_{\perp}),$
(18)

so that, after integrating over parallel variables, we obtain parallel
momentum conservation only. Therefore, we define the resulting one-loop QCD
correction to the vertex matrix as

${\cal
T}=(2\pi)^{2}\delta^{(2)}(q_{\parallel}-p^{\prime}_{\parallel}+p_{\parallel})\,\varepsilon_{\mu}\,\delta{\cal
M}^{\mu}\,,$ (19)

with

$\begin{split}\delta{\cal
M}^{\mu}=-e_{q}g_{s}^{2}\frac{N_{c}^{2}-1}{2N_{c}}\int_{\xi_{\perp}\xi^{\prime}_{\perp}\zeta_{\perp}k}e^{-ik\cdot(\xi-\xi^{\prime})_{\perp}}\,D(k)\\\
\bar{u}^{\prime}_{n}(p_{\parallel}^{\prime})f_{n}^{\dagger}(p_{y}^{\prime},\xi^{\prime})\,\gamma^{\nu}\,G_{n^{\prime}}(k_{\parallel}+p^{\prime}_{\parallel};\zeta_{\perp},\xi^{\prime}_{\perp})\\\
\gamma^{\mu}\,G_{n}(k_{\parallel}+p_{\parallel};\xi_{\perp},\zeta_{\perp})\,\gamma_{\nu}\,f_{n}(p_{y},\xi)u_{n}(p_{\parallel}).\end{split}$
(20)

Figure 2: One-loop QCD correction to the vertex triangular diagram that
contributes to $F_{2}$.

If we restrict our analysis to extreme magnetic fields, we can keep only the
lowest-Landau level (LLL) and set $n=n^{\prime}=l=l^{\prime}=0$. Then, our
expressions simplify to

$f_{0}(p_{y},\xi_{\perp})u_{0}(p_{\parallel})=u_{0}(p_{\parallel})\,\frac{e^{ip_{y}\xi_{y}-\frac{1}{2}(\xi_{x}/\ell_{q}+sp_{y}\ell_{q})^{2}}}{(\pi\ell_{q}^{2})^{1/4}}\,,$
(21)
$\begin{split}G_{0}(k_{\parallel};\zeta_{\perp},\xi_{\perp})=\frac{i}{2\pi\ell_{q}^{2}}S(k_{\parallel})\,{\cal
P}_{+}\qquad\\\
e^{\frac{is}{2\ell_{q}^{2}}(\xi_{x}+\zeta_{x})(\xi_{y}-\zeta_{y})}\,e^{\frac{1}{4\ell_{q}^{2}}(\zeta_{\perp}-\xi_{\perp})^{2}}\,,\end{split}$
(22)

where $s={\rm sgn}(e_{q}B)$, and we have defined the projectors ${\cal
P}_{\pm}=\frac{1}{2}\left[1\pm si\gamma_{1}\gamma_{2}\right]$ or,
equivalently, ${\cal P}_{\pm}=\frac{1}{2}[1\pm\Sigma_{3}]$, with
$\Sigma_{3}=\rm{diag}(1,-1,1,-1)$ being the third component of the spin
projection (cf. also Ref. [8]). The Schwinger phase corresponds to the first
exponential in the second line of Eq. (22).

The LLL matrix can be separated in parallel and perpendicular integrals in the
following form:

$\displaystyle\delta{\cal M}^{\mu}_{\text{\tiny LLL}}$
$\displaystyle=\int_{k_{\perp}}{\cal F}(p_{y},p_{y}^{\prime},k_{\perp})\,{\cal
G}^{\mu}(p_{\parallel},p^{\prime}_{\parallel},k_{\perp}^{2})\,,$ (23)

with

$\displaystyle{\cal F}$
$\displaystyle=\frac{e_{q}g_{s}^{2}}{3\ell_{q}^{5}\pi^{5/2}}\int_{\xi_{\perp}\xi^{\prime}_{\perp}\zeta_{\perp}}\exp\Bigg{\\{}-ik\cdot(\xi-\xi^{\prime})_{\perp}+ip_{y}\xi_{y}-ip^{\prime}_{y}\xi^{\prime}_{y}+\frac{(\zeta-\xi)_{\perp}^{2}}{4\ell_{q}^{2}}+\frac{(\xi^{\prime}-\zeta)_{\perp}^{2}}{4\ell_{q}^{2}}$
$\displaystyle\qquad-\frac{1}{2}(\xi_{x}/\ell_{q}+sp_{y}\ell_{q})^{2}-\frac{1}{2}(\xi_{x}^{\prime}/\ell_{q}+sp_{y}^{\prime}\ell_{q})^{2}-\frac{is}{2\ell_{q}^{2}}(\xi_{x}+\zeta_{x})(\xi_{y}-\zeta_{y})-\frac{is}{2\ell_{q}^{2}}(\zeta_{x}+\xi_{x}^{\prime})(\zeta_{y}-\xi^{\prime}_{y})\Bigg{\\}},$
(24) $\displaystyle{\cal G}^{\mu}$
$\displaystyle=-\int_{k_{\parallel}}\,D(k)\,\bar{u}_{0}(p^{\prime}_{\parallel})\,\gamma^{\nu}S(k_{\parallel}+p^{\prime}_{\parallel})\,{\cal
P}_{+}\gamma^{\mu}S(k_{\parallel}+p_{\parallel})\,{\cal
P}_{+}\gamma_{\nu}\,u_{0}(p_{\parallel}).$ (25)

The integrals in ${\cal F}$ are Gaussian and can be evaluated analytically.
After a long but straightforward calculation, Eq. (23) is reduced to

$\delta{\cal M}^{\mu}_{\text{\tiny
LLL}}=\frac{4}{3}e_{q}g_{s}^{2}\int_{k_{\perp}}e^{\frac{1}{2}k_{\perp}^{2}\ell_{q}^{2}}{\cal
G}^{\mu}(p_{\parallel},p^{\prime}_{\parallel},k_{\perp})\,.$ (26)

Notice that the dependence in $p_{y}$ and $p_{y}^{\prime}$ disappears. Because
of the projectors ${\cal P}_{\pm}$ in Eq. (25), the perpendicular components
of ${\cal G}^{\mu}$ vanish and the form factors are restricted only to the
parallel directions, i.e., $\mu=0$ and $\mu=3$. Moreover, there is only
dependence on parallel momenta, so that all vectors and tensors will be
reduced to their parallel components.

The detailed calculation of ${\cal G}^{\mu}$ is long, but straightforward, and
the full result is not particularly illuminating. In the end, there are two
terms related to $F_{2}$ that are proportional to
$\sigma_{\parallel}^{\mu\nu}q_{\nu}$ and
$\gamma_{1}\gamma_{2}\,\sigma_{\parallel}^{\mu\nu}q_{\nu}$.

###  II.3 Implications to the AMM

The AMM is better described by the polarization projections ${\cal
P}_{+}\,\sigma_{\parallel}^{\mu\nu}q_{\nu}$ and ${\cal
P}_{-}\,\sigma_{\parallel}^{\mu\nu}q_{\nu}$. Defining $\delta{\cal
M}^{\mu}_{\text{\tiny LLL}}=\delta{\cal M}^{\mu}_{+}+\delta{\cal M}^{\mu}_{-}$
and using
$\bar{u}_{0}(p^{\prime}_{\parallel})(\not{p}^{\prime}_{\parallel}-m)=0$,
$(\not{p}_{\parallel}-m)u_{0}(p_{\parallel})=0$ and the Gordon identity, we
obtain the following structure:

$\begin{split}\delta{\cal
M}^{\mu}_{\pm}=e_{q}\bar{u}_{0}(p^{\prime}_{\parallel})\hskip 140.0pt\\\
\Big{[}\delta
F_{1}^{\pm}\gamma^{\mu}_{\parallel}+F_{2}^{\pm}\sigma^{\mu\nu}_{\parallel}\frac{iq_{\nu}}{2m}+F_{3}^{\pm}q_{\parallel}^{\mu}\Big{]}{\cal
P}_{\pm}\,u_{0}(p_{\parallel})\,,\end{split}$ (27)

so that the AMM splits in the two polarization projections.

Here, $\delta F_{1}$ represents a correction to the tree level charge vertex.
All the form factors $F_{i}$ are functions of $q_{\parallel}^{2}/m^{2}$.
$F_{2}$ can be recognized as the AMM contribution, but projected onto parallel
components only. The third term $F_{3}$ should be interpreted as a correction
to the term $\partial\cdot A_{\parallel}$, related to the gauge fixing.

In summary, if we consider the full magnetic moment, we have the following
modification in the presence of a very large magnetic background:

$(1+a)\sigma_{\mu\nu}F^{\mu\nu}\to\sum_{s=\pm}(\sigma_{\mu\nu}F^{\mu\nu}+a^{s}\sigma^{\parallel}_{\mu\nu}F^{\mu\nu}_{\parallel}){\cal
P}_{s}\,,$ (28)

where the parallel-project AMM contributes as $\sim iE_{3}\sigma_{3}$. This
term, however, should not be confused with an anomalous electric dipole moment
[34]. The other AMM contributions to $\bm{B}\cdot\bm{\sigma}$ and
$i\bm{E}\cdot\bm{\sigma}_{\perp}$ are completely suppressed. Indeed, this
becomes clear when the corrections are written in matrix form:

$\sum_{s=\pm}a^{s}\sigma^{\parallel}_{\mu\nu}F^{\mu\nu}_{\parallel}{\cal
P}_{s}=-2iE_{z}\,{\rm diag}(a^{+},-a^{-},-a^{+},a^{-})\,,$ (29)

when compared to the standard anomalous term in vacuum:

$\displaystyle a\,\sigma_{\mu\nu}F^{\mu\nu}$ $\displaystyle=$
$\displaystyle-2a\,\begin{pmatrix}\bm{B}\cdot\bm{\sigma}+i\bm{E}\cdot\bm{\sigma}_{\perp}&0\\\
0&\bm{B}\cdot\bm{\sigma}-i\bm{E}\cdot\bm{\sigma}_{\perp}\end{pmatrix}-$ (30)
$\displaystyle-2iE_{z}\,{\rm diag}(a,-a,-a,a)\,.$

Furthermore, since in general $a^{+}\neq a^{-}$, an asymmetry between the spin
up and spin down corrections develops in the presence of a very large magnetic
background.

Integrating over inner parallel momenta with the use of Feynman parameters,
the AMM component of Eq. (27) reads

$F_{2}^{\pm}(q_{\parallel}^{2})=\int_{xyz}\int_{k_{\perp}}\frac{\frac{4}{3}g_{s}^{2}\,e^{-\frac{1}{2}\bm{k}_{\perp}^{2}\ell_{q}^{2}}\,m^{2}f^{\pm}(z)/4\pi}{\left[m^{2}(1-z)^{2}-q_{\parallel}^{2}xy+(\bm{k}_{\perp}^{2}+m_{g}^{2})z\right]^{2}}$
(31)

where

$\displaystyle f^{+}(z)$ $\displaystyle=\frac{2}{3}(1+6z)$ (32) $\displaystyle
f^{-}(z)$ $\displaystyle=z(4-z)-\frac{2}{3}.$ (33)

To avoid anticipated infrared divergences, we simply add a mass term $m_{g}$
to the gluon propagator, leaving the discussion of possible physical choices
and their consequences to the next subsection.

Since the form factors $F_{i}$ depend on $q_{\parallel}^{2}$, we can choose a
frame where $q_{\parallel}^{2}=0$ if $p_{z}=p_{z}^{\prime}$ or, in the
ultrarelativistic case, for aligned $p_{z}\gg m$ and $p_{z}^{\prime}\gg m$. In
this frame, it is easy to integrate over two of the three Feynman parameters.
On the other hand, the perpendicular momentum integral can be evaluated in
polar coordinates. Changing variables as $|\bm{k}_{\perp}|^{2}=2m^{2}\eta$, we
obtain

$\displaystyle a_{q}^{\pm}$ $\displaystyle\equiv F_{2}^{\pm}(0)$
$\displaystyle=a_{q}^{\mathrm{vac}}\int_{0}^{1}dz\int_{0}^{\infty}d\eta\,\frac{e^{-(m\ell_{q})^{2}\eta}\,f^{\pm}(z)(1-z)}{\left[(1-z)^{2}+(m_{g}/m)^{2}z+2\eta
z\right]^{2}}.$ (34)

One should recall (cf. Eq. (10)) that $a_{q}^{\mathrm{vac}}$ contains the
strong coupling $\alpha_{s}$ which, in principle, runs with the relevant
energy scale.

###  II.4 Effective mass for the gluons

We have introduced a mass term $m_{g}$ for the gluon in the calculation of the
one-loop quark-antiquark-photon vertex in the presence of an extremely large
magnetic background. This scale is necessary to regulate infrared divergences
that appear in this extreme limit111It is important to note that in this
regime, the QED triangle diagram also displays the same infrared divergences
and the full AMM calculation at one-loop would also require an IR
regularization for the photon that we do not discuss here. . We shall analyze
the results in two different scenarios that we motivate in what follows: (i) a
fixed scale and (ii) a magnetically-dressed gluon selfenergy.

First, we consider a fixed infrared scale $m_{g}\sim 0.3$ GeV inspired by
results for correlation functions in the nonperturbative region of QCD (and in
associated pure gauge theories). For Landau and linear covariant gauges,
different Lattice simulations [42, 43, 44, 45], Schwinger-Dyson equations [46,
47, 48] and infrared models [49, 50, 51] seem to indicate the emergence of an
infrared mass scale of this order.

Second, we analyze the consequences of the intense magnetic background on the
gluon polarization tensor to introduce a magnetically-dressed gluon in our
triangle diagram computation. Following [52, 53, 54, 55]. the gluon
polarization tensor in the presence of an intense magnetic field, i.e. within
the LLL approximation, is transverse in the parallel components. Hence:

$\Pi^{\mu\nu}_{g}(k)=\left(g_{\parallel}^{\mu\nu}-\frac{k_{\parallel}^{\mu}k_{\parallel}^{\nu}}{k_{\parallel}^{2}}\right)\Pi_{g}(k_{\parallel}^{2},k_{\perp}^{2})\,,$
(35)

where

$\Pi_{g}(k_{\parallel}^{2},k_{\perp}^{2})=\frac{\alpha_{s}}{\pi}N_{c}\sum_{q}\frac{e^{-\frac{1}{2}\ell_{q}^{2}\bm{k}_{\perp}^{2}}}{\pi^{2}\ell_{q}^{2}}g(m_{q}^{2}/k_{\parallel}^{2})\,.$
(36)

We shall adopt the chiral-limit approximation
$g(m^{2}/k_{\parallel}^{2})\approx g(0)=1$, so that the remaining dependence
is on the perpendicular component of the momentum. Notice that the
polarization factor is exponentially suppressed for $\bm{k}_{\perp}^{2}>|eB|$.
We will consider then $\Pi_{g}$ as a magnetically dressed gluon squared mass.

Figure 3: One-loop QCD correction to the vertex triangular diagram with a
dressed gluon propagator.

##  III QCD correction to the AMM as a function of the magnetic field

Now we proceed to the numerical analysis of the behavior of the one-loop QCD
correction to the AMM in the presence of a very strong magnetic background,
$a_{q}^{\pm}$, compared with its vacuum counterpart, $a_{q}^{\mathrm{vac}}$,
as defined in Eqs. (34) and (10), respectively. In this analysis, we consider
the two scenarios for the effective gluon mass discussed in the previous
section: a constant gluon mass, which we take as $m_{g}=0.3$ GeV, and one
given by the LLL polarization tensor, $m_{g}^{2}=\Pi_{g}$, as defined in Eq.
(36). In both cases we investigate what happens if we set the renormalization
scale to run with the magnetic field,
$\mu_{\mathrm{\tiny\overline{MS}}}=\sqrt{|eB|}$, or to a fixed value,
$\mu_{\mathrm{\tiny\overline{MS}}}=1$ GeV.

Figure 4: Anomalous magnetic moments $a_{u}^{+}$ (black solid), $a_{u}^{-}$
(gray solid), $a_{d}^{+}$ (black dashed), $a_{d}^{-}$ (gray dashed), scaled by
their values in the vacuum, as functions of the external magnetic field. Here,
we consider dressed gluons with $m_{g}=0.3$ GeV and a renormalization scale
$\mu_{\mathrm{\tiny\overline{MS}}}=\sqrt{|eB|}$.

Figure 4 shows the case with a fixed gluon mass, $m_{g}=0.3$ GeV, and
$\mu_{\mathrm{\tiny\overline{MS}}}=\sqrt{|eB|}$. The first thing to be noticed
is the change of sign of the polarization components of the AMM, being
$a_{q}^{+}$ positive and $a_{q}^{-}$ negative. Moreover, we find that
$a_{q}^{-}\approx-a_{q}^{+}$, within $\sim 1\%$. One can also see that
$|a_{u}^{\pm}|>|a_{d}^{\pm}|$. This happens because the quark charge, through
$\ell_{q}$, dominates in the exponential term in Eq. (34). The plots present
an increase in $|a_{q}^{\pm}|$ as a function of $eB$. In fact, $|a_{u}^{\pm}|$
is greater than its corresponding vacuum value for $eB\gtrsim 2$ GeV2, and
$|a_{d}^{\pm}|$ is greater than its vacuum value for $eB\gtrsim 4$ GeV2. This
growth in the ratio $a_{q}^{\pm}/a_{\mathrm{vac}}$ is, however, somewhat
misleading. It basically represents the running of $\alpha_{s}$ with the
magnetic field in the denominator of the ratio. $a^{\pm}_{q}$ is, indeed,
almost constant, as can be seen in Fig. 5.

Figure 5: Anomalous magnetic moments $a_{u}^{+}$ (black solid), $a_{u}^{-}$
(gray solid), $a_{d}^{+}$ (black dashed), $a_{d}^{-}$ (gray dashed), as
functions of the external magnetic field. Here we consider dressed gluons with
$m_{g}=0.3$ GeV and a renormalization scale
$\mu_{\mathrm{\tiny\overline{MS}}}=\sqrt{|eB|}$.

The behavior of the ratio $a_{q}^{\pm}/a_{\mathrm{vac}}$ for a fixed scale
$\mu_{\mathrm{\tiny\overline{MS}}}=1$ GeV is basically the same, since the
main changes are contained in the global factor $\sim\alpha_{s}$. The
modifications in the running quark mass has no significant impact, producing
differences of $\sim 0.1\%$.

Figure 6: Anomalous magnetic moments $a_{u}^{+}$ (black solid), $a_{u}^{-}$
(gray solid), $a_{d}^{+}$ (black dashed), $a_{d}^{-}$ (gray dashed), scaled by
their values in the vacuum, as functions of the external magnetic field. Here
we consider dressed gluons with $m_{g}^{2}=\Pi_{g}$ and a renormalization
scale $\mu_{\mathrm{\tiny\overline{MS}}}=\sqrt{|eB|}$.

A totally different outcome occurs when we consider the dressed gluon mass as
given by the LLL polarization tensor, $m_{g}^{2}=\Pi_{g}$, and a
renormalization scale $\mu_{\mathrm{\tiny\overline{MS}}}=\sqrt{|eB|}$.
Contrary to what happens for a constant gluon mass, the absolute value of the
AMM now diminishes as the magnetic field increases, as illustrated in Fig. 6.
This behavior is expected since $\Pi_{g}$ grows with the magnetic field. If we
consider a fixed renormalization scale, $\mu_{\mathrm{\tiny\overline{MS}}}=1$
GeV, the plot looks very similar. However, the difference is larger than in
the case of a fixed gluon mass. The difference between the two cases increases
$\sim 10\%$ for $eB\sim 5$ GeV2 but stabilizes at higher values.

The strong magnetic background effects on the electron AMM are, of course,
analogous. One has only to substitute the running strong coupling by the
electromagnetic coupling and set the global factor $(N_{c}^{2}-1)/2N_{c}\to 1$
in Eq. (10), besides considering one flavor with charge $e$ and with $N_{c}\to
1$ in Eq. (36).

##  IV Summary and outlook

In this paper we computed the part relevant to the quark AMM of the one-loop
QCD correction to the photon-quark-antiquark vertex in an extremely strong
magnetic background, i.e., one in which $\sqrt{eB}$ is much larger than all
other mass scales. This justified the use of the lowest-Landau level
approximation. The solution exhibits infrared divergences. These can be tamed
by attributing an effective mass to gluons that acts as an infrared regulator,
inspired by what is expected to happen in the infrared sector of strong
interactions [56]. In this vein, we considered two standard scenarios: one
with a constant gluon mass, $m_{g}=0.3$ GeV, which sits in the ballpark of
$\Lambda_{QCD}$, and one in which the gluon mass comes from the one-loop
correction to its polarization tensor in the presence of an external magnetic
field in the LLL approximation, $m_{g}^{2}=\Pi_{g}$. In both cases we
investigated what happens if we set the renormalization scale to run with the
magnetic field, $\mu_{\mathrm{\tiny\overline{MS}}}=\sqrt{|eB|}$, or to a fixed
value, $\mu_{\mathrm{\tiny\overline{MS}}}=1$ GeV.

Considering the the total magnetic moment, we found that it is such that that
the anomalous magnetic moment contributes to the electric part $\sim
iE_{3}\sigma_{3}$, the other contributions to $\bm{B}\cdot\bm{\sigma}$ and
$i\bm{E}\cdot\bm{\sigma}_{\perp}$ being completely suppressed. We found that
the anomalous contribution can be naturally separated into two polarization
projections obtained from ${\cal P}_{\pm}=\frac{1}{2}(1\pm
is\gamma_{1}\gamma_{2})$, or equivalently $\frac{1}{2}(1+\Sigma_{3})$.

We presented results for the behavior of the anomalous magnetic moment of
quarks up and down, scaled by their values in the vacuum, as a function of the
external magnetic field, for the choices of gluon mass and renormalization
scale discussed above. We found that the case with a fixed gluon mass exhibits
an increase of the ratio $a^{\pm}/a_{\mathrm{vac}}$, while this observable in
the case with gluon mass given by $m_{g}^{2}=\Pi_{g}$ decreases as the
magnetic field increases. The latter setup, which seems to be closer to a more
physical description of the gluon mass, shows a strong suppression, i.e., the
quark AMM apparently vanishes for extremely large magnetic fields.

We also found a asymmetry between spin up and spin down components of the
anomalous correction, which might be more promising in searching possible
experimental observables.

A natural extension of the treatment presented here would be considering the
case with other Landau levels in the external lines of the triangular diagram.
This would bring information on the interplay between different levels in the
computation of the quark anomalous magnetic moment.

One could also incorporate effects from a thermal or dense medium to the
framework discussed here. To do that, one would need to consider hot and dense
magnetic QCD [57, 58, 59, 60, 37, 39] in the computation of the photon-quark-
antiquark vertex. The thermal case could be relevant for high-energy heavy-ion
collisions, and even play a role in the chiral magnetic effect scenario, since
the quark AMM could destroy the fermion zero mode in the presence of a
magnetic field. The case at high densities is of interest in magnetars, where
ultra-high magnetic fields can also be achieved. Although in this case the
effect on the equation of state seems to be minor [18], it contributes to
increase the level of pressure anisotropy [16]. Finally, a calculation of the
quark AMM within lattice QCD in the presence of a strong magnetic background
would provide a clean benchmark.

###### Acknowledgements.

This work was partially supported by CAPES (Finance Code 001), CNPq, FAPERJ,
INCT-FNA (Process No. 464898/2014-5) and ANID/FONDECYT under grant 1190192.

## References

  * Schwinger [1948] J. S. Schwinger, On Quantum electrodynamics and the magnetic moment of the electron, Phys. Rev. 73, 416 (1948).
  * Aoyama  _et al._ [2019] T. Aoyama, T. Kinoshita, and M. Nio, Theory of the Anomalous Magnetic Moment of the Electron, Atoms 7, 28 (2019).
  * Schwinger [1949] J. S. Schwinger, Quantum electrodynamics. III: The electromagnetic properties of the electron: Radiative corrections to scattering, Phys. Rev. 76, 790 (1949).
  * Schwinger [1951] J. S. Schwinger, On gauge invariance and vacuum polarization, Phys. Rev. 82, 664 (1951).
  * Baier  _et al._ [1976] V. N. Baier, V. M. Katkov, and V. M. Strakhovenko, Anomalous Magnetic Moment of the electron in the Magnetic Field, Yad. Fiz. 24, 379 (1976).
  * Baier and Katkov [2001] V. N. Baier and V. M. Katkov, Anomalous magnetic moment of the electron in a medium, Phys. Lett. A 280, 275 (2001), arXiv:hep-ph/0011340 .
  * Gusynin  _et al._ [1996] V. P. Gusynin, V. A. Miransky, and I. A. Shovkovy, Dimensional reduction and catalysis of dynamical symmetry breaking by a magnetic field, Nucl. Phys. B 462, 249 (1996), arXiv:hep-ph/9509320 .
  * Shovkovy [2013] I. A. Shovkovy, Magnetic Catalysis: A Review, Lect. Notes Phys. 871, 13 (2013), arXiv:1207.5081 [hep-ph] .
  * Miransky and Shovkovy [2015] V. A. Miransky and I. A. Shovkovy, Quantum field theory in a magnetic field: From quantum chromodynamics to graphene and Dirac semimetals, Phys. Rept. 576, 1 (2015), arXiv:1503.00732 [hep-ph] .
  * Fukushima  _et al._ [2008] K. Fukushima, D. E. Kharzeev, and H. J. Warringa, The Chiral Magnetic Effect, Phys. Rev. D 78, 074033 (2008), arXiv:0808.3382 [hep-ph] .
  * Kharzeev  _et al._ [2008] D. E. Kharzeev, L. D. McLerran, and H. J. Warringa, The Effects of topological charge change in heavy ion collisions: ’Event by event P and CP violation’, Nucl. Phys. A 803, 227 (2008), arXiv:0711.0950 [hep-ph] .
  * Skokov  _et al._ [2009] V. Skokov, A. Y. Illarionov, and V. Toneev, Estimate of the magnetic field strength in heavy-ion collisions, Int. J. Mod. Phys. A 24, 5925 (2009), arXiv:0907.1396 [nucl-th] .
  * Deng and Huang [2012] W.-T. Deng and X.-G. Huang, Event-by-event generation of electromagnetic fields in heavy-ion collisions, Phys. Rev. C 85, 044907 (2012), arXiv:1201.5108 [nucl-th] .
  * Bicudo  _et al._ [1999] P. J. A. Bicudo, J. E. F. T. Ribeiro, and R. Fernandes, The Anomalous magnetic moment of quarks, Phys. Rev. C 59, 1107 (1999), arXiv:hep-ph/9806243 .
  * Ferrer and de la Incera [2010] E. J. Ferrer and V. de la Incera, Dynamically Generated Anomalous Magnetic Moment in Massless QED, Nucl. Phys. B 824, 217 (2010), arXiv:0905.1733 [hep-ph] .
  * Strickland  _et al._ [2012] M. Strickland, V. Dexheimer, and D. P. Menezes, Bulk Properties of a Fermi Gas in a Magnetic Field, Phys. Rev. D 86, 125032 (2012), arXiv:1209.3276 [nucl-th] .
  * Fayazbakhsh and Sadooghi [2014] S. Fayazbakhsh and N. Sadooghi, Anomalous magnetic moment of hot quarks, inverse magnetic catalysis, and reentrance of the chiral symmetry broken phase, Phys. Rev. D 90, 105030 (2014), arXiv:1408.5457 [hep-ph] .
  * Ferrer  _et al._ [2015] E. J. Ferrer, V. de la Incera, D. Manreza Paret, A. Pérez Martínez, and A. Sanchez, Insignificance of the anomalous magnetic moment of charged fermions for the equation of state of a magnetized and dense medium, Phys. Rev. D 91, 085041 (2015), arXiv:1501.06616 [hep-ph] .
  * Mao and Rischke [2019] S. Mao and D. H. Rischke, Dynamically generated magnetic moment in the Wigner-function formalism, Phys. Lett. B 792, 149 (2019), arXiv:1812.06684 [hep-th] .
  * Chaudhuri  _et al._ [2019] N. Chaudhuri, S. Ghosh, S. Sarkar, and P. Roy, Effect of the anomalous magnetic moment of quarks on the phase structure and mesonic properties in the NJL model, Phys. Rev. D 99, 116025 (2019), arXiv:1907.03990 [nucl-th] .
  * Coppola  _et al._ [2019] M. Coppola, D. Gomez Dumm, S. Noguera, and N. N. Scoccola, Neutral and charged pion properties under strong magnetic fields in the NJL model, Phys. Rev. D 100, 054014 (2019), arXiv:1907.05840 [hep-ph] .
  * Mei and Mao [2020] J. Mei and S. Mao, Inverse catalysis effect of the quark anomalous magnetic moment to chiral restoration and deconfinement phase transitions, Phys. Rev. D 102, 114035 (2020), arXiv:2008.12123 [hep-ph] .
  * Chaudhuri  _et al._ [2020] N. Chaudhuri, S. Ghosh, S. Sarkar, and P. Roy, Effects of quark anomalous magnetic moment on the thermodynamical properties and mesonic excitations of magnetized hot and dense matter in PNJL model, Eur. Phys. J. A 56, 213 (2020), arXiv:2003.05692 [nucl-th] .
  * Ghosh  _et al._ [2020] S. Ghosh, N. Chaudhuri, S. Sarkar, and P. Roy, Effects of the anomalous magnetic moment of quarks on the dilepton production from hot and dense magnetized quark matter using the NJL model, Phys. Rev. D 101, 096002 (2020), arXiv:2004.09203 [nucl-th] .
  * Xu  _et al._ [2021] K. Xu, J. Chao, and M. Huang, Effect of the anomalous magnetic moment of quarks on magnetized QCD matter and meson spectra, Phys. Rev. D 103, 076015 (2021), arXiv:2007.13122 [hep-ph] .
  * Ghosh  _et al._ [2021] S. Ghosh, N. Chaudhuri, P. Roy, and S. Sarkar, Thermomagnetic modification of the anomalous magnetic moment of quarks using the NJL model, Phys. Rev. D 103, 116008 (2021), arXiv:2104.14112 [hep-ph] .
  * Farias  _et al._ [2022] R. L. S. Farias, W. R. Tavares, R. M. Nunes, and S. S. Avancini, Effects of the quark anomalous magnetic moment in the chiral symmetry restoration: magnetic catalysis and inverse magnetic catalysis, Eur. Phys. J. C 82, 674 (2022), arXiv:2109.11112 [hep-ph] .
  * Chaudhuri  _et al._ [2021] N. Chaudhuri, S. Ghosh, S. Sarkar, and P. Roy, Dilepton production from magnetized quark matter with an anomalous magnetic moment of the quarks using a three-flavor PNJL model, Phys. Rev. D 103, 096021 (2021), arXiv:2104.11425 [hep-ph] .
  * Kawaguchi and Huang [2023] M. Kawaguchi and M. Huang, Restriction on the form of the quark anomalous magnetic moment from lattice QCD results*, Chin. Phys. C 47, 064103 (2023), arXiv:2205.08169 [hep-ph] .
  * Mao [2022] S. Mao, Inverse catalysis effect of the quark anomalous magnetic moment to chiral restoration and deconfinement phase transitions at finite baryon chemical potential, Phys. Rev. D 106, 034018 (2022), arXiv:2206.12054 [nucl-th] .
  * Chaudhuri  _et al._ [2022] N. Chaudhuri, S. Ghosh, P. Roy, and S. Sarkar, Anisotropic pressure of magnetized quark matter with anomalous magnetic moment, Phys. Rev. D 106, 056020 (2022), arXiv:2209.02248 [hep-ph] .
  * Tavares  _et al._ [2024] W. R. Tavares, S. S. Avancini, R. L. S. Farias, and R. P. Cardoso, Artificial first-order phase transition in a magnetized Nambu–Jona-Lasinio model with a quark anomalous magnetic moment, Phys. Rev. D 109, 016011 (2024), arXiv:2309.04055 [hep-ph] .
  * Lin and Huang [2022] F. Lin and M. Huang, Magnetic correction to the anomalous magnetic moment of electrons, Commun. Theor. Phys. 74, 055202 (2022), arXiv:2112.01051 [hep-ph] .
  * Schwartz [2014] M. D. Schwartz,  _Quantum Field Theory and the Standard Model_ (Cambridge University Press, 2014).
  * Ayala  _et al._ [2020a] A. Ayala, J. L. Hernández, L. A. Hernández, R. L. S. Farias, and R. Zamora, Magnetic corrections to the boson self-coupling and boson-fermion coupling in the linear sigma model with quarks, Phys. Rev. D 102, 114038 (2020a), arXiv:2009.13740 [hep-ph] .
  * Villavicencio [2023] C. Villavicencio, Axial coupling constant in a magnetic background, Phys. Rev. D 107, 076009 (2023), arXiv:2212.04649 [hep-ph] .
  * Dominguez  _et al._ [2023a] C. A. Dominguez, M. Loewe, C. Villavicencio, and R. Zamora, Nucleon axial-vector coupling constant in magnetar environments, Phys. Rev. D 108, 074024 (2023a), arXiv:2308.05663 [hep-ph] .
  * Braghin  _et al._ [2024] F. L. Braghin, M. Loewe, and C. Villavicencio, Yukawa potential under weak magnetic field, Phys. Rev. D 109, 034014 (2024), arXiv:2308.12122 [hep-ph] .
  * Dominguez  _et al._ [2023b] C. A. Dominguez, M. Loewe, C. Villavicencio, and R. Zamora, Magnetic and density effects on the nucleon axial coupling, in  _26th High-Energy Physics International Conference in QCD_ (2023) arXiv:2309.05807 [hep-ph] .
  * Parle [1987] A. J. Parle, Quantum Electrodynamics in Strong Magnetic Fields. 4. Electron Selfenergy, Austral. J. Phys. 40, 1 (1987).
  * Elmfors  _et al._ [1996] P. Elmfors, D. Persson, and B.-S. Skagerstam, Thermal fermionic dispersion relations in a magnetic field, Nucl. Phys. B 464, 153 (1996), arXiv:hep-ph/9509418 .
  * Cucchieri and Mendes [2007] A. Cucchieri and T. Mendes, What’s up with IR gluon and ghost propagators in Landau gauge? A puzzling answer from huge lattices, PoS LATTICE2007, 297 (2007), arXiv:0710.0412 [hep-lat] .
  * Bogolubsky  _et al._ [2007] I. L. Bogolubsky, E. M. Ilgenfritz, M. Muller-Preussker, and A. Sternbeck, The Landau gauge gluon and ghost propagators in 4D SU(3) gluodynamics in large lattice volumes, PoS LATTICE2007, 290 (2007), arXiv:0710.1968 [hep-lat] .
  * Oliveira and Bicudo [2011] O. Oliveira and P. Bicudo, Running Gluon Mass from Landau Gauge Lattice QCD Propagator, J. Phys. G 38, 045003 (2011), arXiv:1002.4151 [hep-lat] .
  * Bicudo  _et al._ [2015] P. Bicudo, D. Binosi, N. Cardoso, O. Oliveira, and P. J. Silva, Lattice gluon propagator in renormalizable $\xi$ gauges, Phys. Rev. D 92, 114514 (2015), arXiv:1505.05897 [hep-lat] .
  * Aguilar  _et al._ [2008] A. C. Aguilar, D. Binosi, and J. Papavassiliou, Gluon and ghost propagators in the Landau gauge: Deriving lattice results from Schwinger-Dyson equations, Phys. Rev. D 78, 025010 (2008), arXiv:0802.1870 [hep-ph] .
  * Mueller  _et al._ [2014] N. Mueller, J. A. Bonnet, and C. S. Fischer, Dynamical quark mass generation in a strong external magnetic field, Phys. Rev. D 89, 094023 (2014), arXiv:1401.1647 [hep-ph] .
  * Aguilar  _et al._ [2016] A. C. Aguilar, D. Binosi, and J. Papavassiliou, The Gluon Mass Generation Mechanism: A Concise Primer, Front. Phys. (Beijing) 11, 111203 (2016), arXiv:1511.08361 [hep-ph] .
  * Dudal  _et al._ [2008] D. Dudal, J. A. Gracey, S. P. Sorella, N. Vandersickel, and H. Verschelde, A Refinement of the Gribov-Zwanziger approach in the Landau gauge: Infrared propagators in harmony with the lattice results, Phys. Rev. D 78, 065047 (2008), arXiv:0806.4348 [hep-th] .
  * Capri  _et al._ [2015] M. A. L. Capri, D. Dudal, D. Fiorentini, M. S. Guimaraes, I. F. Justo, A. D. Pereira, B. W. Mintz, L. F. Palhares, R. F. Sobreiro, and S. P. Sorella, Exact nilpotent nonperturbative BRST symmetry for the Gribov-Zwanziger action in the linear covariant gauge, Phys. Rev. D 92, 045039 (2015), arXiv:1506.06995 [hep-th] .
  * Peláez  _et al._ [2021] M. Peláez, U. Reinosa, J. Serreau, M. Tissier, and N. Wschebor, A window on infrared QCD with small expansion parameters, Rept. Prog. Phys. 84, 124202 (2021), arXiv:2106.04526 [hep-th] .
  * Fukushima [2011] K. Fukushima, Magnetic-field Induced Screening Effect and Collective Excitations, Phys. Rev. D 83, 111501 (2011), arXiv:1103.4430 [hep-ph] .
  * Bandyopadhyay  _et al._ [2016] A. Bandyopadhyay, C. A. Islam, and M. G. Mustafa, Electromagnetic spectral properties and Debye screening of a strongly magnetized hot medium, Phys. Rev. D 94, 114034 (2016), arXiv:1602.06769 [hep-ph] .
  * Ayala  _et al._ [2020b] A. Ayala, J. D. Castaño Yepes, C. A. Dominguez, S. Hernández-Ortiz, L. A. Hernández, M. Loewe, D. Manreza Paret, and R. Zamora, Thermal corrections to the gluon magnetic Debye mass, Rev. Mex. Fis. 66, 446 (2020b), arXiv:1805.07344 [hep-ph] .
  * Ayala  _et al._ [2020c] A. Ayala, J. D. Castaño Yepes, M. Loewe, and E. Muñoz, Gluon polarization tensor in a magnetized medium: Analytic approach starting from the sum over Landau levels, Phys. Rev. D 101, 036016 (2020c), arXiv:1912.07136 [hep-th] .
  * Mena and Palhares [2023] C. Mena and L. F. Palhares, Quark-photon vertex in confining models, (2023), arXiv:2311.14178 [hep-ph] .
  * Blaizot  _et al._ [2013] J.-P. Blaizot, E. S. Fraga, and L. F. Palhares, Effect of quark masses on the QCD pressure in a strong magnetic background, Phys. Lett. B 722, 167 (2013), arXiv:1211.6412 [hep-ph] .
  * Ayala  _et al._ [2021] A. Ayala, L. A. Hernández, M. Loewe, and C. Villavicencio, QCD phase diagram in a magnetized medium from the chiral symmetry perspective: the linear sigma model with quarks and the Nambu–Jona-Lasinio model effective descriptions, Eur. Phys. J. A 57, 234 (2021), arXiv:2104.05854 [hep-ph] .
  * Fraga  _et al._ [2023a] E. S. Fraga, L. F. Palhares, and T. E. Restrepo, Hot perturbative QCD in a very strong magnetic background, Phys. Rev. D 108, 034026 (2023a), arXiv:2303.12140 [hep-ph] .
  * Fraga  _et al._ [2023b] E. S. Fraga, L. F. Palhares, and T. E. Restrepo, Cold and dense perturbative QCD in a very strong magnetic background, (2023b), arXiv:2312.13952 [hep-ph] .


# Modular Growth of Hierarchical Networks: Efficient, General, and Robust
Curriculum Learning

Mani Hamidi1,2,∗, Sina Khajehabdollahi1, Emmanouil Giannakakis1,2, Tim J.
Schäfer 1,2,
Anna Levina1,2, Charley M. Wu1,2
1Department of Computer Science, University of Tübingen, Tübingen, Germany
2Max Planck Institute for Biological Cybernetics, Tübingen, Germany
<EMAIL_ADDRESS>
###### Abstract

Structural modularity is a pervasive feature of biological neural networks,
which have been linked to several functional and computational advantages.
Yet, the use of modular architectures in artificial neural networks has been
relatively limited despite early successes. Here, we explore the performance
and functional dynamics of a modular network trained on a memory task via an
iterative growth curriculum. We find that for a given classical, non-modular
recurrent neural network (RNN), an equivalent modular network will perform
better across multiple metrics, including training time, generalizability, and
robustness to some perturbations. We further examine how different aspects of
a modular network’s connectivity contribute to its computational capability.
We then demonstrate that the inductive bias introduced by the modular topology
is strong enough for the network to perform well even when the connectivity
within modules is fixed and only the connections between modules are trained.
Our findings suggest that gradual modular growth of RNNs could provide
advantages for learning increasingly complex tasks on evolutionary timescales,
and help build more scalable and compressible artificial networks.

## Introduction

From bacteria (Andrews,, 1998), to brains (Sporns and Betzel,, 2016), to man-
made artifacts (Lake et al.,, 2015), many things are composed of modular
components, specialized for different purposes and capable of being recombined
in distinct configurations to solve new problems. In cognitive science,
modularity (Fodor,, 1983) has played an important role in understanding
intelligent behavior as composition of modular, symbolic representations
(Rubino et al.,, 2023; Zhou et al.,, 2024) while evolutionary accounts have
explored how selection pressure towards reducing connection costs may favor
modular solutions (Clune et al.,, 2013). By constraining the search space
(Happel and Murre,, 1994) and favoring more computationally efficient
solutions (Yuan et al.,, 2023), modular architectures (Amer and Maul,, 2019)
offer a complementary and more biologically plausible account of intelligence
(Cosmides and Tooby,, 1997), at a time when current trends in deep learning
have pursued scale at all costs (Shen et al.,, 2023).

Taking inspiration from the modular duplication of entire body parts during
gene duplication events (Garcia-Fernàndez,, 2005), we explore the functional
utility of modular growth in adapting a recurrent network to a memory task of
gradually increasing complexity. We find that coordination of modular growth
together with a learning curriculum facilitates a surprising array of
advantages, both in terms of performance and costs.

Specifically, we focus on the domain of multi-timescale signal processing
(Bathellier et al.,, 2008; Panzeri et al.,, 2010; Safavi et al.,, 2023). In
tasks such as speech recognition (Graves et al.,, 2013), time-series
prediction (Chung et al.,, 2014; Torres et al.,, 2021), and navigation (Morcos
and Harvey,, 2016), both biological and artificial agents commonly need to
represent and remember relatively long timescales in the underlying network
dynamics. In artificial neural networks, such dynamics can arise not only via
the training of recurrent connectivity, but also by explicitly training the
timescales of individual neurons (Perez-Nieves et al.,, 2021; Tallec and
Ollivier,, 2018; Quax et al.,, 2020; Yin et al.,, 2020; Fang et al.,, 2021;
Smith et al.,, 2023). Trainable timescales have been linked with improved
performance for rate- (Tallec and Ollivier,, 2018; Quax et al.,, 2020) and
spiking networks (Yin et al.,, 2020; Fang et al.,, 2021; Perez-Nieves et al.,,
2021). However, recent work has suggested that greater reliance on
connectivity, rather than trainable timescales, is associated with superior
performance and robustness of RNNs in similar memory tasks (Khajehabdollahi et
al.,, 2024). Thus, these current debates suggest that the interaction between
trainable and connectivity-based mechanisms, together with their dependence on
architectural and training decisions remains relatively poorly understood,
even in simple networks performing rather trivial tasks.

One such architectural distinction between artificial and biological networks
is the presence of modular topologies (Litwin-Kumar and Doiron,, 2012;
Chaudhuri et al.,, 2014; Zeraati et al.,, 2023; Shi et al.,, 2023), which
organize various neuronal types into higher level clusters with different
timescales (Greengard,, 2001) and specific local connectivity patterns (Schaub
et al.,, 2015). In particular, hierarchical structures in the cortex are
associated with a gradual increase in intrinsic timescales (Murray et al.,,
2014; Manea et al.,, 2022) that has been linked to the computational
requirements of the tasks performed by different cortical regions (Murray et
al.,, 2014).

Although artificial hierarchical networks have been used for temporal tasks in
the past (El Hihi and Bengio,, 1995), they only explored shallow and static
hierarchies that lacked curriculum-based growing mechanisms. Thus, to our
knowledge, the role of hierarchical structure in accommodating the emergence
of longer timescales has yet to be investigated.

Here, we address this gap by comparing standard RNNs (non-modular) to growing
modular networks, using memory tasks with increasing complexity (Fig. 1).
Building on the success of different curriculum-based training strategies that
support the learning of longer timescales (Khajehabdollahi et al.,, 2024), we
additionally introduce module-level duplicative growth of the network at every
step along the curriculum. This modular strategy is loosely reminiscent of
growth over evolutionary timescales that gradually accommodate adaptation to
more complex tasks and environments (Lui et al.,, 2011). Our networks achieve
superior performance and robustness, in line with other studies that use
neural growth models to train RNNs in a different domain (Najarro et al.,,
2023).

Our findings suggest that the well-calibrated growth of structured networks
can significantly reduce the number of trainable parameters and training steps
required to solve complex temporal tasks.

Figure 1: Network Structures. Non-modular networks consist of a fixed number
of neurons ($M=20,54,91,128$), with connections that are retrained at each
curriculum step. After a given accuracy is reached for a task of length $N$, a
new readout head (consisting of 2 neurons) for $N+1$ is added and the network
is retrained for all previous $N$. In contrast, the modular network adds a
much smaller RNN module ($M_{m}=5,10,15,20$) for every readout head that is
added in the curriculum. Thus, each readout head is attached to a separate RNN
module, rather than one large reservoir as in non-modular networks.

## Methods

We train non-modular and modular networks on the $N$-parity task at increasing
complexity levels $N$ and compare their performance, robustness, and learning
dynamics. Models are trained using back-propagation-through-time with a
stochastic gradient descent optimizer, a cross-entropy loss function, and
curriculum learning.

### Task

An $N$-parity task is a memory task that is commonly used to assess the
capabilities of recurrent network architectures (Stork and Allen,, 1992; Hohil
et al.,, 1999). The task requires the accurate retention of a sequence of
binary digits to perform a modulo 2 summation over the last $N$ digits in the
sequence. Therefore $N$ captures the complexity of the task. The input is a
random binary sequence $S$ with length $L$ chosen uniformly from the interval
$\\{N+2,4N\\}$, with one bit provided to the network at each time step. The
network must output the binary sum ($\mathrm{XOR}$) of the last $N$ digits. To
update the output at every time step, the network must learn to store some
representation of the values and order of the last $N$ digits in memory.

Although simple, this task provides a foundation for testing the
representation learning and memory capacities of artificial neural networks,
while simultaneously allowing us to control the difficulty of the task 1 bit
at a time using $N$. Furthermore, multiple $N$-parity tasks can be computed
for a single sequence, allowing the possibility of training concurrent
tasks/readout layers on the same inputs. This concurrent training is used to
encourage more universal features that can be shared between the different
tasks and to prevent catastrophic forgetting.

The success criterion for completing task $N$ is therefore defined as jointly
satisfying i) $>98\%$ accuracy on task $N$ and ii) maintaining an average of
$>98\%$ accuracy on all previous tasks. Accuracy is reported as an average
over multiple tests, with random chance yielding $50\%$. We report network
performance using $N_{\text{solved}}$ which describes the largest task that it
was able to solve after $60$ epochs.

Figure 2: Performance comparison of modular and non-modular architectures with
different numbers of neurons (color bars) allocated to their recurrent
processing unit. For modular networks, the number of neurons is reported per
module. For non-modular networks, the number indicates the total number of
neurons in the single recurrent core. These numbers were chosen to allow for a
comparable total number of learnable parameters for a given task difficulty
$N$. (a) Learning curves show that modular architectures with sufficient
($>5$) neurons can solve a new task at every epoch while the non-modular
networks plateau in their learning ability. (b) Pareto frontier of performance
after 60 training epochs, shows that the modular architecture always achieves
better performance for the same number of parameters. Note that the saturation
of performance at $N_{\text{solved}}=60$ for the modular architecture is due
to the choice of training epochs. (c) Generalization performance of networks
trained on a task difficulty of $N$ and then tested on tasks of $N+K$.
Accuracy is measured as the percentage of correct trials and error bars
indicate the standard deviation. Results are averages over 4 non-modular and 3
modular networks, and errorbars show STD.

### Networks

We use simple rate neurons with a trainable time constant $\tau$ whose
activity is defined by:

$r_{i}(t)=\left(1-\frac{\Delta t}{\tau_{i}}\right)\cdot r_{i}(t-\Delta
t)+\frac{\Delta t}{\tau_{i}}\cdot\left[C_{t}\right]_{\alpha}$ (1)

$C_{t}$ is the input at each timestep $t$ and the non-linearity
$[\cdot]_{\alpha}$ is the leaky ReLU function with negative slope $\alpha$,
given by:

$[x]_{\alpha}=\left\\{\begin{array}[]{ll}x,&\ x\geq 0\\\ \alpha\cdot x,&\
x<0.\end{array}\right.$ (2)

During learning, we train both the network connectivity (recurrent and
feedforward) as well as the individual neuron timescales $\tau_{i}$. The
trainable timescale $\tau_{i}$ has a minimum value of 1 indicating a neuron
that reacts only to the shortest timescales (i.e., the current input $C_{t}$),
whereas larger values integrate information from longer timescales in the
past.

#### Non-modular networks.

In non-modular networks (Jaeger,, 2002), the input at each timestep $t$ for
neuron $i$ is:

$C_{t}=\sum_{j\neq i}^{M}W^{R}_{ij}\cdot r_{j}(t-\Delta t)+W^{I}_{i}\cdot
S(t)+b_{i}.$ (3)

Here, $W^{R}$ is the recurrent connectivity, $W^{I}$ is the feedforward input
connections, $S(t)$ is the input signal, $b_{i}$ is the neuron bias, and time
discretization $\Delta t=1$.

The non-modular networks follow the multi-head curriculum in (Khajehabdollahi
et al.,, 2024): At the first curriculum step, a single linear readout head is
trained to solve the task for $N=2$. Upon successfully completing task $N$, a
new readout head is added and is trained to solve the $N+1$-parity task.
Previous readout heads continue to be trained after every step in the
curriculum. Thus, at the $m_{th}$ step of the curriculum, the network has $m$
readout heads solving the task for $N=2,\dots m-1$.

#### Modular networks.

In modular networks, each neuron receives input from the other neurons of the
same module, and for modules $m>1$, neurons also receive input from a feed-
forward connection from the previous module $m-1$:

$\displaystyle C_{t}^{m}={}$ $\displaystyle\sum_{j\neq
i}^{M_{m}}W^{R}_{ij}\cdot r_{j}^{m}(t-\Delta t)$ (4) $\displaystyle+$
$\displaystyle\sum_{k}^{M_{m}}W^{FF}_{ik}\cdot r_{k}^{m-1}(t-\Delta t)$
$\displaystyle+$ $\displaystyle\;W^{I}_{i}\cdot S(t)+b_{i}$

where $W^{FF}$ is the feed-forward connectivity from the previous module and
$r_{k}^{m-1}(t-\Delta t)$ is the activity of neuron $k$ in module $m-1$ from
the previous timestep. $M_{m}$ indicates the number of neurons in each module,
which is fixed to be one of $M_{m}\in[5,10,15,20]$ in our experiments.

For training the modular networks, we follow a growing curriculum as follows:
Starting with a single module with a small population (5, 10, 15, or 20
neurons), we train a linear readout for $N=2$. Then, at each curriculum step,
we add a new module and readout head for each new $N$. The connectivity of the
new module is a copy of the corresponding (feedforward and recurrent)
connectivity for the $N-1$ module and receives feedforward input from it. Thus
each module is specialized to solving the task for a single $N$. At each
curriculum step, we freeze all network connections for all modules except the
last and train only the recurrent and feedforward connections (input from
sequence and input from previous module) of the last module.

### Timescale Estimation

Apart from the trainable parameters $\tau_{i}$ that define the timescale of
each neuron explicitly, we also compute the effective timescale of a neuron,
which is influenced by its intrinsic trained $\tau_{i}$ and the extrinsic
modulation of its activity through its connections with the rest of the
network. We determine this network-mediated effective timescale for individual
neurons by calculating the lagged autocorrelations (AC) of their activity
during task performance, following previously developed methods
(Khajehabdollahi et al.,, 2024). The AC of each neuron is then modeled with an
exponential function featuring one or two timescales, selecting the best-fit
model based on the Akaike Information Criterion (Akaike,, 1974). To avoid
estimation bias, we utilize long time series ($10^{5}$ time steps) of activity
(Zeraati et al.,, 2022).

### Robustness

We evaluate robustness by measuring the accuracy of the network after
perturbing one of its three trained parameter groups, $W^{R}$, $W^{FF}$ or
$\tau$. We define the magnitude of the perturbation $\varepsilon_{W}$ as a
function of the magnitude of the weights (Wu et al.,, 2020):

$\widetilde{W}=W+\varepsilon_{W}\frac{\xi_{W}}{||\xi_{W}||}||W||$ (5)

where $\xi_{W}\sim\mathcal{N}(0,\mathbb{I}^{n\times n})$ and $||\cdot||$
represents the Frobenius norm. This normalization allows for comparable
amounts of perturbation across networks of different types and sizes. We also
only perturb $\tau_{i}$ in positive direction to avoid $\tau<1$.

## Results

Our findings are divided into two main sections. In the first part, we
demonstrate that modular networks following an iterative growing curriculum
outperform equivalent non-modular RNNs in terms of task performance, training
speed, generalizability, and robustness to perturbation of their learned
connectivity. In the second part, using perturbation and ablation techniques,
we investigate how different aspects of modular networks (feedforward vs.
recurrent connectivity and trainable vs. effective single-neuron timescales)
contribute to their computational capabilities.

### Modular vs. Non-Modular Networks

#### Performance and Generalization

We first compare a set of modular networks with different module sizes
($M_{m}=5,10,15,20$) against a set of non-modular networks ($M=20,54,91,128$)
with an equivalent number of trainable parameters (at $N=30$). Performance is
calculated as the average task difficulty solved, $N_{\text{solved}}$, both as
a function of training time (Fig. 2a) and the number of parameters (Fig. 2b).

Figure 3:  Weight perturbations degrade performance in both modular and non-
modular networks of different sizes (colours). (a) Modular networks are more
robust to perturbation of connections (modular networks: feedforward and
recurrent weights). (b) Non-modular networks are more robust to the
perturbation of single-neuron trained timescales.

Overall, we find that all modular networks reach a higher task difficulty
($N_{\text{solved}}$) than the corresponding non-modular networks (Fig. 2a).
Except for the smallest modular networks ($M_{m}=5$), all modular networks
maintain a steady, linear progression through the curriculum, solving every
task after a single epoch of training. As an extreme test of the limits of the
modular network, we were able to reach a maximum of $N=200$ in one network
before ending the simulation, with even greater capabilities being
theoretically possible given sufficient training time. The non-modular
networks, on the other hand, plateau much earlier in their progression through
the curriculum. Thus, for a fixed number of trainable parameters, our results
suggest there is a modular architecture that is capable of solving more
complex tasks with the same amount of training.

We also explore the impact of a limited budget for learnable parameters as an
alternative evaluation of the cost-benefit trade-off between two
architectures. Figure 2b shows the Pareto frontiers of the architectures after
60 epochs, indicating that a modular architecture can solve a task with fewer
parameters than a non-modular network.

Finally, we test generalization performance on networks of size $M_{m}=15$ and
$M=91$, by training them to $N_{solved}=10$ and then test their ability to
solve $N=N_{solved}+K$, for $K=1,2,\dots$. For modular networks, the $N=10$
module is now tested for solving $N=10+K$, while the non-modular network is
tested on a new read-out head for $N=10+K$. We only allow 10 epochs of
additional training on the new task. The results are shown in Figure 2c, where
although performance decays with $K$ for both types of networks, modular
networks consistently reach a higher accuracy compared to non-modular
networks, thus demonstrating superior generalization.

#### Robustness to Perturbations.

Next, we tested the robustness of both network architectures to perturbations
of connectivity weights (Fig. 3a) and trained timescale parameters (Fig. 3b).
The perturbations were proportional to the magnitude of the weights of the
layer being targeted (Eq. 5). For different perturbation levels, we measure
the network performance via the number of tasks that the network is still able
to solve with accuracy $>0.9$.

Overall, Figure 3 shows the degradation of performance as a function of the
size of the perturbation on the connections (both recurrent $W^{R}_{ij}$ and
feedforward $W^{FF}_{ij}$ in the case of the modular networks) and the
timescale parameters ($\tau_{i}$) of the network. Modular networks not only
have better initial performance but are also more robust against connectivity
perturbations (Fig. 3a). This is visible in the later inflection point,
reflecting how larger perturbations are required to degrade the performance of
modular networks. However, modular networks are more sensitive to the
perturbation of time-scale parameters (Fig. 3b). A special case is modular
networks with $M_{m}=5$ that appear to be more sensitive to both perturbations
than their larger counterparts. For all other cases, the number of neurons
does not significantly impact the robustness of either network type. So the
inflection point of degrading performance is an inherent feature of network
architecture. In sum, modular networks are significantly more robust against
connectivity perturbations, but are more sensitive to perturbation of time-
scale parameters.

### Functional Analysis of Modular Networks

#### Trained vs. Effective Timescales.

Figure 4: Change in trained and effective timescales for different
$N_{solved}$. (a) The trained timescale of the modular network stays the same
across modules, while the timescales of the non-modular network converge to 1
as the $N_{solved}$ increases. (b) The effective timescales of both networks
increase steadily with $N_{solved}$. A modular network with a module size of
15 neurons and a non-modular network with an equal number of trainable
parameters were used. Results are averages over 4 networks.

To better understand the increased sensitivity of modular networks to
timescale perturbations, we examine the mechanism by which long timescales
emerge in the networks. We do so by examining both the trained timescale
parameterized by $\tau_{i}^{m}$ (Fig. 4a) and also the effective timescale
inferred from the activity of each neuron (Fig. 4b). In the modular case, the
average trained and effective timescales of neurons are reported separately
per module solving $N=N_{\text{solved}}$. In the non-modular case, the average
is taken over all neurons in the entire network at the time when it has just
solved the $N=N_{\text{solved}}$ task.

First, we find that neurons in modular networks maintain stable trained
timescales, $\langle\tau_{i}\rangle_{m}$, throughout the curriculum, but non-
modular networks decay to the fastest possible rate of one (Fig. 4a). In both
cases, these trained timescale values seem too small to account for the long
timescales that are required for the large $N$ tasks they can solve. This
motivated us to examine the effective timescales, which are a product of the
connectivity structure of the whole network.

Figure 4b shows the increase in the average effective timescales of neurons in
both modular and non-modular networks for increasingly complex tasks $N$. This
observation is consistent with recent observations (Khajehabdollahi et al.,,
2024) suggesting a circuit-level implementation of memory by learning
appropriate connection weights ($W$) rather than learning slow timescales
($\tau_{i}$) at the single neuron level. Notably, neurons in the modular
networks appear to harbor memory of longer timescales compared to the non-
modular counterparts at the same stage in the curriculum.

Figure 5:  Feedforward connections are more sensitive than recurrent
connections in modular networks. Here, we use $M_{m}=15$ but achieve
qualitatively similar results for other network sizes. (a) Perturbing
feedforward weights affect downstream modules more strongly than earlier
modules. (b) The recurrent weights exhibit a similar qualitative pattern but
are quantitatively more robust against the same levels of perturbations. (c)
Variability of feedforward versus recurrent connection weights across modules,
where variability is inversely related to the degree of conservation (i.e.,
the amount of shared weights from one module to the next). Accuracy is
averaged over 5 networks, 10 perturbations, and 1000 continuous evaluations.
Error bars show SEM.

#### Recurrent vs Feedforward Connections.

We now turn to how the topological constraints imposed by the modular
architecture suggest different qualitative roles for feedforward and recurrent
connections. For instance, the feedforward connections serve as the only
bottleneck through which the recurrent computations of each module are reused
by subsequent modules (Fig. 1), making them potentially more vulnerable to
noise, whereas the recurrent connections may be more robust. Furthermore, the
dependence of long effective neural timescales on appropriate connectivity
patterns, encourages us to examine the distinct roles that feedforward
($W^{FF}$) and recurrent ($W^{R}$) weights play in performance.

To do so, we performed separate perturbation analyses targeting either
feedforward $W^{FF}$ or recurrent connections $W^{R}$. We evaluate performance
by measuring the accuracy of each module’s prediction for its corresponding
task $N$.

Figure 5a-b shows that the feedforward connections are more sensitive to
perturbations relative to the recurrent weights. Furthermore, modules
corresponding to more difficult tasks (indicated by darker colors; Fig. 5a-b),
suffer more because they are affected by the cumulative effect of all the
perturbations to preceding modules. Note that the effect of a perturbation is
always downstream, due to the sequentially connected architecture of modules,
such that a single perturbation equally affects all downstream modules.

To understand the basis for the sensitivity of $W^{FF}$, we hypothesize that
$W^{FF}$ are subject to more fine-tuning via back-propagation, in order to
facilitate specialization to each new task. On the other hand, we expect
$W^{R}$ to be more conserved, and subject to less modification from one task
to the next. We define the change in the connection weight between nodes $i$
and $j$ across two consecutive modules as $\Delta\overline{W_{ij}}$, where we
normalized each weight by its mean and standard deviation to correct for
systematic differences across module sizes (LeCun et al.,, 2012). The variance
of the normalized weight change, $\sigma^{2}_{\Delta\overline{W}_{ij}}$, is
the change in the normalized weight,
$\Delta\overline{W_{ij}}=\overline{W}_{ij}^{m}-\overline{W}_{ij}^{m-1}$
between corresponding neurons, $i,j$, in consecutive modules, $m,m-1$. This
serves as an inverse proxy for the degree to which the weights are conserved
during training, where lower variance corresponds to greater conservation of
weights.

Figure 5c shows how recurrent weights are more conserved than feedforward
weights (lower variance). This supports the hypothesis that sensitivity of
feedforward connections to perturbations is due to a more general sensitivity
to modification, both from noise and from back-propagated error during
training.

#### Weight Freezing.

The low variability of the recurrent weights motivated us to test the degree
to which a network can continue to learn despite a complete freezing of either
its recurrent or the feedforward weights, after passing the very first step in
the curriculum. A freezing of weights corresponds to a reduction of
computational costs, by converting learnable parameters of a model into fixed
biases. In evolutionary biology, this is known as the Baldwin effect (Gruau
and Whitley,, 1993), where plastic behaviors become fixated in response to
environmental stability. In the case of frozen recurrent weights, we trained
the first module for $N=2$ and then duplicated the learned recurrent weights
on all subsequent modules, training only the feedforward connections between
modules. Conversely, in the case of frozen feedforward weights, we train the
feedforward connections between the $N=2$ and $N=3$ modules and then duplicate
this trained connectivity between all subsequent modules, only training the
internal recurrent connections of each module.

Figure 6a shows that networks with frozen feedforward weights fail
catastrophically to solve the task for any $N\geq 3$. However, equivalent
networks with frozen recurrent connectivity (resembling deep reservoirs;
Gallicchio and Micheli,, 2021), perform much better, comfortably reaching an
$N_{solved}\approx 15$ before stagnating (Fig. 6a) and even outperforming non-
modular networks that taper at $N\approx 10$. Thus, while both network types
underperform the original network for which both recurrent and feedforward
connections are trained, we see that, in line with our previous results, the
learnability of the feedforward connections appears to be much more important
for task performance.

Figure 6: The impact of weight freezing and duplication on feedforward and
recurrent weights. (a) Weight freezing. Solid lines are a modular network with
$M_{m}=5$, selected for the strongest differences. Here, we can see that
freezing recurrent connections (orange) after solving the first task ($N=2$)
still needs to relatively good performance, whereas freezing feedforward
connections (teal) severely impairs learning. In networks of larger module
size ($M_{m}=10$; dotted lines), the difference between frozen recurrent
weights and a fully trained network disappears (both dotted orange), while
frozen feedforward weights remain impaired. (b) Weight duplication. New
modules are initialized with either a duplicated copy (T) of their weights as
they appeared in the previous module, or a random initialization of weights
(F). Duplication of recurrent weights confers a slight advantage (green),
while either duplicative or randomly initialized feedforward weights both
confer a slight impairment to performance. Results are averages over 5
networks and shaded regions indicate SEM.

#### Weight Duplication.

The weight freezing experiments involved the reuse (or duplication) of weights
from the previous module followed by a complete blockage of their training.
Here we consider duplicative versus random initialization of the weights, and
maintain the plasticity of weights during training. Theoretically, duplication
can be beneficial by amortizing the costs of learning through the reuse of
already-trained weights, or it could be detrimental by initializing the
network in a state unsuitable for the new task. Empirically, previous work
using weight-agnostic neural networks (Gaier and Ha,, 2019), compositional
pattern-producing networks (Stanley,, 2007), and hypernetworks (Ha et al.,,
2016), have demonstrated the utility of similar “weight-sharing” schemes that
support high-performing neural networks with a fraction of the parameter
count. Motivated by these past successes, here we explore whether the
duplication of recurrent or feedforward components can accelerate the training
of the network, or if such methods are detrimental by trapping the behavior of
the network in a local optima.

Figure 6b shows the results of these experiments, where weights were either
duplicated (T for true) or not (F for false); in case of the latter, they were
initialized randomly following a uniform distribution. The lines show the
average performance of 6 networks in each of the 4 different conditions, with
and without duplication of feedforward $W^{FF}$ and/or recurrent $W^{R}$
weights. These results show that the duplication of recurrent weights $W^{R}$
appears to account for an unambiguous advantage (Fig. 6b; green curves),
compared to the duplication of feedforward weights, which do not confer any
additional benefits. Indeed, duplicating feedforward weights confers little to
no improvement to when neither weights are duplicated (Fig. 6b, brown curves).
We also experimented with noisy duplication, but given the weak overall
effect, we only show the results of exact duplication, without noise.
Furthermore, different duplication strategies mattered less at larger module
sizes, so we focused only on networks with five nodes per module $M_{m}=5$.

These results are consistent with our previous analysis showing that
feedforward connections must diverge from their past configurations to
specialize for their respective task. Therefore, as long as the feedforward
connections remain plastic after a growth event, they are molded strongly by
the error signal from their new task, making their initial state irrelevant.
In contrast, we have already shown how keeping the same or similar recurrent
weights across multiple tasks is relatively adaptive. Thus, duplication simply
initializes the weights near their last functional state, eliminating the need
to retrain them to that stage at every epoch, and allows for minimal fine-
tuning to improve their adaptiveness over longer timescales.

## Discussion

In sum, we find that modular networks, trained via a growing curriculum to
solve a memory task of adjustable difficulty, outperform standard recurrent
neural networks (RNNs) of fixed size in terms of accuracy, training time,
generalization, and robustness to perturbations of the learned connectivity.

Biological neural networks are highly structured (Hagmann et al.,, 2008;
Sporns and Zwi,, 2004; Sporns and Betzel,, 2016), with multiple neuron types
each following distinct connectivity patterns (Peng et al.,, 2021; Liu et
al.,, 2023) that form complex circuits repeating across brain regions (Douglas
and Martin,, 2004; Shepherd,, 2011). The repetition of similar complex
structures in the brain potentially enables the development of very complex
neural circuits from limited genetic information (Rakic,, 2009; Geschwind and
Rakic,, 2013; Stanley,, 2007), and may also play a role in more efficient
information processing (Bassett and Bullmore,, 2006; Sporns,, 2013). Moreover,
inhomogeneous network topology has been shown to generate non-trivial dynamics
(Litwin-Kumar and Doiron,, 2012) and act as an inductive bias for local
learning (Giannakakis et al.,, 2023). Here, our findings add to this
literature by suggesting relatively simple topological modularity can be
utilized to align network dynamics with task requirements, thus increasing
performance and reducing training time and cost.

Our modular architecture allows us to study the roles played by different
aspects of the network’s structure. In particular, we focus on the distinction
between the recurrent weights within each module and the feedforward weights
connecting adjacent modules. Our analysis indicates that the training and task
specialization of feedforward weights between modules are vastly more
important than the recurrent connections within each module. This suggests
that the function of the recurrent connectivity is largely limited to creating
appropriate and generic dynamical units (and can thus be effectively reused
through duplication), while the feedforward connections control the flow of
information that is necessary for solving the specifics of a given task.

While our network included trainable single-neuron timescales, which have been
utilized in RNNs solving temporal tasks (Tallec and Ollivier,, 2018; Quax et
al.,, 2020; Khajehabdollahi et al.,, 2024), we find that their importance is
minimal in terms of network performance. The distribution of trained
timescales barely changes across modules, suggesting that the network fully
relies on its connectivity (particularly the hierarchical topology) for
developing the long timescales necessary for solving temporal tasks. This
augments previous findings suggesting that network topology can generate long
timescales with relatively random connections (Khajehabdollahi et al.,, 2024).

This understanding of the relative importance of different network components,
allows us to test whether training even fewer parameters can lead to
comparable results. We find that training the networks as deep reservoirs
(Gallicchio and Micheli,, 2021) leads to reasonably good performance despite
reducing the trainable parameters by half. This suggests that functional sub-
networks can be used as building blocks in larger systems with little or no
training, which seems like a promising avenue for future research.

The importance of initializing neural networks in an appropriate dynamical
regime has been widely explored (Zierenberg et al.,, 2020; Khajehabdollahi et
al.,, 2022), and our findings here suggest that hierarchical modular
structures could be used to generate networks with beneficial dynamics for
learning temporal tasks. Recycling functional structures and incorporating
them in complex networks of interacting sub-units is a widespread
characteristic of biological networks (Felleman and Van Essen,, 1991) and has
shown promise in artificial settings (Happel and Murre,, 1994; Sharkey,, 1996;
Amer and Maul,, 2019). Our findings suggest that this approach is not only
more efficient in terms of the number of trainable parameters, but it can also
boost performance as well as the ability to generalize.

#### Limitations & Future Work.

A basic premise of our work was that the power and cost of a neural network
are proportional to the number of its trainable parameters. The number of
parameters, therefore, provided a basis for comparing neural networks of two
different architectures. Our results indicate that one can maintain or even
exceed performance while being more frugal in resource expenditure on
connections. While this result is relevant both for artificial (Amer and
Maul,, 2019) and biological networks (Clune et al.,, 2013), where these
connections carry non-trivial costs, it does overlook the relative cost of
producing and maintaining neurons in addition to their connections. Our
modular architecture incurs a heavy cost in terms of individual neural
processing units, scaling linearly with task complexity in proportion to each
module’s size. Indeed, the ability of modular networks to learn longer
effective timescales points to the importance of large but sparse networks in
solving a memory task. Future work could address the trade-off between these
two resources and their impact on performance.

Finally, even though we found a functional role for duplicative growth in
artificial networks, the observed benefits were modest. In evolutionary
contexts, relaxed selection (Deacon,, 2010) provides a mechanism by which
duplication can accommodate enhanced adaptation by fostering synergistic
interactions between duplicated modules. In our current setup, each duplicated
module is required to solve a new task, subjecting it to heavy selection
pressure that prevents it from exploring the entirety of the weight space.
Relaxing this constraint by introducing redundancy in modules that are
assigned the same task, could amplify the positive effects of duplication in
future research.

In conclusion, modular growth offers a promising avenue for training lighter
networks with more scalable properties. In the era of large models with
increasingly prohibitive costs, there is great interest in developing new
tools to compress models with fewer parameters (Cheung et al.,, 2019) and
fine-tune them to perform novel tasks (Hu et al.,, 2021). Our work, like other
evolutionary-inspired approaches (Akiba et al.,, 2024), can offer new
directions towards this end.

## Acknowledgements

MH, EG, TS and SK thank the International Max Planck Research School for
Intelligent Systems (IMPRS-IS) for their support. MH, EG, and CMW are
supported by the German Federal Ministry of Education and Research (BMBF):
Tübingen AI Center, FKZ: 01IS18039A. MH and CMW are funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy–EXC2064/1–390727645. TS is supported by Else Kröner
Medical Scientist Kolleg “ClinbrAIn: Artificial Intelligence for Clinical
Brain Research”.

## References

  * Akaike, (1974) Akaike, H. (1974).  A new look at the statistical model identification.  IEEE Transactions on Automatic Control, 19(6):716–723.
  * Akiba et al., (2024) Akiba, T., Shing, M., Tang, Y., Sun, Q., and Ha, D. (2024).  Evolutionary optimization of model merging recipes.
  * Amer and Maul, (2019) Amer, M. and Maul, T. (2019).  A review of modularization techniques in artificial neural networks.  Artificial Intelligence Review, 52:527–561.
  * Andrews, (1998) Andrews, J. H. (1998).  Bacteria as modular organisms.  Annual review of microbiology, 52(1):105–126.
  * Bassett and Bullmore, (2006) Bassett, D. S. and Bullmore, E. (2006).  Small-world brain networks.  The neuroscientist, 12(6):512–523.
  * Bathellier et al., (2008) Bathellier, B., Buhl, D. L., Accolla, R., and Carleton, A. (2008).  Dynamic ensemble odor coding in the mammalian olfactory bulb: Sensory information at different timescales.  Neuron, 57(4):586–598.
  * Chaudhuri et al., (2014) Chaudhuri, R., Bernacchia, A., and Wang, X.-J. (2014).  A diversity of localized timescales in network activity.  eLife, 3.
  * Cheung et al., (2019) Cheung, B., Terekhov, A., Chen, Y., Agrawal, P., and Olshausen, B. (2019).  Superposition of many models into one.
  * Chung et al., (2014) Chung, J., Gulcehre, C., Cho, K., and Bengio, Y. (2014).  Empirical evaluation of gated recurrent neural networks on sequence modeling.  arXiv preprint arXiv:1412.3555.
  * Clune et al., (2013) Clune, J., Mouret, J.-B., and Lipson, H. (2013).  The evolutionary origins of modularity.  Proceedings of the Royal Society b: Biological sciences, 280(1755):20122863.
  * Cosmides and Tooby, (1997) Cosmides, L. and Tooby, J. (1997).  The modular nature of human intelligence.  The origin and evolution of intelligence, pages 71–101.
  * Deacon, (2010) Deacon, T. W. (2010).  A role for relaxed selection in the evolution of the language capacity.  Proc. Natl. Acad. Sci. U. S. A., 107 Suppl 2(Suppl 2):9000–9006.
  * Douglas and Martin, (2004) Douglas, R. J. and Martin, K. A. (2004).  Neuronal circuits of the neocortex.  Annu. Rev. Neurosci., 27:419–451.
  * El Hihi and Bengio, (1995) El Hihi, S. and Bengio, Y. (1995).  Hierarchical recurrent neural networks for long-term dependencies.  In Proceedings of the 8th International Conference on Neural Information Processing Systems, NIPS’95, page 493–499, Cambridge, MA, USA. MIT Press.
  * Fang et al., (2021) Fang, W., Yu, Z., Chen, Y., Masquelier, T., Huang, T., and Tian, Y. (2021).  Incorporating Learnable Membrane Time Constant to Enhance Learning of Spiking Neural Networks.  arXiv:2007.05785 [cs].
  * Felleman and Van Essen, (1991) Felleman, D. J. and Van Essen, D. C. (1991).  Distributed hierarchical processing in the primate cerebral cortex.  Cerebral cortex (New York, NY: 1991), 1(1):1–47.
  * Fodor, (1983) Fodor, J. A. (1983).  The modularity of mind.  MIT press.
  * Gaier and Ha, (2019) Gaier, A. and Ha, D. (2019).  Weight agnostic neural networks.  Advances in neural information processing systems, 32.
  * Gallicchio and Micheli, (2021) Gallicchio, C. and Micheli, A. (2021).  Deep Reservoir Computing, pages 77–95.  Springer Singapore, Singapore.
  * Garcia-Fernàndez, (2005) Garcia-Fernàndez, J. (2005).  The genesis and evolution of homeobox gene clusters.  Nat. Rev. Genet., 6(12):881–892.
  * Geschwind and Rakic, (2013) Geschwind, D. H. and Rakic, P. (2013).  Cortical evolution: judge the brain by its cover.  Neuron, 80(3):633–647.
  * Giannakakis et al., (2023) Giannakakis, E., Vinogradov, O., Buendia, V., and Levina, A. (2023).  Recurrent connectivity structure controls the emergence of co-tuned excitation and inhibition.  bioRxiv.
  * Graves et al., (2013) Graves, A., Mohamed, A.-r., and Hinton, G. (2013).  Speech recognition with deep recurrent neural networks.  In 2013 IEEE international conference on acoustics, speech and signal processing, pages 6645–6649. Ieee.
  * Greengard, (2001) Greengard, P. (2001).  The Neurobiology of Slow Synaptic Transmission.  Science, 294(5544):1024–1030.  Publisher: American Association for the Advancement of Science.
  * Gruau and Whitley, (1993) Gruau, F. and Whitley, D. (1993).  Adding learning to the cellular development of neural networks: Evolution and the baldwin effect.  Evolutionary computation, 1(3):213–233.
  * Ha et al., (2016) Ha, D., Dai, A., and Le, Q. V. (2016).  Hypernetworks.  arXiv preprint arXiv:1609.09106.
  * Hagmann et al., (2008) Hagmann, P., Cammoun, L., Gigandet, X., Meuli, R., Honey, C. J., Wedeen, V. J., and Sporns, O. (2008).  Mapping the structural core of human cerebral cortex.  PLoS biology, 6(7):e159.
  * Happel and Murre, (1994) Happel, B. L. and Murre, J. M. (1994).  Design and evolution of modular neural network architectures.  Neural networks, 7(6-7):985–1004.
  * Hohil et al., (1999) Hohil, M. E., Liu, D., and Smith, S. H. (1999).  Solving the n-bit parity problem using neural networks.  Neural Networks, 12(9):1321–1323.
  * Hu et al., (2021) Hu, E. J., Shen, Y., Wallis, P., Allen-Zhu, Z., Li, Y., Wang, S., Wang, L., and Chen, W. (2021).  LoRA: Low-Rank adaptation of large language models.
  * Jaeger, (2002) Jaeger, H. (2002).  Tutorial on training recurrent neural networks, covering bppt, rtrl, ekf and the echo state network approach.
  * Khajehabdollahi et al., (2022) Khajehabdollahi, S., Prosi, J., Giannakakis, E., Martius, G., and Levina, A. (2022).  When to Be Critical? Performance and Evolvability in Different Regimes of Neural Ising Agents.  Artificial Life, pages 1–21.
  * Khajehabdollahi et al., (2024) Khajehabdollahi, S., Zeraati, R., Giannakakis, E., Schäfer, T. J., Martius, G., and Levina, A. (2024).  Emergent mechanisms for long timescales depend on training curriculum and affect performance in memory tasks.  In The Twelfth International Conference on Learning Representations.
  * Lake et al., (2015) Lake, B. M., Salakhutdinov, R., and Tenenbaum, J. B. (2015).  Human-level concept learning through probabilistic program induction.  Science, 350(6266):1332–1338.
  * LeCun et al., (2012) LeCun, Y. A., Bottou, L., Orr, G. B., and Müller, K.-R. (2012).  Efficient BackProp.  In Montavon, G., Orr, G. B., and Müller, K.-R., editors, Neural Networks: Tricks of the Trade: Second Edition, pages 9–48. Springer Berlin Heidelberg, Berlin, Heidelberg.
  * Litwin-Kumar and Doiron, (2012) Litwin-Kumar, A. and Doiron, B. (2012).  Slow dynamics and high variability in balanced cortical networks with clustered connections.  Nature Neuroscience, 15(11):1498–1505.  Number: 11 Publisher: Nature Publishing Group.
  * Liu et al., (2023) Liu, L., Yun, Z., Manubens-Gil, L., Chen, H., Xiong, F., Dong, H., Zeng, H., Hawrylycz, M., Ascoli, G. A., and Peng, H. (2023).  Neuronal connectivity as a determinant of cell types and subtypes.  Research Square.
  * Lui et al., (2011) Lui, J., Hansen, D., and Kriegstein, A. (2011).  Development and evolution of the human neocortex.  Cell, 146:18–36.
  * Manea et al., (2022) Manea, A. M., Zilverstand, A., Ugurbil, K., Heilbronner, S. R., and Zimmermann, J. (2022).  Intrinsic timescales as an organizational principle of neural processing across the whole rhesus macaque brain.  eLife, 11:e75540.
  * Morcos and Harvey, (2016) Morcos, A. and Harvey, C. (2016).  History-dependent variability in population dynamics during evidence accumulation in cortex.  Nature Neuroscience, 19.
  * Murray et al., (2014) Murray, J., Bernacchia, A., Freedman, D., Romo, R., Wallis, J., Cai, X., Padoa Schioppa, C., Pasternak, T., Seo, H., Lee, D., and Wang, X.-J. (2014).  A hierarchy of intrinsic timescales across primate cortex.  Nature neuroscience, 17.
  * Najarro et al., (2023) Najarro, E., Sudhakaran, S., and Risi, S. (2023).  Towards Self-Assembling Artificial Neural Networks through Neural Developmental Programs.  volume ALIFE 2023: Ghost in the Machine: Proceedings of the 2023 Artificial Life Conference, page 80.
  * Panzeri et al., (2010) Panzeri, S., Brunel, N., Logothetis, N. K., and Kayser, C. (2010).  Sensory neural codes using multiplexed temporal scales.  Trends in Neurosciences, 33(3):111–120.
  * Peng et al., (2021) Peng, H., Xie, P., Liu, L., Kuang, X., Wang, Y., Qu, L., Gong, H., Jiang, S., Li, A., Ruan, Z., et al. (2021).  Morphological diversity of single neurons in molecularly defined cell types.  Nature, 598(7879):174–181.
  * Perez-Nieves et al., (2021) Perez-Nieves, N., Leung, V. C. H., Dragotti, P. L., and Goodman, D. F. M. (2021).  Neural heterogeneity promotes robust learning.  Nature Communications, 12(1):5791.
  * Quax et al., (2020) Quax, S. C., D’Asaro, M., and van Gerven, M. A. J. (2020).  Adaptive time scales in recurrent neural networks.  Scientific Reports, 10(1):11360.  Number: 1 Publisher: Nature Publishing Group.
  * Rakic, (2009) Rakic, P. (2009).  Evolution of the neocortex: a perspective from developmental biology.  Nature Reviews Neuroscience, 10(10):724–735.
  * Rubino et al., (2023) Rubino, V., Hamidi, M., Dayan, P., and Wu, C. M. (2023).  Compositionality under time pressure.
  * Safavi et al., (2023) Safavi, S., Chalk, M., Logothetis, N., and Levina, A. (2023).  Signatures of criticality in efficient coding networks.  bioRxiv, pages 2023–02.
  * Schaub et al., (2015) Schaub, M. T., Billeh, Y. N., Anastassiou, C. A., Koch, C., and Barahona, M. (2015).  Emergence of Slow-Switching Assemblies in Structured Neuronal Networks.  PLOS Computational Biology, 11(7):e1004196.  Publisher: Public Library of Science.
  * Sharkey, (1996) Sharkey, A. J. C. (1996).  On combining artificial neural nets.  Connection science, 8(3-4):299–314.
  * Shen et al., (2023) Shen, L., Sun, Y., Yu, Z., Ding, L., Tian, X., and Tao, D. (2023).  On efficient training of large-scale deep learning models: A literature review.
  * Shepherd, (2011) Shepherd, G. M. (2011).  The microcircuit concept applied to cortical evolution: from three-layer to six-layer cortex.  Frontiers in neuroanatomy, 5:30.
  * Shi et al., (2023) Shi, Y.-L., Zeraati, R., Levina, A., and Engel, T. A. (2023).  Spatial and temporal correlations in neural networks with structured connectivity.  Physical Review Research, 5(1):013005.
  * Smith et al., (2023) Smith, J. T., Warrington, A., and Linderman, S. (2023).  Simplified state space layers for sequence modeling.  In The Eleventh International Conference on Learning Representations.
  * Sporns, (2013) Sporns, O. (2013).  Network attributes for segregation and integration in the human brain.  Current opinion in neurobiology, 23.
  * Sporns and Betzel, (2016) Sporns, O. and Betzel, R. F. (2016).  Modular brain networks.  Annual review of psychology, 67:613–640.
  * Sporns and Zwi, (2004) Sporns, O. and Zwi, J. D. (2004).  The small world of the cerebral cortex.  Neuroinformatics, 2:145–162.
  * Stanley, (2007) Stanley, K. O. (2007).  Compositional pattern producing networks: A novel abstraction of development.  Genetic programming and evolvable machines, 8:131–162.
  * Stork and Allen, (1992) Stork, D. G. and Allen, J. D. (1992).  How to solve the n-bit parity problem with two hidden units.  Neural networks, 5(6):923–926.
  * Tallec and Ollivier, (2018) Tallec, C. and Ollivier, Y. (2018).  Can recurrent neural networks warp time?  arXiv preprint arXiv:1804.11188.
  * Torres et al., (2021) Torres, J. F., Hadjout, D., Sebaa, A., Martínez-Álvarez, F., and Troncoso, A. (2021).  Deep learning for time series forecasting: a survey.  Big Data, 9(1):3–21.
  * Wu et al., (2020) Wu, D., Xia, S.-T., and Wang, Y. (2020).  Adversarial weight perturbation helps robust generalization.  Advances in Neural Information Processing Systems, 33:2958–2969.
  * Yin et al., (2020) Yin, B., Corradi, F., and Bohté, S. M. (2020).  Effective and Efficient Computation with Multiple-timescale Spiking Recurrent Neural Networks.  In International Conference on Neuromorphic Systems 2020, ICONS 2020, pages 1–8, New York, NY, USA. Association for Computing Machinery.
  * Yuan et al., (2023) Yuan, X., Savarese, P., and Maire, M. (2023).  Accelerated training via incrementally growing neural networks using variance transfer and learning rate adaptation.  Advances in Neural Information Processing Systems, 36.
  * Zeraati et al., (2022) Zeraati, R., Engel, T. A., and Levina, A. (2022).  A flexible Bayesian framework for unbiased estimation of timescales.  Nature Computational Science, 2(3):193–204.
  * Zeraati et al., (2023) Zeraati, R., Shi, Y.-L., Steinmetz, N. A., Gieselmann, M. A., Thiele, A., Moore, T., Levina, A., and Engel, T. A. (2023).  Intrinsic timescales in the visual cortex change with selective attention and reflect spatial connectivity.  Nature Communications, 14(1):1858.
  * Zhou et al., (2024) Zhou, H., Nagy, D. G., and Wu, C. M. (2024).  Harmonizing program induction with rate-distortion theory.  arXiv preprint arXiv:2405.05294.
  * Zierenberg et al., (2020) Zierenberg, J., Wilting, J., Priesemann, V., and Levina, A. (2020).  Tailored ensembles of neural networks optimize sensitivity to stimulus statistics.  Phys. Rev. Res., 2:013115.