#  DESY 19-104, LTH 1212 August 2019

Soft corrections to inclusive DIS at four loops and beyond

Theory Group, Deutsches Elektronen-Synchrotron (DESY), D-22607 Hamburg,
Germany
E-mail  Sven-Olaf Moch g
II. Institute for Theoretical Physics, Hamburg University, D-22761 Hamburg,
Germany
E-mail<EMAIL_ADDRESS>Andreas Vogt
Department of Mathematical Sciences, University of Liverpool, Liverpool L69
3BX, United Kingdom
E-mail:

<EMAIL_ADDRESS>
###### Abstract:

We study the threshold corrections to the structure functions in deep-
inelastic scattering (DIS) at the fifth logarithmic (N4LL) order of the soft-
gluon exponentiation in massless perturbative QCD. Using recent results for
the splitting functions and the quark form factor, we derive the fourth-order
contribution to the coefficient $f^{\,\rm q}$ of the form factor and from it
the N4LL part of the exponentiation coefficient $B^{\,\rm DIS}$ in the limit
of a large number of colours. An approximation scheme is shown that leads to
sufficiently accurate N4LL results for full QCD. The N4LL corrections are
small and lead to a further stabilization of the perturbative expansion for
the soft-gluon exponent.

##  1 Introduction

The Wilson coefficients (coefficient functions) for the structure functions of
inclusive DIS have been a subject of research since the early days of QCD.
These quantities are not only relevant for determining the parton distribution
functions (PDFs) and the strong coupling constant $\alpha_{\rm s}$ using
structure function data, see, e.g., [1], but also to other processes and less
inclusive observables in DIS, see, e.g., [2]. The main Wilson coefficients for
DIS are presently known up to the third order in $\alpha_{\rm s}$ in massless
perturbative QCD [3]. Their perturbative expansion is well-behaved except
close to the kinematic endpoints $x=0$ and $x=1$ of the Bjorken variable. The
dominant terms $\ln^{\,\ell}(1-x)/(1-x)_{+}$ in the latter (threshold) limit
are resummed by the soft-gluon exponentiation, see, e.g., [4, 5, 6], which is
best formulated in Mellin $N$-space [4]. So far this resummation has been
performed up to the next-to-next-to-next-to-leading logarithmic (N3LL)
accuracy [7].

The threshold resummation coefficients are closely related to the large-$x$
limit of the quark-quark splitting functions for the PDFs and to the quark
form factor [8, 9] which are both fully known to order ${\alpha_{\rm
s}^{\>\\!3}}$ [10, 11, 12, 13]. Recently the computations of these quantities
have been extended to order ${\alpha_{\rm s}^{\>\\!4}}$ in the (L$n_{c}$)
limit of a large number of colours [14, 15]. Together with approximate results
for the ${n_{\\!f}}$-independent [16] and exact expression for the
${n_{\\!f}}$-dependent contributions to the cusp anomalous dimension in full
QCD [17] these results facilitate the effective extension of the threshold
resummation for the DIS Wilson coefficients to the next (N4LL) logarithmic
order. In the following we recall the theoretical framework, present the N4LL
resummation coefficient and briefly address the numerical implications of this
result for the resummation of DIS in QCD.

##  2 Theoretical framework and new fourth-order coefficients

The all-order large-$N$ behaviour of the DIS Wilson coefficients for $F_{1}$,
$F_{2}$ and $F_{3}$ can be written as

$C^{\,N}(Q^{2})\>=\>g_{0}(Q^{2})\cdot\exp\,[G^{\,N}(Q^{2})]\>+\>{\cal
O}(N^{-1}\ln^{\>\\!n}N)\;,$ (2.1)

where the resummation exponent $G^{\,N}$ of the dominant $N^{0}\ln^{\>\\!n}N$
contributions is given by [18]

$\displaystyle
G^{\,N}=\int_{0}^{1}dz\>\frac{z^{\,N-1}-1}{1-z}~{}\bigg{[}\int_{\mu_{f}^{2}}^{(1-z)Q^{2}}\frac{dq^{2}}{q^{2}}\,A^{\rm
q}\big{(}\alpha_{\rm s}(q^{2})\big{)}+B^{\>\\!\rm DIS}\big{(}\alpha_{\rm
s}((1-z)q^{2})\big{)}\bigg{]}\>.$ (2.2)

Here $A^{\rm q}$ is the (light-like) quark cusp anomalous dimension and
$B^{\>\\!\rm DIS}$ is the resummation coefficient for DIS. Both have
perturbative series, $A^{\rm q}=\sum_{i}\,a_{s}^{i}\,A^{\rm q}_{i}$ etc, in
terms of the strong coupling which we normalize as $a_{\rm s}\equiv\alpha_{\rm
s}/4\pi$. Performing the integrations one can organize the exponent as

$\displaystyle G^{\,N}$ $\displaystyle=\ln
N\,g^{(1)}(\lambda)+g^{(2)}(\lambda)+a_{\rm s}g^{(3)}(\lambda)+{a_{\rm
s}^{\>\\!2}}g^{(4)}(\lambda)+{a_{\rm s}^{\>\\!3}}g^{(5)}(\lambda)+\ldots\;,$
(2.3)

where $\lambda=\beta_{0}\>\\!a_{\rm s}\ln\,N$ or $\lambda=\beta_{0}\>\\!a_{\rm
s}\ln\,\widetilde{N}$ with $\ln\,\widetilde{N}=\ln\,N+\gamma_{\,\rm e\,}$. The
first $n\\!+\\!1$ terms in (2.3) are required for the resummation at NnLL
accuracy. The N2LL and N3LL contributions to $G^{\,N}$ have been derived in
[19, 7]; explicit expressions can be found in (3.3) – (3.6) of [7]. The
lengthy new function $g^{(5)}(\lambda)$ entering at N4LL will be presented in
[20]. The $N$-independent prefactor $g_{0}$ is presently known to order
${\alpha_{\rm s}^{\>\\!3}}$ from the all-$N$ calculation in [3], see (4.6) –
(4.8) of [7],

$\displaystyle g_{0}=1+a_{\rm s}g_{01}+{a_{\rm s}^{\>\\!2}}g_{02}+{a_{\rm
s}^{\>\\!3}}g_{03}+{\cal O}({a_{\rm s}^{\>\\!4}})\,.$ (2.4)

The resummation to N4LL requires the terms up to $A_{5}^{\rm q}$ and
$B_{4}^{\>\\!\rm DIS}$ in their corresponding expansions. The impact of the
former quantity, for which a first estimate has been obtained in [21], is very
small. $B^{\>\\!\rm DIS}$ can be calculated from knowledge of the quark form
factor or the DIS Wilson coefficients. The form factor satisfies a
differential equation which follows from the renormalization group and gauge
invariance. Its solution can be found in terms of the cusp anomalous dimension
$A^{\rm q}$ and the function $G^{\rm q}$ containing the quantity $f^{\,\rm q}$
related to a universal eikonal anomalous dimension and the coefficient
$B^{\>\\!\rm q}$ of $\delta(1-x)$ in the quark-quark splitting function. The
four-loop coefficient of $G^{\>\\!\rm q}$ (which appears in the $1/\epsilon$
coefficient in the solution of the form factor) can be written as

$\displaystyle G^{\,\rm q}_{4}=2B_{4}^{\,\rm q}+f_{4}^{\,\rm
q}+\beta_{2}f_{01}^{\,\rm q}+\beta_{1}f_{02}^{\,\rm q}+\beta_{0}f_{03}^{\,\rm
q}+{\cal O}(\epsilon)\>,$ (2.5)

where the quantities $f_{0{\rm n}}^{\,\rm q}$ are (combinations of) known
lower-order coefficients of $G^{\,\rm q}$, see [12] and (20) of [9]. Hence
$f_{4}^{\,\rm q}$ can be determined in the large-$n_{c}$ limit from the
results of [14, 15]. We find

$\displaystyle f_{4}^{\,\rm q}\Bigl{|}_{{\rm L}{n_{c}}}$ $\displaystyle=$
$\displaystyle{C_{F}}\\*{n_{c}^{\,3}}\\*\biggl{(}\frac{9364079}{6561}-\frac{1186735}{729}\\*\zeta_{2}-\frac{837988}{243}\\*\zeta_{3}+\frac{115801}{27}\\*\zeta_{4}+\frac{11896}{9}\\*\zeta_{2}\\*\zeta_{3}+3952\,\\*\zeta_{5}$
(2.6)
$\displaystyle\mbox{}-\frac{4796}{9}\\*\zeta_{3}^{2}-\frac{129547}{54}\\*\zeta_{6}-416\,\\*\zeta_{2}\\*\zeta_{5}-720\,\\*\zeta_{3}\\*\zeta_{4}-1700\,\\*\zeta_{7}\biggr{)}+{C_{F}}\\*{n_{c}^{\,2}}\\*{n_{\\!f}}\\*\biggl{(}-\frac{247315}{432}$
$\displaystyle\mbox{}+\frac{412232}{729}\\*\zeta_{2}+\frac{102205}{243}\\*\zeta_{3}-\frac{7589}{6}\\*\zeta_{4}-\frac{824}{9}\\*\zeta_{2}\\*\zeta_{3}-\frac{740}{9}\\*\zeta_{5}+\frac{2816}{9}\\*\zeta_{3}^{2}+\frac{15611}{27}\\*\zeta_{6}\biggr{)}$
$\displaystyle\mbox{{}}+{C_{F}}\\*{n_{c}}\\*{n^{\,2}_{\\!f}}\\*\left(\frac{329069}{17496}-\frac{22447}{729}\\*\zeta_{2}+\frac{25300}{243}\\*\zeta_{3}+\frac{140}{3}\\*\zeta_{4}-\frac{176}{9}\\*\zeta_{2}\\*\zeta_{3}-\frac{856}{9}\\*\zeta_{5}\right)$
$\displaystyle\mbox{{}}+{C_{F}}\\*{n^{\,3}_{\\!f}}\\*\left(-\frac{16160}{6561}-\frac{16}{81}\\*\zeta_{2}-\frac{400}{243}\\*\zeta_{3}+\frac{128}{27}\\*\zeta_{4}\right)\,.\qquad\quad$

The ${a_{\rm s}^{\>\\!4}}$ contribution to resummation $B^{\,\rm DIS}$ reads,
in terms of the genuine four-loop contributions $f_{4}^{\,\rm q}$ and
$\,B_{4}^{\,\rm q}$, which are exactly known only in the L$n_{c}$ limit for
now, and lower-order coefficients,

$\displaystyle B_{4}^{\>\\!\rm DIS}$ $\displaystyle\\!=\\!$ $\displaystyle-
f_{4}^{\,\rm q}-B_{4}^{\,\rm q}-\beta_{2}\\*\left(f_{01}^{\,\rm
q}+g_{01}-\frac{1}{2}\\*\zeta_{2}\,\\*A_{1}^{\rm
q}\right)+\beta_{0}^{3}\\*\left(3\\*\zeta_{2}\\*f_{01}^{\,\rm
q}+3\\*\zeta_{2}\\*g_{01}+2\\*\zeta_{3}\\*f_{1}^{\,\rm
q}+2\\*\zeta_{3}\\*B_{1}^{\,\rm q}\right.$ (2.7)
$\displaystyle\mbox{{}}\left.+\frac{3}{2}\\*\zeta_{4}\,\\*A_{1}^{\rm
q}-\frac{3}{4}\\*\zeta_{2}^{2}\\*A_{1}^{\rm
q}\right)+\beta_{0}\\*\beta_{1}\\*\left(\frac{5}{2}\\*\zeta_{2}\\*f_{1}^{\,\rm
q}+\frac{5}{2}\\*\zeta_{2}\\*B_{1}^{\,\rm
q}+\frac{5}{3}\\*\zeta_{3}\,\\*A_{1}^{\rm
q}\right)+\beta_{0}^{2}\\*\left(3\\*\zeta_{2}\\*f_{2}^{\,\rm
q}+3\\*\zeta_{2}\\*B_{2}^{\,\rm q}+2\\*\zeta_{3}\,\\*A_{2}^{\rm q}\right)$
$\displaystyle\mbox{{}}-\beta_{1}\\*\left(f_{02}^{\,\rm
q}+2\\*g_{02}-\left(g_{01}\right)^{2}-\zeta_{2}\,\\*A_{2}^{\rm
q}\right)-\beta_{0}\\*\left(f_{03}^{\rm
q}+3\\*g_{03}-3\\*g_{02}\\*g_{01}-\left(g_{01}\right)^{3}-\frac{3}{2}\\*\zeta_{2}\,\\*A_{3}^{\rm
q}\right)\>\>,$

where $g_{0i}$ are to be taken without the $\gamma_{\,\sf e\,}$ terms in (4.6)
– (4.8) of [7]. Its explicit form is given by

$\displaystyle\;B_{4}^{\,\rm DIS}\Bigl{|}_{{\rm L}{n_{c}}}\\!\\!$
$\displaystyle=$
$\displaystyle{C_{F}}\\*{n_{c}^{\,3}}\\*\biggl{(}-\frac{2040092429}{139968}+\frac{23011973}{1944}\\*\zeta_{2}+\frac{517537}{36}\\*\zeta_{3}-\frac{312481}{36}\\*\zeta_{4}-\frac{39838}{9}\\*\zeta_{2}\\*\zeta_{3}$
(2.8)
$\displaystyle\mbox{}-\frac{50680}{9}\\*\zeta_{5}-988\,\\*\zeta_{3}^{2}+\frac{12467}{6}\\*\zeta_{6}+496\,\\*\zeta_{2}\\*\zeta_{5}+688\,\\*\zeta_{3}\\*\zeta_{4}+2260\,\\*\zeta_{7}\biggr{)}$
$\displaystyle\mbox{}+{C_{F}}\\*{n_{c}^{\,2}}\\*{n_{\\!f}}\\*\biggl{(}\frac{83655179}{11664}-\frac{5160215}{972}\\*\zeta_{2}-\frac{639191}{162}\\*\zeta_{3}+\frac{24856}{9}\\*\zeta_{4}+\frac{8624}{9}\\*\zeta_{2}\\*\zeta_{3}$
$\displaystyle\mbox{}+200\,\\*\zeta_{5}-32\,\\*\zeta_{3}^{2}-\frac{1201}{3}\\*\zeta_{6}\biggr{)}+{C_{F}}\\*{n^{\,3}_{\\!f}}\\*\left(\frac{50558}{2187}+\frac{80}{81}\\*\zeta_{3}-\frac{1880}{81}\\*\zeta_{2}+\frac{40}{9}\\*\zeta_{4}\right)$
$\displaystyle\mbox{}+{C_{F}}\\*{n_{c}}\\*{n^{\,2}_{\\!f}}\\*\biggl{(}-\frac{5070943}{5832}+\frac{160903}{243}\\*\zeta_{2}+\frac{14618}{81}\\*\zeta_{3}-\frac{2110}{9}\\*\zeta_{4}-\frac{400}{9}\\*\zeta_{2}\\*\zeta_{3}+\frac{904}{9}\\*\zeta_{5}\biggr{)}.\qquad$

##  3 Numerical implications

The lower-order coefficients $B_{l}^{\,\rm DIS}$ have the same structure as
(2.7), i.e., they contain $-f_{l}^{\,\rm q}-B_{l}^{\,\rm q}$ and lower-order
coefficients. Therefore, by comparing the exact results to an approximation at
NlLL in which the L$n_{c}$ expression for $-f_{l}^{\,\rm q}-B_{l}^{\,\rm q}$
is used together with the exact lower-order coefficients, we can check whether
(2.7) with the L$n_{c}$ results for $f_{4}^{\,\rm q}$ and $B_{4}^{\,\rm q}$
can be expected to provide a good approximation for $B_{4}^{\>\\!\rm DIS}$ and
hence $G^{\,N}$ at the N4LL accuracy of full QCD.

This comparison is carried out in Fig. 1 for $l=2$ and $l=3$ (at $l=1$ there
is no difference between the L$n_{c}$ limit and full QCD). The L$n_{c}$ curves
are off by less than 0.5% at N2LL and 0.25% at N3LL for $G^{\,\rm DIS}$ in the
$N$-range shown and, at $x\leq 0.9$, for the convolution of its exponential
with a schematic but sufficiently realistic form for a quark PDF. Therefore we
can safely expect that the L$n_{c}$ numbers will deviate from (presumably
exceed) the exact QCD results by well below 1%.

Figure 1:  The ratio of the large-${n_{c}}$ approximation, defined as above,
and the exact results at N2LL and N3LL for the DIS resummation exponent
$G^{\,N}$ (left) and for the convolution of the exponential with the schematic
quark PDF shape $xf=x^{0.5}(1-x)^{3}$ (right) for $\alpha_{\rm s}=0.2$ and
${n_{\\!f}}=3$ flavours.

The cumulative effect, relative to the NLL results, of the exact N2LL and N3LL
contributions and our new N4LL corrections, as above determined using the
L$n_{c}$ limit of $-f_{l}^{\,\rm q}-B_{l}^{\,\rm q}$ in (2.7), is illustrated
in Fig. 2. Unlike the N3LL contribution, the N4LL correction is almost
negligible at $N\leq 15$ and $x\leq 0.9$. Even at $N=40$, the functions
$g^{(n)}(\lambda)$ add only 6%, 1.6% and 1% to the NLL result, respectively,
for $n=2$, $n=3$ and $n=4$, where the latter L$n_{c}$ result is presumably a
slight overestimate. The corresponding N2LL, N3LL and N4LL percentages for the
convolution of $\exp G^{\,N}$ with $xf=x^{0.5}(1-x)^{3}$ at $x=0.95$ read
9.5%, 1.5% and 0.5%, where we have performed the Mellin inversion using a
standard contour, see, e.g., [22], which constitutes a ‘minimal prescription’
contour [4] in the context of the present exponentiation. It appears that the
expansion of $G^{\,N}$ to N4LL for the structure functions in inclusive DIS is
sufficient for all practical purposes.

Figure 2:  Left: The DIS resummation exponent $G^{\,N}$ in (2.3) up to N4LL
accuracy, normalized to the NLL result at the standard reference point
$\alpha_{\rm s}=0.2$ for $n_{f}=3$. Right: corresponding $x$-space results for
$\exp G^{\,N}$ after convoluted with a schematic form of a quark PDF of the
proton.

##  4 Summary and outlook

We have studied the soft-gluon exponentiation (SGE) of inclusive DIS at the
fifth logarithmic (N4LL) order. Recent four-loop results on splitting
functions and the quark form factor [14, 15] facilitate the exact
determination of the form-factor coefficient $f^{\,\rm q}$ and the SGE
coefficient $B^{\,\rm DIS}$ at order ${\alpha_{\rm s}^{\>\\!4}}$ in the
large-$n_{c}$ (L$n_{c}$) limit. Both coefficients are relevant beyond the
context of DIS: Like the lightlike quark and gluon cusp anomalous dimensions
$A^{\rm q,g}$ [10], the quantities $f^{\,\rm q,g}$ are maximally non-Abelian
and related by a simple Casimir scaling up to three loops. We expect that the
generalized Casimir scaling of [16] also applies to $f^{\,\rm q,g}$, hence our
result (2.6) fixes also $f^{\,\rm g}$ at large $n_{c}$. The coefficient here
called $B^{\,\rm DIS}$ is due to the outgoing unobserved quark; hence it
contributes to the SGE for many other processes including, e.g., direct photon
production [6].

The L$n_{c}$ approximation to the N4LL resummation exponent $G^{\,N}$ for
inclusive DIS, defined as discussed above, is sufficiently accurate to
demonstrate that the N4LL corrections are small: they contribute well below 1%
over a wide range in $N$ and $x$. As shown in [23], the $1/N\ln^{\>\\!\ell}N$
non-SGE contributions are larger; the highest four of these logarithms are
currently known to all orders [23, 24] – recall the parameter $\xi_{\rm
DIS_{4}}$ unspecified in [23] was fixed in [24]. We have considered the case
of ${n_{\\!f}}=3$ light flavours. In electromagnetic and neutral-current DIS,
also charm production close to threshold needs to be taken into account beyond
the threshold for $c\bar{c}$ production, see [25].

## References

  * [1] A. Accardi et al.,  Eur. Phys. J. C76 (2016) 471, arXiv:1603.08906.
  * [2] P. Bolzoni, F. Maltoni, S. Moch and M. Zaro,  Phys. Rev. Lett. 105 (2010) 011801, arXiv:1003.4451,
M. Cacciari et al.,  Phys. Rev. Lett. 115 (2015) 082002, arXiv:1506.02660,
F.A. Dreyer and A. Karlberg,  Phys. Rev. Lett. 117 (2016) 072001,
arXiv:1606.00840,
J. Currie et al.,  JHEP 05 (2018) 209, 1803.09973.

  * [3] J.A.M. Vermaseren, A. Vogt and S. Moch,  Nucl. Phys. B724 (2005) 3, hep-ph/0504242,
S. Moch, J.A.M. Vermaseren and A. Vogt,  Nucl. Phys. B813 (2009) 220,
arXiv:0812.4168.

  * [4] S. Catani, M.L. Mangano, P. Nason and L. Trentadue,  Nucl. Phys. B478 (1996) 273, hep-ph/9604351.
  * [5] H. Contopanagos, E. Laenen and G.F. Sterman,  Nucl. Phys. B484 (1997) 303, hep-ph/9604313.
  * [6] S. Catani, M.L. Mangano and P. Nason,  JHEP 07 (1998) 024, hep-ph/9806484.
  * [7] S. Moch, J.A.M. Vermaseren and A. Vogt,  Nucl. Phys. B726 (2005) 317, hep-ph/0506288.
  * [8] S. Moch and A. Vogt,  Phys. Lett. B631 (2005) 48, hep-ph/0508265,
E. Laenen and L. Magnea,  Phys. Lett. B632 (2006) 270, hep-ph/0508284.

  * [9] V. Ravindran,  Nucl. Phys. B752 (2006) 173, hep-ph/0603041.
  * [10] S. Moch, J.A.M. Vermaseren and A. Vogt,  Nucl. Phys. B688 (2004) 101, hep-ph/0403192,
A. Vogt, S. Moch and J.A.M. Vermaseren,  Nucl. Phys. B691 (2004) 129, hep-
ph/0404111

  * [11] S. Moch, J.A.M. Vermaseren and A. Vogt,  JHEP 08 (2005) 049, hep-ph/0507039.
  * [12] S. Moch, J.A.M. Vermaseren and A. Vogt,  Phys. Lett. B625 (2005) 245, hep-ph/0508055.
  * [13] P.A. Baikov et al.,  Phys. Rev. Lett. 102 (2009) 212002, arXiv:0902.3519,
R.N. Lee, A.V. Smirnov and V.A. Smirnov,  JHEP 04 (2010) 020, arXiv:1001.2887,
T. Gehrmann et al.,  JHEP 06 (2010) 094, arXiv:1004.3653,
T. Gehrmann et al.,  JHEP 11 (2010) 102, arXiv:1010.4478.

  * [14] J. Henn et al.,  JHEP 03 (2017) 139, arXiv:1612.04389.
  * [15] S. Moch et al.,  JHEP 10 (2017) 041, arXiv:1707.08315.
  * [16] S. Moch et al.,  Phys. Lett. B782 (2018) 627, arXiv:1805.09638.
  * [17] J.M. Henn, A.V. Smirnov, V.A. Smirnov and M. Steinhauser,  JHEP 05 (2016) 066, arXiv:1604.03126, A. Grozin,  PoS LL2016 (2016) 053, arXiv:1605.03886,
J. Davies et al.,  Nucl. Phys. B915 (2017) 335, arXiv:1610.07477,
R.N. Lee, A.V. Smirnov, V.A. Smirnov and M. Steinhauser,  Phys. Rev. D96
(2017) 014008, arXiv:1705.06862, A. Grozin,  JHEP 06 (2018) 073,
arXiv:1805.05050,
R.N. Lee, A.V. Smirnov, V.A. Smirnov and M. Steinhauser,  JHEP 02 (2019) 172,
arXiv:1901.02898,
J.M. Henn, T. Peraro, M. Stahlhofen and P. Wasser,  (2019), arXiv:1901.03693.

  * [18] S. Moch, J.A.M. Vermaseren and A. Vogt,  Nucl. Phys. B646 (2002) 181, hep-ph/0209100,
S. Forte and G. Ridolfi,  Nucl. Phys. B650 (2003) 229, hep-ph/0209154,
E. Gardi and R.G. Roberts,  Nucl. Phys. B653 (2003) 227, hep-ph/0210429.

  * [19] A. Vogt,  Phys. Lett. B497 (2001) 228, hep-ph/0010146,
S. Catani, D. de Florian, M. Grazzini and P. Nason,  JHEP 07 (2003) 028, hep-
ph/0306211.

  * [20] G. Das, S. Moch and A. Vogt,  DESY 19-088, LTH 1205, to appear.
  * [21] F. Herzog et al.,  Phys. Lett. B790 (2019) 436, 1812.11818.
  * [22] A. Vogt,  Comput. Phys. Commun. 170 (2005) 65, hep-ph/0408244.
  * [23] S. Moch and A. Vogt,  JHEP 11 (2009) 099, arXiv:0909.2124.
  * [24] G.Grunberg,  Phys. Lett. B687 (2010) 405, arXiv:0911.4471v5,
A.A. Almasy, G. Soar and A. Vogt,  JHEP 03 (2011) 030, arXiv:1012.3352.

  * [25] H. Kawamura, N.A. Lo Presti, S. Moch and A. Vogt,  Nucl. Phys. B864 (2012) 399, arXiv:1205.5727


# Reply to comment on ”Exact Solution for Three-Dimensional Ising Model” added
by Reply to the Reply by Jacques H.H. Perk

Degang Zhang  College of Physics and Electronic Engineering, Sichuan Normal
University, Chengdu 610101, China  Institute of Solid State Physics, Sichuan
Normal University, Chengdu 610101, China

###### Abstract

In this reply, I point out that the comment by Jacques H.H. Perk on my paper
”Exact Solution for Three-Dimensional Ising Model” [Symmetry 2021, 13, 1837]
is wrong because the author seems not to understand the Onsager operators
along both X and Y directions and the relation between them.

It is known that the three-dimensional (3D) classical Ising model can be
exactly solved by employing the transfer matrix method. When this approach is
applied, the Hamiltonian of 3D classical Ising model becomes the 2D transverse
quantum Ising model on the transfer matrix plane. For the spin configurations
on a plane, there are two kinds of the boundary conditions: (I) the periodic
boundary conditions along both X and Y directions, and (II) the screw boundary
condition along X direction and the periodic boundary condition along Y
direction. Obviously, the Hamiltonian of the 2D quantum Ising model with (II)
is composed of that with (I) and an alternating Ising chain with the periodic
boundary condition along Y direction. I note that such a classical Ising chain
could not affect the thermodynamic quantities of the 3D Ising model in the
bulk. However, it is nothing but the screw boundary condition along X
direction that makes the local Lie algebra produced by the Onsager operators
closed.

In Ref. [1], I have solved exactly the 3D Ising model with (II) by an
algebraic approach, similar to that in solving the 2D Ising model [2]. The
analytical expression of the partition function per atom for the infinite
crystal is presented. The critical temperatures associated with three
magnetization directions are determined. Some thermodynamic quantities, e.g.
the free energy, the internal energy, and the specific heat, and the high-
temperature expansions of the partition function per atom with
$J_{1}=J_{2}=J$, are also calculated. In addition, the critical exponents
$\alpha_{3D}=0,\nu_{3D}=\frac{2}{3}$ and $\mu_{3D}=\frac{4}{3}$ are obtained.
Here I emphasize that for (I), the Onsager operators along X or Y direction do
not construct a Lie algebra, and hence the 3D Ising model with (I) cannot be
solved exactly. Although the boundary conditions play a crucial role in
solving 3D Ising model, the physical properties are irrelevant to them in the
thermodynamic limit.

In the comment [3], Jacques H.H. Perk does not understand the operators ${\cal
L}^{p}_{a,a+1}$, $L^{p}_{a,a+1}$, ${\cal A}_{p,1}$ and $A_{p,1}$ and the
relation between them. Here $p$ ($p=1,2,\cdots,m)$ denotes the $p$th subspace
with $2^{2n}$ dimensions along Y direction in Fig. 1 [1]. ${\cal
A}_{p,1}\equiv A_{p,1}$ comes from the fact that $[{\cal A}_{p,1},A_{p,1}]=0$
and both ${\cal A}_{p,1}$ and $A_{p,1}$ satisfy the same local Lie algebra [1]
and can be easily proved by defining new spin variables. Due to ${\cal
A}_{p,1}=\sum_{a=1}^{n}{\cal L}^{p}_{a,a+1}$ and
$A_{p,1}=\sum_{a=1}^{n}L^{p}_{a,a+1}$, where ${\cal L}^{p}_{a,a+1}$ in
$2^{2n}$ subspace and $L^{p}_{a,a+1}$ in $2^{2mn}$ space are the Onsager
operators along Y and X directions, respectively (see Fig. 2 in Ref. [1]),
then
$A_{p,1}=\sum_{a=1}^{n}L^{p}_{a,a+1}=\sum_{a=1}^{n}L_{p+(a-1)m,p+am}=\sum_{a=1}^{n}\sigma^{z}_{p+(a-1)m}(\prod_{i=p+(a-1)m+1}^{p+am-1}\sigma^{x}_{i})\sigma^{z}_{p+am}=\sum_{a=1}^{n}\sigma^{z}_{p+(a-1)m}{\sigma^{\prime}}^{z}_{p+am}=\sum_{a=1}^{n}{\cal
L}^{p}_{a,a+1}={\cal A}_{p,1}$, where
${\sigma^{\prime}}^{z}_{p+am}=(\prod_{i=p+(a-1)m+1}^{p+am-1}\sigma^{x}_{i})\sigma^{z}_{p+am}$.
Because $[\prod_{i=p+(a-1)m+1}^{p+am-1}\sigma^{x}_{i},{\cal
A}_{p,1}]=[\prod_{i=p+(a-1)m+1}^{p+am-1}\sigma^{x}_{i},A_{p,1}]=0$,
$\prod_{i=p+(a-1)m+1}^{p+am-1}\sigma^{x}_{i}$ are conserved quanities.
Compared to the equation above, Eq. (2) in Ref. [3] is obviously wrong! In
addition, the operators $\sigma^{z}_{p}\sigma^{z}_{p+m}$ and
$\sigma^{z}_{p}\sigma^{x}_{p+1}\cdots\sigma^{x}_{p+m-1}\sigma^{z}_{p+m}$ in
Eq. (2) in Ref. [3] are not used in my paper [1]. Furthermore, because the
eigenvalue of $\sigma^{x}_{p+1}\cdots\sigma^{x}_{p+m-1}$ is 1 or $-1$,
$\sigma^{z}_{p}\sigma^{z}_{p+m}\not=\sigma^{z}_{p}\sigma^{x}_{p+1}\cdots\sigma^{x}_{p+m-1}\sigma^{z}_{p+m}$.
Therefore, such a fatal error in Eq. (2) in Ref. [3] should lead to the
ridiculous results.

It is noted that the exact partition function of 3D Ising model in the
thermodynamic limit oscillates with the size $m$ along X direction [1].
However, all the previous numerical calculations for the physical quantities,
including the critical temperatures, did not take the oscillatory effect into
account [4]. Therefore, the conjectured values extrapolating to the infinite
system in the numerical calculations seem to be inaccurate, and the 3D finite-
size scaling theory must be modified. Here I emphasize that the exact critical
temperatures at the ferromagnetic and paramagnetic phase transitions
determined by Eq. (32) in Ref. [1] coincide completely with the results found
in the domain wall analysis [5] and the asymptotically exact values in the
anisotropic limit [6,7].

The specific heat $C$ of 2D Ising model with the quadratic symmetry (i.e.
$H_{1}=H_{2D}$) have been calculated analytically by Onsager and can be
expressed in terms of the complete elliptic integrals [2]. The critical
exponent associated with the specific heat $\alpha_{2D}=0$ due to the log
divergence at the critical point [8]. Because 3D Ising model with the simple
cubic symmetry (i.e. $H_{1}=H_{2}=H$) can be mapped exactly into 2D one by Eq.
(41) in Ref. [1], the expression of $C$ in three dimensions has similar form
with that in two dimensions. So the critical exponent $\alpha_{3D}$ of the 3D
Ising model is identical to $\alpha_{2D}$, i.e. $\alpha_{3D}=0$. According to
the scaling laws $d\nu=2-\alpha$ and $\mu+\nu=2-\alpha$ [8], we have
$\nu_{3D}=\frac{2}{3}$ and $\mu_{3D}=\frac{4}{3}$.

In conclusion, there is no question that Jacques H.H. Perk′s comment [3] is
definitely wrong due to misunderstanding my work.

This work was supported by the Sichuan Normal University and the ”Thousand
Talents Program” of Sichuan Province, China.

Reply to ”Reply to Symmetry 2023, 15, 375”

Note that
$\sum_{a=1}^{n}\sigma^{z}_{p+(a-1)m}(\prod_{i=p+(a-1)m+1}^{p+am-1}\sigma^{x}_{i})\sigma^{z}_{p+am}=\sum_{a=1}^{n}\sigma^{z}_{p+(a-1)m}{\sigma^{\prime}}^{z}_{p+am}\equiv\sum_{a=1}^{n}\sigma^{z}_{p+(a-1)m}\sigma^{z}_{p+am}$,
which cannot lead to
$\sigma^{z}_{p+(a-1)m}{\sigma^{\prime}}^{z}_{p+am}=\sigma^{z}_{p+(a-1)m}\sigma^{z}_{p+am}$.
Obviously, such a person misunderstands the summation over $a$, which implies
the sum of spin configurations $(\pm 1)$ on $a$ [1,2,9]. Jacques H.H. Perk
needs to carefully read the Onsager′s paper [2] and my work [1,9].

## References

  * (1) Degang Zhang, Exact Solution for Three-Dimensional Ising Model. Symmetry 2021, 13, 1837.
  * (2) Onsager, L. Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition. Phys. Rev. 1944, 65, 117-149.
  * (3) Jacques H.H. Perk, arXiv:2202.03136.
  * (4) Ferrenberg, A. M.; Xu, J.; Landau, D. P. Pushing the limits of Monte Carlo simulations for the three-dimensional Ising model. Phys. Rev. E 2018, 97, 043301.
  * (5) Zandvliet, H.J.W.; Hoede, C. Boundary tension of 2D and 3D Ising models. Memorandum 2008, 1880. Available online: https://ris.utwente.nl/ws/portalfiles/portal
/5100205/memo1880.pdf (accessed on 1 September 2021).

  * (6) Weng, C.-Y.; Griffiths, R.B.; Fisher, M.E. Critical Temperatures of Anisotropic Ising Lattices. I. Lower Bounds. Phys. Rev. 1967, 162, 475-479.
  * (7) Fisher, M.E. Critical Temperatures of Anisotropic Ising Lattices. II. General Upper Bounds. Phys. Rev. 1967, 162, 480-485.
  * (8) Baxter, R.J. Exactly Solved Models in Statistical Mechanics; Academic Press: London, UK, 1982.
  * (9) Degang Zhang, Reply to Perk, J.H.H. Comment on ”Zhang, D. Exact Solution for Three-Dimensional Ising Model. Symmetry 2021, 13, 1837”. Symmetry 2023, 15, 375.


# Birth Quota of Non-Generic Degeneracy Points

Gergő Pintér  Department of Theoretical Physics, Institute of Physics,
Budapest University of Technology and Economics, Műegyetem rkp. 3., H-1111
Budapest, Hungary  György Frank  Department of Theoretical Physics, Institute
of Physics, Budapest University of Technology and Economics, Műegyetem rkp.
3., H-1111 Budapest, Hungary  Dániel Varjas  Department of Physics, Stockholm
University, AlbaNova University Center, 106 91 Stockholm, Sweden  Max Planck
Institute for the Physics of Complex Systems, Nöthnitzer Strasse 38, 01187
Dresden, Germany  András Pályi  Department of Theoretical Physics, Institute
of Physics, Budapest University of Technology and Economics, Műegyetem rkp.
3., H-1111 Budapest, Hungary  MTA-BME Quantum Dynamics and Correlations
Research Group, Műegyetem rkp. 3., H-1111 Budapest, Hungary

###### Abstract

Weyl points are generic and stable features in the energy spectrum of
Hamiltonians that depend on a three-dimensional parameter space. Non-generic
isolated two-fold degeneracy points, such as multi-Weyl points, split into
Weyl points upon a generic perturbation that removes the fine-tuning or
protecting symmetry. The number of the resulting Weyl points is at least
$|Q|$, where $Q$ is the topological charge associated to the non-generic
degeneracy point. Here, we show that such a non-generic degeneracy point also
has a birth quota, i.e., a maximum number of Weyl points that can be born from
it upon any perturbation. The birth quota is a local multiplicity associated
to the non-generic degeneracy point, an invariant of map germs known from
singularity theory. This holds not only for the case of a three-dimensional
parameter space with a Hermitian Hamiltonian, but also for the case of a two-
dimensional parameter space with a chiral-symmetric Hamiltonian. We illustrate
the power of this result for band structures of two- and three-dimensional
crystals. Our work establishes a strong and powerful connection between
singularity theory and topological band structures, and more broadly,
parameter-dependent quantum systems.

_Introduction._ Weyl semimetals are a class of topological materials whose
electronic band structure exhibits pointlike linear band-touchings [1, 2, 3].
These Weyl points are stable and generic, requiring no fine-tuning or
symmetries. Crystal symmetries, however, can stabilize isolated non-generic
two-fold degeneracy points, such as multi-Weyl points, in the electronic band
structure of three-dimensional solids [4, 5, 6, 7]. If the symmetry is broken,
e.g., by changing an external magnetic (Zeeman) field or applying mechanical
strain, the non-generic degeneracy point splits into multiple Weyl points [4].
As the symmetry-breaking perturbation is switched on gradually, these newly
born Weyl points follow continuous trajectories in the Brillouin zone,
originating from the original degeneracy point.

For a given non-generic degeneracy point, how many newborn Weyl points are
allowed? From topological charge conservation and the generic character of
Weyl points, it follows that the minimum number of newborn Weyl points upon a
generic perturbation is the absolute value $|Q|$ of the topological charge $Q$
associated to the non-generic degeneracy point (e.g. the Chern number in 3D).
The same consideration implies that the number of newborn Weyl points may also
be $|Q|+2M$ with $M$ being a positive integer, in such a way that $M$ of the
excess Weyl points have unit positive charge and $M$ have unit negative
charge, hence the sum of the charges of the newborn Weyl points equals $Q$. Is
there also a ‘birth quota’, i.e., an upper limit of the number of newborn Weyl
points? This is a fundamental question that, to our knowledge, has not been
addressed in the context of band structure theory.

In this work, we answer this question positively, and show that such a non-
generic degeneracy point does have a birth quota, which turns out to be the
so-called _local multiplicity_ known in singularity theory. Our result is not
exclusive to Weyl points in three dimensions (3D), we also obtain an analogous
result for two-dimensional (2D) crystals with chiral symmetry illustrated on
the example of bilayer graphene, as well as a minimal example in 1D (Fig. 1).
We showcase the power of these results on quasiparticle (electronic, photonic,
phononic) band structures: we compute the birth quota of all four types of
two-fold degeneracies stabilised by the 230 crystalline space groups [7]
(Table 1). However, the notion of the birth quota is more generic: it is
applicable to quantum systems controlled by external parameters, such as
interacting spin systems [8, 9, 10, 11, 12] or quantum circuits [13, 14];more
generally, it is applicable to any physical system that is described by a
matrix, e.g., linearly coupled mechanical oscillators or electromagnetic
modes.

_Example: Bilayer graphene._ We exemplify the birth of Weyl points from a non-
generic degeneracy point with the well-known example of electrons in A-B-
stacked bilayer graphene [15]. For more details and illustrations, see Sec. I
of SI.

Bilayer graphene is a two-dimensional crystal, whose simplest band-structure
models have chiral (a.k.a. sublattice) symmetry, which protects the generic
degeneracies (2D Weyl points) at zero energy in the two-dimensional Brillouin
zone. A simple tight-binding model of delocalized electrons in graphene yields
an envelope-function Hamiltonian
$H(k_{x},k_{y})\propto\left(k_{x}^{2}-k_{y}^{2}\right)\sigma_{x}+2k_{x}k_{y}\sigma_{y}$,
valid in the vicinity of the high-symmetry point $K$ (and $K^{\prime}$) of the
Brillouin zone [15, 16]. Here, $k_{x,y}$ are momentum components measured from
$K$, $\sigma_{x,y}$ are Pauli matrices, the matrix structure represents a
combined layer-sublattice degree of freedom, and chiral symmetry forbids the
$\sigma_{z}$ term [17].

This Hamiltonian $H$ exemplifies a non-generic degeneracy point in 2D, located
at $(k_{x},k_{y})=0$ with a quadratic dispersion $E_{\pm}\propto\pm k^{2}$.
Furthermore, the Hamiltonian $H$ can be reinterpreted as the dimensionless
_effective Hamiltonian map_
$h:\mathbb{R}^{2}\to\mathbb{R}^{2},(k_{x},k_{y})\mapsto(k_{x}^{2}-k_{y}^{2},\
2k_{x}k_{y})$. This map has a ‘winding number’ or ‘topological charge’
$Q_{\text{2D}}=+2$, as opposed to 2D Weyl points whose topological charge is
$\pm 1$.

One type of perturbation to this degeneracy point is mechanical strain:
applying it along $x$ adds a perturbation term $H_{\text{s}}\propto\sigma_{x}$
to $H$ [15]. The effect of the strain is that the degeneracy point splits into
two 2D Weyl points along the $k_{x}$ axis, both having charge
$Q_{\text{2D}}=+1$, i.e., the total charge is conserved upon perturbation.

Another type of perturbation appears when the skew interlayer hopping is taken
into account, causing _trigonal warping_ [15, 16]. This perturbation term is
described by a Hamiltonian
$H_{\text{tw}}\propto\left(k_{x}\sigma_{x}-k_{y}\sigma_{y}\right)$. The effect
of this perturbation is that the non-generic degeneracy point of $H$ splits
into four 2D Weyl points: 3 having charge $Q_{\text{2D}}=+1$ and 1 having
$Q_{\text{2D}}=-1$. Again, the charge of the degeneracy point of $H$ is
conserved by the perturbation.

The observation that the non-generic degeneracy point at $K$ can split to two
or four 2D Weyl points triggers the question: is there a perturbation of $H$
such that the number of newborn 2D Weyl points is different from 2 or 4? For
example, is it possible to have 4 newborn Weyl points with $Q_{\text{2D}}=+1$
and 2 newborn Weyl points with $Q_{\text{2D}}=-1$? One implication of our
analysis below is that the answer is ‘no’: The lower bound (2) is governed by
topological charge conservation, whereas the upper bound (4), i.e., the ‘birth
quota’, follows from our analysis below.

_Minimal 1D model._ First, we consider a minimal mathematical model of the
birth quota effect: the birth of generic roots (‘1D Weyl points’) of the
polynomial $f:\mathbb{R}\to\mathbb{R},\,f(x)=x^{3}$, from its non-generic root
at the origin $x=0$. Any generic perturbation of the form
$f_{t}(x)=x^{3}+t_{1}x^{2}+t_{2}x+t_{3}$ has one or three real roots, and all
of them converge to $0$ as the vector of control parameters
$t=(t_{1},t_{2},t_{3})$ tends to 0. We draw such perturbations in Fig. 1a,
where, for simplicity, we set $(t_{1},t_{2},t_{3})=(0,p,q)$. In the case of
three roots (green region), two of them have positive signs (i.e. the function
is increasing, charge $+1$, red dot), and one has negative sign (the function
is decreasing, charge $-1$, blue dot). In the case of one root (orange
region), it has positive sign (charge $+1$). Hence the sum of the charges is
$+1$, same as the charge of the root of the unperturbed function $f$,
exemplifying charge conservation.

Figure 1:  Number of generic roots (‘1D Weyl points’) born from the non-
generic root of $f(x)=x^{3}$ is lower-bounded by the charge of the root (which
is 1) and upper-bounded by the local multiplicity of the root
($\text{mult}_{0}f=3$). a) Phase diagram of perturbations of $f(x)=x^{3}$ of
the form $f_{p,q}(x)=x^{3}+px+q$. Orange: perturbation has one generic root.
Green: perturbation has three generic roots. Purple line: perturbation has one
generic root and one non-generic root. b) Complexified perturbation
$f_{\mathbb{C},p,q}$ with $p=1$, $q=0$. c) Complexified map
$f_{\mathbb{C}}(x)=x^{3}$, exhibiting a single root with a winding number or
local multiplicity $\text{mult}_{0}f=3$. d) Complexified perturbation
$f_{\mathbb{C},p,q}$ with $p=-1$, $q=0$. The complex values of the function
are represented as arrows, so phase windings can be read off. b), d) exhibits
three complex roots with winding number $1$ each, even though the real
function $f_{p,q}$ has b) one real root and d) three real roots.

The number of the real roots changes when the control vector $t$ steps through
the zero locus $\mathcal{D}=\\{(t_{1},t_{2},t_{3})\in\mathbb{R}^{3}\ |\
D(t_{1},t_{2},t_{3})=0\\}$ of the discriminant $D$ of $f$ (see Sec. II. of
SI.). This _discriminant set_ $\mathcal{D}$ in the control space is
illustrated as the purple solid line in Fig. 1a.

Importantly, when we consider the complexification 111For a real analytic
function $f:\mathbb{R}^{m}\to\mathbb{R}^{n}$, its complexification is defined
as the function $f_{\mathbb{C}}:\mathbb{C}^{m}\to\mathbb{C}^{n}$ that has the
same power series as $f$. $f_{\mathbb{C},t}:\mathbb{C}\to\mathbb{C}$, the
_complex discriminant set_ $\mathcal{D}_{\mathbb{C}}$ does not separate
different root configurations: one can always find a path between different
configurations such that the roots avoid collision in the complex plane. In
fact, $f_{\mathbb{C},t}$ has three distinct complex roots for any
$t\in\mathbb{C}^{3}\setminus\mathcal{D}_{\mathbb{C}}$. This is illustrated by
Fig. 1b,c,d, displaying the (c) unperturbed map $f_{\mathbb{C}}$ and (b,d) two
perturbed maps, and the roots of those $\mathbb{C}\to\mathbb{C}$ maps as red
points. The number of complex roots born upon perturbation from the root of
$f_{\mathbb{C}}(x)=x^{3}$ is called the _local multiplicity_
$\text{mult}_{0}f_{\mathbb{C}}=3$ of $f_{\mathbb{C}}$ at the root $x=0$. As
illustrated in Fig. 1c, this local multiplicity is the winding number of
$f_{\mathbb{C}}$ evaluated on a small circle (green) around the root; the
value 3 can be read off.

This winding number is also called the _local degree_
$\text{deg}_{0}f_{\mathbb{C}}$, which is defined as the local degree
$\text{deg}_{0}f_{\mathbb{C},\mathbb{R}}$ of the real-imaginary decomposition
(‘realification’)

$\displaystyle f_{\mathbb{C},\mathbb{R}}$ $\displaystyle:$
$\displaystyle\mathbb{R}^{2}\to\mathbb{R}^{2},$ $\displaystyle(u,v)$
$\displaystyle\mapsto$
$\displaystyle(\text{Re}(f_{\mathbb{C}}(u+iv)),\text{Im}(f_{\mathbb{C}}(u+iv))$
(1)

of the complexification $f_{\mathbb{C}}$. Here, the explicit form of
$f_{\mathbb{C},\mathbb{R}}$ reads
$f_{\mathbb{C},\mathbb{R}}(u,v)=(u^{3}-3uv^{2},3u^{2}v-v^{3})$. Since every
generic root of $f_{t}$ is also a generic root of $f_{\mathbb{C},t}$, we
conclude that the number of roots born from the original root of $f$ at $x=0$
upon perturbation is upper-bounded by
$\text{deg}_{0}f_{\mathbb{C},\mathbb{R}}$.

Note that higher-degree perturbations of $f$ can have more than three roots,
but these extra roots are ‘born at infinity’, not from the original root at
$x=0$. For example, the root $x=1/t$ of the perturbed map
$f_{t}(x)=x^{3}-tx^{4}=-tx^{3}(x-1/t)$ goes to infinity as $t$ tends to 0.

In the rest of this paper, we outline the generalization of this minimal
model, and use it to derive the birth quota of non-generic isolated two-fold
degeneracy points.

_Effective Hamiltonian at a twofold degeneracy point._ The models we study,
e.g., tight-binding models of the electronic band structure, are described by
a Hamiltonian map $H:M^{m}\to\mbox{Herm}(N)$. Here, $M^{m}$ is an
$m$-manifold, and $\mbox{Herm}(N)$ is the real vector space of $N\times N$
Hermitian matrices. We focus on the example when $M^{m}=\text{BZ}$ is the
Brillouin zone of a 3D crystal ($m=3$), but will also comment on the 2D case
($m=2$), relevant for bilayer graphene described above.

A key quantity in our analysis of the birth quota, analogous to the function
$f$ in the minimal model above, is the _effective Hamiltonian map_ $h$
associated to the non-generic isolated two-fold degeneracy point
$P\in\text{BZ}$. Denote the ordinal number of the degenerate levels by $i$ and
$i+1$. Then, for the eigenvalues of $H(P)$, it holds that
$E_{i-1}<E_{i}=E_{i+1}<E_{i+2}$, and $E_{i}<E_{i+1}$ holds in a neighborhood
$U_{P}$ of $P\in\text{BZ}$. The effective Hamiltonian map $h$ is obtained by
an exact Schrieffer–Wolff transformation [19] at $P$, which provides a map
from the neighborhood $U_{P}$ of the degeneracy point $P$ into the space of
traceless Hermitian $2\times 2$ matrices. The latter matrix space is
identified with $\mathbb{R}^{3}$ via the standard Pauli-matrix decomposition,
$X\sigma_{x}+Y\sigma_{y}+Z\sigma_{z}\equiv(X,Y,Z)$, hence the Schrieffer–Wolff
transformation yields the effective Hamiltonian map
$h:U_{P}\to\mathbb{R}^{3}$.

The fact that $P$ is an isolated degeneracy point can be reformulated as
$h^{-1}(0)=\\{P\\}$, i.e. the pre-image set of the origin $0\in\mathbb{R}^{3}$
contains only one point. We assume that in the natural wave-vector coordinates
of the BZ, measured from $P$ as the reference point, $h$ is an analytic map,
that is, its Taylor series is convergent and it produces $h$ 222Note that our
analysis also works in $\mathcal{C}^{\infty}$ (i.e. smooth) category, using
the notion of _map germs_ , however the analytic condition simplifies the
discussion.. This results in an $h:(\mathbb{R}^{3},0)\to(\mathbb{R}^{3},0)$
map, consisting of 3 locally convergent power series $h_{1}$, $h_{2}$ and
$h_{3}$ of 3 variables, fulfilling $h(0)=0$. Our following analysis works for
$m=3$ in general, and also for $m=2$ assuming chiral symmetry, when the
degeneracy is at zero energy. The latter conditions imply $h_{3}=0$ and
provide a map $h:(\mathbb{R}^{2},0)\to(\mathbb{R}^{2},0)$ 333Similarly,
spinless $P\mathcal{T}$ symmetry restricts the Hamiltonian to be real valued,
yielding $h_{2}=0$. Hamiltonians with both chiral and $P\mathcal{T}$ symmetry
correspond to the case $m=1$ with effective Hamiltonian map
$h:(\mathbb{R},0)\to(\mathbb{R},0)$. Such symmetric 1D systems can provide a
physical realization of the minimal example with $H(k)=f(k)\sigma_{x}$..

_Birth quota._ Up to now, we converted band-structure features to an
$h:\mathbb{R}^{3}\to\mathbb{R}^{3}$ analytic map. This enables us to establish
the birth quota of a non-generic degeneracy point using concepts and relations
from singularity theory. The _local degree_ $\deg_{0}h$ of the effective
Hamiltonian map $h:(\mathbb{R}^{m},0)\to(\mathbb{R}^{m},0)$ at $0$ is the
_global degree_ of the normalised map (‘pseudospin texture’ [5])
$\tilde{h}=\frac{h}{|h|}:S^{m-1}_{\epsilon}\to S^{m-1}$ defined on a
sufficiently small sphere $S^{m-1}_{\epsilon}$ around the origin, see [22, 23,
24, 25, 26] for details. For $m=2$ the local degree is a winding number, and
for $m=3$ it agrees with the first Chern number of the eigenstate
corresponding to the $i$-th eigenvalue, see [17]. For both cases, the local
degree is often referred to as the (topological) charge of the degeneracy
point. We will also refer to the complexification
$h_{\mathbb{C}}:(\mathbb{C}^{m},0)\to(\mathbb{C}^{m},0)$ of $h$, and its local
degree $\deg_{0}h_{\mathbb{C}}$.

The effect of a physical perturbation on the crystal electrons, e.g.,
mechanical strain or a change of the magnetic (Zeeman) field, is described in
our framework as a _deformation_ of the effective Hamiltonian map $h$. An
unfolding of $h$ with $k$ _control parameters_ is an analytic map
$\mathcal{H}:\mathbb{R}^{3}\times\mathbb{R}^{k}\to\mathbb{R}^{3}\times\mathbb{R}^{k}$
of the special form $\mathcal{H}(x,y,z,t)=(h_{t}(x,y,z),t)$, where
$t\in\mathbb{R}^{k}$, such that $h_{0}=h$. For fixed control parameters $t$,
we call $h_{t}$ an analytic deformation, or simply deformation, of $h$. To
derive the birth quota of $h$, we will use the complex generalizations of
these concepts: a holomorphic unfolding $\mathcal{H}_{\mathbb{C}}$ and the
corresponding complex deformation $h_{\mathbb{C},t}$ of the complexification
$h_{\mathbb{C},0}=h_{\mathbb{C}}$ are defined similarly to the real case above
[26, 23].

We call the points of $h_{t}^{-1}(0)\subset\mathbb{R}^{3}$ _degeneracy points_
of $h_{t}$, and the points of $h_{\mathbb{C},t}^{-1}(0)\subset\mathbb{C}^{3}$
_complex degeneracy points_. A complex degeneracy point $p$ is generic, and we
call it a _complex Weyl point_ , if the Jacobian of $h_{\mathbb{C},t}$ at $p$
has maximal rank, i.e., rank 3. A real degeneracy point $p$ is generic, and we
call it a (real) Weyl point, if the Jacobian of $h_{t}$ at $p$ has maximal
rank. This happens if and only if $p$ is generic as a complex degeneracy
point. At a real Weyl point $p$, the local degree is $\deg_{p}h_{t}=\pm 1$,
determined by the sign of the Jacobian determinant, cf. [4]. At a complex Weyl
point $p$, the local degree of $h_{\mathbb{C},t}$ is always
$\deg_{p}h_{\mathbb{C},t}=1$.

Our goal is to characterize the birth of (real) Weyl points from a non-generic
degeneracy point. Hence, we need to distinguish between Weyl points born from
the original degeneracy point, and all other Weyl points. This we already
noted in the minimal model above, by pointing out that 1D Weyl points can be
born at infinity. To make this distinction here, we consider the complex case
first. We take a spherical boundary $S^{5}_{\epsilon}\subset\mathbb{C}^{3}$ in
the configuration space, centered at the origin, which we call the _separator_
, such that the only degeneracy point of $h_{\mathbb{C}}$ inside this sphere
is the origin. The separator bounds the closed ball
$B^{6}_{\epsilon}\subset\mathbb{C}^{3}$. When a deformation is applied
continuously, the degeneracy points $h^{-1}_{\mathbb{C},t}(0)$ follow
continuous trajectories in the complex configuration space. Therefore, there
is a neighborhood $\mathcal{U}_{\mathbb{C}}$ of the origin of the control
space such that for all $t\in\mathcal{U}_{\mathbb{C}}$, (i) the degeneracy
points born from the original degeneracy point do not reach the separator, and
(ii) the degeneracy points of the undeformed map $h_{\mathbb{C}}$ that are
outside of the separator do stay outside. For convenience, we take a
neighborhood $\mathcal{U}_{\mathbb{C}}$ that is an open ball centered at the
origin of the control space. Furthermore, from now on, we restrict the
unfolding $\mathcal{H}_{\mathbb{C}}$ onto
$B^{6}_{\epsilon}\times\mathcal{U}_{\mathbb{C}}$.

We define the _complex discriminant set_
$\mathcal{D}_{\mathbb{C}}\subset\mathcal{U}_{\mathbb{C}}$ as those control
vectors $t$ for which $h_{\mathbb{C},t}$ has non-generic complex degeneracy
points in $B^{6}_{\epsilon}$. That is, for any control vector
$t\in\mathcal{U}_{\mathbb{C}}\setminus\mathcal{D}_{\mathbb{C}}$, all
degeneracy points of $h_{\mathbb{C},t}$ are complex Weyl points. A key
observation is that the number of complex Weyl points within
$B^{6}_{\epsilon}$, denoted as $\sharp h_{\mathbb{C},t}^{-1}(0)$, is the same
for any control vector
$t\in\mathcal{U}_{\mathbb{C}}\setminus\mathcal{D}_{\mathbb{C}}$. This is
related to the fact that the complex codimension of the discriminant set
$\mathcal{D}_{\mathbb{C}}$ is at least 1, and therefore the real codimension
is at least 2. This implies that
$\mathcal{U}_{\mathbb{C}}\setminus\mathcal{D}_{\mathbb{C}}$ is path-connected,
and hence the number $\sharp h_{\mathbb{C},t}^{-1}(0)$ of complex Weyl points
cannot change along any control trajectory in
$\mathcal{U}_{\mathbb{C}}\setminus\mathcal{D}_{\mathbb{C}}$. For details, we
refer to Sec. II.A of SI and [26]. In conclusion, the number of preimages of
$0$, $\sharp h_{\mathbb{C},t}^{-1}(0)$ is a property of $h_{\mathbb{C}}$, and
hence a property of $h$. It is called the _local multiplicity_ of
$h_{\mathbb{C}}$ ($\text{mult}_{0}h_{\mathbb{C}}$), and also called the _local
multiplicity_ of $h$ ($\text{mult}_{0}h$).

In the band-structure context, the effective Hamiltonian map $h$ is a real
map. We denote the set of real control vectors $t\in\mathbb{R}^{k}$ in
$\mathcal{U}_{\mathbb{C}}$ and $\mathcal{D}_{\mathbb{C}}$ by $\mathcal{U}$ and
$\mathcal{D}$, respectively. With these, we can express the key message of
this work: the local multiplicity $\text{mult}_{0}h$ is the birth quota of the
original non-generic degeneracy point, i.e, $\text{mult}_{0}h$ is the upper
bound of the number of Weyl points born from the original degeneracy point:

$\sharp h^{-1}_{t}(0)\leq\text{mult}_{0}h.$ (2)

Here, $h_{t}$ is restricted to $B^{3}_{\epsilon}$ and
$t\in\mathcal{U}\setminus\mathcal{D}$. Equation (2) is a consequence of the
previous paragraph, and the fact that a real Weyl point is also a complex Weyl
point.

_Methods to calculate the birth quota._ We mention three different methods to
compute the birth quota $\mbox{mult}_{0}h_{\mathbb{C}}$. The first one is the
direct application of the definition for the _constant deformation_
$h_{\mathbb{C},q}(x)=h_{\mathbb{C}}-q$ with $q\in\mathbb{C}^{3}$ close to 0.
Note that by Sard’s lemma [22], for almost all values $q$, all degeneracy
points of $h_{\mathbb{C},q}$ are generic. By definition, the number of the
roots of $h_{\mathbb{C},q}$ inside the ball $B^{6}_{\epsilon}$ is the local
multiplicity. Hence the computation of the birth quota
$\text{mult}_{0}h_{\mathbb{C}}$ reduces to solving the equation
$h_{\mathbb{C}}=q$, which, in general, can be done numerically.

Another method is based on the fact that the local multiplicity is equal to
the local degree $\deg_{0}h_{\mathbb{C}}$ of the complexification
$h_{\mathbb{C}}$, see [26, E.3]. Here, the local degree is understood in the
sense described for the minimal model, i.e. it is the local degree of the
real-imaginary decomposition (realification)
$h_{\mathbb{C},\mathbb{R}}:\mathbb{R}^{6}\to\mathbb{R}^{6}$. The computation
of the local degree also allows the use of integral formulas, see e.g. in [27,
17].

We note that an alternative, commonly used definition of the local
multiplicity uses algebraic methods: the local multiplicity is defined as the
dimension of the so-called _local algebra_ of $h$ at $0$, see Sec. III. of SI
and Refs. [23, 26]. This definition also provides a practical computational
method of the birth quota. In the physical applications below, we computed the
birth quotas using deformations, as well as this algebraic method.

_Applications._ In what follows, we derive the local multiplicity for two
different families of isolated twofold degeneracy points: for chiral symmetric
band-structure models of few-layer graphene, and for all stable isolated
twofold degeneracy points that arise in time-reversal-symmetric crystals. The
results are summarised in Table 1.

For bilayer graphene, discussed above, we find that the local multiplicity of
the unperturbed effective Hamiltonian map $h$ is 4, confirming that the number
of newborn 2D Weyl points upon a generic deformation is either 2 or 4. For
chiral-symmetric models of trilayer graphene with ABC stacking, or multi-layer
($n$-layer) graphene with ABCA… stacking, the local degree of the degeneracy
point at the $K$ point is $n$, whereas the local multiplicity is $n^{2}$, as
shown in the top panel of Table 1. (For derivation, see Example 2 in Sec. IV.
of SI.)

Our methods provide the birth quota for non-generic degeneracy points
appearing in 3D crystals as well. We focus on time-reversal symmetric
crystals, which are classified in 230 space groups. It is known [7] that there
are four types of isolated twofold degeneracy points in the quasiparticle band
structures of such crystals, listed in the bottom panel of Table 1.

The charge-1 Weyl point is the generic degeneracy point (called Weyl point
throughout this paper), whereas the charge-2, -3, -4 Weyl points are non-
generic degeneracy points showing nonlinear dispersion in certain directions.
Charge-2 and -3 Weyl points have been proposed in Ref. [4], with effective
Hamiltonians in the form
$H=ak_{z}\sigma_{z}+\left(bk_{+}^{n}+ck_{-}^{n}\right)\sigma_{+}+\textnormal{h.c.}$
with $n\in\\{2,3\\}$, $|b|\neq|c|$, $k_{\pm}=k_{x}\pm ik_{y}$ and
$\sigma_{\pm}=\sigma_{x}\pm i\sigma_{y}$. The charge-4 Weyl point has been
proposed in Ref. [5], and is described by the effective Hamiltonian
$H=Ak_{x}k_{y}k_{z}\sigma_{z}+B\left(k_{x}^{2}+\omega
k_{y}^{2}+\omega^{2}k_{z}^{2}\right)\sigma_{+}+\textnormal{h.c.}$ with
$\omega=\exp(-2\pi i/3)$. Note that this map corresponds to the map studied in
Ref. [25, Pg. 24], and also that such charge-4 Weyl points have been observed
experimentally [28, 29].

We have computed the local multiplicities of these degeneracy points, and list
the results in Table 1. (Derivations are shown in Example 2 and Example 3 in
Sec. IV. of SI.) The table also indicates the number of space groups where the
corresponding non-generic degeneracy points are stabilised by symmetries.

$\mathbb{R}^{2}\to\mathbb{R}^{2}$ |  local degree (winding number, top. charge) |  local multiplicity (birth quota)
---|---|---
monolayer graphene | 1 | 1
bilayer graphene (AB) | 2 | 4
$n$-layer graphene (ABCA…) | $n$ | $n^{2}$

$\mathbb{R}^{3}\to\mathbb{R}^{3}$ |  local degree (Chern number, top. charge) |  local multiplicity (birth quota) |  #SG SO |  #SG nSO
---|---|---|---|---
Charge-1 Weyl point | 1 | 1 |  |
Charge-2 Weyl point | 2 | 4 | 36 | 48
Charge-3 Weyl point | 3 | 9 | 26 | 12
Charge-4 Weyl point | $4$ | $12$ | - | 13

Table 1: Local multiplicities of isolated twofold degeneracy points. The upper
panel ($\mathbb{R}^{2}\to\mathbb{R}^{2}$) lists the absolute value of the
local degree, and the local multiplicity, associated to the electronic
quasiparticles that emerge in the simplest chiral-symmetric tight-binding
models of mono- or multilayer graphene. The lower panel
($\mathbb{R}^{3}\to\mathbb{R}^{3}$) lists the same invariants, for all four
types of isolated two-fold degeneracy points in crystals. The absolute value
of the local degree is the minimum number of newborn Weyl points upon a
generic deformation. The local multiplicity is the birth quota, i.e., the
maximum number of newborn Weyl points. The last two columns of the lower panel
indicate the number of space groups where symmetries stabilise non-generic
degeneracy points, based on Ref. [7]. The label SO (nSO) denotes band
structures with (without) spin-orbit coupling, e.g. electronic (phononic,
photonic) band structures.

_Conclusions._ We have shown that any isolated two-fold non-generic degeneracy
point in a 3D configuration space has a birth quota, a maximum number of Weyl
points that can be born from the degeneracy point. This birth quota can be
computed as the local multiplicity of the effective Hamiltonian map associated
to the degeneracy point. We also extended our result to the 2D case with
chiral symmetry. We have computed the birth quota for the isolated twofold
degeneracy points of quasiparticles in 3D crystals, and for the $K$-point
electronic degeneracies of multilayer graphene. An important open question
triggered is whether 3-fold or higher-fold isolated degeneracy points exhibit
a similar birth quota effect. Our result establishes a strong connection
between singularity theory and topological semimetals, and also applies to
phononic, magnonic, and photonic band structures. Furthermore, the concepts
and methods discussed here are applicable to parameter-dependent quantum
systems, e.g. multi-terminal Josephson circuits or magnetically controlled
quantum spin systems, but also more generally, to physical systems described
by matrices, e.g., linearly coupled mechanical oscillators, or linear
electronic circuits.

###### Acknowledgements.

We acknowledge helpful discussions and correspondence with J. Asbóth, A.
Bernevig, L. Fehér, R. Giménez Conejero, M. Hirschmann, L. Oroszlány, and A.
Schnyder. This research was supported by the Ministry of Innovation and
Technology (MIT) and the National Research, Development and Innovation Office
(NKFIH) within the Quantum Information National Laboratory of Hungary and the
Quantum Technology National Excellence Program (Project No.
2017-1.2.1-NKP-2017-00001), by the NKFIH fund TKP2020 IES (Grant No. BME-IE-
NAT) under the auspices of the MIT, by the NKFIH through the OTKA Grants FK
124723 and FK 132146, and by the EU Horizon Europe grants IGNITE and ONCHIPS.
D. V. was supported by the Swedish Research Council (VR) and the Knut and
Alice Wallenberg Foundation.

## References

  * von Neumann and Wigner [1929] J. von Neumann and E. P. Wigner, Über das Verhalten von Eigenwerten bei adiabatischen Prozessen, Physikalische Zeitschrift 30, 467 (1929).
  * Herring [1937] C. Herring, Accidental degeneracy in the energy bands of crystals, Phys. Rev. 52, 365 (1937).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * Fang  _et al._ [2012] C. Fang, M. J. Gilbert, X. Dai, and B. A. Bernevig, Multi-Weyl topological semimetals stabilized by point group symmetry, Phys. Rev. Lett. 108, 266802 (2012).
  * Zhang  _et al._ [2020] T. Zhang, R. Takahashi, C. Fang, and S. Murakami, Twofold quadruple weyl nodes in chiral cubic crystals, Phys. Rev. B 102, 125148 (2020).
  * Hirschmann  _et al._ [2021] M. M. Hirschmann, A. Leonhardt, B. Kilic, D. H. Fabini, and A. P. Schnyder, Symmetry-enforced band crossings in tetragonal materials: Dirac and weyl degeneracies on points, lines, and planes, Phys. Rev. Materials 5, 054202 (2021).
  * Yu  _et al._ [2022] Z.-M. Yu, Z. Zhang, G.-B. Liu, W. Wu, X.-P. Li, R.-W. Zhang, S. A. Yang, and Y. Yao, Encyclopedia of emergent particles in three-dimensional crystals, Science Bulletin 67, 375 (2022).
  * Wernsdorfer and Sessoli [1999] W. Wernsdorfer and R. Sessoli, Quantum phase interference and parity effects in magnetic molecular clusters, Science 284, 133 (1999).
  * Bruno [2006] P. Bruno, Berry phase, topology, and degeneracies in quantum nanomagnets, Phys. Rev. Lett. 96, 117208 (2006).
  * Scherübl  _et al._ [2019] Z. Scherübl, A. Pályi, G. Frank, I. E. Lukács, G. Fülöp, B. Fülöp, J. Nygård, K. Watanabe, T. Taniguchi, G. Zaránd, and S. Csonka, Observation of spin–orbit coupling induced Weyl points in a two-electron double quantum dot, Communications Physics 2, 108 (2019).
  * Frank  _et al._ [2020] G. Frank, Z. Scherübl, S. Csonka, G. Zaránd, and A. Pályi, Magnetic degeneracy points in interacting two-spin systems: Geometrical patterns, topological charge distributions, and their stability, Physical Review B 101, 245409 (2020).
  * Frank  _et al._ [2022] G. Frank, D. Varjas, P. Vrana, G. Pintér, and A. Pályi, Topological charge distributions of an interacting two-spin system, Physical Review B 105, 035414 (2022).
  * Fatemi  _et al._ [2021] V. Fatemi, A. R. Akhmerov, and L. Bretheau, Weyl josephson circuits, Phys. Rev. Research 3, 013288 (2021).
  * Frank  _et al._ [2021] G. Frank, D. Varjas, G. Pintér, and A. Pályi, Weyl-point teleportation, arXiv:2112.14556 (2021).
  * McCann and Koshino [2013] E. McCann and M. Koshino, The electronic properties of bilayer graphene, Reports on Progress in physics 76, 056503 (2013).
  * Péterfalvi  _et al._ [2012] C. G. Péterfalvi, L. Oroszlány, C. J. Lambert, and J. Cserti, Intraband electron focusing in bilayer graphene, New Journal of Physics 14, 063028 (2012).
  * Asbóth  _et al._ [2016] J. K. Asbóth, L. Oroszlány, and A. Pályi,  _A Short Course on Topological Insulators_ (Springer, Heidelberg, 2016).
  * Note [1] For a real analytic function $f:\mathbb{R}^{m}\to\mathbb{R}^{n}$, its complexification is defined as the function $f_{\mathbb{C}}:\mathbb{C}^{m}\to\mathbb{C}^{n}$ that has the same power series as $f$.
  * Bravyi  _et al._ [2011] S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer–Wolff transformation for quantum many-body systems, Annals of physics 326, 2793 (2011).
  * Note [2] Note that our analysis also works in $\mathcal{C}^{\infty}$ (i.e. smooth) category, using the notion of _map germs_ , however the analytic condition simplifies the discussion.
  * Note [3] Similarly, spinless $P\mathcal{T}$ symmetry restricts the Hamiltonian to be real valued, yielding $h_{2}=0$. Hamiltonians with both chiral and $P\mathcal{T}$ symmetry correspond to the case $m=1$ with effective Hamiltonian map $h:(\mathbb{R},0)\to(\mathbb{R},0)$. Such symmetric 1D systems can provide a physical realization of the minimal example with $H(k)=f(k)\sigma_{x}$.
  * Milnor and Weaver [1997] J. Milnor and D. W. Weaver,  _Topology from the differentiable viewpoint_ (Princeton university press, 1997).
  * Arnold  _et al._ [2012] V. I. Arnold, A. N. Varchenko, and S. M. Gusein-Zade,  _Singularities of Differentiable Maps: Volume II Monodromy and Asymptotic Integrals_ , Vol. 83 (Springer Science & Business Media, 2012).
  * Eisenbud [1978] D. Eisenbud, An algebraic approach to the topological degree of a smooth map, Bulletin of the American Mathematical Society 84, 751 (1978).
  * Eisenbud  _et al._ [1977] D. Eisenbud, H. I. Levine, and B. Teissier, An algebraic formula for the degree of a $C^{\infty}$ map germ / Sur une inégalité à la Minkowski pour les multiplicités, Annals of Mathematics 106, 19 (1977).
  * Mond and Nuño-Ballesteros [2020] D. Mond and J. J. Nuño-Ballesteros, Singularities of mappings (2020).
  * Bott  _et al._ [1982] R. Bott, L. W. Tu, _et al._ ,  _Differential forms in algebraic topology_ , Vol. 82 (Springer, 1982).
  * [28] L. Luo, W. Deng, Y. Yang, M. Yan, J. Lu, X. Huang, and Z. Liu, Observation of quadruple Weyl point in hybrid-weyl phononic crystals, arXiv:2203.14600 .
  * [29] Q. Chen, F. Chen, Q. Yan, L. Zhang, Z. Gao, S. A. Yang, Z.-M. Yu, H. Chen, B. Zhang, and Y. Yang, Discovery of a maximally charged Weyl point, arXiv:2203.10722 .