DontShowMe[1][1]

# Supplementary materials for:
Long-term memory and synapse-like dynamics in two-dimensional nanofluidic
channels

P. Robin1,†, T. Emmerich1,†, A. Ismail,2,3,†,
A. Niguès1, Y. You2,3, G.-H. Nam2,3, A. Keerthi2,4,
A. Siria1, A.K. Geim2,3, B. Radha2,3,∗, L. Bocquet1,∗
1Laboratoire de Physique de l’Ecole normale Supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université Paris Cité, 75005 Paris, France
2 National Graphene Institute, The University of Manchester, Manchester, UK
3 Department of Physics and Astronomy, The University of Manchester,
Manchester, UK
∗To whom correspondence should be addressed;
4Department of Chemistry, University of Manchester, Manchester, UK
E-mail: lyderic.bocquet<EMAIL_ADDRESS>
† Equal contributions

Supplementary Text
Figs. S1 to S17
References (?, ?, ?, ?, ?, ?, ?)

## Supplementary Materials

###### Contents

  1. 1 Materials and methods
    1. 1.1 Nanofabrication of pristine MoS2 channels
    2. 1.2 Nanofabrication of activated  carbon channels
    3. 1.3 Current measurements
  2. 2 Bio-inspired algorithms
    1. 2.1 Long-term potentiation
    2. 2.2 Hebbian learing
  3. 3 Theory of the nanofluidic memristor
    1. 3.1 Unipolar memristors – Wien effect under 2D confinement
    2. 3.2 Bipolar memristors – Ionic rectification
    3. 3.3 Minimal model of a nanochannel with long-term memory
      1. 3.3.1 What does it mean for a nanochannel to have memory?
      2. 3.3.2 Adsorption-desorption model
      3. 3.3.3 Wien effect
      4. 3.3.4 A simple ansatz - Determination of the memory time
    4. 3.4 Mixed mechanisms and shape of the IV curve
  4. 4 Additional experimental data
    1. 4.1 Pristine MoS2 channels
    2. 4.2 Activated  carbon channels
  5. 5 Detailed list of supplementary figures

##  1 Materials and methods

###  1.1 Nanofabrication of pristine MoS2 channels

We fabricated the pristine MoS2 channels via van der Waals assembly following
the protocol reported in our previous work (?). Briefly, the process has two
major parts: I) the preparation of the top-spacer layers, II) the assembly of
the resulting top-spacer layers with the bottom layer to form tri-crystal
(top-spacer-bottom) stack. Graphene and MoS2 flakes were prepared by the
mechanical exfoliation of their bulk layered forms, Graphenium graphite and
natural MoS2 crystals (purchased from Manchester Nanomaterials). On a Si/SiO2
substrate, an exfoliated graphene flake with a specific thickness between 0.68
and 86 nm, was searched for using an optical microscope, and the flake
thickness was confirmed by atomic force microscopy (AFM) (Fig. S1A). Parallel
strips (width of $\sim$120 nm, spacing $\sim$150 nm) were made from this
graphene layer using e-beam lithography (EBL) and dry etching using oxygen
plasma (Fig. S1B). The spacing between the graphene strips would become the
channel width $w$ in the final device, while the thickness of the graphene
would become the channel height $h$. Fig. S1B (bottom panel) shows the AFM
height profile of a three-layer thin graphene spacer, with channel width of
$\sim$150 nm and a channel height of $\sim 1.2\pm 0.1$ nm. Following this, a
MoS2 crystal (thickness of 150 to 200 nm) was transferred on top of the
graphene spacer (Fig. S1C). This MoS2 layer served as the top wall of the
channel. The nanofluidic chip fabrication process began by drilling a
microhole ($\sim 3\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$\times
50\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$) in a SiNx membrane (thickness of
500 nm) on a silicon wafer, using photolithography and dry etching (mixture of
SF6 and CHF3 gases). A thin MoS2 layer (thickness of 20 to 40 nm) was then
transferred on top of the microhole on SiNx membrane to act as the bottom
layer of the channel, using polymethylmethacrylate (PMMA) based wet transfer
method (Fig. S1D). Then, the bottom layer was etched from the back of the SiNx
membrane via dry etching with CHF3 and O2 gases, to protrude the microhole
onto MoS2. Following this, the previously prepared stack of top MoS2-graphene
spacer was wet-transferred onto the bottom MoS2 layer on SiNx (Fig. S1E).
During this transfer, the spacer was oriented in such a way that the channels
are perpendicular to the rectangular microhole. Then, the graphene spacer was
dry etched (O2 plasma) from the back of SiNx membrane to further open the
microhole into the channels. At this stage, the resulting channels had
variable lengths determined by the shape of the top MoS2 flake, and any thin
edges of the top layer could lead to its sagging into the channels thus
blocking the channel entries. To address this, a Cr/Au strip (thickness 5
nm/70 nm, respectively) was deposited on the tri-crystal stack after
photolithography to open the channels (Fig. S1E). The gold layer aided the
device stability by minimizing the lifting of the MoS2 top layer during
measurements at high voltages. Moreover, the gold strip served as a mask to
create uniform and desired channel length for all channels across the device.
Regions of the tri-crystal stack not masked by the Au strip were etched away,
hence only the channel region underneath the gold strip remained (Fig. S1H).
Throughout the device fabrication process, after each flake transfer, the
device was placed in a furnace under H2:Ar (1:10) gas for annealing
(300$\mathrm{\SIUnitSymbolCelsius}$ for 3 hours and
400$\mathrm{\SIUnitSymbolCelsius}$ for 4 hours) to clean the polymer
contamination. The optical images of the final channel devices on SiNx are
shown in Fig. S1I, both in reflection mode and transmission mode, with the
channel length $L$ (from the microhole to the end of the Au strip) indicated.

###  1.2 Nanofabrication of activated  carbon channels

We briefly recall here the nanofabrication process of activated carbon
nanochannels. A more detailed description can be found in Ref. (?).
Bidimensional graphite crystals were obtained from commercially available
graphite (GRAPHENIUM) by mechanical exfoliation on a Si/SiO2 substrate using
cleanroom tape. A first graphite flake (‘bottom layer’) was pierced using
electron beam induced direct etching (EBIE) in water vapor inside a scanning
electron microscope, and then several trenches were dug from the hole using
the same technique (Fig. S2, step 1). A second graphite crystal (‘top layer’)
of roughly 50 nm in thickness was deposited above the hole and covered
partially the trenches, closing them to form channels (Fig. S2, step 2). This
first transfer was realized using the dry-transfer techniques with a droplet-
shaped polydimethylsiloxane (PDMS) stamp spin-coated with polypropylene
carbonate (PPC). Finally, this bi-layer heterostructure was transferred above
a Si/SiN window with a circular aperture in its middle, by making sure that
the hole in the bottom layer landed above the aperture (Fig. S2, step 3). This
second transfer was realized by wet transfer using a polymethylmethacrylate
(PMMA) sacrificial layer.

In typical cases, activated  carbon channels are
$5-10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ long,
$100\,$\mathrm{n}\mathrm{m}$$ wide, with a height ranging from 5 to
$15\,$\mathrm{n}\mathrm{m}$$.

###  1.3 Current measurements

Devices were placed into nanofluidic measurement cells separating two
reservoirs filled with electrolyte solutions of various salt concentrations
and ionic species (KCl, CaCl2, AlCl3, NiSO4, NaCl, LiCl). We used Ag/AgCl
electrodes to apply a potential drop across the channels and measure the
resulting ionic current. Our electrodes were connected to KEITHLEY amplifiers
(models 2636B and 2401). We used AC voltage of various frequencies ($0.1$ to
$200\,$\mathrm{m}\mathrm{H}\mathrm{z}$$) with a sampling rate of
$2\,$\mathrm{H}\mathrm{z}$$. Due to experimental technicalities, triangular
waveforms were used to probe pristine MoS2 channels and sinusoidal for
activated carbon channels; it was later checked, however, that our
nanochannels react identically to both types of waveforms.

Even in the thinnest channels, the measured current was found to be several
orders of magnitude higher than what was measured on a control system obtained
by following the “pristine” protocol, but without spacers to create the
channels. This shows that no leakage through the channel’s walls is not
possible.

##  2 Bio-inspired algorithms

In this section, we detail how we performed basic neuromorphic operations with
activated  carbon channels. In all cases, input voltage was generated via
MATLAB and exported as a text file, and then applied on the system by a
LABVIEW program with a sampling rate of $2\,$\mathrm{H}\mathrm{z}$$.

###  2.1 Long-term potentiation

Reversible long-term modification of a nanochannel’s conductance was achieved
by applying positive “write” pulses ($+1\,$\mathrm{V}$$ during
$10\,$\mathrm{s}$$) or negative “erase” pulses ($-1\,$\mathrm{V}$$ during
$10\,$\mathrm{s}$$). After each pulse, the conductance relaxation was tracked
by applying 10 weak “read” pulses ($\Delta V_{\text{read}}=+0.1\,$\mathrm{V}$$
during $5\,$\mathrm{s}$$) separated by $5\,$\mathrm{s}$$ (see Fig. S2A-B), and
computing for each of these pulses the conductance from:

$G=\frac{I}{\Delta V_{\text{read}}}$ (S1)

The conductance was found to stabilize after roughly two minutes of
relaxation.

We then checked that these modifications were incremental and reversible by
applying 30 write pulses followed by 30 erase pulses, see Fig. 5C from main
text. After each pulse, the system was let to relax for two minutes to let the
conductance stabilize.

###  2.2 Hebbian learing

The algorithm used to implement Hebbian learning is detailed in main text but
we recall it here for the sake of clarity. A computer generated a voltage time
series that emulated the behaviour of two neurons. This time series was then
applied on a nanofluidic channel. The activation of the presynaptic neuron A
was modeled by a weak positive voltage pulse. Whenever it activated, a flip-
flop mechanism was triggered, connecting the channel to a generator $E_{-}$
that applied negative voltage spikes. This behaviour lasted until the
postsynaptic neuron B activated and the system was branched on another
generator $E_{+}$ applying positive spikes instead. The opposite chain of
events occured if neuron B activated first: in that case, the channel first
received positive spikes from $E_{+}$, followed by negative spikes from
$E_{-}$ once neuron A activated. In both cases, the flip-flop reseted if a
given total amount of time $t_{0}$ passed since its activation, allowing the
process to start over.

To mimick Hebbian learning, we assumed that neuron B always activated with a
delay $\Delta t$ compared to neuron A. Note that $\Delta t$ could be negative,
if B actually fired before A. Then, for a given value of $\Delta t$, we
applied the above procedure 8 times, and measured the conductance change of
the channel at the end (see Fig. S4).

If A activates just before B, the system will be subject to a few negative
spikes followed by many positive ones, increasing its overall conductance (see
Fig. S4D). Likewise, If A activates just after B, the system receives a few
positive spikes and many negative spikes, and its conductance is lowered. On
the other hand, if $|\Delta t|$ is comparable to $t_{0}$, the channel will
receive an almost equal amount of positive and negative spikes, leaving its
conductance unchanged.

Additional data measured on different systems than the ones presented in Fig.
5-6 from main text are displayed in Fig. S5, demonstrating the robustness of
the observed phenomena.

In all cases, we prepared the channel in an intermediate conductance state so
that saturation to the state of maximum or minimum conductance is not a
problem during the learning process.

##  3 Theory of the nanofluidic memristor

In this section, we detail an analytical model of the nanofluidic memristor,
highlighting two distinct mechanisms. The most illustrative difference between
the two is the shape of the IV and GV curves under a periodic excitation. In
some cases, the GV curve displayed a self-crossing point (unipolar memristor,
proposed mechanism based on Wien effect, section 3.1); in others, it had the
shape of an open loop and the IV curve had a self-crossing point instead
(bipolar memristor, proposed mechanism based on ionic rectification, section
3.2).

If considered in a vacuum, both these effects only yield a memory time on the
timescales of milliseconds at best. In section 2.3, we will therefore provide
a minimal model explaining the emergence of long-term memory from the coupling
of surface processes to bulk transport between two reservoirs.

We start by recalling that a memristor is a resistor with a hysteretic
conductance. In terms of elementary electronics, it is decribed by a set of
two equations:

$\displaystyle I$ $\displaystyle=G\left[n(t)\right]\Delta V(t),$ (S2)
$\displaystyle\dot{n}$ $\displaystyle=f(n,\Delta V(t))$ (S3)

where $I$ is the electrical current flowing through the device under a time-
varying voltage drop $\Delta V(t)$, $G$ is the conductance that depends on an
internal parameter $n$, which can be seen as the system’s memory. It evolves
according to a dynamical equation (S3), where the dot represents the time
derivative.

The goal of this section is to detail, for the two mechanisms, what $n$
physically represents and to derive its evolution equation from the underlying
physics. In other words, to explain the memristor effect observed in confined
electrolytes, we need to first show that they possess several internal
conductance states, and then that they are able to retain such states over
long periods.

###  3.1 Unipolar memristors – Wien effect under 2D confinement

In this section, we discuss equation $(1)$ from main text, which describes
unipolar memristors (with a self-crossing GV curve). The full derivation can
be found in Ref. (?): herein, we only state and discuss the various results
for the sake of clarity. In this mechanism, $n$ represents the proportion of
ions which are able to move under an electric field. All other ions form
neutral (and therefore non-conducting) ion pairs, also refered to as Bjerrum
pairs (?). If the field is strong enough, it will tear some pairs apart,
increasing $n$ and the global conductance, in a process known as the (second)
Wien effect. If the field is turned off, pairs will eventually form again.
However, this pairing and unpairing process takes some time, allowing for a
memristor effect, akin to an electric arc in a discharge tube: an external
voltage is required to ionize the gas and make it conduct current, but the gas
will stay conducting for a short time after the voltage is removed (?).

The second Wien effect is a well-known phenomenon observed in weak
electrolytes. It was extensively studied by Onsager for bulk electrolytes (?),
resulting in an approximate law for the conductance $G$ under an external
field $E$:

$G(E)\simeq G(0)\left[1+\frac{\beta
e^{2}}{4\pi\epsilon\ell_{E}}+\dots\right]\simeq
G(0)\left[1+\frac{\ell_{B}}{\ell_{E}}+\dots\right]$ (S4)

with $\beta=1/k_{B}T$, $\ell_{E}=k_{B}T/eE$ a lengthscale defined by the
external field and $\ell_{B}$ the Bjerrum length. Here, $G(0)$ corresponds to
the conductance of ions that are already free at thermal equilibrium (i.e.
with no external field). This results in a sligthly non-linear IV curve:

$I=G(E)\Delta V\simeq G(0)\Delta V\left[1+\alpha\Delta V\right]$ (S5)

$\Delta V$ being the voltage drop associated with the field $E$. Without
detailing the (mathematically involved) exact derivation by Onsager, the
critical point is that ion pairs form according to a chemical equilibrium
given by:

$\dot{n}_{p}=\frac{1}{\tau_{a}}n_{f}^{2}-\frac{1}{\tau_{d}}n_{p}$ (S6)

with $n_{p}$ the proportion of ion pairs, $n_{f}$ the proportion of free ions,
$\tau_{a}$ and $\tau_{d}$ the ion pair association and dissociation times,
respectively. One has $n_{f}+n_{p}=1$ and the conductance is simply given by:

$G[n_{f}]=n_{f}G_{\infty}$ (S7)

where $G_{\infty}$ is the conductance in the fully dissociated case. Overall,
the system is able to remember the application of an electric field in its
recent past over a timescale $\tau_{a}$, and Bjerrum pairs can be used to
create memristors.

The above process cannot be achieved in bulk water, which fully dissociates
all commonly used salts. However, as noted by Ref. (?), ionic interactions are
greatly increased under confinement. This results in the formation of non-
conducting Bjerrum pairs, sometimes to the point that there are no ‘free’ ions
left at thermal equilibrium. In this case, the conductance vanishes in absence
of an electric field, $G(0)=0$, and it can be shown that:

$G(E)\propto\ell_{E}^{-\alpha}$ (S8)

with $\alpha>1$ scaling like the strength of ionic interactions. In typical
cases, one has $\alpha\sim 2$. In many experimental examples, however, a non-
zero conductance remains even in the absence of voltage. Taking into account
the fact that there may actually be a few free ions left at equilibrium, we
then write:

$G(E)=G_{0}+G_{1}\left(\frac{|E|}{E_{0}}\right)^{\alpha}$ (S9)

with $G_{0}\ll G_{1}$. This corresponds to equation $(1)$ of main text. The
number of ion pairs play the role of an internal state variable governing the
system’s conductance. It should be noted, however, that the value of $G_{0}$
was found to vary from device to device (even with similar dimensions),
although the general shape of the IV curves was preserved, see Fig. S11.

###  3.2 Bipolar memristors – Ionic rectification

The above mechanism, based on Wien effect, can only provide an explanation for
the unipolar memristor effect, as only the absolute magnitude of voltage, and
not its sign, governs the dynamics of ion pairs. Consequently, bipolar
memristors must have a different origin, as they display different conductance
states depending on the sign of applied voltage

This observation is remindful of ionic rectification and nanofluidic diode
(Fig. S6). This phenomenon typically appears in systems that are spatially
asymmetric, e.g. conical pores (geometrical asymmetry) (?), channels
connecting reservoirs of different salt concentrations (chemical asymmetry)
(?) or with an inhomogeneous surface charge (electrostatic asymmetry) (?). In
Ref. (?), the authors show that in all those cases the underlying mechanism is
the same: the local value of the conductivity is increased where surface
effects are the strongest (e.g., where the channel is the thinnest, salt
concentration the lowest or surface charge the highest). Globally, the system
then behaves like a PN junction: when current flow from the side of highest
conductivity to that of lowest conductivity, ions accumulate inside the
channel and conductance increases. Conversely, conductance is lowered if
current flows from the side of lowest conductivity.

For the sake of example, let us consider a 2D nanochannel with varying
thickness $h(x)$ or surface charge $\Sigma(x)$, $x$ being the direction along
the channel. The local value of salt concentration $c(x)$ will depend on
$h(x)$, $\Sigma(x)$ and the salt concentration in the reservoirs $c_{0}$. The
exact relation will depend on wether or not the system is in the Debye overlap
regime. If the Debye length is large compared to the channel’s thickness, then
conductivity reads:

$g(x)=2w\frac{e^{2}D}{k_{B}T}\sqrt{(ec_{0}h(x))^{2}+\Sigma(x)^{2}}$ (S10)

whereas in the opposite regime (Debye length small compared to $h$) it reads:

$g(x)=2w\frac{e^{2}D}{k_{B}T}(ec_{0}h(x)+|\Sigma(x)|)$ (S11)

with $w$ the channel’s width, $D$ the diffusion coefficient of ions (assumed
to be the same for both cations and anions), $e$ the elementary charge and $T$
temperature.  In particular, in the case of strong surface charges ($\Sigma\gg
hec$), conductance is governed by the surface charge alone:

$g\sim 2w\frac{e^{2}D}{k_{B}T}\Sigma$ (S12)

In addition, carbon-based nanochannels like activated carbon channels are
known for their high hydrodynamic slip-length (?). Due to the differential
friction of ions and water molecules on the channel’s walls, hydrodynamic
slippage yields a correction to the conductance that scales like:

$\Delta g\propto\frac{b\Sigma^{2}}{\eta}$ (S13)

where $b$ is the slip length and $\eta$ the viscosity of water. In what
follows, we will not discuss further the impact of slippage, but the two above
equations show how even thin channels can reach a conductance of the order of
$100\,$nS thanks to their strong surface charge.

In all cases, conductivity can always be expressed in terms of the Dukhin
number defined as:

$\mathrm{Du}(x)=\frac{\Sigma(x)}{ec_{0}h(x)}$ (S14)

and it increases with $\mathrm{Du}$ in all cases. Several regimes can now be
identified:

  1. 1.

$\mathrm{Du}\ll 1$ in the entire channel: bulk conduction dominates and the
system conducts linearly; there is no ionic rectification.

  2. 2.

$\mathrm{Du}\gg 1$ in the entire channel: surface conduction dominates, and
conductance is essentially given by counterions of the surface charge. As the
number of counterions is fixed, there can be no accumulation nor depletion of
ions in the channel, and no ionic rectification.

  3. 3.

$\mathrm{Du}\sim 1$ in at least parts of the channel, and is inhomogeneous:
ionic rectification can occur.

Quantitatively, Ref. (?) computes the global conductance of a channel with a
constant gradient of Dukhin number, in the absence of Debye overlap. The
obtained ratio of the maximum conductance of the channel (for a strongly
positive voltage) to its minimum conductance (for a strongly negative voltage)
is given by:

$\beta=\frac{G_{\text{max}}}{G_{\text{min}}}=\frac{1+2\mathrm{Du}_{\min}}{1+2\mathrm{Du}_{\max}}\frac{\mathrm{Du}_{\max}}{\mathrm{Du}_{\min}}$
(S15)

where $\mathrm{Du}_{\min}$ and $\mathrm{Du}_{\max}$ are the minimum and
maximum values of the Dukhin number, respectively.

The bipolar memristor effect was typically observed in activated  carbon
channels at low salt concentration. Typical parameters are $\Sigma\sim
0.1\,$\mathrm{C}\mathrm{/}\mathrm{m}^{2}$$, $c_{0}\sim 1\,$mM and $h\sim
10\,$nm, corresponding to $\mathrm{Du}\sim 10^{2}$. Yet, we observed ionic
rectification with $\beta\sim 1-5$, corresponding to $\mathrm{Du}_{\min}\sim
1$ in the above equation (assuming $\mathrm{Du}_{\max}\gg 1$). Such variations
of the Dukhin number over several orders of magnitude across the channel
cannot be accounted for by fluctuations of the channel’s height or surface
charge, especially in the case of pristine MoS2 channels, which are atomically
smooth throughout. Therefore, it seems more reasonable to attribute the ionic
rectification observed in 2D nanochannels to entry effects: the nanofluidic
device is deposited on a SiNx membrane, which is typically much thicker than
the device itself (500 nm versus $\sim$50 nm, respectively). The hole in the
bottom layer thus could act as an additional channel in series with the
nanofluidic device, with similar surface charge but much larger spatial
extension ($\sim 1\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$), yielding
$\mathrm{Du}\sim 1$.

This asymmetry in entry effects could be the source of an accumulation or a
depletion of ions inside the channel, depending on the sign of voltage.

###  3.3 Minimal model of a nanochannel with long-term memory

In this section, we complement the above analysis by showing how long-term
memory can emerge from interfacial processes in 2D nanochannels. We first show
how the rectification mechanism can be slowed down considerably if there is
adsorption of ions on the channel’s walls. This allows us to derive a minimal
model with analytical solutions, which can then be extended to more complex
cases, like that of the second Wien effect. However, we start by taking a step
back to analyze the nature of nanofluidic memory.

####  3.3.1 What does it mean for a nanochannel to have memory?

In this section, we discuss what we really measure in a memristor experiment,
and how we can quantify a system’s memory, with the help of a simplistic
advection-diffusion model.

We first focus on channels exhibiting ionic rectification. Rather than taking
the model presented in previous section, with geometrical asymmetry, we
simplify the problem even further. We assume that the channel has a
concentration $c$ of ions, and is connected to two reservoirs. We introduce an
ad hoc asymmetry by assuming that the reservoir on the left of the channel has
concentration $c_{L}$ and the reservoir on the right $c_{R}$, to mimick ionic
rectification while keeping mathematical complexity at a minimal level.
Lastly, we assume that there is an external forcing $f$ driving particles from
the left to the right. For the sake of simplicity, we work with units such
that the diffusion time across the channel, $L^{2}/D$, is equal to $1/2$.

In a continuous problem, the particle flux would be:

$j(x,t)=-\frac{1}{2}\partial_{x}c+f(t)c(x,t)=j_{\text{diff}}+j_{\text{adv}}$
(S16)

Note that the diffusive part $j_{\text{diff}}$ is an artefact: it exists
because we replaced a spatial asymmetry by a concentration gradient, resulting
in a global diffusive flux; it is irrelevant in actual experiments. The real
physical quantity of interest is therefore $j_{\text{adv}}$. In practice, we
measure this flux at the electrodes, i.e. at both ends of the channel, without
access to the full spatial dependence of $j$. In a time-independent problem
with $f(t)=f_{0}$ at all times, flux is conserved and this does not matter; we
obtain $j_{\text{measured}}=f_{0}\times\left\langle
c\right\rangle_{\infty,f_{0}}$, where $\left\langle
c\right\rangle_{\infty,f_{0}}$ is the spatial average of $c$ when the channel
is subject for a forcing $f_{0}$ for a very long time. In the quasistatic
limit, one would replace $f_{0}$ by a slow-varying forcing $f(t)$:

$j_{\text{quasistatic}}=f(t)\times\left\langle c\right\rangle_{\infty,f(t)}$
(S17)

where $\left\langle c\right\rangle_{\infty,f(t)}$ is now the spatial average
of $c$ when the channel is subject for a forcing “frozen” at a specific value
$f(t)$. A memory effect is any deviation from the above equation; it
translates the fact that the forcing is varying faster than the time needed to
equilibrate the system quasistatically. In other words, memory is stored in
the number of particles that can contribute to conduction.

However, when the system is not in the quasistatic limit, flux is not
conserved across the system. It is then hard to make exact sense of what is
being actually measured at the electrodes; in the following, we admit we may
still assume we measure a spatial average of the advection flux, but that this
average is now instantaneous:

$j_{\text{measured}}\simeq f(t)\left\langle
c(x,t)\right\rangle_{x}=f(t)\frac{1}{L}\int c(x,t)\,\mathrm{d}x$ (S18)

Assuming the forcing is sinusoidal $f(t)=f_{0}\cos\omega t$, one can then
compute the area of the hysteretic loop in the conductance-voltage curve:

$\mathcal{A}(\omega)=\left|\oint\left\langle
c(x,t)\right\rangle_{x}\,\mathrm{d}f\right|=\int_{0}^{2\pi/\omega}\omega
f_{0}\left\langle c(x,t)\right\rangle_{x}\sin\omega t\,\mathrm{d}t$ (S19)

The bigger $\mathcal{A}$ is, the more memory the system has of its recent
past. We thus define the memory timescale as $\tau_{m}=2\pi/\omega_{m}$ such
that $\mathcal{A}(\omega_{m})$ is maximum.

There is, however, no simple way to solve the advection-diffusion equation
under periodic forcing, even in 1D, so we simplify the problem further by
writing an approximate equation for the mean $\left\langle c\right\rangle$
only, see Fig. S7A:

$\dot{\left\langle c\right\rangle}+\left\langle
c\right\rangle\simeq\frac{c_{L}+c_{R}}{2}+f(t)\frac{c_{L}-c_{R}}{2}$ (S20)

This model is summed up in Fig. S7. In all what follows, we drop the
$\left\langle\cdot\right\rangle$ for the sake of simplicity. This yields:

$c(t)=\frac{c_{L}+c_{R}}{2}+f_{0}\frac{c_{L}-c_{R}}{2}\frac{\cos\omega
t+\omega\sin\omega t}{1+\omega^{2}}$ (S21)

so that the loop area is:

$\mathcal{A}(\omega)=f_{0}\frac{c_{L}-c_{R}}{2}\frac{\omega}{\omega^{2}+1}$
(S22)

This yields $\tau_{m}=2\pi$, or, in dimensional terms:

$\tau_{m}=\pi\frac{L^{2}}{D}$ (S23)

This is perfectly intuitive: since information is encoded in the particle
number, it cannot be retained for more than the diffusion time. However,
nanochannels have $L<10\,$\mathrm{\SIUnitSymbolMicro m}$$, meaning that the
diffusion time cannot exceed a second, contrary to what is observed in
experiments ($\tau_{m}\sim 1$ hour).

Before we move on to a slightly modified version of this model to account for
this, let us make the following remark. The quasistatic solution to the above
problem is:

$c_{\text{qs}}[f]=\frac{c_{L}+c_{R}}{2}+f\frac{c_{L}-c_{R}}{2}$ (S24)

The instantaneous solution we obtained can be rewritten into the following
form:

$c(t)=\int_{0}^{+\infty}c_{\text{qs}}[f(t-s)]\frac{e^{-s/\tau}}{\tau}\,\mathrm{d}s$
(S25)

where $\tau=1$, equal to $\tau_{m}$ up to a factor of order unity. This
equation will allow us to model more complex situation, where analytical
solutions do not exist.

####  3.3.2 Adsorption-desorption model

To complement the above model, which does not account for long memory times
observed in experiments, we consider the possibility of ion adsorbing on the
channel’s walls, and denote the number of adsorbed particles by $\sigma$.
Introducing $k$ and $\lambda$, the adsorption and desorption rates,
respectively, we obtain (see Fig. S7B):

$\displaystyle\dot{c}+c$
$\displaystyle=\frac{c_{L}+c_{R}}{2}+f(t)\frac{c_{L}-c_{R}}{2}-kc+\lambda\sigma$
(S26) $\displaystyle\dot{\sigma}$ $\displaystyle=kc-\lambda\sigma$ (S27)

It again can be solved analytically:

$\displaystyle c(t)$
$\displaystyle=\frac{c_{L}+c_{R}}{2}+f_{0}\frac{c_{L}-c_{R}}{2}\frac{\left[\lambda^{2}+\left(1+k\right)\omega^{2}\right]\cos\omega
t+\omega\left[\lambda\left(k+\lambda\right)+\omega^{2}\right]\sin\omega
t}{\lambda^{2}+\left[\left(1+k\right)^{2}+2k\lambda+\lambda^{2}\right]\omega^{2}+\omega^{4}}$
(S28) $\displaystyle\mathcal{A}(\omega)$
$\displaystyle=f_{0}\frac{c_{L}-c_{R}}{2}\frac{\omega\left[\lambda\left(k+\lambda\right)+\omega^{2}\right]}{\lambda^{2}+\left[\left(1+k\right)^{2}+2k\lambda+\lambda^{2}\right]\omega^{2}+\omega^{4}}$
(S29)

The memory timescale is then given by:

$\displaystyle\left\\{\lambda\left(k+\lambda\right)+3\omega^{2}\right\\}\left[\lambda^{2}+\left\\{\left(1+k\right)^{2}+2k\lambda+\lambda^{2}\right\\}\omega^{2}+\omega^{4}\right]$
(S30)
$\displaystyle...=\omega\left\\{\lambda\left(k+\lambda\right)+\omega^{2}\right\\}\left[2\left\\{\left(1+k\right)^{2}+k\lambda+\lambda^{2}\right\\}\omega+4\omega^{3}\right]\
$ (S31)

There is no closed-form solution to this last equation, but we can extract an
approximate solution when surface effects strongly dominate ($k\gg\lambda$ and
$k\gg 1$):

$\omega_{m}\sim\frac{\lambda}{k}\ll 1$ (S32)

In other words, the memory time reads in this case:

$\tau_{m}=2\pi\frac{\sigma_{\infty}}{c_{\infty}}\frac{L^{2}}{D}$ (S33)

where $\sigma_{\infty}$ and $c_{\infty}$ are the values of $\sigma$ and $c$ at
chemical equilibrium, respectively. The ratio $\sigma_{\infty}/{c_{\infty}}$
appearing in the above equation is the analogue of the Dukhin number (?),
which measures the importance of the surface charge of a channel with respect
to the bulk concentration in ions. This number can reach several hundreds, so
$\tau_{m}$ can be of the order of several minutes.

Qualitatively, the above equation can be recovered from a semi-quantitative
argument as follows. The memory time is given by the maximum time a particle
can stay within the channel. It will reach one of the reservoirs if left free
for more than $L^{2}/D$, by randomly diffusing along the channel. However,
every $1/k\ll 1$, the particle is adsorbed and stops moving, only liberated
after a time $1/\lambda$. Along its course through the channel, there will
therefore be $kL^{2}/D\gg 1$ such events of duration $1/\lambda$. The total
time spent inside the channel is the sum of the “travelling time” and “resting
time”:

$\tau_{m}\sim\frac{1}{\lambda}\frac{kL^{2}}{D}+\frac{L^{2}}{D}\sim\frac{k}{\lambda}\frac{L^{2}}{D}$
(S34)

The last approximation holds since we assumed that $k\gg\lambda$, such that
the total resting time is much larger than the travelling time.

To complement the above minimal model, we also computed a numerical solution
of the 1D advection-diffusion under periodic forcing and with adsorption,
yielding similar results (not reported here).

####  3.3.3 Wien effect

In the previous model, long-term memory emerges from a stop-and-go mechanism
of ions being adsorbed and desorbed. However, this process is more general
than the specific physics of adsorption, and we can write a similar system
taking into account a Wien effect mechanism:

$\displaystyle\dot{c}+c$
$\displaystyle=c_{0}-\frac{1}{\tau_{a}}c^{2}+\frac{1}{\tau_{d}[f(t)]}p$ (S35)
$\displaystyle\dot{p}$
$\displaystyle=\frac{1}{\tau_{a}}c^{2}-\frac{1}{\tau_{d}[f(t)]}p$ (S36)

Here, $c$ is the concentration of ions that can contribute to conduction (free
ions or polyelectrolytes), while $p$ represents the concentration of pairs.
This system is, however, non-linear in both $c$ and $f$, and as such admits no
analytical solution, but is qualitatively similar to the previous linear case.
Rather than looking for an approximate solution of this already simplistic
model, we instead use the ansatz derived above.

####  3.3.4 A simple ansatz - Determination of the memory time

Here, we detail how we fitted the curves presented in Figure 3 of main text.
As per sections 3.1 and 3.2, we have two models describing the conductance of
a 2D nanochannel under a constant electrical forcing. As shown by equation
(S25), this conductance becomes in the timevarying-case:

$G(t)=\int_{0}^{+\infty}G_{\text{qs}}[\Delta
V(t-s)]\frac{e^{-s/\tau}}{\tau}\,\mathrm{d}s$ (S37)

where $G_{\text{qs}}$ is the conductance in the stationary case (i.e. if the
voltage was ‘frozen’ at the value $\Delta V(t-s)$ for an infinite amount of
time). Depending on cases, we use the following expressions:

$\displaystyle G_{\text{qs}}(\Delta V)$ $\displaystyle=G_{1}\Delta
V^{\alpha},\quad\alpha=2-3\quad\text{(Wien effect)}$ (S38) $\displaystyle
G_{\text{qs}}(\Delta V)$ $\displaystyle=\begin{cases}G_{+},&\text{ if }\Delta
V>0,\\\ G_{-},&\text{ if }\Delta V<0\end{cases}\quad\text{(Ionic
rectification)}$ (S39)

We first use these models to extract the quasistatic limit of experimental
curves, and we then “turn on” memory effects by plugging the chosen model into
equation (S25). To assess the robustness of our approach, we can determine the
memory time through two methods:

  * •

By using $\tau$ as a fitting parameter in equation (S25).

  * •

By noticing that, upon correct renormalization (see below), the area of the
loop in the IV curve should take the form:

$\mathcal{A}^{*}(\omega)\simeq
K\frac{\omega\tau_{m}}{1+\omega^{2}\tau_{m}^{2}}$ (S40)

with a theoretical value $K=1$ in the adsorption-desorption model, and using
$\tau_{m}$ as a fitting parameter in this last equation. In the case of the
pairing-unpairing model, we cannot derive the exact expression of
$\mathcal{A}^{*}$, but we use the above expression as an approximative ansatz,
with $K$ as a fitting parameter.

Both quantities $\tau$ and $\tau_{m}$ can be interpreted as memory timescales,
being equal up to a factor $\pi$ in the minimal model. In practice, both
methods yield similar results ($\tau\sim\tau_{m}\sim 100\,$\mathrm{s}$$, see
Figure 4 from main text), but we believe the second one ($\tau_{m}$) to be
more robust as it is determined using more experimental data points obtained
using several frequencies.

Let us now detail how we normalize the loop area. Memristive devices cycle
between different conductance states. The largest loop that could possibly be
observed would be in the case where the device switches abruptly between the
lowest and the highest conductance states, $G_{\text{off}}$ and
$G_{\text{on}}$, whenever it reaches $\Delta V=\pm V_{0}$ (for bipolar
memristors) or $\Delta V=0$ and $|\Delta V|=V_{0}$ (for unipolar memristors).
The conductance-voltage curve then takes the shape of a rectangle of size
$2V_{0}\times(G_{\text{on}}-G_{\text{off}})$. The IV curve then takes the
shape of two triangles of total area
$V_{0}^{2}\times(G_{\text{on}}-G_{\text{off}})$, which we therefore use for
normalization. For each case, we determine $G_{\text{off}}$ and
$G_{\text{on}}$ graphically using data with the lowest frequency for a given
device and salt concentration. In the case of devices with mixed unipolar-
bipolar behavior, we determined both extremal conductance for positive and
negative voltage, and then used
$0.5V_{0}^{2}\times(G_{\text{on},+}-G_{\text{off},+})+0.5V_{0}^{2}\times(G_{\text{on},-}-G_{\text{off},-})$
as normalization. This normalization process is summed up in Fig. S8.

###  3.4 Mixed mechanisms and shape of the IV curve

The two mechanisms detailed above are not mutually exclusive, as some
nanofluidic systems can display both ionic rectification and allow the
formation of ionic pairs. In particular, pristine MoS2 were found to display
both behaviour depending on experimental conditions (Fig. S10). Additionally,
some of our experimental data do not fall in either of the two
phenomenologies: they displayed two crossing points in their IV curve, for
example.

To shed light on these ‘mixed’ cases, we designed a simple model of a
memristor that is intermediate between unipolar and bipolar cases:

$\displaystyle I(t)$ $\displaystyle=G(t)\times V(t)$ (S41) $\displaystyle
G(t)$
$\displaystyle=\int_{0}^{+\infty}G_{\text{qs}}(V(t-s))\frac{e^{-s/\tau}}{\tau}$
(S42) $\displaystyle G_{\text{qs}}(v)$
$\displaystyle=G_{\text{unipolar}}(1+\alpha v^{2})+G_{\text{bipolar}}(1+\beta
v)$ (S43)

with $I$ the ionic current, $V$ voltage, $G$ the instantaneous conductance,
$G_{\text{qs}}$ the quasistatic conductance and $\tau$, $G_{\text{unipolar}}$,
$G_{\text{bipolar}}$, $\alpha$ and $\beta$ some constants. For different
values of the parameters, this model can be purely unipolar, purely bipolar,
or intermediate between the two. Various examples of IV and GV curves for
different values of parameters are shown on Fig. S17, along an example of
experimental data displaying this mixed behaviour. We note that there exists a
“critical point” such that the GV curve has no crossing point; however, it
develops a cusp, which transforms into a crossing point when parameters are
varied.

##  4 Additional experimental data

###  4.1 Pristine MoS2 channels

We provide additional data to characterize the memristor effect in pristine
MoS2 channels, which is found to be robust when we vary salt concentration,
channel height and the chemical nature of the electrolyte.

On Fig. S9, we show that the memristor effect can be observed in pristine MoS2
channels regardless of the electrolyte used (KCl, NaCl, LiCl, CaCl2, NiSO4),
asserting its robustness. Despite variations in conductance, all these curves
display the same non-linear general shape remindful of the Wien effect.

On Fig. S10, we provide additional data for the memristor effect in bilayer
pristine MoS2 channels with potassium chloride solutions of various
concentrations. In particular, we notice that the system behaves like a
bipolar memristor at low salt concentration.

On Fig. S11, we show the evolution of loops in IV curve at fixed frequency and
salt concentration, but with increasing channel height. We find that the
effect is most visible in thin channels, and that the loop collapses to a
straight line in larger channels. This shows that memory effects can only be
observed if confinement is sufficiently strong so that ion-ion interactions
are enhanced, allowing the formation of ion pairs.

Lastly, Fig. S12 shows the determination of the memory time of three different
devices (and two salt concentrations) from the evolution of the loop area with
frequency. We obtain values between $50$ and $400\,$\mathrm{s}$$.

###  4.2 Activated  carbon channels

Additional data characterizing the influence of salt concentration on
memristive effects in activated  carbon channels are presented in Fig. S13.
They notably include raw data for Fig. 2A of main text. We observe that salt
concentration variations have little effect on the memory of thinner channels
(Fig. S13A-B), while this influence is more visible for devices above
$10\,$\mathrm{n}\mathrm{m}$$ in thickness. This notably shows that, in more
confined systems, interfacial processes are stronger than bulk effects.

In addition, the memristor effect can also be observed with salts other than
CaCl2, as shown in Fig. S14 for KCl and AlCl3.

We also provide additional data regarding the effect of voltage frequency
(Fig. S15), corresponding to Fig. 4C of main text. These data allow us to
compute the memory time of each activated system as the inverse of the
frequency such that the loop in the IV curve is the largest. We obtain values
spanning from $50$ to $400\,$\mathrm{s}$$, in a similar range as pristine MoS2
channels. Such variations may be explained by the variability of the surface
state of activated carbon channels following the etching in low pressure water
vapor.

##  5 Detailed list of supplementary figures

We summarize below the content of all supplementary figures:

  * •

Fig. S1: Nanofabrication of pristine MoS2 channels.

  * •

Fig. S2: Nanofabrication of activated  carbon channels.

  * •

Fig. S3: Details of the long-term potentation algorithm.

  * •

Fig. S4: Details of the Hebbian learning algorithm

  * •

Fig. S5: Additional data for long-term potentiation and Hebbian learning.

  * •

Fig. S6: Sources of ionic rectification in nanofluidics.

  * •

Fig. S7: Minimal model of nanofluidic memory.

  * •

Fig. S8: Normalization process of loop area in memristive GV curves.

  * •

Fig. S9: Additional data: pristine MoS2 channels with different salt types

  * •

Fig. S10: Additional data: pristine MoS2 channels with different salt
concentrations.

  * •

Fig. S11: Additional data: pristine MoS2 channels of different heights.

  * •

Fig. S12: Determination of the memory time of pristine MoS2 channels.

  * •

Fig. S13: Additional data: activated carbon channels with different salt
concentrations and different heights.

  * •

Fig. S14: Activated carbon channels with different salt types.

  * •

Fig. S15: Determination of the memory time of activated carbon channels.

  * •

Fig. S16: Influence of pH.

  * •

Fig. S17: ‘Mixed’ memristor types (theoretical model and experimental
example).

Fig. S1. Schematic flow-chart in cross-sectional view and corresponding
optical and AFM images. Step I: Preparation of the top-spacer layers on
silicon/silicon oxide (SiO2) wafer. A Graphene is mechanically exfoliated; the
thickness of this graphene flake will determine the height ($h$) of the
channel. Bottom panel: Optical image of a 3-layer graphene spacer. B
Patterning of the spacer using electron beam lithography (EBL) and etching
into parallel strips. Bottom panel: atomic force microscopy (AFM) image and
height profile of the patterned graphene-spacer ($h\sim$ 1.2 nm and $w\sim$
150 nm). C Transfer of MoS2 flake as a top layer over graphene spacer. Bottom
panel: Optical image of the top-spacer stack. Step II: Assembly of the tri-
crystal (top-spacer-bottom) stack on silicon/silicon nitride (SiNx) wafer. D
Transfer of a MoS2 flake onto SiNx membrane bearing a hole ($\sim 3\times
50\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$), to serve as bottom wall of the
channel. E Dry etching of the MoS2 bottom layer from the back of the SiNx. F
Transfer of the top (MoS2)-spacer (graphene) stack prepared in C over the
bottom MoS$2$ prepared in E. G Patterned gold (Au) deposition over the tri-
crystal stack. H With the Au strip as a mask to protect the underneath
channels, the surrounding regions are etched away. The gold strip thus
determines the channel length ($L$). I Optical images (left: reflection mode,
right: transmission mode) of the final channel device. The tricrystal device
is underneath the gold strip on the SiNx membrane. Scale bar in all images
represents $20\,$\mathrm{\SIUnitSymbolMicro m}$$. Fig. S2. Fabrication of
activated carbon channels. Step 1: Patterning of the bottom layer. Up: 2D side
view of a graphite bottom layer crystal after the patterning on a Si/SiO2
substrate. Middle: 3D view of the patterning process. A square-shaped hole and
four trenches connected to the hole are represented. Bottom layer graphite is
dark grey and the Si/SiO2 substrate is light pink. The electron flux,
represented as a green tip, enables selective removal of matter. Down: AFM
image of the bottom layer after etching. Scale bar represents
$1\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$. Step 2: Dry transfer of the top
layer. Up: 2D side view. A top layer crystal is added above the bottom layer.
Middle: 3D view of the device after the transfer of the top layer, represented
in glassy transparent grey. Down: SEM image of a device at that stage. Four
channels are visible in white. The bottom layer hole remains visible through
the top layer. Scale bar represents
$5\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$. Step 3: Wet transfer on the
Si/SiNx membrane. Up: 2D side view. Middle: 3D view, with the SiN membrane in
green. The circular aperture in the SiNx membrane is visible by transparency.
Down: Optical microscope image of a finished device. Scale bar represents
$10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$.  Fig. S3. Long-term
modification of a nanochannel conductance using voltage pulses. (A) A ‘write’
pulse ($+1\,$\mathrm{V}$$, $10\,$\mathrm{s}$$), followed by ten ‘read’ pulses
($+0.1\,$\mathrm{V}$$, $10\,$\mathrm{s}$$) to study the relaxation of the
conductance. (B) An ‘erase’ pulse ($-1\,$\mathrm{V}$$, $10\,$\mathrm{s}$$),
followed by ten ‘read’ pulses ($+0.1\,$\mathrm{V}$$, $10\,$\mathrm{s}$$). Fig.
S4. Algorithm for the implementation of Hebb’s rule with activated carbon
channels. (A) Voltage input emulating 8 successive activation of a pre-
synaptic and a post-synaptic neuron with a relative spike timing $\Delta
t=10\,$\mathrm{s}$$. The conductance is read before and after with low
amplitude ‘read’ pulses, highlighted by red circles. (B to E) Examples of
input voltage for various relative spike timings of the two neurons. Fig. S5.
Additional data for the implementation of neuromorphic computing with
activated carbon channels. The salt used is CaCl2,
$1\,$\mathrm{m}\mathrm{M}$$. (A) Evolution of the ionic current (red) under
voltage pulses of constant polarity (blue) (see Fig. 5A of main text). (B)
Conductance change following a positive voltage pulse of $1\,$\mathrm{V}$$ in
amplitude and a duration of $20\,$\mathrm{s}$$ (see Fig. 5B of main text). The
conductance is read with an alternating square voltage of $0.1\,$\mathrm{V}$$
in amplitude and a period of $20\,$\mathrm{s}$$. (C) Revsersible, long-term
modification of a nanochannel’s conductance. 30 write spikes (+1 V, 10 s) are
applied, followed by 30 erase spikes (-1 V, 10 s) which bring back the system
to its initial state. Between each spike, the conductance is let to stabilize
during two minutes and is then measured with a read pulse (0.1 V, 5 s), see
Fig. 5C of main text. (D) Conductance change after 8 successive activations of
the two neurons, in percentage of the initial conductance and as function of
the relative activation timing of the simulated neurons (see Fig. 6D of main
text). Fig. S6. Possible sources of ionic rectification. Nanochannel with
variable (A) height, (B) surface charge or (C) concentration gradient. In all
three cases, a gradient of Dukhin number is established across the channel,
resulting in a conductivity gradient and ionic rectification. Fig. S7. Minimal
model of nanofluidic memory. (A) Minimal model of purely diffusive memory. The
nanochannel is modeled by a single point that exchanges particles with two
reservoirs. The concentration difference between reservoirs plays the role of
geometrical asymmetry. Memory is stored in the concentration inside the
channel $c$, and the memory time is the diffusion timescale. (B) Same minimal
model, but with particle adsorption on the channel’s walls. Transport is now
limited by a stop-and-go mechanism of particles adsorbing and desorbing from
the walls, giving rise to a memory time orders of magnitude larger than
diffusion. (C) Memristor effect in the minimal model, as shown by the loop in
the IV curve, in dimensionless units. The blue curve corresponds to the model
described in panel A (section 3.3.1), and the red one in panel B (section
3.3.2). Memory effects are visible at low frequency only when surface effects
are taken into account. Fig. S8. Normalization of loop area. (A) Example of an
experimental conductance-voltage curve for a bipolar memristor (see Figure 1C
from main text), corresponding to an adsorption-desorption memory mechanism,
with the graphical determination of $G_{\text{on}}$ and $G_{\text{off}}$. (B)
Idealized loop in the conductance-voltage curve. Its area is used as a
normalization factor for a bipolar memristor. (C) Example of an experimental
conductance-voltage curve for a unipolar memristor (see Fig. 1B from main
text), corresponding to an Wien effect memory mechanism, with the graphical
determination of $G_{\text{on},\pm}$ and $G_{\text{off},\pm}$. (D) Idealized
loop in the conductance-voltage curve, defining the normalization factor for a
unipolar memristor. Note that the difference between $G_{\text{off},-}$ and
$G_{\text{off},+}$ is exaggerated compared to experimental data, to allow
easier visualization. Fig. S9. Evolution of the memristor effect with various
electrolytes (pristine channels). (A to E) Current-voltage characteristics of
an activated carbon channel with height under a voltage sweep of amplitude
$0.75\,$\mathrm{V}$$. Other parameters (salt type, salt concentration, channel
height and frequency) are specified on each panel. Fig. S10. Evolution of the
memristor effect with salt concentration (pristine channels). (A to D)
Current-voltage characteristics of a pristine MoS2 channel with height
$h=0.68\,$\mathrm{n}\mathrm{m}$$, filled with potassium chloride at various
concentrations, under a voltage sweep of amplitude $0.75\,$\mathrm{V}$$ and
frequency $3\,$\mathrm{m}\mathrm{H}\mathrm{z}$$. Orange curve indicates a
self-crossing loop, while blue curves do not self-intersect. Fig. S11.
Evolution of memristive effects with channel height (pristine channels).
Current-voltage characteristics of pristine MoS2 channels of different heights
filled with potassium chloride, under a voltage sweep (frequency
$f=3\,$\mathrm{m}\mathrm{H}\mathrm{z}$$). (A and B) Two different devices with
height $h=0.68\,$\mathrm{n}\mathrm{m}$$ (with salt concentration
$1\,$\mathrm{M}$$, voltage amplitude $0.75\,$\mathrm{V}$$). (C to E) Devices
with height $h=$ 2.8, 8.5 and $86\,$\mathrm{n}\mathrm{m}$$, respectively (salt
concentration $1\,$\mathrm{M}$$, voltage amplitude $0.75\,$\mathrm{V}$$ or
$1\,$\mathrm{V}$$). Fig. S12. Evolution of the memristor effect with voltage
frequency (pristine channels). Current-voltage characteristics of different
pristine MoS2 channels with channel heights $h=0.68\,$\mathrm{n}\mathrm{m}$$
(A and B) and $7\,$\mathrm{n}\mathrm{m}$$ (C and D), the electrolyte is 1M KCl
for A, C and 3M KCl for B, D; applied voltage is sinusoidal with frequency
ranging from 0.6 mHz to 200 mHz. The normalized loop area vs frequency of A is
presented in main Fig. 4C of main text. Insets represent the corresponding
normalized areas of the different current-voltage characteristics for that
device. The error bars represent the area variation between three successive
voltage sweeps. Fig. S13. Evolution of the memristor effect with salt
concentration (activated channels). (A to D) Current-voltage characteristics
of four different activated carbon channels, with CaCl2 and AC voltage
oscillating between $\pm 0.8\,$\mathrm{V}$$ at
$1\,$\mathrm{m}\mathrm{H}\mathrm{z}$$. The current is normalized by its
maximum value for each salt concentration. In each case, current is normalized
by its maximum absolute value, to allow easier comparison between different
datasets. Fig. S14. Evolution of the memristor effect with various
electrolytes (activated channels). (A to C) Current-voltage characteristics of
an activated carbon channel with height $h=13\,$\mathrm{n}\mathrm{m}$$, under
a voltage sweep of frequency $1\,$\mathrm{m}\mathrm{H}\mathrm{z}$$ and
amplitude $0.8\,$\mathrm{V}$$. Fig. S15. Evolution of the memristor effect
with voltage frequency (activated channels). (A to D) Current-voltage
characteristics of four different devices filled with 100 mM or 1 mM CaCl2 (as
indicated on each panel). Inset: Normalized area versus frequency. Squares are
experimental values and solid lines are theoretical model with memory time
parameter, $\tau_{m}$ equals to $50\,$\mathrm{s}$$ (B), $100\,$\mathrm{s}$$
(C) and $400\,$\mathrm{s}$$ (D). The normalized area vs frequency of the
$4\,$\mathrm{n}\mathrm{m}$$ device (A) is presented in main text Fig 4.C. Fig.
S16. Evolution of the memristor effect with pH. (A) Current-voltage
characteristics of a pristine MoS2 channel (height $0.68\,$nm) filled with 3 M
KCl under a voltage sweep at 10 mHz, for different pH values. (B) IV curve of
an activated carbon channel (height $10\,$nm) filled with AlCl3 1 mM under a
voltage sweep at 1 mHz, for different pH values. Fig. S17. Transition between
bipolar and unipolar behavior. (A to D) IV and GV (inset) curves for the
minimal model of a ”mixed” memristor (equations S47-S49). All units are
arbitrary. Parameters: $G_{\text{unipolar}}=3$, $G_{\text{bipolar}}=2$,
$\tau=1$, $V(t)=1.5\sin t$. (A) $\alpha=0.8$, $\beta=0.2$ (bipolar case). The
IV curve self-intersects once, but the GV curve does not. (B) $\alpha=0$,
$\beta=0.8$ (unipolar case). The IV curve does not self-intersect; the GV
curve self-intersects on the $y$ axis. (C) $\alpha=\beta=0.8$ (intermediate
case). The IV curve self-intersects twice and the GV curves self-intersects
outside of the $y$ axis. (D) $\alpha=0.8$, $\beta=0.45$ (criticality). The IV
curve self-intersects once and the GV curve develops a cusp. (E) Left panel:
Example of experimental IV curve showing a ”mixed” memristor behavior,
obtained for a pristine MoS2 channel (height $0.68\,$nm) filled with CaCl2
$3\,$M at voltage frequency $3\,$mHz. Red arrows indicate the two crossing
points. Inset is the corresponding GV curve. Right panel: Zoom-in on the two
crossing points of the IV curve.

## References

  * 1. N. Kavokine, R. R. Netz, L. Bocquet, Annu. Rev. Fluid Mech. 53, 377–410 (2021).
  * 2. L. Bocquet, Nat. Mater. 19, 254–256 (2020).
  * 3. M. Wang, Y. Hou, L. Yu, X. Hou, Nano Lett. 20, 6937–6946 (2020).
  * 4. K. Celebi, et al., Science 344, 289–292 (2014).
  * 5. A. Marcotte, T. Mouterde, A. Niguès, A. Siria, L. Bocquet, Nat. Mater. 19, 1057–1061 (2020).
  * 6. S. Garaj, et al., Nature 467, 190–193 (2010).
  * 7. C. A. Merchant, et al., Nano Lett. 10, 2915–2921 (2010).
  * 8. G. F. Schneider, et al., Nano Lett. 10, 3163–3167 (2010).
  * 9. J. Feng, et al., Nature 536, 197–200 (2016).
  * 10. E. Secchi, et al., Nature 537, 210–213 (2016).
  * 11. R. H. Tunuguntla, et al., Science 357, 792–796 (2017).
  * 12. N. Kavokine, M.-L. Bocquet, L. Bocquet, Nature 602, 84–90 (2022).
  * 13. B. Radha, et al., Nature 538, 222–225 (2016).
  * 14. A. Esfandiar, et al., Science 358, 511–513 (2017).
  * 15. T. Mouterde, et al., Nature 567, 87–90 (2019).
  * 16. T. Emmerich, et al., Nat. Mater. 21, 696–702 (2022).
  * 17. P. Robin, N. Kavokine, L. Bocquet, Science 373, 687–691 (2021).
  * 18. L. Chua, IEEE Trans. Circuit Theory 18, 507–519 (1971).
  * 19. D. B. Strukov, G. S. Snider, D. R. Stewart, R. S. Williams, Nature 453, 80–83 (2008).
  * 20. A. Sebastian, M. Le Gallo, R. Khaddam-Aljameh, E. Eleftheriou, Nat. Nanotechnol. 15, 529–544 (2020).
  * 21. R. Ge, et al., Nano Lett. 18, 434–441 (2018).
  * 22. Y. Bu, Z. Ahmed, L. Yobas, Analyst 144, 7168–7172 (2019).
  * 23. Q. Sheng, Y. Xie, J. Li, X. Wang, J. Xue, Chem. Commun. 53, 6125–6127 (2017).
  * 24. P. Zhang, et al., Nano Lett. 19, 4279–4286 (2019).
  * 25. B. Hille, Biophys. J. 22, 283–294 (1978).
  * 26. W. Gerstner, W. M. Kistler, Spiking Neuron Models: Single Neurons, Populations, Plasticity (Cambridge Univ. Press, 2002).
  * 27. Y. V. Pershin, M. Di Ventra, Adv. Phys. 60, 145–227 (2011).
  * 28. F. H. Van der Heyden, D. Stein, K. Besteman, S. G. Lemay, C. Dekker, Phys. Rev. Lett. 96, 224502 (2006).
  * 29. L. Fumagalli, et al., Science 360, 1339–1342 (2018).
  * 30. L. Bocquet, E. Charlaix, Chem. Soc. Rev. 39, 1073–1095 (2010).
  * 31. J. Comtet, et al., Nature Nanotechnol. 15, 598–604 (2020).
  * 32. S. Gravelle, R. R. Netz, L. Bocquet, Nano Lett. 19, 7265–7272 (2019).
  * 33. B. L. Werkhoven, et al., Phys. Rev. Lett. 120, 264502 (2018).
  * 34. P. Ober, et al., Nat. Commun. 12, 1–11 (2021).
  * 35. J.-X. Bao, E. R. Kandel, R. D. Hawkins, Science 275, 969–973 (1997).
  * 36. R. S. Zucker, W. G. Regehr, Annu. Rev. Physiol. 64, 355–405 (2002).
  * 37. M. Bear, B. Connors, M. A. Paradiso, Neuroscience: Exploring the Brain, Enhanced Edition: Exploring the Brain (Jones & Bartlett Learning, 2020).
  * 38. T. V. P. Bliss, G. L. Collingridge, Nature 361, 31–39 (1993).
  * 39. W. Gerstner, W. M. Kistler, R. Naud, L. Paninski, Neuronal dynamics: From single neurons to networks and models of cognition (Cambridge Univ. Press, 2014).
  * 40. G.-q. Bi, M.-m. Poo, J. Neurosci. 18, 10464–10472 (1998).
  * 41. P. Robin, et al., Experimental data for: Long-term memory and synapse-like dynamics in two-dimensional nanofluidic channels, Zenodo (2022); https://doi.org/10.5281/zenodo.7085645
  * 42. N. J. Bjerrum, Untersuchungen über Ionenassoziation (A. F. Høst,1926).
  * 43. T. Prodromakis, C. Tomazou, L. Chua, Nat. Mater. 11, 478–481 (2012).
  * 44. L. Onsager, J. Chem. Phys. 2, 599–615 (1934).
  * 45. Z. Siwy, et al., Europhys. Lett. 60, 349 (2002).
  * 46. C. B. Picallo, S. Gravelle, L. Joly, E. Charlaix, L. Bocquet, Phys. Rev. Lett. 111, 244501 (2013).
  * 47. R. Karnik, C. Duan, K. Castelino, H. Daiguji, A. Majumdar, Nano Lett. 7, 547–551 (2007).
  * 48. A. R. Poggioli, A. Siria, L. Bocquet, J. Phys. Chem. B 123, 1171–1185 (2019).


# Field theory of the Fermi function

Richard J. Hill<EMAIL_ADDRESS>University of Kentucky, Department of
Physics and Astronomy, Lexington, KY 40506, USA  Fermilab, Theoretical Physics
Department, Batavia, IL 60510, USA  Ryan Plestid<EMAIL_ADDRESS>University of Kentucky, Department of Physics and Astronomy, Lexington, KY
40506, USA  Fermilab, Theoretical Physics Department, Batavia, IL 60510, USA
Walter Burke Institute for Theoretical Physics, California Institute of
Technology, Pasadena, CA, 91125 USA

###### Abstract

The Fermi function $F(Z,E)$ accounts for QED corrections to beta decays that
are enhanced at either small electron velocity $\beta$ or large nuclear charge
$Z$. For precision applications, the Fermi function must be combined with
other radiative corrections and with scale- and scheme-dependent hadronic
matrix elements. We formulate the Fermi function as a field theory object and
present a new factorization formula for QED radiative corrections to beta
decays. We provide new results for the anomalous dimension of the
corresponding effective operator complete through three loops, and resum
perturbative logarithms and $\pi$-enhancements with renormalization group
methods. Our results are important for tests of fundamental physics with
precision beta decay and related processes.

††preprint: FERMILAB-PUB-23-453-T††preprint: CALT-TH/2023-029

Introduction. Many precision measurements and new physics searches involve
charged leptons interacting with nucleons or nuclei. Examples include neutrino
scattering to obtain fundamental neutrino parameters Nunokawa  _et al._
(2008); Abe  _et al._ (2014); Diwan  _et al._ (2016); Acero  _et al._ (2022);
Abe  _et al._ (2021); Abi  _et al._ (2020); muon-to-electron conversion to
search for charged lepton flavor violation de Gouvea and Vogel (2013);
Bernstein and Cooper (2013); Lee (2018); Bernstein (2019); and beta decay to
measure fundamental constants Bopp  _et al._ (1986); Ando  _et al._ (2004);
Darius  _et al._ (2017); Seng  _et al._ (2018, 2019); Fry  _et al._ (2019);
Czarnecki  _et al._ (2019); Hayen (2021); Seng  _et al._ (2020); Gorchtein and
Seng (2021); Hardy and Towner (2020); Gonzalez  _et al._ (2021); Shiells  _et
al._ (2021) and search for new physics Cirigliano  _et al._ (2013); Glick-
Magid  _et al._ (2017); González-Alonso  _et al._ (2019); Glick-Magid  _et
al._ (2022); Falkowski  _et al._ (2021); Brodeur  _et al._ (2023). It is
important to control radiative corrections to these processes Sirlin (1967);
Jaus (1972); Wilkinson (1982); Sirlin and Zucchini (1986); Czarnecki  _et al._
(2004); Cirigliano  _et al._ (2022). QED corrections are enhanced relative to
naive power counting in the fine structure constant $\alpha\approx 1/137$ for
large-$Z$ nuclei and for small-$\beta$ leptons ($Z$ denotes the nuclear
charge, and $\beta$ the lepton velocity). In this Letter we present new
results for long-distance QED corrections to beta decay Hill and Plestid
(2023a, b); Vander Griend  _et al._ (2023).

The Fermi function in beta decay is an extremely large QED effect Fermi
(1934), representing an $O(1)$ correction to beta decay rates Hayen  _et al._
(2018). It describes the enhancement (suppression) for negatively (positively)
charged leptons propagating in a nuclear Coulomb field. For a nuclear charge
$Z$ and electron energy $E$ it is defined as Fermi (1934); Hayen  _et al._
(2018):

$\displaystyle F(Z,E)$
$\displaystyle={2(1+\eta)\over|\Gamma(2\eta+1)|^{2}}[\Gamma(\eta+{\rm
i}\xi)]^{2}{\rm e}^{\pi\xi}(2pr)^{2(\eta-1)}\,,$ (1)

where $\eta\equiv\sqrt{1-(Z\alpha)^{2}}$, $\xi=Z\alpha/\beta$,
$p=\sqrt{E^{2}-m^{2}}$ and $m$ is the electron mass. The quantity $r$ denotes
a short distance regulator identified approximately as the nuclear size 111The
ambiguity in defining $r$ is often resolved by including an explicit nuclear
charge density distribution Hayen  _et al._ (2018). This procedure does not
address how to interface the Fermi function with EFT matching conditions, or
how to properly interface with higher-order radiative corrections.. Several
questions arise in the application of $F(Z,E)$ to physical processes: 1) What
is the scale $r^{-1}$ and how does it relate to conventional renormalization
in quantum field theory? 2) How can other radiative corrections be included
systematically? 3) What is the relation between the Fermi function with $Z=1$
and the radiative correction to neutron beta decay? Answers to these questions
impact the interpretation of modern beta decay experiments. For example,
corrections at order $\alpha(Z\alpha)^{2}\log(\Lambda_{\rm nuc.}/E)$ impact
beta decay rates of moderate $Z$ nuclei at permille level and must be included
at the current precision ($\sim 3\times 10^{-4}$) of $|V_{ud}|$ extractions
Hardy and Towner (2020). To answer these questions, we re-formulate the Fermi
function in effective field theory (EFT), and study its interplay with
subleading radiative corrections.

Factorization and all-orders matching. Factorization arises from the
separation of different energy scales involved in a physical process Collins
_et al._ (1989); Bodwin  _et al._ (1995); Bauer  _et al._ (2002). In a
sequence of EFTs, the components of a factorization formula are identified
with a corresponding sequence of matching coefficients, and a final low-energy
matrix element.

Consider the corrections to a tree-level contact interaction with a
relativistic electron in the final state. Ladder diagrams from a Coulomb
potential with source charge $+Ze$ correct the tree level amplitude, ${\cal
M}_{\rm tree}$, with explicit loop integrals given by

$\displaystyle\bar{u}(p){\cal M}$
$\displaystyle=\sum_{n=0}^{\infty}(Ze^{2})^{n}\int{\differential^{D}L_{1}\over(2\pi)^{D}}\int{\differential^{D}L_{2}\over(2\pi)^{D}}\cdots\int{\differential^{D}L_{n}\over(2\pi)^{D}}{1\over\mathbf{L}_{1}^{2}+\lambda^{2}}{1\over(\mathbf{L}_{1}-\mathbf{p})^{2}-\mathbf{p}^{2}-{\rm
i}0}$ $\displaystyle\hskip
32.52017pt\times{1\over(\mathbf{L}_{1}-\mathbf{L}_{2})^{2}+\lambda^{2}}{1\over(\mathbf{L}_{2}-\mathbf{p})^{2}-\mathbf{p}^{2}-{\rm
i}0}\cdots{1\over(\mathbf{L}_{n-1}-\mathbf{L}_{n})^{2}+\lambda^{2}}{1\over(\mathbf{L}_{n}-\mathbf{p})^{2}-\mathbf{p}^{2}-{\rm
i}0}$ $\displaystyle\hskip
130.08731pt\times\bar{u}(p)\gamma^{0}(\not{p}-\not{L}_{1}+m)\gamma^{0}(\not{p}-\not{L}_{2}+m)\cdots\gamma^{0}(\not{p}-\not{L}_{n}+m){\cal
M}_{\rm tree}\,.$ (2)

Integrals are evaluated in dimensional regularization with $D=3-2\epsilon$
dimensions, and we have included a photon mass, $\lambda$, to regulate
infrared divergences 222 Parameter $\lambda$ may be interpreted physically as
the scale of atomic screening, which influences decay rates at sub-leading
power Hayen  _et al._ (2018)..

In contrast with the analogous non-relativistic problem 333The classic
Sommerfeld enhancement is obtained in a perturbative expansion by replacing
the Dirac structures in the last line of (Field theory of the Fermi function)
with $2m$ in every propagator. This renders the problem UV finite. See Ref.
Hill and Plestid (2023b) for a more detailed discussion., the relativistic
expression (Field theory of the Fermi function) is UV divergent beginning at
two-loop order, indicating sensitivity to short-distance structure. The
factorization theorem reads Hill and Plestid (2023b)

$\displaystyle{\cal M}={\cal M}_{S}(\lambda/\mu_{S}){\cal
M}_{H}(p/\mu_{S},p/\mu_{H}){\cal M}_{\rm UV}(\Lambda/\mu_{H})\,,$ (3)

counting $p\sim m\sim E$ and where $\Lambda$ denotes the scale of hadronic and
nuclear structure. After ${\overline{\rm MS}}$ renormalization, to all orders
in $Z\alpha$, the soft function is given by ${\cal
M}_{S}=\exp\left(i\xi\log{\mu\over\lambda}\right)$ Yennie  _et al._ (1961);
Weinberg (1965). Our result for the hard function is new Hill and Plestid
(2023b), and is given (again to all orders in $Z\alpha$) by 444Covariant
expressions are obtained by replacing $\beta\rightarrow\sqrt{1-E^{2}/m^{2}}$,
$\gamma_{0}\rightarrow v_{\mu}\gamma^{\mu}$, and $E\rightarrow v_{\mu}p^{\mu}$
where $v_{\mu}$ is the reference vector introduced in Eqs. (6) and (7).

$\displaystyle\mathcal{M}_{H}$ $\displaystyle=e^{\tfrac{\pi}{2}\xi+{\rm
i}\phi_{H}}{2\Gamma(\eta-{\rm i}\xi)\over\Gamma(2\eta+1)}\sqrt{\eta-{\rm
i}\xi\over 1-{\rm i}\xi{m\over E}}\sqrt{E+\eta m\over E+m}{\sqrt{2\eta\over
1+\eta}}\left(2p\over e^{\gamma_{\rm
E}}\mu_{H}\right)^{\eta-1}\bigg{[}{1+\gamma^{0}\over 2}+{E+m\over E+\eta
m}\left(1-{\rm i}\xi{m\over E}\right){1-\gamma^{0}\over 2}\bigg{]}\,,$ (4)

where $\phi_{H}=\xi\left(\log{2p\over\mu_{S}}-\gamma_{\rm
E}\right)-(\eta-1){\pi\over 2}$, $\gamma^{0}$ is a Dirac matrix, and
$\gamma_{\rm E}\approx 0.577$ is the Euler constant.

The leading-in-$Z$ radiative correction to unpolarized observables from the
soft and hard functions is given by

$\displaystyle\big{\langle}|{\cal
M}_{H}|^{2}\big{\rangle}=F(Z,E)\big{|}_{r_{H}}\times\frac{4\eta}{(1+\eta)^{2}}\,,$
(5)

where we define $r_{H}^{-1}=\mu_{H}{\rm e}^{\gamma_{\rm E}}$. The angle
brackets denote contraction with lepton spinors, ${\cal M}_{H}\to\bar{e}{\cal
M}_{H}\gamma^{0}\nu_{L}$, sum over final state spins, and division by the same
expression in the absence of Coulomb corrections. Note that there is a finite
multiplicative correction relating the $\overline{\rm MS}$ hard function to
$F(Z,E)$.

Effective operators and anomalous dimension. The structure-dependent factor
${\cal M}_{\rm UV}$ appearing in Eq. (3) depends on the process of interest.
Important examples are beta decay transitions
$[A,Z]\to[A,Z+1]e^{-}\bar{\nu}_{e}$ or $[A,Z+1]\to[A,Z]e^{+}\nu_{e}$.
Superallowed beta decays are governed by an EFT consisting of QED for
electrons, and heavy charged scalar fields Caswell and Lepage (1986); Georgi
(1990); Manohar and Wise (2000); Paz (2015),

$\displaystyle{\cal L}_{\rm eff}$ $\displaystyle=-{\cal
C}(\phi_{v}^{[A,Z+1]})^{*}\phi_{v}^{[A,Z]}\bar{e}\not{v}(1-\gamma_{5})\nu_{e}+{\rm
H.c.}\,,$ (6)

where $\phi_{v}^{[A,Z]}$ denotes a heavy scalar with electric charge $Z$ whose
momentum fluctuations are expanded about $p^{\mu}=M_{[A,Z]}v^{\mu}$, with
$v^{\mu}=(1,0,0,0)$ in the nuclear rest frame. For neutron decay, the EFT
involves spin-1/2 heavy fields Caswell and Lepage (1986); Georgi (1990);
Manohar and Wise (2000); Paz (2015),

$\displaystyle{\cal L}_{\rm eff}$ $\displaystyle=-\bar{h}_{v}^{(p)}\left({\cal
C}_{V}\gamma^{\mu}+{\cal
C}_{A}\gamma^{\mu}\gamma_{5}\right)h_{v}^{(n)}\bar{e}\gamma^{\mu}(1-\gamma_{5})\nu_{e}+{\rm
H.c.}\,,$ (7)

where $h_{v}^{(p)}$ and $h_{v}^{(n)}$ denote spin-$1/2$ heavy fields with
electric charge 1 and 0, respectively. Matching to the EFT represented by Eqs.
(6) or (7), we identify the components of (3) in terms of operator
coefficients and matrix elements: ${\cal M}_{\rm UV}$ is proportional to (a
linear combination of) ${\cal C}_{i}$, while ${\cal M}_{H}$ and ${\cal M}_{S}$
give the hard and soft contributions to the EFT matrix element. In ${\cal
M}_{H}$, at each order in $\alpha$, the leading power of $Z$ is given by the
explicit expression (4).

We may proceed to analyze the renormalization group properties of weak-current
operators in the EFT. Radiative corrections enhanced by large logarithms,
$L\sim\log(\Lambda_{\rm nuc.}/E)$, are determined by the anomalous dimensions
of the operators in (6) and (7), which are spin-structure independent, i.e.,
$\gamma_{A}=\gamma_{V}=\gamma_{\mathcal{O}}$. Writing

$\begin{split}\gamma_{\cal O}={d\log{\cal C}\over
d\log\mu}&=\sum_{n=0}^{\infty}\sum_{i=0}^{n+1}\left(\alpha\over
4\pi\right)^{n+1}\gamma_{n}^{(i)}Z^{n+1-i}\\\
&\equiv\gamma^{(0)}(Z\alpha)+\alpha\gamma^{(1)}(Z\alpha)+\dots\,,\end{split}$
(8)

we note several interesting all-orders properties:

  * •

Powers of $Z$ greater than the power of $\alpha$ do not appear 555This can be
shown explicitly using heavy particle Feynman rules in Coulomb gauge, and
confirms an old theorem due to Sirlin and Bég Beg  _et al._ (1969)..

  * •

The leading series involving $(Z\alpha)^{n}$ sums to

$\displaystyle\gamma^{(0)}$ $\displaystyle=\sqrt{1-(Z\alpha)^{2}}-1\,.$ (9)

This result is obtained by differentiating Eq. (4) with respect to $\mu_{H}$.

  * •

At each order in perturbation theory, the leading and first subleading powers
of $Z$ are related 666More generally, the anomalous dimension at a given order
$n$ in $\alpha$ is a linear combination of powers $Z^{i}(1+Z)^{i}$, where
$2i\leq n$. This is a consequence of symmetries that emerge in the limit of a
massless lepton Hill and Plestid (2023a).,

$\displaystyle\gamma_{2n-1}^{(1)}=n\gamma_{2n-1}^{(0)}~{}~{},~{}~{}\gamma_{2n}^{(2)}=n\gamma_{2n}^{(1)}\quad(n\geq
1)\,.$ (10)

When $Z=0$, the problem reduces to a heavy-light current operator. Using our
new result for $\gamma_{2}^{(1)}=16\pi^{2}(6-\pi^{2}/3)$ Hill and Plestid
(2023a) and property (10), the complete result through three-loop order at
arbitrary $Z$ is

$\displaystyle\gamma_{\cal O}$ $\displaystyle={\alpha\over
4\pi}\gamma_{0}^{(1)}+\left({\alpha\over
4\pi}\right)^{2}\left[-8\pi^{2}Z(Z+1)+\gamma_{1}^{(2)}\right]$ (11)
$\displaystyle\hskip 21.68231pt+\left({\alpha\over
4\pi}\right)^{3}\left[16\pi^{2}\,Z(Z+1)\left(6-{\pi^{2}\over
3}\right)+\gamma_{2}^{(3)}\right]\,,$

where $\gamma_{n-1}^{(n)}$, $n=1,2,3$, are known from the heavy quark
literature 777The relevant results are $\gamma_{0}^{(1)}=-3$ Sirlin (1967) at
one-loop, $\gamma_{1}^{(2)}=-16\zeta_{2}+\frac{5}{2}+\frac{10}{3}n_{e}$ at
two-loops Ji and Musolf (1991), and
$\gamma_{2}^{(3)}=-80\zeta_{4}-36\zeta_{3}+64\zeta_{2}-\frac{37}{2}+n_{e}\left(-\frac{176}{3}\zeta_{3}+\frac{448}{9}\zeta_{2}+\frac{470}{9}\right)+\frac{140}{27}n_{e}^{2}$
at three-loops Chetyrkin and Grozin (2003) with $n_{e}=1$ for the number of
dynamical charged leptons. The application of $\gamma_{0}^{(1)}$,
$\gamma_{1}^{(2)}$ and $\gamma_{2}^{(3)}$ to neutron beta decay has been
discussed in Ref. Cirigliano  _et al._ (2023).. Our result for
$\gamma_{2}^{(1)}$ disagrees with Ref. Jaus (1972). Note that properties (9)
and (10) also determine the anomalous dimension at order $Z^{4}\alpha^{4}$ and
$Z^{3}\alpha^{4}$.

Figure 1: Radiative correction to the beta decay rate as a function of nuclear
charge, normalized to leading Fermi function. Black, red, blue and green
curves show results correct through resummed order $\alpha^{0}$,
$\alpha^{\frac{1}{2}}$, $\alpha$ and $\alpha^{\frac{3}{2}}$ respectively.
Illustrative values $E_{0}=2\,{\rm MeV}$, $\Delta_{0}=5\,{\rm MeV}$,
$\Lambda_{0}=100\,{\rm MeV}$ are used for the electron energy, nuclear mass
difference (which enters the one-loop matrix element Sirlin (1967)), and fixed
renormalization scale $\mu_{H}=\Lambda_{0}$. The width of the curves is given
by varying $m_{e}<\mu_{L}<\Delta_{0}$.

Renormalization group analysis. Consider the solution to the renormalization
group equation

$\displaystyle\differential\log{\cal
C}={\gamma(\alpha)\over\beta(\alpha)}\differential\alpha\,,$ (12)

where $\alpha$ is the ${\overline{\rm MS}}$ QED coupling (for one dynamical
electron flavor) and
$\beta=\differential\alpha/\differential\log\mu=-2\alpha[\beta_{0}\alpha/(4\pi)+\beta_{1}\alpha^{2}/(4\pi)^{2}+\dots]$
888Analytic results for the QED beta function are available through five loops
Baikov  _et al._ (2012); Herzog  _et al._ (2017). For our purposes we require
$\beta_{0}=-(4/3)n_{f}$ and $\beta_{1}=-4n_{f}$ (for $n_{f}=1$), cf. Ref. Hill
(2017). . Expanding $\gamma$ and $\beta$ in powers of $\alpha$ and $Z$, then
integrating, we obtain a systematic expansion for the ratio of the
renormalized operator coefficient at different scales,
$C(\mu_{H})/C(\mu_{L})$. Setting $\mu_{H}\sim\Lambda$ and $\mu_{L}\sim m$, we
thus resum large logarithms $\log(\Lambda/m)$. Let us consider several regimes
of $Z$:

  * •

Large $Z$ asymptotics.

Consider a large $Z$ nucleus, counting $\log^{2}(\Lambda/m)\sim\alpha^{-1}$
and $Z\sim\alpha^{-1}$. Through $O(\alpha^{1/2})$,

$\displaystyle\log\quantity({C(\mu_{L})\over C(\mu_{H})})=$
$\displaystyle\bigg{[}-\gamma^{(0)}(Z\alpha_{L})L\bigg{]}+\bigg{[}b_{0}\alpha_{L}L^{2}{(Z\alpha_{L})^{2}\over
2\sqrt{1-(Z\alpha_{L})^{2}}}\bigg{]}~{}$ (13)
$\displaystyle+\bigg{[}b_{0}^{2}\alpha_{L}^{2}L^{3}{(Z\alpha_{L})^{2}(3-2(Z\alpha_{L})^{2})\over
6(1-(Z\alpha_{L})^{2})^{\frac{3}{2}}}-\alpha_{L}L\gamma^{(1)}(Z\alpha_{L})\bigg{]}\,,$

where $\alpha_{H,L}\equiv\alpha(\mu_{H,L})$, $L=\log(\mu_{H}/\mu_{L})$, and
$b_{0}=-\beta_{0}/(2\pi)$. Resummation in $Z\alpha$ is an important effect for
large $Z$ nuclei Hayen  _et al._ (2018). Consider separately the terms in
$\gamma^{(1)}$ with odd and even powers of $(Z\alpha)$. Using property (10),

$\displaystyle\gamma^{(1)}_{\rm
odd}={\partial\over\partial(Z\alpha)}\gamma^{(0)}=\frac{-Z\alpha}{2\sqrt{1-(Z\alpha)^{2}}}\,.$
(14)

The corresponding decay rate corrections involve (less the known $Z\alpha^{2}$
correction)

$\begin{split}&\delta{|C(\mu_{L})|^{2}\over|C(\mu_{H})|^{2}}-\alpha(Z\alpha)\log{\Lambda\over
E}\\\ &=\alpha\log{\Lambda\over
E}\bigg{[}\frac{1}{2}(Z\alpha)^{3}+\frac{3}{8}(Z\alpha)^{5}+\dots\bigg{]}\,.\end{split}$
(15)

This exact result replaces the ansatz of Wilkinson Wilkinson (1997); Hayen
_et al._ (2018), with differences beginning at order $Z^{3}\alpha^{4}$. The
even series, $\gamma^{(1)}_{\rm even}$, is determined through three loop order
by Eq. (11).

  * •

Intermediate $Z$.

Consider a medium $Z$ nucleus, counting $\log^{2}(\Lambda/m)\sim
Z^{2}\sim\alpha^{-1}$. Through $O(\alpha^{3/2})$, the scale dependence is

$\displaystyle\log\quantity({C(\mu_{L})\over C(\mu_{H})})$
$\displaystyle={\gamma_{0}^{(1)}\over
2\beta_{0}}\Bigg{\\{}\bigg{[}\log{a_{H}\over
a_{L}}+{Z^{2}\gamma_{1}^{(0)}\over\gamma_{0}^{(1)}}\left(a_{H}-a_{L}\right)\bigg{]}+\bigg{[}{Z\gamma_{1}^{(1)}\over\gamma_{0}^{(1)}}(a_{H}-a_{L})\bigg{]}$
(16)
$\displaystyle\quad+\bigg{[}\left({\gamma_{1}^{(2)}\over\gamma_{0}^{(1)}}-{\beta_{1}\over\beta_{0}}\right)(a_{H}-a_{L})+\left({Z^{2}\gamma_{2}^{(1)}\over\gamma_{0}^{(1)}}-{\beta_{1}\over\beta_{0}}{Z^{2}\gamma_{1}^{(0)}\over\gamma_{0}^{(1)}}\right)\frac{1}{2}(a_{H}^{2}-a_{L}^{2})+{Z^{4}\gamma_{3}^{(0)}\over\gamma_{0}^{(1)}}\frac{1}{3}(a_{H}^{3}-a_{L}^{3})\bigg{]}\Bigg{\\}}\,,$

where $a_{H,L}=\alpha(\mu_{H,L})/(4\pi)$ and the square brackets account for
effects at order $\alpha^{\frac{1}{2}}$, $\alpha^{1}$, $\alpha^{\frac{3}{2}}$,
etc.

This “intermediate” regime of $Z$ applies to typical nuclei employed for
$V_{ud}$ determinations using super-allowed nuclear beta decay, with
$\alpha^{-\frac{1}{2}}\sim Z\sim\log(\Lambda^{2}/m^{2})\sim 10$. Achieving
permille precision thus demands proper treatment of terms through resummed
order $\alpha^{\frac{3}{2}}$ in Eq. (16). This result replaces (and disagrees
with) logarithmically enhanced contributions at order $Z^{2}\alpha^{3}$ in the
“heuristic estimate” of Sirlin and Zucchini 999The heuristic estimate Sirlin
and Zucchini (1986) uses results from Ref. Jaus (1972) and is defined relative
to the contributions contained in $F(Z,E)(1+\delta_{1}+\delta_{2})$, where
$F(Z,E)$ is the classical Fermi function, $\delta_{1}$ determines the order
$\alpha$ correction, and $\delta_{2}$ determines the $Z\alpha^{2}$ correction.
. Using our new result for $\gamma_{2}^{(1)}$ Hill and Plestid (2023a) we
investigate the convergence of perturbation theory in Fig. 1. Here we fix
$\mu_{H}$, and plot the product of $|C(\mu_{L})/C(\mu_{H})|^{2}$ and the
squared operator matrix element at $\mu_{L}$, varying $\mu_{L}$ as an estimate
of perturbative uncertainty 101010We include the full ${\cal O}(\alpha)$
matrix element Sirlin (1967) (translated to $\overline{\rm MS}$ in heavy
particle effective theory). For this comparison we include only the scale-
dependent logarithm, $\log(\mu/m_{e})$, in the $\alpha(Z\alpha)$ matrix
element Sirlin and Zucchini (1986).. We note that the results presented here
are in fact sufficient for a resummation valid through $O(\alpha^{2})$,
although for practical applications this would demand currently unknown
operator matrix elements.

  * •

Neutron beta decay.

Neutron beta decay corresponds to the case $Z=0$; we therefore define
$\gamma_{n-1}\equiv\gamma_{n-1}^{(n)}$. Again counting
$\log^{2}(\Lambda/m)\sim\alpha^{-1}$, the resummation is 111111 The two-loop
anomalous dimension $\gamma_{1}$ and other subleading two-loop corrections
were omitted in the leading-logarithm analysis of Ref. Czarnecki  _et al._
(2004). Reference Cirigliano  _et al._ (2023) included $\gamma_{1}$ in the RGE
evolution between nucleon and electron-mass scales.

$\begin{split}~{}~{}~{}~{}~{}&\log\quantity({C(\mu_{L})\over C(\mu_{H})})=\\\
&~{}~{}~{}{\gamma_{0}\over 2\beta_{0}}\Bigg{\\{}\log{a_{H}\over
a_{L}}+\left({\gamma_{1}\over\gamma_{0}}-{\beta_{1}\over\beta_{0}}\right)(a_{H}-a_{L})\Bigg{\\}}\,,\end{split}$
(17)

where the first term is of order $\alpha^{\frac{1}{2}}$, and the second term
is of order $\alpha^{\frac{3}{2}}$. The complete result, correct through order
$\alpha^{\frac{3}{2}}$, is obtained using (LABEL:eq:resummed_neutron) together
with the one-loop low-energy matrix element

Even after resumming logarithms in the ratio of hadronic and electron mass
scales, $\log(\Lambda/m)$, large coefficients remain in the perturbative
expansion of the hard matrix element. While the class of amplitudes summed in
the Fermi function are enhanced at small $\beta$ and large $Z$, neither limit
holds for neutron beta decay 121212At tree level, less than $0.1\%$ of the
total decay rate involves electron velocity $\beta<0.1$, and less than $10\%$
of the rate involves $\beta<0.5$. . The large coefficients can instead be
traced to an analytic continuation of the decay amplitude from spacelike to
timelike values of momentum transfers. The enhancements are systematically
resummed by renormalization of the hard factor ${\cal M}_{H}$ in the
factorization formula (3) from negative to positive values of $\mu_{S}^{2}$
Vander Griend  _et al._ (2023) (cf. Refs. Ahrens  _et al._ (2009a, b)). The
dependence of ${\cal M}_{H}$ on $\mu_{S}$ is governed (to all orders) by the
cusp anomalous dimension,

$\displaystyle{\cal M}_{H}(\mu_{S+})$ $\displaystyle={\rm e}^{{\rm
i}\alpha\phi_{\beta}}\exp[{\pi\alpha\over 2\beta}\bigg{]}{\cal
M}_{H}(\mu_{S-}^{2})\,,$ (18)

with $\mu_{S\pm}^{2}=\pm 4p^{2}-{\rm i}0$ and the phase
$\phi_{\beta}=\tfrac{1}{2}\big{(}-1+{1\over 2\beta}\log{1+\beta\over
1-\beta}\big{)}$. The hard contribution to the physical matrix element is
given on the right hand side of Eq. (18) in terms of an irrelevant phase
factor, an enhancement factor $\exp[\pi\alpha/(2\beta)]$, and
$\mathcal{M}_{H}$ evaluated at negative $\mu_{S}^{2}$. The latter matrix
element is free of $\pi$-enhancements. This analysis systematically resums
$\pi$-enhanced contributions, and does not rely on a non-relativistic
approximation. The result (18) differs from the nonrelativistic Fermi function
ansatz Wilkinson (1982); Cirigliano  _et al._ (2023) beginning at two loop
order.

Discussion. We have studied the all-orders factorization of long-distance QED
corrections in beta decay. This includes leading-$Z$ resummation traditionally
treated with the Fermi function, and subleading corrections in $\alpha$. Our
EFT analysis allows us to systematically resum large perturbative logarithms,
and to incorporate corrections that are suppressed by $1/Z$ or $E/\Lambda$.
New results include:

  1. 1.

New coefficients in the expansion of the anomalous dimension for beta decay
operators. We have computed the order $Z^{2}\alpha^{3}$ coefficient for the
first time 131313Existing calculations Jaus (1972) for order $Z^{2}\alpha^{3}$
terms included only a subset of diagrams and are therefore incomplete
estimates., and found a new symmetry linking leading-$Z$ and subleading-$Z$
terms in the perturbative expansion. Using our new result, and the existing
HQET literature, we show that the first unknown coefficient occurs at four
loops, at order $Z^{2}\alpha^{4}$ Hill and Plestid (2023a).

  2. 2.

New results for the large-$Z$ asymptotics of QED radiative corrections to beta
decay. We supply the infinite series of terms of order
$\alpha(Z\alpha)^{2n+1}\log(\Lambda/E)$, replacing Wilkinson’s ansatz, and
present a new result for the term of order
$\alpha(Z\alpha)^{2}\log(\Lambda/E)$, replacing Sirlin’s heuristic estimate.
We provide the EFT matrix element to all orders in $Z\alpha$ and clarify its
relation to the historically employed Fermi function Hill and Plestid (2023b).

  3. 3.

An all-orders resummation of “$\pi$-enhanced” terms in neutron beta decay,
replacing the Fermi function ansatz. This substantially improves the
convergence of perturbation theory, and is important for modern applications
to neutron beta decay Vander Griend  _et al._ (2023).

Each of these results has important phenomenological implications for ongoing
and near-term precision beta decay programs Brown  _et al._ (2018); Darius
_et al._ (2017); Fry  _et al._ (2019); Gonzalez  _et al._ (2021); Shidling
_et al._ (2014); Eibach  _et al._ (2015); Sternberg  _et al._ (2015); Gulyuz
_et al._ (2016); Fenker  _et al._ (2016); Long  _et al._ (2017); Fenker  _et
al._ (2018); Brodeur  _et al._ (2016); Shidling, P. D. and others (2018);
Valverde  _et al._ (2018); O’Malley  _et al._ (2020); Burdette  _et al._
(2020); Long  _et al._ (2020); Müller  _et al._ (2022); Long  _et al._ (2022).
Detailed computations are presented elsewhere Hill and Plestid (2023a, b);
Vander Griend  _et al._ (2023). Related work on new eikonal identities for
charged current processes is presented in Ref. Plestid (2023). The same
formalism applies to any situation where a charged lepton appears in a
reaction with a nucleus, provided its energy is small compared to the inverse
nuclear radius. Future work will address factorization at subleading power,
and investigate the impact on phenomenology including hadronic Ando  _et al._
(2004); Cirigliano  _et al._ (2022) and nuclear Seng  _et al._ (2019); Gennari
_et al._ (2023) matching uncertainties.

Acknowledgments. We thank Susan Gardner and Oleksandr Tomalak for useful
discussions regarding radiative corrections for beta decays, and Peter Vander
Griend for collaboration on Ref. Vander Griend  _et al._ (2023). RP thanks the
Institute for Nuclear Theory at the University of Washington for its kind
hospitality and stimulating research environment during program INT 23-1B.
This research was supported in part by the INT’s U.S. Department of Energy
grant No. DE-FG02-00ER41132. This work was supported by the U.S. Department of
Energy, Office of Science, Office of High Energy Physics, under Award DE-
SC0019095. Fermilab is operated by Fermi Research Alliance, LLC under Contract
No. DE-AC02-07CH11359 with the United States Department of Energy. RP is
supported by the U.S. Department of Energy, Office of Science, Office of High
Energy Physics, under Award Number DE-SC0011632 and by the Walter Burke
Institute for Theoretical Physics. RP acknowledges support from the U.S.
Department of Energy, Office of Science, Office of High Energy Physics, under
Award Number DE-SC0011632 and the Neutrino Theory Network Program Grant under
Award Number DEAC02-07CHI11359 and the US DOE under Award Number DE-SC0020250.
RJH gratefully acknowledges support from the Institute for Advanced Study,
where a part of this work was completed.

## References

  * Nunokawa  _et al._ (2008) H. Nunokawa, S. J. Parke, and J. W. F. Valle, Prog. Part. Nucl. Phys. 60, 338 (2008), arXiv:0710.0554 [hep-ph] .
  * Abe  _et al._ (2014) K. Abe _et al._ (Hyper-Kamiokande Working Group) (2014) arXiv:1412.4673 [physics.ins-det] .
  * Diwan  _et al._ (2016) M. V. Diwan, V. Galymov, X. Qian, and A. Rubbia, Ann. Rev. Nucl. Part. Sci. 66, 47 (2016), arXiv:1608.06237 [hep-ex] .
  * Acero  _et al._ (2022) M. A. Acero _et al._ (NOvA), Phys. Rev. D 106, 032004 (2022), arXiv:2108.08219 [hep-ex] .
  * Abe  _et al._ (2021) K. Abe _et al._ (T2K), Phys. Rev. D 103, 112008 (2021), arXiv:2101.03779 [hep-ex] .
  * Abi  _et al._ (2020) B. Abi _et al._ (DUNE), Eur. Phys. J. C 80, 978 (2020), arXiv:2006.16043 [hep-ex] .
  * de Gouvea and Vogel (2013) A. de Gouvea and P. Vogel, Prog. Part. Nucl. Phys. 71, 75 (2013), arXiv:1303.4097 [hep-ph] .
  * Bernstein and Cooper (2013) R. H. Bernstein and P. S. Cooper, Phys. Rept. 532, 27 (2013), arXiv:1307.5787 [hep-ex] .
  * Lee (2018) M. Lee, Front. in Phys. 6 (2018), 10.3389/fphy.2018.00133.
  * Bernstein (2019) R. H. Bernstein (Mu2e), Front. in Phys. 7, 1 (2019), arXiv:1901.11099 [physics.ins-det] .
  * Bopp  _et al._ (1986) P. Bopp, D. Dubbers, L. Hornig, E. Klemt, J. Last, H. Schutze, S. J. Freedman, and O. Scharpf, Phys. Rev. Lett. 56, 919 (1986), [Erratum: Phys.Rev.Lett. 57, 1192 (1986)].
  * Ando  _et al._ (2004) S. Ando, H. W. Fearing, V. P. Gudkov, K. Kubodera, F. Myhrer, S. Nakamura, and T. Sato, Phys. Lett. B 595, 250 (2004), arXiv:nucl-th/0402100 .
  * Darius  _et al._ (2017) G. Darius _et al._ , Phys. Rev. Lett. 119, 042502 (2017).
  * Seng  _et al._ (2018) C.-Y. Seng, M. Gorchtein, H. H. Patel, and M. J. Ramsey-Musolf, Phys. Rev. Lett. 121, 241804 (2018), arXiv:1807.10197 [hep-ph] .
  * Seng  _et al._ (2019) C. Y. Seng, M. Gorchtein, and M. J. Ramsey-Musolf, Phys. Rev. D 100, 013001 (2019), arXiv:1812.03352 [nucl-th] .
  * Fry  _et al._ (2019) J. Fry _et al._ , EPJ Web Conf. 219, 04002 (2019), arXiv:1811.10047 [nucl-ex] .
  * Czarnecki  _et al._ (2019) A. Czarnecki, W. J. Marciano, and A. Sirlin, Phys. Rev. D 100, 073008 (2019), arXiv:1907.06737 [hep-ph] .
  * Hayen (2021) L. Hayen, Phys. Rev. D 103, 113001 (2021), arXiv:2010.07262 [hep-ph] .
  * Seng  _et al._ (2020) C.-Y. Seng, X. Feng, M. Gorchtein, and L.-C. Jin, Phys. Rev. D 101, 111301 (2020), arXiv:2003.11264 [hep-ph] .
  * Gorchtein and Seng (2021) M. Gorchtein and C.-Y. Seng, JHEP 10, 053 (2021), arXiv:2106.09185 [hep-ph] .
  * Hardy and Towner (2020) J. C. Hardy and I. S. Towner, Phys. Rev. C 102, 045501 (2020).
  * Gonzalez  _et al._ (2021) F. M. Gonzalez _et al._ (UCN$\tau$), Phys. Rev. Lett. 127, 162501 (2021), arXiv:2106.10375 [nucl-ex] .
  * Shiells  _et al._ (2021) K. Shiells, P. G. Blunden, and W. Melnitchouk, Phys. Rev. D 104, 033003 (2021), arXiv:2012.01580 [hep-ph] .
  * Cirigliano  _et al._ (2013) V. Cirigliano, S. Gardner, and B. Holstein, Prog. Part. Nucl. Phys. 71, 93 (2013), arXiv:1303.6953 [hep-ph] .
  * Glick-Magid  _et al._ (2017) A. Glick-Magid, Y. Mishnayot, I. Mukul, M. Hass, S. Vaintraub, G. Ron, and D. Gazit, Phys. Lett. B 767, 285 (2017), arXiv:1609.03268 [nucl-ex] .
  * González-Alonso  _et al._ (2019) M. González-Alonso, O. Naviliat-Cuncic, and N. Severijns, Prog. Part. Nucl. Phys. 104, 165 (2019), arXiv:1803.08732 [hep-ph] .
  * Glick-Magid  _et al._ (2022) A. Glick-Magid, C. Forssén, D. Gazda, D. Gazit, P. Gysbers, and P. Navrátil, Phys. Lett. B 832, 137259 (2022), arXiv:2107.10212 [nucl-th] .
  * Falkowski  _et al._ (2021) A. Falkowski, M. González-Alonso, A. Palavrić, and A. Rodríguez-Sánchez, (2021), arXiv:2112.07688 [hep-ph] .
  * Brodeur  _et al._ (2023) M. Brodeur _et al._ (2023) arXiv:2301.03975 [nucl-ex] .
  * Sirlin (1967) A. Sirlin, Phys. Rev. 164, 1767 (1967).
  * Jaus (1972) W. Jaus, Phys. Lett. B 40, 616 (1972).
  * Wilkinson (1982) D. H. Wilkinson, Nucl. Phys. A 377, 474 (1982).
  * Sirlin and Zucchini (1986) A. Sirlin and R. Zucchini, Phys. Rev. Lett. 57, 1994 (1986).
  * Czarnecki  _et al._ (2004) A. Czarnecki, W. J. Marciano, and A. Sirlin, Phys. Rev. D70, 093006 (2004), arXiv:hep-ph/0406324 [hep-ph] .
  * Cirigliano  _et al._ (2022) V. Cirigliano, J. de Vries, L. Hayen, E. Mereghetti, and A. Walker-Loud, Phys. Rev. Lett. 129, 121801 (2022), arXiv:2202.10439 [nucl-th] .
  * Hill and Plestid (2023a) R. J. Hill and R. Plestid, “Renormalization of beta decay at three and four loops,” In preparation (2023a).
  * Hill and Plestid (2023b) R. J. Hill and R. Plestid, “All orders factorization and the Coulomb problem,” To appear (2023b).
  * Vander Griend  _et al._ (2023) P. Vander Griend, R. J. Hill, and R. Plestid, “Resummation of infrared $\pi$-enhancements and neutron beta decay,” In preparation (2023).
  * Fermi (1934) E. Fermi, Z. Phys. 88, 161 (1934).
  * Hayen  _et al._ (2018) L. Hayen, N. Severijns, K. Bodek, D. Rozpedzik, and X. Mougeot, Rev. Mod. Phys. 90, 015008 (2018), arXiv:1709.07530 [nucl-th] .
  * Note (1) The ambiguity in defining $r$ is often resolved by including an explicit nuclear charge density distribution Hayen  _et al._ (2018). This procedure does not address how to interface the Fermi function with EFT matching conditions, or how to properly interface with higher-order radiative corrections.
  * Collins  _et al._ (1989) J. C. Collins, D. E. Soper, and G. F. Sterman, Adv. Ser. Direct. High Energy Phys. 5, 1 (1989), arXiv:hep-ph/0409313 .
  * Bodwin  _et al._ (1995) G. T. Bodwin, E. Braaten, and G. P. Lepage, Phys. Rev. D 51, 1125 (1995), [Erratum: Phys.Rev.D 55, 5853 (1997)], arXiv:hep-ph/9407339 .
  * Bauer  _et al._ (2002) C. W. Bauer, S. Fleming, D. Pirjol, I. Z. Rothstein, and I. W. Stewart, Phys. Rev. D 66, 014017 (2002), arXiv:hep-ph/0202088 .
  * Note (2) Parameter $\lambda$ may be interpreted physically as the scale of atomic screening, which influences decay rates at sub-leading power Hayen  _et al._ (2018).
  * Note (3) The classic Sommerfeld enhancement is obtained in a perturbative expansion by replacing the Dirac structures in the last line of (Field theory of the Fermi function) with $2m$ in every propagator. This renders the problem UV finite. See Ref. Hill and Plestid (2023b) for a more detailed discussion.
  * Yennie  _et al._ (1961) D. Yennie, S. C. Frautschi, and H. Suura, Annals Phys. 13, 379 (1961).
  * Weinberg (1965) S. Weinberg, Phys. Rev. 140, B516 (1965).
  * Note (4) Covariant expressions are obtained by replacing $\beta\rightarrow\sqrt{1-E^{2}/m^{2}}$, $\gamma_{0}\rightarrow v_{\mu}\gamma^{\mu}$, and $E\rightarrow v_{\mu}p^{\mu}$ where $v_{\mu}$ is the reference vector introduced in Eqs. (6) and (7).
  * Caswell and Lepage (1986) W. E. Caswell and G. P. Lepage, Phys. Lett. B 167, 437 (1986).
  * Georgi (1990) H. Georgi, Phys. Lett. B 240, 447 (1990).
  * Manohar and Wise (2000) A. V. Manohar and M. B. Wise,  _Heavy quark physics_ , Vol. 10 (2000).
  * Paz (2015) G. Paz, Mod. Phys. Lett. A 30, 1550128 (2015), arXiv:1503.07216 [hep-ph] .
  * Note (5) This can be shown explicitly using heavy particle Feynman rules in Coulomb gauge, and confirms an old theorem due to Sirlin and Bég Beg  _et al._ (1969).
  * Note (6) More generally, the anomalous dimension at a given order $n$ in $\alpha$ is a linear combination of powers $Z^{i}(1+Z)^{i}$, where $2i\leq n$. This is a consequence of symmetries that emerge in the limit of a massless lepton Hill and Plestid (2023a).
  * Note (7) The relevant results are $\gamma_{0}^{(1)}=-3$ Sirlin (1967) at one-loop, $\gamma_{1}^{(2)}=-16\zeta_{2}+\frac{5}{2}+\frac{10}{3}n_{e}$ at two-loops Ji and Musolf (1991), and $\gamma_{2}^{(3)}=-80\zeta_{4}-36\zeta_{3}+64\zeta_{2}-\frac{37}{2}+n_{e}\left(-\frac{176}{3}\zeta_{3}+\frac{448}{9}\zeta_{2}+\frac{470}{9}\right)+\frac{140}{27}n_{e}^{2}$ at three-loops Chetyrkin and Grozin (2003) with $n_{e}=1$ for the number of dynamical charged leptons. The application of $\gamma_{0}^{(1)}$, $\gamma_{1}^{(2)}$ and $\gamma_{2}^{(3)}$ to neutron beta decay has been discussed in Ref. Cirigliano  _et al._ (2023).
  * Note (8) Analytic results for the QED beta function are available through five loops Baikov  _et al._ (2012); Herzog  _et al._ (2017). For our purposes we require $\beta_{0}=-(4/3)n_{f}$ and $\beta_{1}=-4n_{f}$ (for $n_{f}=1$), cf. Ref. Hill (2017).
  * Wilkinson (1997) D. Wilkinson, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 401, 275 (1997).
  * Note (9) The heuristic estimate Sirlin and Zucchini (1986) uses results from Ref. Jaus (1972) and is defined relative to the contributions contained in $F(Z,E)(1+\delta_{1}+\delta_{2})$, where $F(Z,E)$ is the classical Fermi function, $\delta_{1}$ determines the order $\alpha$ correction, and $\delta_{2}$ determines the $Z\alpha^{2}$ correction.
  * Note (10) We include the full ${\cal O}(\alpha)$ matrix element Sirlin (1967) (translated to $\overline{\rm MS}$ in heavy particle effective theory). For this comparison we include only the scale-dependent logarithm, $\log(\mu/m_{e})$, in the $\alpha(Z\alpha)$ matrix element Sirlin and Zucchini (1986).
  * Note (11) The two-loop anomalous dimension $\gamma_{1}$ and other subleading two-loop corrections were omitted in the leading-logarithm analysis of Ref. Czarnecki  _et al._ (2004). Reference Cirigliano  _et al._ (2023) included $\gamma_{1}$ in the RGE evolution between nucleon and electron-mass scales.
  * Note (12) At tree level, less than $0.1\%$ of the total decay rate involves electron velocity $\beta<0.1$, and less than $10\%$ of the rate involves $\beta<0.5$.
  * Ahrens  _et al._ (2009a) V. Ahrens, T. Becher, M. Neubert, and L. L. Yang, Phys. Rev. D 79, 033013 (2009a), arXiv:0808.3008 [hep-ph] .
  * Ahrens  _et al._ (2009b) V. Ahrens, T. Becher, M. Neubert, and L. L. Yang, Eur. Phys. J. C 62, 333 (2009b), arXiv:0809.4283 [hep-ph] .
  * Cirigliano  _et al._ (2023) V. Cirigliano, W. Dekens, E. Mereghetti, and O. Tomalak, (2023), arXiv:2306.03138 [hep-ph] .
  * Note (13) Existing calculations Jaus (1972) for order $Z^{2}\alpha^{3}$ terms included only a subset of diagrams and are therefore incomplete estimates.
  * Brown  _et al._ (2018) M. A. P. Brown _et al._ (UCNA), Phys. Rev. C 97, 035505 (2018), arXiv:1712.00884 [nucl-ex] .
  * Shidling  _et al._ (2014) P. D. Shidling _et al._ , Phys. Rev. C 90, 032501 (2014), arXiv:1407.1742 [nucl-ex] .
  * Eibach  _et al._ (2015) M. Eibach _et al._ , Phys. Rev. C 92, 045502 (2015).
  * Sternberg  _et al._ (2015) M. G. Sternberg _et al._ , Phys. Rev. Lett. 115, 182501 (2015).
  * Gulyuz  _et al._ (2016) K. Gulyuz _et al._ , Phys. Rev. Lett. 116, 012501 (2016).
  * Fenker  _et al._ (2016) B. Fenker _et al._ , New J. Phys. 18, 073028 (2016), arXiv:1602.04526 [physics.atom-ph] .
  * Long  _et al._ (2017) J. Long _et al._ , Phys. Rev. C 96, 015502 (2017).
  * Fenker  _et al._ (2018) B. Fenker _et al._ , Phys. Rev. Lett. 120, 062502 (2018), arXiv:1706.00414 [nucl-ex] .
  * Brodeur  _et al._ (2016) M. Brodeur, J. Kelly, J. Long, C. Nicoloff, and B. Schultz, Nucl. Instrum. Meth. B 376, 281 (2016).
  * Shidling, P. D. and others (2018) Shidling, P. D. and others, Phys. Rev. C 98, 015502 (2018).
  * Valverde  _et al._ (2018) A. A. Valverde _et al._ , Phys. Rev. C 97, 035503 (2018).
  * O’Malley  _et al._ (2020) P. D. O’Malley _et al._ , Nucl. Instrum. Meth. B 463, 488 (2020).
  * Burdette  _et al._ (2020) D. P. Burdette _et al._ , Phys. Rev. C 101, 055504 (2020).
  * Long  _et al._ (2020) J. Long _et al._ , Phys. Rev. C 101, 015501 (2020).
  * Müller  _et al._ (2022) P. Müller _et al._ , Phys. Rev. Lett. 129, 182502 (2022), arXiv:2206.00742 [nucl-ex] .
  * Long  _et al._ (2022) J. Long _et al._ , Phys. Rev. C 106, 045501 (2022).
  * Plestid (2023) R. Plestid, “Generalized eikonal identities for charged currents,” To appear (2023).
  * Gennari  _et al._ (2023) M. Gennari _et al._ , “Standard model corrections to fermi transitions in light nuclei,”  (2023), INT PROGRAM 23-1B, New physics searches at the precision frontier .
  * Beg  _et al._ (1969) M. A. B. Beg, J. Bernstein, and A. Sirlin, Phys. Rev. Lett. 23, 270 (1969).
  * Ji and Musolf (1991) X.-D. Ji and M. J. Musolf, Phys. Lett. B 257, 409 (1991).
  * Chetyrkin and Grozin (2003) K. G. Chetyrkin and A. G. Grozin, Nucl. Phys. B 666, 289 (2003), arXiv:hep-ph/0303113 .
  * Baikov  _et al._ (2012) P. A. Baikov, K. G. Chetyrkin, J. H. Kuhn, and J. Rittinger, JHEP 07, 017 (2012), arXiv:1206.1284 [hep-ph] .
  * Herzog  _et al._ (2017) F. Herzog, B. Ruijl, T. Ueda, J. A. M. Vermaseren, and A. Vogt, JHEP 02, 090 (2017), arXiv:1701.01404 [hep-ph] .
  * Hill (2017) R. J. Hill, Phys. Rev. D 95, 013001 (2017), arXiv:1605.02613 [hep-ph] .