All-order renormalization of electric charge
in the Standard Model and beyond

Stefan Dittmaier

Albert-Ludwigs-Universität Freiburg, Physikalisches Institut, D-79104
Freiburg, Germany

*<EMAIL_ADDRESS>
15th International Symposium on Radiative Corrections:
Applications of Quantum Field Theory to Phenomenology,
FSU, Tallahasse, FL, USA, 17-21 May 2021 10.21468/SciPostPhysProc.?

## Abstract

Electric charge, as defined in the Thomson limit of the electron–photon
interaction vertex, is renormalized to all orders both in the Standard Model
and in any spontaneously broken gauge theory with gauge group $G\times$U(1)
with a group factor U(1) that mixes with electromagnetic gauge symmetry. In
the framework of the background-field method the charge renormalization
constant $Z_{e}$ is directly obtained from the photon wave-function
renormalization constant, similar to the situation in QED, which proves charge
universality as a byproduct. Exploiting charge universality in arbitrary
$R_{\xi}$ gauge by formulating the charge renormalization condition for a
“fake fermion” that couples only via an infinitesimal electric charge, $Z_{e}$
can be expressed in terms of renormalization constants that are obtained
solely from gauge-boson self-energies.

###### Contents

  1. 1 Introduction
  2. 2 Charge renormalization and charge universality in the background-field method
  3. 3 Charge renormalization in arbitrary $R_{\xi}$-gauge
  4. 4 Generalization to non-standard gauge groups
  5. 5 Conclusion

##  1 Introduction

The definition of electric unit charge is carried over from classical
electrodynamics to quantum electrodynamics (QED) and to more comprehensive
quantum field theories such as the Standard Model (SM) upon imposing the
Thomson renormalization condition which demands that the electron–photon
interaction vertex for physical (on-shell) electrons does not receive any
radiative corrections in the limit of low-energy photons. The charge
renormalization constant $Z_{e}$, which relates the bare charge $e_{0}=Z_{e}e$
to the renormalized unit charge $e$, can, thus, be determined from the
ee$\gamma$ vertex correction in this limit by direct calculation of this
correction in a given perturbative order. Exploiting electromagnetic gauge
invariance in the form of the famous Ward identity for the ee$\gamma$ vertex,
in QED it is possible to calculate $Z_{e}$ from the wave-function
renormalization constant of the photon, i.e. merely from a self-energy. This
is not only a welcome technical simplification, but also an important field-
theoretical result that is helpful in proving structural theorems such as
charge universality or Thirring’s theorem for low-energy Compton scattering
[1, 2]. Charge universality is the non-trivial statement that the renormalized
unit charge does not depend on the fermion species $f$ (actually not even on
the type of charged particle) that is used in the formulation of the Thomson
renormalization condition on the $ff\gamma$ vertex.

In the SM and most of its extensions, the derivation of $Z_{e}$ from self-
energies is very complicated (and not known beyond the one-loop level) if the
derivation is directly based on gauge invariance, which manifests itself via
Slavnov–Taylor identities for Green functions and Lee identities for (one-
particle irreducible, 1PI) vertex functions. This is due to the fact that the
unbroken electromagnetic U(1)em gauge symmetry is not a mere factor of the
gauge group, but mixes with a non-abelian subgroup in a non-trivial way. Other
U(1) group factors, such as the weak hypercharge group in the SM, which mix
with electromagnetic gauge symmetry cannot fully play the role of U(1)em of
QED, because they are spontaneously broken. Owing to this difficulty, in the
early pioneering proposals for electroweak one-loop renormalization [3, 4, 5,
6, 7, 8, 9, 10] the derivation of $Z_{e}$ was either based on explicit
calculations of the vertex correction or on Ward identities that were verified
by explicit calculations. A derivation fully based on Lee identities at one
loop without any explicit loop calculation has only been given quite recently
[11];111An alternative proof based on Slavnov–Taylor identities is sketched in
the slides of the corresponding talk given at the conference. owing to its
complexity a generalization of this approach beyond the one-loop level seems
rather complicated, if not infeasible.

The first all-order result for $Z_{e}$ in the SM, again expressed in terms of
the photon wave-function renormalization constant, was given in the framework
of the background-field method (BFM) in Ref. [12]. This result, in particular,
proves charge universality in the SM to all orders. In Ref. [13] it was shown
by explicit calculation that the BFM result for $Z_{e}$ indeed is in line with
the Thomson condition for the electron–photon vertex at the two-loop level in
arbitrary $R_{\xi}$ gauge. Similarly, in Refs. [14, 15, 16] it was shown by
explicit two-loop calculation in conventional ’t Hooft–Feynman gauge that the
sum of all genuine vertex corrections and fermionic wave-function corrections
to the $ff\gamma$ vertex vanishes in the Thomson limit; this is exactly the
part in the calculation of $Z_{e}$ that is ruled by gauge invariance but not
yet proven on the basis of Slavnov–Taylor or Lee identities to all orders. In
fact a general result on $Z_{e}$ expressed in terms of other renormalization
constants that are related to gauge-boson self-energies, was suggested in Ref.
[17]—although correct, this result was actually more conjectured than
derived.222Detailed comments on the arguments given in Ref. [17] can be found
in the appendix of Ref. [18]. This form of $Z_{e}$ was subsequently used in
the few existing explicit electroweak two-loop calculations, such as for the
muon decay [19, 20] or Z-boson decay widths [21].

After briefly sketching the BFM derivation of $Z_{e}$ within the SM to all
orders, in the following we review the all-order derivation [18] of $Z_{e}$
within arbitrary $R_{\xi}$ gauge, exploiting charge universality by imposing
the Thomson renormalization condition on the electromagnetic interaction
vertex of a “fake fermion” $\eta$ with infinitesimal electric charge
$Q_{\eta}$ but no other charges or couplings, so that $\eta$ fully decouples
from all other particles for $Q_{\eta}\to 0$. We first describe the derivation
within the SM, where the result confirms the earlier “conjecture” of Ref.
[17], and then generalize it to a more general gauge group $G\times$U(1),
where $G$ is any compact Lie group and U(1) has an admixture of
electromagnetic gauge symmetry. This article is just a brief summary of Ref.
[18], where a much more complete treatment of the subject can be found.

##  2 Charge renormalization and charge universality in the background-field
method

Denoting the relative charge and mass of the fermion $f$ by $Q_{f}$ and
$m_{f}$, respectively, the Thomson renormalization condition reads

$\displaystyle\left.\bar{u}(p)\,\hat{\Gamma}^{\hat{A}\bar{f}f}_{{\mathrm{R}},\mu}(0,-p,p)\,u(p)\right|_{p^{2}=m_{f}^{2}}=-Q_{f}e\,\bar{u}(p)\gamma_{\mu}u(p),$
(1)

where $\hat{A}$ is the background photon field and
$\hat{\Gamma}^{\hat{A}\bar{f}f}_{{\mathrm{R}},\mu}$ the renormalized
$\hat{A}\bar{f}f$ vertex function in the BFM. Here $\bar{u}(p)$ and $u(p)$ are
Dirac spinors of the fermion $f$ with momentum $p$ fulfilling
$p^{2}=m_{f}^{2}$ with the renormalized on-shell mass $m_{f}$. In the notation
and conventions for all field-theoretical quantities we follow Ref. [11]
throughout.

The needed low-energy limit of the $\hat{A}\bar{f}f$ for on-shell fermions can
be obtained from its BFM Ward identity, which follows from the background-
field gauge invariance of the BFM effective action (see Refs. [12, 11] and
references therein). This Ward identity for the unrenormalized
$\hat{A}\bar{f}f$ vertex function reads [12]333We ignore fermion generation
mixing here; for its restoration see Ref. [18].

$\displaystyle
k^{\mu}\hat{\Gamma}^{\hat{A}\bar{f}f}_{\mu}(k,\bar{p},p)={}-e_{0}Q_{f}\left[\hat{\Gamma}^{\bar{f}f}(\bar{p},-\bar{p})-\hat{\Gamma}^{\bar{f}f}(-p,p)\right],$
(2)

where $\hat{\Gamma}^{\bar{f}f}$ is the unrenormalized two-point vertex
function of the fermions. To exploit this identity in condition (1), we have
to replace unrenormalized by renormalized quantities. Indicating bare
quantities consistently by subscripts 0, the relevant parts of this
renormalization transformation reads

$\displaystyle
f_{0}^{\sigma}=\left(Z^{f,\sigma}\right)^{1/2}\,f^{\sigma},\quad\bar{f}_{0}^{\sigma}=\left(Z^{f,\sigma\,*}\right)^{1/2}\,\bar{f}^{\sigma},\quad\quad\begin{pmatrix}\hat{Z}_{0}\\\
\hat{A}_{0}\end{pmatrix}={}\begin{pmatrix}Z_{\hat{Z}\hat{Z}}^{1/2}&Z_{\hat{Z}\hat{A}}^{1/2}\\\\[4.30554pt]
Z_{\hat{A}\hat{Z}}^{1/2}&Z_{\hat{A}\hat{A}}^{1/2}\end{pmatrix}\begin{pmatrix}\hat{Z}\\\
\hat{A}\end{pmatrix},$ (3)

where $\sigma={\mathrm{R}},{\mathrm{L}}$ refers to the right- and left-handed
parts of the fermion field $f$ and $\hat{Z}$ is the background Z-boson field.
The resulting relation between renormalized and unrenormalized vertex
functions reads

$\displaystyle\hat{\Gamma}^{\bar{f}f}_{{\mathrm{R}},\mu}(-p,p)={}$
$\displaystyle\left({Z^{f,\sigma}}^{*}\right)^{1/2}\,\left(Z^{f,\sigma}\right)^{1/2}\,\hat{\Gamma}^{\bar{f}f}_{\mu}(-p,p),$
(4)
$\displaystyle\hat{\Gamma}^{\hat{A}\bar{f}f}_{{\mathrm{R}},\mu}(k,\bar{p},p)={}$
$\displaystyle\sum_{\hat{V}=\hat{A},\hat{Z}}Z_{\hat{V}\hat{A}}^{1/2}\,\left({Z^{f,\sigma}}^{*}\right)^{1/2}\,\left(Z^{f,\sigma}\right)^{1/2}\,\hat{\Gamma}^{\hat{V}\bar{f}f}_{\mu}(k,\bar{p},p).$
(5)

The introduced field-renormalization constants $Z^{f,\sigma}$ for the fermions
and $Z_{\hat{V}\hat{V}^{\prime}}$ for the photon–Z-boson system are fixed by
on-shell (OS) renormalization conditions, which require canonical
normalization of the residues of particle propagators and eliminates mixing
between different fields for on-shell momenta. Background-field gauge
invariance automatically implies $Z_{\hat{Z}\hat{A}}=0$ (see, e.g., Refs. [12,
11] for more details). Using additionally $e_{0}=Z_{e}e$, turns the Ward
identity (2) into an analogous identity for renormalized quantities,

$\displaystyle
k^{\mu}\hat{\Gamma}^{\hat{A}\bar{f}f}_{{\mathrm{R}},\mu}(k,\bar{p},p)={}-eQ_{f}\,Z_{e}\,Z_{\hat{A}\hat{A}}^{1/2}\,\left[\hat{\Gamma}^{\bar{f}f}_{\mathrm{R}}(\bar{p},-\bar{p})-\hat{\Gamma}^{\bar{f}f}_{\mathrm{R}}(-p,p)\right],$
(6)

which is valid for arbitrary momenta $k$, $\bar{p}$, $p$ obeying
$k+\bar{p}+p=0$. Taking $k\to 0$ for fixed $p$, the terms linear in $k$ obey
the relation

$\displaystyle\hat{\Gamma}^{\hat{A}\bar{f}f}_{{\mathrm{R}},\mu}(0,-p,p)={}$
$\displaystyle-
eQ_{f}\,Z_{e}\,Z_{\hat{A}\hat{A}}^{1/2}\,\frac{\partial\hat{\Gamma}^{\bar{f}f}_{\mathrm{R}}(-p,p)}{\partial
p^{\mu}}.$ (7)

Sandwiching this relation between Dirac spinors, its l.h.s. becomes identical
to the one of (7), and its r.h.s. can be simplified with the OS
renormalization condition for the fermion field, which can be written as

$\displaystyle\bar{u}(p)\,\frac{\partial\hat{\Gamma}^{\bar{f}f}_{\mathrm{R}}(-p,p)}{\partial
p^{\mu}}\,u(p)={}$ $\displaystyle\bar{u}(p)\,\gamma_{\mu}\,u(p).$ (8)

Thus, combining the charge renormalization condition (1) with (7) and (8)
leads to the simple equation [12]

$\displaystyle Z_{e}=Z_{\hat{A}\hat{A}}^{-1/2}$ (9)

in the BFM, which is formally identical to the well-known relation in QED.
Note also that Eq. (9) shows that $e_{0}\hat{A}_{0,\mu}(x)=e\hat{A}_{\mu}(x)$,
i.e. that the product of electromagnetic coupling and background photon field
is not renormalized, again in analogy to a QED relation.

##  3 Charge renormalization in arbitrary $R_{\xi}$-gauge

We now extend the SM by adding a fermion field $\eta$ with vanishing weak
isospin, $I_{{\mathrm{w}},\eta}^{a}=0$, and weak hypercharge
$Y_{{\mathrm{w}},\eta}=2Q_{\eta}$. Taking the limit of vanishing electric
charge, $Q_{\eta}\to 0$, the fermion $\eta$ decouples from all other
particles, and we recover the original theory—therefore the terminology “fake
fermion”. The Lagrangian ${\cal L}$ of the SM is, thus, modified by adding

$\displaystyle{\cal L}_{\eta}=\bar{\eta}\Bigl{(}{\mathrm{i}}\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-6.97566pt}\hbox{$\displaystyle\partial$}}{\hbox
to0.0pt{\hss$\textstyle/$\hss\kern-6.97566pt}\hbox{$\textstyle\partial$}}{\hbox
to0.0pt{\hss$\scriptstyle/$\hss\kern-5.38295pt}\hbox{$\scriptstyle\partial$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-4.32114pt}\hbox{$\scriptscriptstyle\partial$}}-{\textstyle\frac{1}{2}}g_{1}Y_{{\mathrm{w}},\eta}\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-8.0868pt}\hbox{$\displaystyle B$}}{\hbox
to0.0pt{\hss$\textstyle/$\hss\kern-8.0868pt}\hbox{$\textstyle B$}}{\hbox
to0.0pt{\hss$\scriptstyle/$\hss\kern-5.66075pt}\hbox{$\scriptstyle B$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-4.0434pt}\hbox{$\scriptscriptstyle
B$}}-m_{\eta}\Bigr{)}\eta=\overline{\eta}\left[{\mathrm{i}}\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-6.97566pt}\hbox{$\displaystyle\partial$}}{\hbox
to0.0pt{\hss$\textstyle/$\hss\kern-6.97566pt}\hbox{$\textstyle\partial$}}{\hbox
to0.0pt{\hss$\scriptstyle/$\hss\kern-5.38295pt}\hbox{$\scriptstyle\partial$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-4.32114pt}\hbox{$\scriptscriptstyle\partial$}}-Q_{\eta}e\left(\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-7.50002pt}\hbox{$\displaystyle
A$}}{\hbox to0.0pt{\hss$\textstyle/$\hss\kern-7.50002pt}\hbox{$\textstyle
A$}}{\hbox to0.0pt{\hss$\scriptstyle/$\hss\kern-5.25pt}\hbox{$\scriptstyle
A$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-3.75pt}\hbox{$\scriptscriptstyle
A$}}+\frac{s_{\scriptscriptstyle{\mathrm{W}}}}{c_{\scriptscriptstyle{\mathrm{W}}}}\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-7.54167pt}\hbox{$\displaystyle
Z$}}{\hbox to0.0pt{\hss$\textstyle/$\hss\kern-7.54167pt}\hbox{$\textstyle
Z$}}{\hbox to0.0pt{\hss$\scriptstyle/$\hss\kern-5.27916pt}\hbox{$\scriptstyle
Z$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-3.77083pt}\hbox{$\scriptscriptstyle
Z$}}\right)-m_{\eta}\right]\eta,$ (10)

with $m_{\eta}$ denoting the arbitrary mass of the fermion $\eta$. We note
that the mass term for $\eta$ is gauge invariant, that $\eta$ is stable, and
that no anomalies are introduced owing to the non-chirality of $\eta$. As in
the SM, $g_{1}$ is the U(1)Y gauge coupling, $B^{\mu}$ the U(1)Y gauge field,
and $s_{\scriptscriptstyle{\mathrm{W}}}=\sin\theta_{\mathrm{w}}$ and
$c_{\scriptscriptstyle{\mathrm{W}}}=\cos\theta_{\mathrm{w}}$ the sine and
cosine of the weak mixing angle $\theta_{\mathrm{w}}$. Employing charge
universality, we can take the Thomson limit of the $A\bar{\eta}\eta$ vertex to
define the renormalized electric unit charge $e$,

$\displaystyle\left.\bar{u}(p)\,\Gamma^{A\bar{\eta}\eta}_{{\mathrm{R}},\mu}(0,-p,p)\,u(p)\right|_{p^{2}=m_{\eta}^{2}}=-Q_{\eta}e\,\bar{u}(p)\gamma_{\mu}u(p).$
(11)

The relation between $\Gamma^{A\bar{\eta}\eta}_{{\mathrm{R}},\mu}$ and its
bare counterpart $\Gamma^{A\bar{\eta}\eta}_{\mu}$ follows from the field
renormalization transformation $\eta_{0}=Z^{1/2}_{\eta}\eta$ and the analog of
(3) for the photon–Z-boson system and reads

$\displaystyle\Gamma^{A\bar{\eta}\eta}_{{\mathrm{R}},\mu}(k,\bar{p},p)=Z_{\eta}Z_{AA}^{1/2}\,\Gamma^{A\bar{\eta}\eta}_{\mu}(k,\bar{p},p)+Z_{\eta}Z_{ZA}^{1/2}\,\Gamma^{Z\bar{\eta}\eta}_{\mu}(k,\bar{p},p).$
(12)

The bare vertex functions $\Gamma^{V\bar{\eta}\eta}_{\mu}$ ($V=A,Z$) receive
lowest-order contributions and bare vertex corrections
$\Lambda^{V\bar{\eta}\eta}_{\mu}$, which consist of 1PI loop diagrams and
tadpole corrections,

$\displaystyle\Gamma^{A\bar{\eta}\eta}_{\mu}(k,\bar{p},p)={}-Q_{\eta}e_{0}\gamma_{\mu}+e_{0}\Lambda^{A\bar{\eta}\eta}_{\mu},\qquad\Gamma^{Z\bar{\eta}\eta}_{\mu}(k,\bar{p},p)={}-Q_{\eta}e_{0}\frac{s_{{\mathrm{w}},0}}{c_{{\mathrm{w}},0}}\gamma_{\mu}+e_{0}\Lambda^{Z\bar{\eta}\eta}_{\mu}.$
(13)

The important observation is now that all diagrammatic contributions to
$\Lambda^{V\bar{\eta}\eta}_{\mu}$ and $Z_{\eta}$ involve at least two
couplings of photons or Z bosons to the $\eta$ line that passes through the
whole diagram. Some sample diagrams are shown in Fig. 1.

(a) (b) (c)

(d) (e) (f)

Figure 1: Some higher-order diagrams contributing to
$\Gamma^{V\bar{\eta}\eta}_{\mu}$ (a–c) and $\Gamma^{\bar{\eta}\eta}$ (d–f),
which receive contributions from 1PI diagrams (a, b, d, e), from explicit
tadpole diagrams (c), and from diagrams involving tadpole counterterms $\delta
t$ (f).

For 1PI diagrams it is obvious that at least two couplings to the $\eta$ line
exist, for diagrams with tadpole loops or tadpole counterterms the same holds
true, because the Higgs field $H$ does not couple to $\eta$. Since both the
photon and the Z boson couple to $\eta$ proportional to $Q_{\eta}$, this means
that $\Lambda^{V\bar{\eta}\eta}_{\mu}={\cal O}(Q_{\eta}^{2})$ and
$Z_{\eta}=1+{\cal O}(Q_{\eta}^{2})$. Inserting, thus,
$\Gamma^{A\bar{\eta}\eta}_{{\mathrm{R}},\mu}$ from (12) into condition (11)
and keeping only terms linear in $Q_{\eta}$ for $Q_{\eta}\to 0$, we get

$\displaystyle-Q_{\eta}e\,\bar{u}(p)\gamma_{\mu}u(p)={}$
$\displaystyle\left.\bar{u}(p)\,Z_{\eta}\left[Z_{AA}^{1/2}\,\Gamma^{A\bar{\eta}\eta}_{\mu}(0,-p,p)+Z_{ZA}^{1/2}\,\Gamma^{Z\bar{\eta}\eta}_{\mu}(0,-p,p)\right]u(p)\right|_{p^{2}=m_{\eta}^{2}}$
$\displaystyle={}$ $\displaystyle-
Q_{\eta}e_{0}\left[Z_{AA}^{1/2}+Z_{ZA}^{1/2}\frac{s_{{\mathrm{w}},0}}{c_{{\mathrm{w}},0}}\right]\bar{u}(p)\gamma_{\mu}u(p)\,+\,{\cal
O}(Q_{\eta}^{2}).$ (14)

This relation implies

$\displaystyle
e=e_{0}\left[Z_{AA}^{1/2}+Z_{ZA}^{1/2}\frac{s_{{\mathrm{w}},0}}{c_{{\mathrm{w}},0}}\right],$
(15)

which is the desired relation between $e$ and $e_{0}$. Introducing the
renormalization constant $\delta c_{\scriptscriptstyle{\mathrm{W}}}^{2}$
according to

$\displaystyle
c_{{\mathrm{w}},0}^{2}=1-s_{{\mathrm{w}},0}^{2}=c_{\scriptscriptstyle{\mathrm{W}}}^{2}+\delta
c_{\scriptscriptstyle{\mathrm{W}}}^{2}=1-s_{\scriptscriptstyle{\mathrm{W}}}^{2}+\delta
c_{\scriptscriptstyle{\mathrm{W}}}^{2},$ (16)

we can determine $Z_{e}$ from (15),

$\displaystyle
Z_{e}=\left[Z_{AA}^{1/2}+Z_{ZA}^{1/2}\sqrt{\frac{s_{\scriptscriptstyle{\mathrm{W}}}^{2}-\delta
c_{\scriptscriptstyle{\mathrm{W}}}^{2}}{c_{\scriptscriptstyle{\mathrm{W}}}^{2}+\delta
c_{\scriptscriptstyle{\mathrm{W}}}^{2}}}\right]^{-1}.$ (17)

This is fully equivalent to the result quoted and used in Refs. [17, 19, 20].
Since $Z_{ZA}$ vanishes at tree level, the renormalization constant $\delta
c_{\scriptscriptstyle{\mathrm{W}}}^{2}$ is only required to $(\ell-1)$-loop
order in the $\ell$-loop calculation of $Z_{e}$.

##  4 Generalization to non-standard gauge groups

The concept of charge universality and charge renormalization outlined above
for the SM can be generalized easily to gauge groups of the type
$G{\times}$U(1), where $G$ is any compact Lie group of rank $r$ and the U(1)
group factor plays the analogous role of weak hypercharge in the SM. More
precisely, we mean by this that the U(1)em subgroup of electromagnetic gauge
symmetry mixes transformations of U(1) and $G$, so that the photon field
$A^{\mu}$ is a non-trivial linear combination of the U(1) gauge field
$B^{\mu}$ and the gauge fields $C_{k}^{\mu}$ ($k=1,\dots,r$) of $G$
corresponding to the diagonal group generators in the Lie algebra of $G$. For
the mechanism of electroweak symmetry breaking we only assume that
electromagnetic gauge invariance is unbroken.

The original gauge fields $B^{\mu}$ and $\\{C_{k}^{\mu}\\}$ can be transformed
into fields that correspond to mass eigenstates,

$\displaystyle\begin{pmatrix}B^{\mu}\\\ C_{1}^{\mu}\\\ \vdots\\\
C_{r}^{\mu}\end{pmatrix}=R\begin{pmatrix}A^{\mu}\\\ Z_{1}^{\mu}\\\ \vdots\\\
Z_{r}^{\mu}\end{pmatrix},\qquad
R=\begin{pmatrix}R_{BA}&R_{BZ_{1}}&\cdots&R_{BZ_{r}}\\\
R_{C_{1}A}&R_{C_{1}Z_{1}}&\cdots&R_{C_{1}Z_{r}}\\\
\vdots&\vdots&\ddots&\vdots\\\
R_{C_{r}A}&R_{C_{r}Z_{1}}&\cdots&R_{C_{r}Z_{r}}\end{pmatrix},$ (18)

where $Z_{k}^{\mu}$ ($k=1,\dots,r$) describe neutral massive gauge bosons
similar to the Z boson of the SM. The matrix $R$ is a generalization of the SM
rotation matrix parametrized by the weak mixing angle, but is not necessarily
orthogonal or unitary. The OS renormalization of the gauge fields proceeds
exactly as in the SM, with an obvious generalization of the $2\times 2$ matrix
of renormalization constants $Z_{VV^{\prime}}$ of (3) to a $(r+1)\times(r+1)$
matrix. The constants $Z_{VV^{\prime}}$ can be computed from the $VV^{\prime}$
gauge-boson self-energies order by order.

The BFM derivation of the charge renormalization constant generalizes to the
more general gauge group without any problems. Background-field gauge
invariance implies $Z_{\hat{Z}_{k}\hat{A}}=0$ for all $k=1,\dots,r$, i.e.
there is no mixing of on-shell photons with any of the Zk bosons. The Ward
identities (2) and (6), thus, carry over without modification. As a result,
the charge renormalization constant $Z_{e}$ is given by
$Z_{e}=Z_{\hat{A}\hat{A}}^{-1/2}$ as in (9), proving charge universality as in
the SM.

To exploit charge universality in the determination of $Z_{e}$ in arbitrary
$R_{\xi}$-gauge, we again introduce a fake fermion $\eta$ with the same
properties as above, i.e. $\eta$ only carries infinitesimal U(1) charge
$Y_{{\mathrm{w}},\eta}=2Q_{\eta}$, but no non-trivial quantum number of $G$,
so that the unit charge is given by $e=g_{1}R_{BA}$. If the model contains
singlet scalars $S_{i}$, the scalars $S_{i}$ may couple to $\eta$ via Yukawa
couplings. The corresponding couplings $y_{i}$ are free parameters of the
model and can be taken to be infinitesimally small in analogy to $Q_{\eta}\to
0$, so that decoupling of $\eta$ is guaranteed. The Lagrangian ${\cal
L}_{\eta}$ reads

$\displaystyle{\cal L}_{\eta}={}$
$\displaystyle\bar{\eta}\Bigl{(}{\mathrm{i}}\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-6.97566pt}\hbox{$\displaystyle\partial$}}{\hbox
to0.0pt{\hss$\textstyle/$\hss\kern-6.97566pt}\hbox{$\textstyle\partial$}}{\hbox
to0.0pt{\hss$\scriptstyle/$\hss\kern-5.38295pt}\hbox{$\scriptstyle\partial$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-4.32114pt}\hbox{$\scriptscriptstyle\partial$}}-{\textstyle\frac{1}{2}}g_{1}Y_{{\mathrm{w}},\eta}\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-8.0868pt}\hbox{$\displaystyle B$}}{\hbox
to0.0pt{\hss$\textstyle/$\hss\kern-8.0868pt}\hbox{$\textstyle B$}}{\hbox
to0.0pt{\hss$\scriptstyle/$\hss\kern-5.66075pt}\hbox{$\scriptstyle B$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-4.0434pt}\hbox{$\scriptscriptstyle
B$}}-m_{\eta}-\sum_{i}y_{i}S_{i}\Bigr{)}\eta$ $\displaystyle={}$
$\displaystyle\overline{\eta}\bigg{[}{\mathrm{i}}\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-6.97566pt}\hbox{$\displaystyle\partial$}}{\hbox
to0.0pt{\hss$\textstyle/$\hss\kern-6.97566pt}\hbox{$\textstyle\partial$}}{\hbox
to0.0pt{\hss$\scriptstyle/$\hss\kern-5.38295pt}\hbox{$\scriptstyle\partial$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-4.32114pt}\hbox{$\scriptscriptstyle\partial$}}-Q_{\eta}e\bigg{(}\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-7.50002pt}\hbox{$\displaystyle
A$}}{\hbox to0.0pt{\hss$\textstyle/$\hss\kern-7.50002pt}\hbox{$\textstyle
A$}}{\hbox to0.0pt{\hss$\scriptstyle/$\hss\kern-5.25pt}\hbox{$\scriptstyle
A$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-3.75pt}\hbox{$\scriptscriptstyle
A$}}+\sum_{k}\frac{R_{BZ_{k}}}{R_{BA}}\mathchoice{\hbox
to0.0pt{\hss$\displaystyle/$\hss\kern-7.54167pt}\hbox{$\displaystyle
Z$}}{\hbox to0.0pt{\hss$\textstyle/$\hss\kern-7.54167pt}\hbox{$\textstyle
Z$}}{\hbox to0.0pt{\hss$\scriptstyle/$\hss\kern-5.27916pt}\hbox{$\scriptstyle
Z$}}{\hbox
to0.0pt{\hss$\scriptscriptstyle/$\hss\kern-3.77083pt}\hbox{$\scriptscriptstyle
Z$}}_{k}\bigg{)}-m_{\eta}-\sum_{i}y_{i}S_{i}\bigg{]}\eta.$ (19)

Following the same reasoning as for the SM above, the renormalized
$A\bar{\eta}\eta$ vertex function is given by

$\displaystyle\Gamma^{A\bar{\eta}\eta}_{{\mathrm{R}},\mu}(k,\bar{p},p)=Z_{\eta}Z_{AA}^{1/2}\,\Gamma^{A\bar{\eta}\eta}_{\mu}(k,\bar{p},p)+\sum_{k}Z_{\eta}Z_{Z_{k}A}^{1/2}\,\Gamma^{Z_{k}\bar{\eta}\eta}_{\mu}(k,\bar{p},p),$
(20)

with the unrenormalized $V\bar{\eta}\eta$ vertex functions

$\displaystyle\Gamma^{A\bar{\eta}\eta}_{\mu}(k,\bar{p},p)={}-Q_{\eta}e_{0}\gamma_{\mu}+e_{0}\Lambda^{A\bar{\eta}\eta}_{\mu},\qquad\Gamma^{Z_{k}\bar{\eta}\eta}_{\mu}(k,\bar{p},p)={}-Q_{\eta}e_{0}\,\frac{R_{0,BZ_{k}}}{R_{0,BA}}\,\gamma_{\mu}+e_{0}\Lambda^{Z_{k}\bar{\eta}\eta}_{\mu}.$
(21)

Again the vertex corrections $\Lambda^{A\bar{\eta}\eta}_{\mu}$ and
$\Lambda^{Z_{k}\bar{\eta}\eta}_{\mu}$ as well as the field renormalization
constant $Z_{\eta}$ receive only corrections that are suppressed at least by
quadratic factors in the new couplings, such as $Q_{\eta}^{2}$ or
$Q_{\eta}y_{i}$. Typical diagrams contributing to those corrections at the
order ${\cal O}(Q_{\eta}^{2})$ (or higher in $Q_{\eta}$) can be obtained from
the graphs shown in Fig. 1 upon interpreting the field $Z$ as any of the
$Z_{k}$ and taking the Higgs field $H$ as a any Higgs field of the model.
Equation (3) then generalizes to the considered model in an obvious way, and
we obtain the final result for the charge renormalization constant:

$\displaystyle
Z_{e}=\left[Z_{AA}^{1/2}+\sum_{k}Z_{Z_{k}A}^{1/2}\,\frac{R_{BZ_{k}}+\delta
R_{BZ_{k}}}{R_{BA}+\delta R_{BA}}\right]^{-1},$ (22)

where $\delta R_{BA}$ and $\delta R_{BZ_{k}}$ are renormalization constants
for the matrix elements of $R$, i.e. $R_{0}=R+\delta R$. Recalling that the
constants $Z_{Z_{k}A}$ vanish at tree level, we see that the constants $\delta
R_{BA}$ and $\delta R_{BZ_{k}}$ are only required to $(\ell-1)$-loop order in
the $\ell$-loop calculation of $Z_{e}$.

The case of the SM is trivially recovered from the results of this section
upon identifying $G={}$SU(2)w, $r=1$, $Z_{1}^{\mu}=Z^{\mu}$,
$R_{BZ_{1}}=s_{\scriptscriptstyle{\mathrm{W}}}$, and
$R_{BA}=c_{\scriptscriptstyle{\mathrm{W}}}$.

##  5 Conclusion

Employing the property of charge universality, the determination of the charge
renormalization constant $Z_{e}$ in arbitrary $R_{\xi}$ gauge can be greatly
simplified upon applying the Thomson renormalization condition to a “fake
fermion” that has infinitesimal electric charge but no other charges or
couplings. Both in the SM and in the wider class of gauge theories with gauge
group $G\times$U(1), where the U(1) subgroup contains some component of
electromagnetic gauge transformations, $Z_{e}$ can be deduced from gauge-boson
wave-function and parameter renormalization constants which can be calculated
from gauge-boson self-energies only. Using Slavnov–Taylor or Lee identities to
derive this result is already very cumbersome at the one-loop level, and a
consistent generalization beyond one loop is not known, if not infeasible.

The assumed property of charge universality can for instance proven in the
framework of the background-field method from which it is known that $Z_{e}$
can be directly obtained from the photon wave-function renormalization
constant.

For even more general gauge groups without a U(1) factor mixing with
electromagnetic gauge symmetry, the strategy with the decoupling “fake
fermion” seems not to be applicable. In this case the consistent use of the
background-field method, however, still bears the possibility to calculate the
charge renormalization constant in terms of gauge-boson wave-function
renormalization constants in contrast to conventional $R_{\xi}$-like gauges.

## References

  * [1] W. E. Thirring,  _Radiative corrections in the nonrelativistic limit_ ,  Phil. Mag. Ser. 7 41, 1193 (1950).
  * [2] S. Dittmaier,  _Thirring’s low-energy theorem and its generalizations in the electroweak standard model_ ,  Phys. Lett. B 409, 509 (1997),  10.1016/S0370-2693(97)00888-5,  hep-ph/9704368.
  * [3] D. A. Ross and J. C. Taylor,  _Renormalization of a unified theory of weak and electromagnetic interactions_ ,  Nucl. Phys. B 51, 125 (1973),  10.1016/0550-3213(73)90505-1,  [Erratum: Nucl.Phys.B 58, 643–643 (1973)].
  * [4] A. Sirlin,  _Radiative Corrections in the SU(2)-L x U(1) Theory: A Simple Renormalization Framework_ ,  Phys. Rev. D 22, 971 (1980),  10.1103/PhysRevD.22.971.
  * [5] D. Y. Bardin, P. K. Khristova and O. M. Fedorenko,  _On the Lowest Order Electroweak Corrections to Spin 1/2 Fermion Scattering. 1. The One Loop Diagrammar_ ,  Nucl. Phys. B 175, 435 (1980),  10.1016/0550-3213(80)90021-8.
  * [6] J. Fleischer and F. Jegerlehner,  _Radiative Corrections to Higgs Decays in the Extended Weinberg-Salam Model_ ,  Phys. Rev. D 23, 2001 (1981),  10.1103/PhysRevD.23.2001.
  * [7] K. I. Aoki, Z. Hioki, M. Konuma, R. Kawabe and T. Muta,  _Electroweak Theory. Framework of On-Shell Renormalization and Study of Higher Order Effects_ ,  Prog. Theor. Phys. Suppl. 73, 1 (1982),  10.1143/PTPS.73.1.
  * [8] M. Böhm, H. Spiesberger and W. Hollik,  _On the One Loop Renormalization of the Electroweak Standard Model and Its Application to Leptonic Processes_ ,  Fortsch. Phys. 34, 687 (1986),  10.1002/prop.19860341102.
  * [9] W. F. L. Hollik,  _Radiative Corrections in the Standard Model and their Role for Precision Tests of the Electroweak Theory_ ,  Fortsch. Phys. 38, 165 (1990),  10.1002/prop.2190380302.
  * [10] A. Denner,  _Techniques for calculation of electroweak radiative corrections at the one loop level and results for W physics at LEP-200_ ,  Fortsch. Phys. 41, 307 (1993),  10.1002/prop.2190410402,  0709.1075.
  * [11] A. Denner and S. Dittmaier,  _Electroweak Radiative Corrections for Collider Physics_ ,  Phys. Rept. 864, 1 (2020),  10.1016/j.physrep.2020.04.001,  1912.06823.
  * [12] A. Denner, G. Weiglein and S. Dittmaier,  _Application of the background field method to the electroweak standard model_ ,  Nucl. Phys. B 440, 95 (1995),  10.1016/0550-3213(95)00037-S,  hep-ph/9410338.
  * [13] G. Degrassi and A. Vicini,  _Two loop renormalization of the electric charge in the standard model_ ,  Phys. Rev. D 69, 073007 (2004),  10.1103/PhysRevD.69.073007,  hep-ph/0307122.
  * [14] S. Actis, A. Ferroglia, M. Passera and G. Passarino,  _Two-Loop Renormalization in the Standard Model. Part I: Prolegomena_ ,  Nucl. Phys. B 777, 1 (2007),  10.1016/j.nuclphysb.2007.04.021,  hep-ph/0612122.
  * [15] S. Actis and G. Passarino,  _Two-Loop Renormalization in the Standard Model Part II: Renormalization Procedures and Computational Techniques_ ,  Nucl. Phys. B 777, 35 (2007),  10.1016/j.nuclphysb.2007.03.043,  hep-ph/0612123.
  * [16] S. Actis and G. Passarino,  _Two-Loop Renormalization in the Standard Model Part III: Renormalization Equations and their Solutions_ ,  Nucl. Phys. B 777, 100 (2007),  10.1016/j.nuclphysb.2007.04.027,  hep-ph/0612124.
  * [17] S. Bauberger,  _Two-loop contributions to muon decay_ ,  Phd thesis, Würzburg U. (1997).
  * [18] S. Dittmaier,  _Electric charge renormalization to all orders_ ,  Phys. Rev. D 103(5), 053006 (2021),  10.1103/PhysRevD.103.053006,  2101.05154.
  * [19] A. Freitas, W. Hollik, W. Walter and G. Weiglein,  _Electroweak two loop corrections to the $M_{W}-M_{Z}$ mass correlation in the standard model_,  Nucl. Phys. B 632, 189 (2002),  10.1016/S0550-3213(02)00243-2,  [Erratum: Nucl.Phys.B 666, 305–307 (2003)],  hep-ph/0202131.
  * [20] M. Awramik, M. Czakon, A. Onishchenko and O. Veretin,  _Bosonic corrections to Delta r at the two loop level_ ,  Phys. Rev. D 68, 053004 (2003),  10.1103/PhysRevD.68.053004,  hep-ph/0209084.
  * [21] I. Dubovyk, A. Freitas, J. Gluza, T. Riemann and J. Usovitsch,  _Complete electroweak two-loop corrections to Z boson production and decay_ ,  Phys. Lett. B 783, 86 (2018),  10.1016/j.physletb.2018.06.037,  1804.10236.


# Quarkonium production as a function of charged-particle multiplicity in pp
and p–Pb collisions measured by ALICE at the LHC

Discipline of Physics, School of Basic Sciences, Indian Institute of
Technology Indore, Indore- 453552, INDIA
E-mail

###### Abstract:

Quarkonium production as a function of the charged-particle multiplicity could
provide an insight into particle production processes at the partonic level in
hadronic collisions. It is believed that multiple partonic interactions play
an important role in particle production and affect both soft and hard
processes. The study of correlations between quarkonia and charged-particle
multiplicity may provide information about this. In this contribution, ALICE
measurements of J$/\psi$ and $\Upsilon$ production as a function of charged-
particle multiplicity are presented for pp collisions at center-of-mass
energies $\sqrt{s}$ = 5.02 and 13 TeV. A similar measurement performed in p–Pb
collisions at $\sqrt{s_{\rm{NN}}}$ = 8.16 TeV at both forward and backward
rapidity is also discussed.

##  1 Introduction

In pp collisions at the LHC, where the center-of-mass energy is very high, the
measurement of the inclusive production of charged hadrons and their
correlation to quarkonia is an important tool to characterize particle
production. Quarkonium (bound state of charm and anti-charm quarks or beauty
and anti-beauty quarks) production in pp collisions allows us to test
perturbative QCD calculations and is used as a reference for studies in heavy-
ion collisions. The interplay between hard and soft QCD processes in hadronic
collisions can be studied by measuring quarkonium production as a function of
charged-particle multiplicity. This study is believed to bring information
about the structure of the collision in terms of Multiple Partonic Interaction
(MPI), in which several interactions occur in parallel in a single pp
collision. The role of MPI in heavy-flavor production has already been
investigated by measuring the relative J$/\psi$ and D-meson yields as a
function of charged-particle multiplicity in pp collisions at $\sqrt{s}$ = 7
TeV [1]. In both cases, an increasing trend has been observed. With the help
of new data at $\sqrt{s}$ = 5.02 and 13 TeV, one can study even larger
relative multiplicities and learn about the energy dependence. Similarly, the
study of the multiplicity dependence of J$/\psi$ production in p–Pb collisions
at $\sqrt{s_{\rm{NN}}}$ = 5.02 TeV is extended to $\sqrt{s_{\rm{NN}}}$ = 8.16
TeV, for which higher relative multiplicity can be reached, compared to
previous measurements.

##  2 Experimental setup and analysis procedure

In ALICE [2] , J$/\psi$ mesons are reconstructed via the
J$/\psi\rightarrow\rm{e}^{+}\rm{e}^{-}$ and
J$/\psi\rightarrow\rm{\mu}^{+}\rm{\mu}^{-}$ decays using two different
spectrometers. The central barrel detector covers the pseudo-rapidity range
$|\eta|<0.9$ and includes the Inner Tracking System (ITS) and Time Projection
Chamber (TPC). It is used for the reconstruction of J$/\psi$ via the
$\rm{e}^{+}\rm{e}^{-}$ decay channel in the rapidity range $|y|<0.9$. In
contrast, J$/\psi$ are reconstructed in the rapidity range $-4.0<y<-2.5$ via
the $\rm{\mu}^{+}\rm{\mu}^{-}$ decay channel using the forward muon
spectrometer. Two V0 detectors, V0A and V0C, consisting of two scintillator
arrays located at $-3.7<\eta<-1.7$ and at 2.8 $<\eta<$ 5.1 are used for
triggering. The V0 detectors are also used as a high-multiplicity trigger.
Tracklets, i.e. track segments reconstructed in the Silicon Pixel Detector
(SPD) at $|\eta|<1$, are used for the determination of charged-particle
pseudo-rapidity density ($\rm{d}\it{N}_{\rm{ch}}/\rm{d}\eta$). Several cuts
are applied to determine the accurate position of the $\it{z}$-coordinate of
the collision vertex. The $z$ positions of the collision vertex is restricted
to the range $|z_{\rm{vtx}}|<$ 10 cm, to have sufficient acceptance of the SPD
for tracklets to be measured in the range $|\eta|<1$. The actual charged-
particle density in the $i^{\rm{th}}$ multiplicity interval is calculated as:

$\frac{{\langle d\it{N}_{\rm{ch}}/d\eta\rangle}_{i}}{{\langle
d\it{N}_{\rm{ch}}/d\eta\rangle}}=\frac{{f(\langle
N_{\rm{trk}}^{\rm{corr}}\rangle}_{i})}{\langle
dN_{\rm{ch}}/d\eta\rangle}_{\scalebox{0.5}{\rm{INEL }> 0}}$,

(1)

where ${\langle N_{\rm{trk}}^{\rm{corr}}\rangle}_{i}$ is the mean number of
corrected SPD tracklets in $i^{\rm{th}}$ multiplicity interval. The
correlation function $f$ is calculated using Monte Carlo (MC) simulations. The
measurement is performed for inelastic pp collisions with at least one charged
particle in $|\eta|<1$. The relative J$/\psi$ yield in the $i^{\rm{th}}$
charged-particle multiplicity interval is calculated as:

$\frac{\frac{\rm{d}\it{N}_{\rm{J}/\psi}^{\rm{~{}i}}}{\rm{d}y}}{\langle\frac{\rm{d}\it{N}_{\rm{J}/\psi}}{\rm{d}y}\rangle}=\frac{\it{N}_{\rm{J}/\psi}^{~{}\it{i}}}{\it{N}_{\rm{J}/\psi}}\times\frac{\it{N}_{\rm{MB}}}{\it{N}_{\rm{MB}}^{\it{i}}}\times\epsilon,$
(2)

where ($N_{J/\psi}^{\it{i}}$, $N_{\rm{J/\psi}}$) and ($N_{\rm{MB}}^{i}$,
$N_{\rm{MB}}$) are the number of J$/\psi$ and number of minimum bias events
(V0-triggered events) in the $i^{\rm{th}}$ multiplicity interval and
integrated over all multiplicity intervals, respectively. The factor
$\epsilon$ corrects for the efficiencies concerning trigger selection, event
cuts, pile up rejection and correction factor for INEL¿ 0 event selection.

For the analysis of p–Pb collisions, correction factors are applied on both
the measured J$/\psi$ relative yields and the relative charged-particle
density so that they correspond to Non-Single-Diffractive (NSD) [3] events.
Additionally, the J$/\psi$ mean transverse momentum ($\langle
p_{\rm{T}}^{\rm{J}/\psi}\rangle$) is measured as a function of the charged-
particle multiplicity by fitting the $\langle p_{\rm{T}}\rangle$ of unlike-
sign muon pairs as a function of the dimuon invariant mass.

##  3 Results

The top left panel of Fig. 1 shows the $p_{\rm{T}}$-integrated relative
J$/\psi$ yield at mid and forward rapidity as a function of charged-particle
multiplicity in pp collisions at $\sqrt{s}=$13 TeV, compared to that at
$\sqrt{s}=$ 5.02 TeV. In all cases the J$/\psi$ yield increases with
increasing charged-particle multiplicity. A similar trend has been observed
for pp collisions at $\sqrt{s}=$ 7 TeV and it is believed that this behavior
is due to the role of MPI on the production of J$/\psi$ [4, 5]. It can also be
observed that the relative J$/\psi$ yield increases as a function of relative
charged-particle density faster than linear when the J$/\psi$ signal and the
multiplicity are measured in the same rapidity window (mid-rapidity). On the
contrary, with the introduction of a rapidity gap between the signal
($2.5<y<4.0$) and the charged-particle density estimator (—$\eta|<1$), an
approximately linear increase of the relative J$/\psi$ yield with increasing
relative charged-particle density is observed. The different trends observed
with the introduction of a rapidity gap might hint to autocorrelations or jet
biases. One can also observe that the multiplicity dependence of the J$/\psi$
production at forward rapidity is the same in pp collisions at $\sqrt{s}$ =13
TeV and 5 TeV, indicating that particle production is dominantly driven by the
relative event activity independently of $\sqrt{s}$.

In the top-right panel of Fig. 1, the mid-rapidity measurement in pp collision
at $\sqrt{s}$ = 13 TeV is compared to different theoretical models in which
MPI is an important ingredient. They are: Ferreiro et al. [6], EPOS3 [7],
PYTHIA 8 [8] and Kopeliovich et al. [9]. All the models qualitatively
reproduce the trend and EPOS3 reproduces the behavior quantitatively well.
EPOS3 includes MPI and hydrodynamical expansion of the system and describes
the azimuthal correlation of D mesons with charged particles well [10]. The
EPOS3 calculation shown in this work is for D mesons in a $p_{\rm{T}}$ range
that is supposed to serve as good proxy of inclusive J$/\psi$ production. The
good agreement of the EPOS3 model with the data might show that the energy
density reached in pp collisions at the LHC might be high enough to be
described by a hydrodynamical evolution, but other explanations are also
possible. Additionally, the multiplicity dependence of J$/\psi$ production has
also been studied for several J$/\psi$ $p_{\rm{T}}$ intervals. In that case,
the observed trends are qualitatively well reproduced by PYTHIA 8 simulations.

To investigate the possible dependence of the observed increase on the quark
content of the quarkonium state, and on its mass, a similar study of the
relative $\Upsilon$(1S) and $\Upsilon$(2S) yields as a function of the
relative charged-particle density has been performed in pp collisions at
$\sqrt{s}$ = 13 TeV. A linear increase of the relative $\Upsilon$(1S) and
$\Upsilon$(2S) yields as a function of multiplicity is observed as well. To
better quantify the similarities between the different particle species, we
have measured the $\Upsilon(1\rm{S})$ to J$/\psi$ and the $\Upsilon(2\rm{S})$
to $\Upsilon(1\rm{S})$ relative yield ratios, as a function of multiplicity.
It can be seen in the bottom panels of Fig.1 that these ratios are consistent
with unity for all multiplicities. This study reveals that the multiplicity
dependence is the same for the various quarkonium states despite their
different masses and binding energies.

Figure 1: Top left: relative J$/\psi$ yield as a function of the relative
charged-particle density measured at mid and forward rapidity in pp collisions
at $\sqrt{s}=13$ TeV and at forward rapidity in pp collisions at
$\sqrt{s}=5.02$ TeV. The dotted line corresponds to a linear increase. Top
right: model comparison for the mid-rapidity result at $\sqrt{s}=13$ TeV.
Bottom panels: ratio of relative yield of $\Upsilon(1\rm{S})$ to J$/\psi$
(left) and $\Upsilon(2\rm{S})$ to $\Upsilon(1\rm{S})$ (right) as a function of
multiplicity.

The multiplicity dependence of J$/\psi$ production as a function of charged-
particle multiplicity has also been studied in p–Pb collisions at
$\sqrt{s_{\rm{NN}}}$ = 8.16 TeV. The result is compared to that obtained in
p–Pb collision at $\sqrt{s_{\rm{NN}}}$ = 5.02 TeV and is shown in the left
panel of Fig. 2. The measurements have been performed in different rapidity
ranges by inverting the directions of the lead and proton beams. The backward
rapidity result ($-4.46<y_{\rm{cms}}<-2.96$) corresponds to J$/\psi$ produced
in the direction of the Pb beam, whereas the forward rapidity result
($2.03<y_{\rm{cms}}<3.53$) corresponds to J$/\psi$ produced in the direction
of the proton beam. We observe a colliding energy independent but rapidity
dependent behavior for p–Pb collisions. A saturation towards higher
multiplicities is observed at forward rapidity. In this kinematic region, the
proton probes the small Bjorken-$x$ region of the Pb-nucleus, where cold
nuclear matter effects, e.g. gluon saturation and/or gluon shadowing, are
expected.

Similarly, the J$/\psi$ $\langle{p_{\rm{T}}}\rangle$ has been measured as a
function of the relative charged-particle density at $\sqrt{s_{\rm{NN}}}$ =
8.16 TeV for p–Pb collisions and compared to that measured at
$\sqrt{s_{\rm{NN}}}$ = 5.02 TeV. We observed that the relative J$/\psi$
$\langle{p_{\rm{T}}}\rangle$ is independent of rapidity as well as center-of-
mass energy within the current uncertainties. It first increases with
multiplicity, and then saturates to a similar value for both rapidity ranges.

Figure 2: Relative yields (left) and mean transverse momentum (right) of
J$/\psi$ as a function of relative charged-particle multiplicity for p–Pb
collisions at $\sqrt{s_{\rm{NN}}}$ = 8.16 TeV and $\sqrt{s_{\rm{NN}}}$ = 5.02
TeV. In the left panel the results are shown at both forward
($2.03<y_{\rm{cms}}<3.53$) and backward ($-4.46<y_{\rm{cms}}<-2.96$)
rapidities. In the right panel only the backward rapidity result is shown, but
similar trends are also observed at forward rapidity.

##  4 Summary

In this contribution, recent ALICE measurements of the multiplicity dependence
of quarkonium production have been presented for pp collisions at $\sqrt{s}$ =
5.02 and 13 TeV at forward rapidity. A similar study, performed in p–Pb
collisions at $\sqrt{s_{\rm{NN}}}$ = 8.16, along with a measurement of the
multiplicity dependence of the J$/\psi$ mean transverse momentum have also
been presented. The comparison of these new results to that obtained at mid-
rapidity hints for possible autocorrelations and jet bias effects when the
J$/\psi$ and the charged-particle multiplicity are measured in the same
rapidity range.

## References

  * [1] J. Adam et al. [ALICE Collaboration], JHEP 1509, 148 (2015).
  * [2] K. Aamodt et al. [ALICE Collaboration], JINST 3, S08002 (2008).
  * [3] D. Adamov et al. [ALICE Collaboration], Phys. Lett. B 776, 91 (2018)
  * [4] B. Abelev et al. [ALICE Collaboration], Phys. Lett. B 712, 165 (2012).
  * [5] D. Thakur et al., Phys. Rev. D 97, 094002 (2018).
  * [6] E. G. Ferreiro and C. Pajares, Phys. Rev. C 86, 034903 (2012).
  * [7] H. J. Drescher, M. Hladik, S. Ostapchenko, T. Pierog and K. Werner, Phys. Rept. 350, 93 (2001).
  * [8] T. Sjostrand, S. Mrenna and P. Z. Skands, Comput. Phys. Commun. 178, 852 (2008).
  * [9] B. Z. Kopeliovich, H. J. Pirner, I. K. Potashnikova, K. Reygers and I. Schmidt, Phys. Rev. D 88, 116002 (2013).
  * [10] J. Adam et al. [ALICE Collaboration], Eur. Phys. J. C 77, 245 (2017).


# Predicting the pseudocapacitive windows for MXene electrodes with voltage-
dependent cluster expansion models

James M. Goff  Department of Materials Science and Engineering, Materials
Research Institute, The Pennsylvania State University, University Park, PA
16802, USA<EMAIL_ADDRESS>Francisco Marques dos Santos Vieira  Department of
Materials Science and Engineering, Materials Research Institute, The
Pennsylvania State University, University Park, PA 16802, USA  Nathan D.
Keilbart  Lawrence Livermore National Laboratory, Livermore, California 94550,
USA  Yasuaki Okada  Murata Manufacturing Co., Ltd., 1-10-1 Higashikotari,
Nagaokakyo-shi, Kyoto 617-8555, Japan  Ismaila Dabo  Department of Materials
Science and Engineering, Materials Research Institute, The Pennsylvania State
University, University Park, PA 16802, USA

###### Abstract

MXene transition-metal carbides and nitrides are of growing interest for
energy storage applications. These compounds are especially promising for use
as pseudocapacitive electrodes due to their ability to convert energy
electrochemically at fast rates. Using voltage-dependent cluster expansion
models, we predict the charge storage performance of MXene pseudocapacitors
for a range of electrode compositions. $M_{3}C_{2}O_{2}$ electrodes based on
group-VI transition metals have up to 80 percent larger areal energy densities
than prototypical titanium-based ( e.g. $Ti_{3}C_{2}O_{2}$ ) MXene electrodes.
We attribute this high pseudocapacitance to the Faradaic voltage windows of
group-VI MXene electrodes, which are predicted to be 1.2 to 1.8 times larger
than those of titanium-based MXenes. The size of the pseudocapacitive voltage
window increases with the range of oxidation states that is accessible to the
MXene transition metals. By similar mechanisms, the presence of multiple ions
in the solvent ( $Li^{+}$ and $H^{+}$ ) leads to sharp changes in the
transition-metal oxidation states and can significantly increase the charge
capacity of MXene pseudocapacitors.

Department of Materials Science and Engineering, Materials Research Institute,
The Pennsylvania State University, University Park, PA 16802, USA Penn State
Institutes of Energy and the Environment, The Pennsylvania State University,
University Park, PA 16802, USA

Keywords: Cluster expansion, Monte Carlo sampling, Electrochemistry, Energy
storage, Pseudocapacitor, MXene

Pseudocapacitive voltage windows predicted with voltage-dependent cluster
expansions of MXene electrode surfaces

##  1 Introduction

MXenes have garnered attention for a range of applications, including
catalysis, electronics, gas sensing, and energy storage.1, 2, 3 Their layered
structures (Fig. 1a) and their electronic properties make them excellent
candidates for pseudocapacitor and battery electrodes.4, 5, 6, 7
Pseudocapacitors are energy storage devices exhibiting higher energy densities
than supercapacitors and higher power densities than batteries, which can be
used in hybrid electrochemical systems to achieve electrical characteristics
beyond the nominal performance of devices operating via a single charge
storage mechanism (e.g. bulk intercalation or double-layer polarization).8
They are of practical relevance to high-power storage technologies for
regenerative braking on trains, power electronics, and heavy machinery, among
other potential uses.9, 10

Figure 1: (a) Hexagonal crystal structure of the prototypical
$\rm{Ti_{3}C_{2}O_{2}}$ MXene. (b) Idealized pseudocapacitive current-voltage
response with corresponding ion adsorption states with representative
structures (viewed along the $c$ axis) shown as insets (top). As the ion
coverage decreases with the electrode voltage, the oxidation state of the
transition metal increases (bottom). The nearly continuous set of redox
reactions allowed by MXenes is responsible for the pseudocapacitive response.

Pseudocapacitors store charge by reversible redox reactions occurring near the
surface of the material, and have characteristic capacitor-like
current–voltage responses, as depicted schematically in Fig. 1b. They differ
from electrochemical double-layer (EDL) capacitors in their charge storage
mechanisms; the primary process of energy storage in pseudocapacitors is
Faradaic charging—i.e charge accumulation due to near-surface chemical
reactions. As evidenced by this and other studies, pseudocapacitors are
characterized by a nearly continuous set of closely related redox states along
the pseudocapacitive voltage window of the material.1, 9, 11, 12 These
different redox states can be induced in the transition metals by the
adsorption of ions (H, Li, and Na, for example). Hence, MXene compositions
that exhibit larger ranges of transition-metal oxidation states have wider
Faradaic voltage windows. In this work we present predictions of the voltage
windows for MXene pseudocapacitor electrodes as a function of composition. We
highlight the effects of transition-metal valence and oxidation state on the
energy capacity.

While ruthenia is a pseudocapacitive material with promising energy density,
it has found limited commercial applications due to its high cost.13, 14 MXene
materials, Mn+1XnT2, are exempt from this constraint—being composed of earth-
abundant carbon (X = C, N) and oxygen (T = O, OH, F), and of relatively
inexpensive transition-metal elements (M). These transition metals provide a
wide compositional space, justifying the development of computational
descriptors to rapidly assess pseudocapacitive performance.15, 16 Electronic
and thermodynamic descriptors have been used to predict trends in energy
storage for MXene materials and have often been motivated by experimental
observations. The oxidation state of the transition metal is closely related
to the pseudocapacitive performance of MXenes.9, 17 Related quantities such as
the Bader charge of the transition metals, and the Fermi energy shift upon
adsorption of hydrogen have been used to describe charge storage trends.18
Here, we aim to make connections between charge transfer descriptors (such as
the transition-metal valence or Fermi energy shifts) and the active Faradaic
voltage windows of MXene electrodes.

A number of models have been used to predict the energy storage performance of
pseudocapacitive electrodes.19 Trends in Faradaic capacitance can be estimated
for different MXene compositions using density-functional theory (DFT).17
Local averaging methods can also be used to predict equilibrium
charge–discharge curves.20 These DFT models can provide electronic
descriptors, such as shifts in the potential of zero charge or metal oxidation
states, which reflect the ability to transfer electronic charge between the
adsorbed species and the pseudocapacitive surface.21 While the calculation of
thermodynamic quantities over a range of adsorbate coverages can be generally
accomplished with high accuracy in few-atom unit cells, the stability of
larger surfaces with a high degree of configurational disorder at finite
temperatures is rarely addressed. Additionally, DFT-based models generally
provide an incomplete representations of lateral adsorbate interactions due to
the limited length scales accessible to electronic-structure calculations.
Adsorption phenomena at low and high voltages are key to predicting the active
voltage windows, but are not accurately captured using DFT calculations alone.

A model that describes extended MXene electrodes at finite temperatures is
needed to predict charge–discharge characteristics in the low- and high-
voltage regimes. This can be accomplished with a two-dimensional surface
cluster expansion model on the adsorbate sublattice(s) that contribute to the
Faradaic charging. In the case of MXene pseudocapacitors, protic and aprotic
lithium-containing solvents can be used to achieve excellent pseudocapacitive
performance. By training a cluster expansion model with a dataset of DFT
calculations that include an implicit solvent representation combined with
double-layer charging, the effects of solvent and voltage can be accounted for
along extended surfaces. The inclusion of these environmental conditions
allows for a more realistic representation of the solid–liquid interface that
approaches the accuracy of DFT calculations in some cases. This computational
methodology is described and applied below.

The pseudocapacitive response of MXene electrodes is highly dependent on the
solvent environment. For MXenes in aqueous electrolytes, intercalated water
facilitates fast proton transfer to the electrode resulting in Faradaic
charging for proton-driven pseudocapacitance.22 For other ions such as Li+ and
Na+, it was demonstrated that whether or not the ions were solvated leads to
an electrochemical double-layer or pseudocapacitive response.23 In fact,
different solvents facilitate the desolvation of participating ions.24 Wang et
al. have shown that pseudocapacitance at MXene electrodes in propylene
carbonate solvents is dominated by charging processes involving desolvated
ions. Using representative implicit solvent calculations, a clear distinction
between electrochemical double-layer and Faradaic charging may be drawn.25 We
generate models in the regime that the response is dominated by near-surface
charging processes; these also serve as first-order approximations of
intercalation-based charging responses in the absence of inter-layer
interactions. In reality, the overall electrochemical response will be the
result of some combination of all charging processes. The predicted
pseudocapacitve responses of MXene surfaces in this work can help determine
which mechanisms contribute to charging, and may guide further experimentation
and modeling efforts to understand and optimize MXene pseudocapacitors.

##  2 Computational theory

To model the adsorption mechanisms that contribute to charge storage, we begin
by examining elementary surface processes. The charge stored per formula unit
due to Faradaic charging can be calculated directly based on the coverage of a
given species assuming integer charge transfer from the adsorbates:

$Q^{\rm{Faradaic}}=\sum_{i}N_{i}z_{i}$ (1)

where the sum runs over all adsorbing species, $i$. The product of the number
of an adsorbing ion, $N_{i}$, and the respective valency of the adsorbing ion,
$z_{i}$, gives the overall Faradaic charge stored. The equilibrium coverages
and reversible charge-discharge curves can be determined via thermodynamic
analysis.

In this approach, we first consider the adsorption of cations from the
solution. Specifically, we aim to model electrodes in contact with lithium-
and/or proton-containing solvents:

$\displaystyle*+{\rm H}^{+}+e^{-}\longrightarrow{\rm H}^{*}$ (2)
$\displaystyle*+{\rm Li}^{+}+e^{-}\longrightarrow{\rm Li}^{*}$

where $*$ indicates an adsorption site. The associated adsorption energies are
given as

$\displaystyle\Delta F({\rm H}^{*})=E({\rm H}^{*})-E_{\rm MXene}-N_{\rm
H}\frac{E_{\rm{H_{2}}}}{2}$ (3) $\displaystyle\Delta F({\rm Li}^{*})=E({\rm
Li}^{*})-E_{\rm MXene}-N_{{\rm Li}}E_{{\rm Li}}$

and can be calculated from first principles. These adsorption energies vary
with electrochemical conditions. We simulate these conditions using continuum
embedding methods, which will be discussed below. In simulated electrochemical
conditions, the free energies in Eq. (3) should be modified by including
(electro)chemical potentials:

$\displaystyle\Delta F({\rm H}^{*})^{\rm solv}=E({\rm H}^{*})-E_{\rm
MXene}-N_{\rm H}\big{(}\mu_{{\rm H}^{+}}+\mu_{e^{-}}\big{)}$ (4)
$\displaystyle\Delta F({\rm Li}^{*})^{\rm solv}=E({\rm Li}^{*})-E_{\rm
MXene}-N_{{\rm Li}}\big{(}\mu_{{\rm Li}^{+}}+\mu_{e^{-}}\big{)}$

where $\mu_{{\rm H}^{+}}$, $\mu_{{\rm Li}^{+}}$, and $\mu_{e^{-}}$ are the
electrochemical potentials of hydrogen ions, lithium ions, and the electrons
in the electrode, respectively. These chemical potentials are evaluated as

$\displaystyle\mu_{{\rm H}^{+}}=\frac{1}{2}E({\rm H}_{2})+k_{\rm{B}}T{\rm
ln}(a_{\rm H})+e\phi_{{\rm H}|{\rm H}^{+}}$ (5) $\displaystyle\mu_{{\rm
Li}^{+}}=E({\rm Li})+k_{\rm{B}}T{\rm ln}(a_{\rm Li})+e\phi_{\rm Li|Li^{+}}$
$\displaystyle\mu_{e^{-}}=E_{\rm{F}}$

using the free energies of hydrogen and lithium in their gaseous and metallic
states, $E({\rm{H}}_{2})$ and $E$(Li), respectively. The activity of the
hydrogen and lithium in the solution, $a_{\rm{H}}$ and $a_{\rm{Li}}$, are
taken to be 1 (under standard conditions). The standard reduction potentials
of hydrogen and lithium, $\phi_{{\rm{H}}|{\rm{H}}^{+}}$ and
$\phi_{{\rm{Li}}|{\rm{Li}}^{+}}$ on the absolute scale, multiplied by the
fundamental charge $e$, is the chemical work to oxidize the adsorbates. The
reduction potentials of hydrogen and lithium in water under standard
conditions are used; these are 3.04 and 0.0 V on the lithium scale,
respectively. The absolute electrochemical potential of the electron is equal
to the Fermi energy, $E_{\rm{F}}$ (by definition).

Depending on the electrode voltage, the equilibrium surface coverage—and the
stored Faradaic charge—on the MXene surface will change. We account for this
effect by expanding the adsorption energies in Eq. (4) with respect to the
charge

$\Delta F^{\rm{ads}}(q)=\Delta F^{\rm{ads}}(0)+q\phi_{\rm
pzc}+\frac{q^{2}}{2{\mathscr{A}}C_{\rm{dl}}}+...$ (6)

where the charge on the surface is $q$, $\mathscr{A}$ is the surface area, and
$C_{\rm dl}$ is the differential double layer capacitance of the MXene
electrode. The terms of this expression can be calculated using continuum
embedding models. The potential at which there is zero excess charge on the
surface is the potential of zero charge, $\phi_{\rm pzc}$, which is given by
$-E_{\rm{F}}/e$. A voltage-dependent adsorption energy can be obtained by
truncating this expansion at second order and taking the Legendre transform
with respect to charge:

$\Delta\mathscr{F}^{\rm{ads}}(\phi)=\Delta F^{\rm{ads}}(q)-q\phi$ (7)

where $\mathscr{F}^{\rm{ads}}$ is the transformed free energy and $\phi$ is
the electrode voltage. In the calculation of this voltage-dependent adsorption
energy, we note that the electron chemical potential in Eq. (4) is set by the
electrode voltage $\mu_{e^{-}}=-e\phi$.

The prediction of voltage-dependent adsorption energies in Eq. (7) using DFT
models is limited to the length scales accessible. This often results in
coarse samplings of ion coverages, and models produced from these do not often
capture longer range lateral interactions between adsorbates. Additionally,
the convex hull is densely populated with adsorption configurations near the
lowest energy ground state, as seen in Fig. S2 of the Supporting Information.
In order to accurately predict adsorption phenomena as a function of voltage,
statistical ensemble models of extended surfaces are required. Cluster
expansions and Monte Carlo simulations can be used to build such models.

Figure 2: Construction of the cluster basis on the hexagonal MXene surface
(left) for the hydrogen sublattice (center) and lithium sublattice (right).

The energies in Eq. (4) vary as a function of the coverage of the respective
ions, and the nearly continuous set of adsorption states on the electrode
surface contributes to the pseudocapacitive response described in Eq. (1).
Therefore statistical ensemble representations of the electrodes and inclusion
of lateral adsorbate interactions are important to accurately predict the
pseudocapacitance. It is possible to represent scalar extensive functions of
configuration as a linear expansion with the basis functions, following
Sanchez and coworkers.26 The adsorption energy on a single sublattice can be
represented as an expansion in this spectral basis,

$F(\boldsymbol{\sigma})=\sum_{\alpha}j_{\alpha}\Phi_{\alpha}(\boldsymbol{\sigma})$
(8)

where $\boldsymbol{\sigma}$ is a vector of occupational variables giving the
occupancy of each site on the sublattice (e.g. –1 for vacancy and +1 for
lithium ions). The sum in the expansion runs over clusters, $\alpha$, which
are illustrated in Fig. 2. There is an expansion coefficient associated with
each cluster, $j_{\alpha}$, and the spectral basis function for this cluster,
$\Phi_{\alpha}(\boldsymbol{\sigma})$. This expansion can be used to predict
the adsorption energy of configurations for extended surfaces that are more
representative of electrode surfaces. In practice, the cluster basis functions
in Eq. 8 are averaged over the crystal and the same expansion coefficient is
used for symmetrically equivalent clusters.27

We can also expand the voltage-dependent adsorption energy to obtain an energy
expression for arbitrary adsorption configurations on a single adsorption
sublattice as a function of voltage:28

$\mathscr{F}(\boldsymbol{\sigma},\phi)=\sum_{\alpha}J_{\alpha}(\phi)\,\Phi_{\alpha}(\boldsymbol{\sigma})$
(9)

In this voltage-dependent cluster expansion, the coefficients are defined at
each voltage. These can be obtained by training against a dataset of
calculated voltage-dependent adsorption energies, as described in Eq. (7).

The co-adsorption of multiple ions cannot be modeled with Eq. (8) and Eq. (9)
if the adsorption sites are on different, non-exchanging sublattices.
Adsorption configurations on one sublattice need to be designated with a
vector of spin variables, $\boldsymbol{\sigma}$. The adsorption configurations
on the sublattice corresponding to the other sites need to be specified with
another vector of spin variables, $\boldsymbol{\delta}$, and the
multicomponent cluster expansion formalism developed by Tepesch, Garbulsky,
and Ceder should be used.29 We adapt this approach to also include voltage
dependence:

$\mathscr{F}(\boldsymbol{\sigma},\boldsymbol{\delta},\phi)=\sum_{\alpha\beta}V_{\alpha\beta}(\phi)\,\Theta_{\alpha\beta}(\boldsymbol{\sigma},\boldsymbol{\delta})$
(10)

This expansion is facilitated by the construction of basis functions of the
coupled sublattices. These are constructed by taking the tensor product of the
basis functions {$\Phi_{\alpha}$} for each sublattice forming the basis
functions for the coupled sublattices,
$\Theta_{\alpha\beta}(\boldsymbol{\sigma},\boldsymbol{\delta})$ (e.g. all
possible products of clusters in Fig. 2). Thus there are two cluster indices,
$\alpha$ and $\beta$. The expansion coefficients, $V_{\alpha\beta}(\phi)$,
include effective cluster interactions (ECIs) from single-sublattice clusters
as well as multiple-sublattice clusters. Using this formalism, a method for
modeling the simultaneous adsorption of hydrogen and lithium is achieved.

The expansion coefficients include the effects of solvent polarization and the
diffuse double layer when fit to DFT calculations of surface cells in implicit
solvents. The MXene electrode systems are embedded in a dielectric cavity and
a smooth Poisson–Boltzmann density of countercharges set at a fixed distance
of 3.5 Å away from the MXene surface using a planar interface.30 This
interface function is chosen to reduce excessive screening between adsorbates.
At the boundary, a planar cavity function for the continuum dielectric and
implicit ion charges, $s(\mathbf{r})$, is used with the dielectric function,
$\epsilon(\mathbf{r})=\exp({\rm ln}(\epsilon_{0})[1-s(\mathbf{r})])$, where
$\epsilon_{0}$ is the dielectric constant of the bulk solvent.31 The
dielectric was specified here as 12 for an approximate representation of
solvent polarization in simple organic solvents. The charge density for the
$i^{\rm th}$ ion, is given by the concentration in the solvent, $c_{i}$,
multiplied by a Boltzmann weight. A sum over implicit ion types gives the
total density,
$\rho_{\rm{PB}}(\mathbf{r})=\sum_{i}(1-s(\mathbf{r}))c_{i}\exp\left({-({z_{i}e\phi(\mathbf{r})})/({k_{\rm
B}T})}\right)$, and contributions to the DFT energy are solved for self-
consistently.31 Examples of the resulting implicit Poisson–Boltzmann
distribution of ions and bound polarization charges for the solvent are
provided in Fig. S1 of the Supporting Information.

##  3 Computational methods

Cluster expansion models were parameterized for $\rm{M_{2}CO_{2}}$ and
$\rm{M_{3}C_{2}O_{2}}$ MXene electrodes where M is an early transition metal:
Cr, Hf, La, Mo, Nb, Sc, Ta, Ti, V, W, Y, Zr. The adsorption of lithium and
hydrogen on the MXene surfaces were considered independently to simulate
lithium-containing and protic solvents. Single-ion adsorption models were made
using Eq. (9). In protic, lithium-containing solvents, it is possible for
lithium and hydrogen to co-adsorb. Hydrogen prefers the oxygen-top sites and
lithium generally prefers the hexagonal closed-packed (hollow) site.

The cluster expansion models were fit using the Clusters Approach for
Statistical Mechanics (CASM) code 27, with automated Python wrappers to
incorporate the voltage dependence into the cluster expansions. All of the
symmetrically distinct pair clusters within 3 adsorption sites, triplets
within 3 adsorption sites, and quadruplet clusters within 2 adsorption sites
were used in the cluster expansion models along with the identity and single
site clusters. Discrete Chebychev site basis functions were used to construct
the cluster basis functions. Optimal sets of expansion coefficients are
determined at each voltage using $l_{1}$-regularized least squares fit against
DFT data. The test errors of the fits are estimated with $k$-fold ($k$ = 10)
cross validation along with cluster elimination as implemented in CASM. The
test errors are reported in the Supporting Information in Fig. S6.

The DFT training sets were comprised of 184 and 53 symmetrically distinct
adsorption configurations for multiple-ion and single-ion adsorption systems,
respectively. Semi-local DFT calculations were performed at the GGA level in
Quantum ESPRESSO to calculate the energies of the structures in the fitting
set.32 Plane wave basis sets were used with a kinetic energy cutoff of 120 Ry,
and norm-conserving pseudopotentials were employed to represent atomic
cores.33 The Brillouin zone was sampled in all cells using a (11/$n$ $\times$
11/$m$ $\times$ 1) Monkhorst-pack k-point mesh where $n$ and $m$ are the
multiples of the primitive surface cell vectors. A cold smearing value of 0.01
Ry was used to smooth electronic occupations.34 All fitting cells were
embedded in a 15 Å continuum dielectric and implicit ion distribution to
represent the solvent (cf. the previous section).

The implicit ion distribution in our continuum solvent model allows for a DFT-
level approximation of the differential diffuse layer capacitance. Explicit
charge applied to the MXene electrode is screened by the ion distribution in
the solvent region. The energy as a function of charge can be used to fit the
quadratic form of Eq. (6), where the second-order term is related to the
diffuse double-layer capacitance, $C_{\rm dl}$. This double layer capacitance
varies with adsorption configuration with an average value of 8.6 $\mu$F/cm2
in a test MXene composition, Mo3C2O2. It was found that using a constant
diffuse double layer capacitance, 10 $\mu$F/cm2, close to the average yields
similar results. A constant capacitance of 10 $\mu$F/cm2 was used for all
MXene compositions to evaluate voltage-dependent adsorption energies in Eq.
(9).

Monte Carlo simulations were performed using the CASM software package.
Simulations for hydrogen or lithium ion adsorption on MXene electrodes were
performed in the grand-canonical ensemble on 30 $\times$ 30 surface cells. In
these simulations, the steady-state coverage of lithium or hydrogen adsorbates
was determined at each voltage in 10 separate simulations and averaged; the
resulting adsorption isotherms are reported in Figs. S7-S10 of the Supporting
Information with an example highlighted for $\rm{Ti}_{2}\rm{C}\rm{O}_{2}$ in
Fig. 3. To obtain these isotherms, each Monte Carlo simulation was ran for
1,000 passes after the system had been equilibrated, and were initialized with
a random adsorption configuration. The steady-state coverage of lithium ions
as a function of voltage allows us to determine the active pseudocapacitive
voltage window. This is the voltage window where the slope of the adsorption
isotherm is non-zero and the ion coverage is between 0 and 1 monolayers.
Additionally, the slopes of the isotherms are directly related to the
pseudocapacitance of the MXene electrodes.

Figure 3: Predicted charge–voltage response of the a
$\rm{Ti}_{2}\rm{C}\rm{O}_{2}$ MXene electrode from Monte Carlo simulations.
The blue curve is the charge stored due to diffuse electrochemical double-
layer charging, the red curve is the Faradaic charge stored, and the black
curve is the net charge stored. The decoupled pseudocapactive charge storage
response is used to predict the active voltage window of the electrode
indicated by the horizontal lines.

In lithium-containing protic solvents, the effect of hydrogen adsorption at
small coverages can help extend the active voltage window for lithium. To
predict the extent of lithium stabilization in the presence of other ions, we
simulate lithium and hydrogen adsorption on MXene electrodes in low coverage
regimes. These simulations were also performed on 30 $\times$ 30 surface cells
to represent the electrode with a single Monte Carlo simulation (10,000
passes) at each voltage. The cluster expansions fit at each voltage include
small cells with both hydrogen and lithium on the surface. Using the
multilattice cluster expansion models in Eq. (10) for fixed composition Monte
Carlo simulations, we predict the stabilization of lithium adsorption in the
upper voltage limit. Solving for the voltages with zero formation energy in
Eq. (7) in regimes with small lithium and small hydrogen ion adsorption can
reveal the extent to which the introduction of hydrogen may favor lithium
adsorption.

##  4 Results and discussion

The equilibrium charge–discharge curves (Fig. 3) were obtained for all M3C2O2
and M2CO2 early-transition-metal compositions. From these, the voltage windows
with pseudocapacitive activity are calculated as the range of voltages where
the derivative of the Faradaic charge with respect to voltage is non-zero and
the coverage of the adsorbing ion is between but not equal to 0% and 100%.
These voltage windows are compiled in Fig. 4a for the different MXene
compositions. The Faradaic capacitances, were evaluated by taking the
derivative of the Faradaic charge with respect to voltage over the active
window and averaged. We note that an instantaneous capacitance could be
obtained by a simple numerical derivative of this charge–voltage curve, but to
be accurate, this would require a much finer numerical grid than is sampled.
The average capacitances are reported in Table S1 of the Supporting
Information along with the adsorption isotherms used to calculate them. The
predicted isotherms are consistent with other first-principles based models of
MXene electrode surfaces.17, 21 Relatively linear decreases in ion coverages
as a function of voltage correspond to Faradaic charging and the
characteristic rectangular cyclic voltammetry response for capacitors. The
equilibrium charge/discharge curves can be contrasted with similar simulations
of lithium-ion battery materials, which are associated with rapid decrease in
Li content over a short voltage range, to show that many MXene compositions
yield pseudocapacitive charge/discharge behavior.35 Pseudocapacitive responses
and voltage windows were highlighted for the group-VI MXene electrode surfaces
in simulated organic solvent environments. Improved accuracy is expected over
existing models with the voltage-dependent cluster expansion method at high
and low ion coverages because of the third-nearest neighbor lateral adsorbate
interactions trained in the models.

Figure 4:  (a) Pseudocapacitive voltage windows for oxygen terminated MXenes
determined from adsorption isotherms at fixed voltage. Lithium voltage windows
are given in blue and hydrogen voltage windows in red. (b) Average voltage
window extension calculated by evaluating the average $\phi(\mathscr{F}=0)$
for a range of low H and Li coverages, giving the upper voltage limit. The
lower voltage limit is $\phi(\mathscr{F}=0)$ predicted by Monte Carlo
simulations without H present.

The size of the voltage windows increases from left to right in the periodic
table. The voltage windows tend to decrease with increasing period, except
with the group-VI M3C$2$ MXenes. An ideal pseudocapacitor has a constant
Faradaic capacitance over a large voltage window while desorption over a small
voltage range is more characteristic of a battery electrode. Thus, Zr and Hf
are not suitable lithium-based pseudocapacitor candidates, but may perform
better as lithium-ion battery electrodes. The active voltage windows for
hydrogen are higher and closer to the standard redox potential of hydrogen.
For some transition-metal compositions, the active voltage windows of hydrogen
and lithium overlap, indicating that in protic, lithium-containing solvents
there may be voltage ranges where both ions contribute to energy storage. We
draw attention to the group-VI MXene compositions in Fig. 4a; these
compositions have large active windows for both hydrogen and lithium
pseudocapacitance.

The pseudocapacitive windows of lithium and hydrogen overlap, suggesting that,
at some potentials, there is simultaneous adsorption of hydrogen and lithium.
While adsorption of large amounts of lithium and hydrogen were predicted to
lead to significant reconstructions for all MXenes, small amounts of lithium
and hydrogen, corresponding to the upper voltage lithium voltage window, are
stable. The hydrogen partially stabilizes the adsorption of lithium at higher
voltages. This observation is supported by fixed composition Monte Carlo
simulations with small amounts of hydrogen $\theta_{\rm H}=$
$\nicefrac{{1}}{{9}}$, $\nicefrac{{2}}{{9}}$, and $\nicefrac{{3}}{{9}}$ and
small amounts of lithium, $\theta_{\rm Li}=$ $\nicefrac{{1}}{{10}}$. With
these simulations, the voltage at which adsorption is stable can be obtained
with Eq. (7) and the average extension of the Li voltage window can be
evaluated for both M3C2 and M2C MXenes. The average voltage window extension
upon introduction of hydrogen is given for each MXene composition in Fig. 4b.
The effect is larger for M2C MXenes and the voltage window tends to be larger
for earlier transition metals. For the W3C2O2 MXene, the voltage window
extension due to multiple ion adsorption is negligible. Lithium adsorption by
itself leads to reconstructions on lanthanum-based MXenes. Lithium adsorption
is stabilized to a small extent with the introduction of hydrogen, and is
indicated by the red bars in Fig. 4b. The stabilization of Li ion adsorption
at high voltages is due to increased interactions between Li and the
transition metals in the presence of H adsorbates. This observation can be
explained by examining the density of states reported the Supporting
Information, which reveal increased charge transfer to the MXene transition
metal in the presence of multiple ion species.

Figure 5: (a) Pseudocapacitive energy density for lithium adsorption on M3C2O2
MXene as a function of Lowdin d-orbital population. The d-orbital filling is
calculated on MXene slabs and is proportional to the transition-metal valence.
The energy density increases with d-orbital filling and follows row and column
periodic table trends. (b) The pseudocapacitive voltage window size increases
with the scaled change in oxidation states during adsorption induced redox on
MXene surfaces with the M3C2O2 compositions.

The pseudocapacitive energy storage is calculated using the average Faradaic
capacitances $C_{{\rm{Faradaic}}}$ obtained in Table S1 of the Supporting
Information and the active voltage windows $\Delta V$ in Fig. 3 with the
equation $E_{\rm{pseudo}}=\frac{1}{2}C_{{\rm Faradaic}}\Delta V^{2}$. We
report the lithium pseudocapacitive energy densities for
$\rm{M_{3}C_{2}O_{2}}$ MXenes in Fig. 5a in the presence of lithium-containing
solvents. The trends for $\rm{M_{2}CO_{2}}$ MXenes and hydrogen-containing
solvents are similar and are included in the Supporting Information. We note
that these pseudocapacitive energies are decoupled from the electrochemical
double-layer energy storage. These predictions of MXene energy storage
capabilities show that the pseudocapacitive energy stored correlates with the
valence of the transition metal. This trend can be explained by the fact that
redox reactions involving the transition metal are primarily responsible for
the pseudocapacitive response in MXenes. There are generally more oxidation
states available for late-transition-metal (group-V and group-VI) MXenes, and
this facilitates Faradaic charging from adsorbed ions and larger
pseudocapacitive windows. We confirm this with relative ranges of oxidation
states for M3C2O2 MXenes in Fig. 5b. Larger ranges of oxidation states
generally indicate larger pseudocapactive windows, but magnetic MXenes, V and
Nb-based, follow more irregular trends. The V- and Nb-based MXenes deviate
most from linear increases in voltage window size with oxidation state range;
a thorough study of the effects of magnetism on pseudocapacitive charge
storage would be beneficial for these systems.

It has been shown that the Bader charge is linearly proportional to
experimental oxidation states in MXenes, but the constants of proportionality
may be different from composition to composition.17 For this reason, the
oxidation state ranges in Fig. 5b were predicted as relative changes in Bader
charge upon redox with Li+. The changes in Bader charge for the surface
transition metals, reported in the Supporting Information, were multiplied by
a scaling factor, $f=\Delta\phi_{\rm{pzc}}({\rm M})/\Delta\phi_{\rm{pzc}}({\rm
Ti})$. Other scaling factors could be used to describe relative ranges of
oxidation states, but the ratio $f$ bears valuable physical meaning. It
describes the change in the affinity of the MXene surface for charged species
relative to another and is an important, reliable descriptor for the Faradaic
capacity.18 The combination of these descriptors help describe the trend found
with group-VI MXenes that when the MXene transition metals access larger
ranges of oxidation states, larger windows of Faradaic charging are observed.
We note that the group-VI $\rm{M_{3}C_{2}O_{2}}$ MXenes are not typically
stable due to the transition-metal layer stacking preference. Group-VI metals
prefer the ABAB stacking but possess ABCABC stacking in the
$\rm{M_{3}C_{2}O_{2}}$ compositions.16 Late transition metals can be
stabilized in $\rm{M_{3}C_{2}O_{2}}$ MXenes by forming ordered alloys with
other transition metals such as Ti. We show in Fig. 5b that one such MXene
composition, $\rm{W_{2}TiC_{2}O_{2}}$, accesses a wide range of oxidation
states as well, and is expected to have a large pseudocapacitive voltage
window.

The search for new MXene pseudocapacitor electrodes should consider MXene
compositions that allow for large ranges of transition-metal oxidation states,
and solvent environments that facilitate them. For different solvent,
electrolyte, and electrode combinations, different charging processes will
dominate the electrochemical response.23 The voltage-dependent cluster
expansion models used in this work correctly represent MXene-electrode systems
where surface charging dominates, highly exfoliated MXenes and systems where
intercalation charging is expected to contribute less. Though our results
serve as an approximate model of intercalation charging in the absence of
inter-layer interactions, they also highlight systems where Li+ intercalation
charging effects are important and would benefit from further direct
investigation, such as with heavy group-V $\rm{M_{2}C}$ MXenes given the small
predicted voltage windows.

Considering alloy MXenes as pseudocapacitive electrodes beyond the
${\rm{W_{2}TiC_{2}O_{2}}}$ composition will be important. The presence of
other transition metals or X elements could induce varied oxidation states,
depending on chemical ordering or lack thereof, that enhance pseudocapacitive
energy storage through a combination of larger voltage windows or capacity.15
Increased pseudocapacitive voltage windows have been observed experimentally.
It was confirmed by Vahidmohammadi et al. that $\rm{V_{2}CO_{2}}$ MXenes have
larger pseudocapacitive voltage windows compared to the early-transition-metal
or Ti MXene electrodes.36, 1 There are promising results for group VI (Mo-
based) MXenes in similar conditions.37 The key effect of transition-metal
valence and the range of induced oxidation states in the MXene transition
metals on the MXene pseudocapacitance could be confirmed with x-ray adsorption
spectroscopy experiments, but are qualitatively anticipated given the wider
range of oxidation states allowed for vanadium relative to titanium.17

##  5 Conclusion

Voltage-dependent cluster expansion models were used to predict Faradaic
charge storage in MXenes across various transition-metal compositions. We
found that group-VI MXenes may be promising pseudocapacitor candidates due to
large pseudocapacitive voltage windows that lead to high energy storage
densities. The predicted pseudocapacitive voltage windows only describe
voltages where the MXene electrodes will be active with Faradaic charging, but
do not provide information about the breakdown voltage of the MXene electrode.
These breakdown voltages will depend on the solvent, the ionic species, and
the composition of the electrodes. A study of the electrochemical durability
of MXene electrodes would show how much of the Faradaically active window is
accessible for different solvent–electrode combinations. The Faradaic voltage
windows for lithium adsorption on these MXenes may be accessible with
propylene carbonate (PC) or other carbonate solvents. The voltage windows of
MXenes can possibly be extended by introducing multiple cation species as
demonstrated with model protic, lithium-containing solvents. The voltage
window extension is the most pronounced for $\rm{M_{2}CO_{2}}$ early-
transition-metal MXenes. Though these MXene compositions are often not
considered for lithium-ion battery electrode applications due to relatively
slow ion intercalation kinetics, these results suggest that some of the early-
transition-metal MXenes such as $\rm{Ti_{2}CO_{2}}$ would still be excellent
pseudocapacitor electrodes. The relative positions and sizes of voltage
windows for the MXene compositions could assist the selection of MXenes for
hybrid and asymmetric supercapacitor devices. The Zr and Hf M2C MXenes may be
more suitable in hybrid supercapacitor devices, while compositions with more
psuedocapacitor-like responses may be applicable to asymmetric supercapacitor
devices.

We have determined that the valence of the transition metals in MXenes
correlates with the energy storage. Faradaic ion adsorption induces larger
ranges of oxidation states of group-VI MXenes compared to early-transition-
metal compositions such as $\rm{Ti_{3}C_{2}O_{2}}$. We confirm that the
oxidation state of the transition metals in MXene electrodes is strongly
related to the pseudocapacitance, but we also show that the range of oxidation
states induced in the transition metals yields larger pseudocapacitive windows
as predicted by voltage-dependent cluster expansion models. The relative
ranges in Bader charges serve as a computationally tractable and readily
accessible descriptor for predicting MXene compositions with large
pseudocapacitive voltage windows. On the basis of these trends, group-VI alloy
MXenes should have large pseudocapacitive (lithium- or hydrogen-mediated)
energy storage densities. Additionally, the introduction of hydrogen
adsorbates to lithiated electrodes surfaces leads to an increased
hybridization of the lithium s-orbitals with the transition-metal d-orbitals.
Not only is the active lithium adsorption window extended by the introduction
of hydrogen ions, but the transition-metal oxidation state is modulated more
strongly compared to aprotic environments. These predictions suggest that
improved pseudocapacitance can be achieved with MXene electrodes in solvents
containing multiple ions that participate in Faradaic charging.

##  6 Supporting Information

Isocontours of bound charge densities and electrolyte densities near MXene
electrodes, ion adsorption convex hull as a function of voltage,
pseudocapacitive energy densities for other MXene compositions, Bader charge
ranges for M3C2O2 MXenes, electronic-structure changes of MXenes upon
adsorption of multiple ions, cross-validation scores for cluster expansion
models, predicted areal pseudocapacitances, predicted adsorption isotherms for
all MXene compositions

##  7 Acknowledgements

J. G. and I. D. acknowledge primary financial support from the U.S. Department
of Energy, Office of Science, Basic Energy Sciences, CPIMS Program, under
Award No. DE-SC0018646. F. V. and N. K. acknowledge financial support from
Murata Manufacturing Co., Ltd. and the Center for Dielectrics and
Piezoelectrics. A portion of this work was performed by N. K. under the
auspices of the DOE by Lawrence Livermore National Laboratory (LLNL) under
Contract No. DE-AC52-07NA27344. J. G. acknowledges support and training
provided by the Computational Materials Education and Training (CoMET) NSF
Research Training program under grant number DGE-1449785. Calculations for
this research were performed on Roar supercomputer of the Institute for
Computational and Data Sciences at the Pennsylvania State University.

## References

  * Okubo et al. 2018 Okubo, M.; Sugahara, A.; Kajiyama, S.; Yamada, A. MXene as a Charge Storage Host. _Accounts of Chemical Research_ 2018, _51_ , 591–599
  * Seh et al. 2016 Seh, Z. W.; Fredrickson, K. D.; Anasori, B.; Kibsgaard, J.; Strickler, A. L.; Lukatskaya, M. R.; Gogotsi, Y.; Jaramillo, T. F.; Vojvodic, A. Two-Dimensional Molybdenum Carbide (MXene) as an Efficient Electrocatalyst for Hydrogen Evolution. _ACS Energy Letters_ 2016, _1_ , 589–594
  * Gogotsi and Anasori 2019 Gogotsi, Y.; Anasori, B. The Rise of MXenes. _ACS Nano_ 2019, _13_ , 8491–8494
  * Lukatskaya et al. 2013 Lukatskaya, M. R.; Mashtalir, O.; Ren, C. E.; Dall’Agnese, Y.; Rozier, P.; Taberna, P. L.; Naguib, M.; Simon, P.; Barsoum, M. W.; Gogotsi, Y. Cation Intercalation and High Volumetric Capacitance of Two-Dimensional Titanium Carbide. _Science_ 2013, _341_ , 1502–1505
  * Naguib et al. 2011 Naguib, M.; Kurtoglu, M.; Presser, V.; Lu, J.; Niu, J.; Heon, M.; Hultman, L.; Gogotsi, Y.; Barsoum, M. W. Two-Dimensional Nanocrystals Produced by Exfoliation of Ti3AlC2. _Advanced Materials_ 2011, _23_ , 4248–4253
  * Naguib et al. 2012 Naguib, M.; Come, J.; Dyatkin, B.; Presser, V.; Taberna, P.-L.; Simon, P.; Barsoum, M. W.; Gogotsi, Y. MXene: a promising transition metal carbide anode for lithium-ion batteries. _Electrochemistry Communications_ 2012, _16_ , 61–64
  * Er et al. 2014 Er, D.; Li, J.; Naguib, M.; Gogotsi, Y.; Shenoy, V. B. Ti3C2 MXene as a High Capacity Electrode Material for Metal (Li, Na, K, Ca) Ion Batteries. _ACS Applied Materials & Interfaces_ 2014, _6_ , 11173–11179
  * Choudhary et al. 2017 Choudhary, N.; Li, C.; Moore, J.; Nagaiah, N.; Zhai, L.; Jung, Y.; Thomas, J. Asymmetric Supercapacitor Electrodes and Devices. _Advanced Materials_ 2017, _29_ , 1605336
  * Lukatskaya et al. 2015 Lukatskaya, M. R.; Bak, S.-M.; Yu, X.; Yang, X.-Q.; Barsoum, M. W.; Gogotsi, Y. Probing the Mechanism of High Capacitance in 2D Titanium Carbide Using In Situ X-Ray Absorption Spectroscopy. _Advanced Energy Materials_ 2015, _5_ , 1500589
  * Wang et al. 2015 Wang, X.; Kajiyama, S.; Iinuma, H.; Hosono, E.; Oro, S.; Moriguchi, I.; Okubo, M.; Yamada, A. Pseudocapacitance of MXene nanosheets for high-power sodium-ion hybrid capacitors. _Nature Communications_ 2015, _6_ , 6544
  * Li et al. 2017 Li, J.; Yuan, X.; Lin, C.; Yang, Y.; Xu, L.; Du, X.; Xie, J.; Lin, J.; Sun, J. Achieving High Pseudocapacitance of 2D Titanium Carbide (MXene) by Cation Intercalation and Surface Modification. _Advanced Energy Materials_ 2017, _7_ , 1602725
  * Fleischmann et al. 2020 Fleischmann, S.; Mitchell, J. B.; Wang, R.; Zhan, C.; Jiang, D.-e.; Presser, V.; Augustyn, V. Pseudocapacitance: from fundamental understanding to high power energy storage materials. _Chemical Reviews_ 2020, _120_ , 6738–6782
  * Zheng et al. 1995 Zheng, J. P.; Cygan, P. J.; Jow, T. R. Hydrous Ruthenium Oxide as an Electrode Material for Electrochemical Capacitors. _Journal of The Electrochemical Society_ 1995, _142_ , 2699
  * Emsley 2011 Emsley, J. Nature’s building blocks: an AZ guide to the elements. 2011
  * Tan et al. 2017 Tan, T. L.; Jin, H. M.; Sullivan, M. B.; Anasori, B.; Gogotsi, Y. High-Throughput Survey of Ordering Configurations in MXene Alloys Across Compositions and Temperatures. _ACS Nano_ 2017, _11_ , 4407–4418
  * Anasori et al. 2017 Anasori, B.; Lukatskaya, M. R.; Gogotsi, Y. 2D metal carbides and nitrides (MXenes) for energy storage. _Nature Reviews Materials_ 2017, _2_ , 1–17
  * Zhan et al. 2018 Zhan, C.; Naguib, M.; Lukatskaya, M.; Kent, P. R. C.; Gogotsi, Y.; Jiang, D.-e. Understanding the MXene Pseudocapacitance. _The Journal of Physical Chemistry Letters_ 2018, _9_ , 1223–1228
  * Zhan et al. 2019 Zhan, C.; Sun, W.; Kent, P. R. C.; Naguib, M.; Gogotsi, Y.; Jiang, D.-e. Computational Screening of MXene Electrodes for Pseudocapacitive Energy Storage. _The Journal of Physical Chemistry C_ 2019, _123_ , 315–321
  * Ozoliņš et al. 2013 Ozoliņš, V.; Zhou, F.; Asta, M. Ruthenia-Based Electrochemical Supercapacitors: Insights from First-Principles Calculations. _Accounts of Chemical Research_ 2013, _46_ , 1084–1093
  * Keilbart et al. 2017 Keilbart, N.; Okada, Y.; Feehan, A.; Higai, S.; Dabo, I. Quantum-continuum simulation of the electrochemical response of pseudocapacitor electrodes under realistic conditions. _Physical Review B_ 2017, _95_ , 115423
  * Okada et al. 2019 Okada, Y.; Keilbart, N.; Goff, J. M.; Higai, S.; Shiratsuyu, K.; Dabo, I. MXene Electrode Materials for Electrochemical Energy Storage: First-Principles and Grand Canonical Monte Carlo Simulations. _MRS Advances_ 2019, _4_ , 1833–1841
  * Lukatskaya et al. 2017 Lukatskaya, M. R.; Kota, S.; Lin, Z.; Zhao, M.-Q.; Shpigel, N.; Levi, M. D.; Halim, J.; Taberna, P.-L.; Barsoum, M. W.; Simon, P.; Gogotsi, Y. Ultra-high-rate pseudocapacitive energy storage in two-dimensional transition metal carbides. _Nature Energy_ 2017, _2_ , 1–6
  * Gao et al. 2020 Gao, Q.; Sun, W.; Ilani-Kashkouli, P.; Tselev, A.; Kent, P. R. C.; Kabengi, N.; Naguib, M.; Alhabeb, M.; Tsai, W.-Y.; Baddorf, A. P.; Huang, J.; Jesse, S.; Gogotsi, Y.; Balke, N. Tracking ion intercalation into layered Ti3C2 MXene films across length scales. _Energy Environ. Sci._ 2020, _13_ , 2549–2558
  * Wang et al. 2019 Wang, X.; Mathis, T. S.; Li, K.; Lin, Z.; Vlcek, L.; Torita, T.; Osti, N. C.; Hatter, C.; Urbankowski, P.; Sarycheva, A.; Tyagi, M.; Mamontov, E.; Simon, P.; Gogotsi, Y. Influences from solvents on charge storage in titanium carbide MXenes. _Nature Energy_ 2019, _4_ , 241–248
  * Ando et al. 2020 Ando, Y.; Okubo, M.; Yamada, A.; Otani, M. Capacitive versus Pseudocapacitive Storage in MXene. _Advanced Functional Materials_ 2020, 2000820
  * Sanchez et al. 1984 Sanchez, J. M.; Ducastelle, F.; Gratias, D. Generalized cluster description of multicomponent systems. _Physica A: Statistical Mechanics and its Applications_ 1984, _128_ , 334–350
  * Puchala and Van der Ven 2013 Puchala, B.; Van der Ven, A. Thermodynamics of the Zr-O system from first-principles calculations. _Physical Review B_ 2013, _88_ , 094108
  * Weitzner and Dabo 2017 Weitzner, S. E.; Dabo, I. Voltage-dependent cluster expansion for electrified solid-liquid interfaces: Application to the electrochemical deposition of transition metals. _Physical Review B_ 2017, _96_ , 205134
  * Tepesch et al. 1995 Tepesch, P. D.; Garbulsky, G. D.; Ceder, G. Model for Configurational Thermodynamics in Ionic Systems. _Physical Review Letters_ 1995, _74_ , 2272–2275
  * Andreussi et al. 2019 Andreussi, O.; Hörmann, N. G.; Nattino, F.; Fisicaro, G.; Goedecker, S.; Marzari, N. Solvent-Aware Interfaces in Continuum Solvation. _Journal of Chemical Theory and Computation_ 2019, _15_ , 1996–2009
  * Nattino et al. 2018 Nattino, F.; Truscott, M.; Marzari, N.; Andreussi, O. Continuum models of the electrochemical diffuse layer in electronic-structure calculations. _The Journal of Chemical Physics_ 2018, _150_ , 041722
  * Giannozzi et al. 2009 Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Corso, A. D.; Gironcoli, S. d.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. _J. Phys.: Condens. Matter_ 2009, _21_ , 395502
  * van Setten et al. 2018 van Setten, M. J.; Giantomassi, M.; Bousquet, E.; Verstraete, M. J.; Hamann, D. R.; Gonze, X.; Rignanese, G. M. The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table. _Computer Physics Communications_ 2018, _226_ , 39–54
  * Marzari et al. 1999 Marzari, N.; Vanderbilt, D.; De Vita, A.; Payne, M. C. Thermal Contraction and Disordering of the Al(110) Surface. _Physical Review Letters_ 1999, _82_ , 3296–3299
  * Wolverton and Zunger 1998 Wolverton, C.; Zunger, A. First-principles prediction of vacancy order-disorder and intercalation battery voltages in Li x CoO 2. _Physical review letters_ 1998, _81_ , 606
  * VahidMohammadi et al. 2019 VahidMohammadi, A.; Mojtabavi, M.; Caffrey, N. M.; Wanunu, M.; Beidaghi, M. Assembling 2D MXenes into highly stable pseudocapacitive electrodes with high power and energy densities. _Advanced Materials_ 2019, _31_ , 1806931
  * Tao et al. 2017 Tao, Q.; Dahlqvist, M.; Lu, J.; Kota, S.; Meshkian, R.; Halim, J.; Palisaitis, J.; Hultman, L.; Barsoum, M. W.; Persson, P. O.; Rosen, J. Two-dimensional Mo 1.33 C MXene with divacancy ordering prepared from parent 3D laminate with in-plane chemical ordering. _Nature communications_ 2017, _8_ , 1–7