anisotropic Heisenberg spin-$1/2$ chain in the easy-axis regime with weak
anisotropies. Using the framework of algebraic integrability, we have
established that these semi-classical eigenstates emerge classically as
interacting nonlinear spin waves governed by the Landau–Lifshitz field theory
with uniaxial anisotropy. Firstly, we have expressed the asymptotic Bethe
equations in the form of a singular integral equation for the spectral
resolvent, which we subsequently recast in the form of the Riemann–Hilbert
problem, providing jump discontinuity conditions for a double-valued complex
function called quasi-momentum. The latter encodes the moduli of hyperelliptic
Riemann surfaces which provide complete information about the spectrum of
nonlinear eigenmodes for a class of finite-gap solutions of the anisotropic
Landau–Lifshitz ferromagnet.

We have outlined the main ingredients of the algebro-geometric integration
technique. The starting point of this approach is the usual Lax representation
which realises an auxiliary linear problem of parallel transport on a smooth
manifold with a flat connection. In our formulation we made use of the adjoint
representation, enabling us to parametrise the solutions in terms of squared
polynomial eigenfunctions whose zeros contain information about the dynamical
phases evolving on a finite-genus Riemann surface. We have demonstrated how
their dynamics can be mapped to a linear evolution on the Liouville hypertorus
using the Abel–Jacobi transform. Finally, individual components of the
physical spin field can be retrieved with aid of reconstruction formulae.

We have implemented the finite-gap integration procedure for two simplest
classes of solutions: (i) the single-mode (one-cut) solutions, describing
precessional motion around the anisotropy axis, and (ii) the two-mode (two-
cut) solutions which take the shape of elliptic waves. Amongst the two-cut
solutions, there are special elliptic solutions that describe bions, a bound
state of kink and antikink. In a particular singular limit, the bion solution
degenerates into the static kink.

One central result of our work is an algorithm for performing semi-classical
quantisation of classical finite-gap solutions. In this respect, the key
object is the density resolvent associated to the classical quasi-momentum.
The spectral resolvent is supported on a union of one-dimensional segments in
the complex plane which may be adopted as branch cuts of a Riemann surface of
finite genus. By following the programme of Ref. [13], we described and
implemented a numerical algorithm for determining the locations of physical
cuts (associated with the density contours). Each branch cut is a magnon
condensate that represents a nonlinear mode in the spectrum of the effective
classical equation of motion. In this view, semi-classical quantisation
amounts to dissolve each branch cut of a finite-gap quasi-momentum into a
large but finite number of magnetisation quanta (carried by magnons) with a
prescribed accuracy; the resolvent density along each contour specifies a
local density of magnon excitations. When the local density of magnons exceeds
a critical threshold value, the physical contours experience a certain ‘non-
perturbative effect’ which leads to the formation of special condensates with
uniform unit density of Bethe roots. A proper resolution of such situations
necessitates to take into account quantum corrections.

In this work we devote most attention to various formal properties of semi-
classical eigenstates and other related mathematical underpinnings. We hope
this can provide a foundation for further developments which would ultimately
pave the way to physical applications. Particularly in the domain of out-of-
equilibrium dynamics there has been tremendous progress recently in employing
integrability techniques that enabled us, among others, to study late-time
relaxation dynamics from highly-excited many-body initial states which goes
commonly under the name of ‘quantum quenches’ [3], see also [81, 82, 83]). A
particularly useful tool in this regard is the functional integral
representation, dubbed the Quench Action [84, 85], which exploits exact
knowledge of thermodynamic overlap coefficients. Our hope is to obtain its
semi-classical counterpart. Despite that general expressions for the overlap
coefficients between a semi-classical eigenstates and an off-shell state are
explicitly known due to Ref. [32, 33, 71], we have not been successful in
employing them in practice yet. We have nonetheless been able to provide
several benchmarks for the computation of Gaudin norms and Slavnov overlaps
for a few simple finite-gap solutions and found good convergence.

The main difficulty when dealing with the overlap formulae was to satisfy the
requirements for the contour prescription. Since avoiding all the branch cuts
of the integrand does not appear to be easily overcome, it seems that an
alternative formula based solely on the resolvent densities (similarly to that
for the semi-classical limit of the Gaudin norm [31]) might be preferable.
Overcoming this issue would be a stepping stone for formulating a quench
problem at the level of semi-classical states, a prominent example of which
would be the semi-classical version of the domain wall melting which has
recently been solved analytically by the authors in [16] using the inverse
scattering transformation. This could help to solidify the classical–quantum
correspondence also brought forward in [16] and corroborated in [17].

Lastly, we shortly examined the structure of correlation functions in the
semi-classical eigenstates of the Heisenberg XXZ chain and compared them to
their classical counterparts, namely correlators of classical fields as
finite-gap solutions. We found empirical evidence for a classical–quantum
correspondence between static multipoint correlators on both sides, in
alignment with the earlier results of Ref. [79]. Importantly, the semi-
classical correlators can only be compared to correlators of classical fields
after computing phase-space averages, as demonstrated on a few basic examples.
While there are strong indications that such a correspondence should hold
generally in quantum integrable models that possess (integrable) classical
limits, a proof is still lacking at the moment. We believe that it would be
fruitful to investigate this matter in the framework of quantum separated
variables, see e.g. [53, 86, 87], to learn how classical separated variables
on Riemann surfaces [36] emerge from the microscopic quantum model that sits
underneath. In our opinion, quantum integrable lattice models and spin chains
provide paradigmatic examples to address these aspects.

## Acknowledgements

We are very grateful to Jean-Sébastien Caux for his involvement at the
starting stage of the project and numerous discussions. We thank Filippo
Colomo, Andrea De Luca, Andrii Liashyk, Grégoire Misguich, Vincent Pasquier,
and Dmytro Volin for valuable discussions. We are greatly indebted to Nikolay
Gromov and Ivan Kostov for sharing their insights on various facets of the
problem. Y.M. thanks Vincenzo Alba and Alvise Bastianello for useful
suggestions on numerical implementations.

Y.M., and O.G. acknowledge the support from the European Research Council
under ERC Advanced grant 743032 DYNAMINT. E.I. is supported by the research
programme P1-0402 of Slovenian Research Agency.

Appendices

##  Appendix A Riemann-Hilbert problem in $\zeta$-plane

In order to study the formation of condensates, the Riemann–Hilbert problem is
most conveniently written in terms of spectral parameter $\zeta=1/\mu$, namely

$\mathfrak{p}(\zeta+{\rm i}0)+\mathfrak{p}(\zeta-{\rm i}0)=2\pi
n_{j},\quad\zeta\in\mathcal{C}_{j},$ (A.1)

where $\mathcal{C}_{j}$ denotes the $j$-th branch cut in $\zeta$ plane,
whereas quasi-momentum $\mathfrak{p}(\zeta)$ is defined as

$\mathfrak{p}(\zeta)=G(\zeta)-\frac{\ell}{2\zeta}=\ell\int{\rm
d}\xi\tilde{\mathcal{K}}_{\delta}(\zeta,\xi)\rho(\xi)-\frac{\ell}{2\zeta},$
(A.2)

with integration kernel

$\tilde{\mathcal{K}}_{\delta}(\zeta,\xi)=\frac{1+\delta\xi\zeta}{\zeta-\xi}.$
(A.3)

The density (of the Bethe roots) is accordingly given by

$\rho(\zeta)=\frac{\mathfrak{p}(\zeta+{\rm i}0)-\mathfrak{p}(\zeta-{\rm
i}0)}{2{\rm i}\pi\ell(1+\delta\zeta^{2})},\quad\zeta\in\mathcal{C}_{j}.$ (A.4)

Note that the orientation of integration along $\mathcal{C}_{j}$ is now in the
opposite direction as previously, i.e. it goes from the branch point with
negative imaginary part to the one with positive imaginary part.

##  Appendix B Finite size corrections to Riemann-Hilbert problem

Here we outline how to take the semi-classical limit of the logarithm of
$Q_{j}^{[\pm 2]}$. The first step is to split the term into the anomalous part
and normal part [55, 88], i.e.

$\begin{split}\log Q_{j}^{[\pm 2]}(\vartheta_{j})&=\sum_{k\neq
j}^{M}\log\sin(\vartheta_{j}-\vartheta_{k}\pm{\rm i}\eta)\\\
&=\sum_{0<|k-j|\leq K}\log\sin(\vartheta_{j}-\vartheta_{k}\pm{\rm
i}\eta)+\sum_{|k-j|>K}\log\sin(\vartheta_{j}-\vartheta_{k}\pm
i\eta),\end{split}$ (B.1)

where parameter $K$ is a cut-off with the following properties,

$\vartheta_{j}-\vartheta_{k}\sim\begin{cases}\mathcal{O}\left(1/L\right),\quad|k-j|\leq
K\\\ \mathcal{O}\left(1\right),\quad|k-j|>K\\\ \end{cases}.$ (B.2)

We denote the anomalous part as

$\log Q_{j}^{a}(\vartheta_{j}\pm{\rm i}\eta)=\sum_{0<|k-j|\leq
K}\log\sin(\vartheta_{j}-\vartheta_{k}\pm{\rm i}\eta),$ (B.3)

while the normal part is

$\log Q_{j}^{n}(\vartheta_{j}\pm{\rm
i}\eta)=\sum_{|k-j|>K}\log\sin(\vartheta_{j}-\vartheta_{k}\pm{\rm i}\eta).$
(B.4)

For the normal part, we can perform the same expansion as in Eq. (3.12),
namely

$\begin{split}\log Q_{j}^{n}(\vartheta_{j}\pm{\rm i}\eta)=\log
Q_{j}^{n}(\vartheta_{j})\pm{\rm i}\eta\frac{d}{d\vartheta}\log
Q_{j}^{n}(\vartheta)\rvert_{\vartheta=\vartheta_{j}}\\\
-\frac{\eta^{2}}{2}\frac{d^{2}}{d\vartheta^{2}}\log
Q_{j}^{n}(\vartheta)\rvert_{\vartheta=\vartheta_{j}}+\mathcal{O}\left(\frac{1}{L^{2}}\right),\end{split}$
(B.5)

and

${\rm i}\eta\frac{d}{d\vartheta}\log
Q_{j}^{n}(\vartheta)\rvert_{\vartheta=\vartheta_{j}}=\frac{\epsilon\ell}{L}\sum_{|k-j|>K}\frac{1}{\tan(\vartheta_{j}-\vartheta_{k})}=\frac{\ell}{L}\sum_{|k-j|>K}\frac{\mu_{j}\mu_{k}+\delta}{\mu_{j}-\mu_{k}}.$
(B.6)

Combining the two parts, we obtain

$\log Q_{j}^{n}(\vartheta_{j}+{\rm i}\eta)-\log Q_{j}^{n}(\vartheta_{j}-{\rm
i}\eta)=\frac{2\ell}{L}\sum_{|k-j|>K}\frac{\mu_{j}\mu_{k}+\delta}{\mu_{j}-\mu_{k}}+\mathcal{O}\left(\frac{1}{L^{2}}\right).$
(B.7)

Meanwhile, for the anomalous part, denoting $m=k-j$, we have

$\begin{split}\log Q_{j}^{n}(\vartheta_{j}+{\rm i}\eta)-\log
Q_{j}^{n}(\vartheta_{j}-{\rm
i}\eta)&=\sum_{0<|m|<K}\log\frac{\sin(\vartheta_{j}-\vartheta_{j+m}+{\rm
i}\eta)}{\sin(\vartheta_{j}-\vartheta_{j+m}+{\rm i}\eta)}\\\
&=\sum_{0<|m|<K}\log\frac{L(\mu_{j}-\mu_{j+m})+{\rm
i}\ell(\mu_{j}^{2}+\delta)}{L(\mu_{j}-\mu_{j+m})-{\rm
i}\ell(\mu_{j}^{2}+\delta)}.\end{split}$ (B.8)

We can develop an expansion

$L\mu_{j+m}\sim
c_{1}L+c_{2}m+\frac{1}{2}\frac{c_{3}m^{2}}{L}+\mathcal{O}\left(\frac{1}{L^{2}}\right),\quad|m|\leq
K,$ (B.9)

where all the “constants” can be expressed in terms of density $\rho(\mu)$,
i.e.

$c_{1}=\mu_{j},\quad c_{2}=\frac{1}{\rho(\mu_{j})},\quad
c_{3}=-\frac{\rho^{\prime}(\mu_{j})}{\rho(\mu_{j})^{3}},$ (B.10)

and

$\rho(\mu)=\frac{1}{L}\sum_{j=1}^{M}\delta(\mu-\mu_{j}),\quad\rho(\mu)\simeq\frac{{\rm
d}j}{{\rm d}\mu}.$ (B.11)

By combining the $m$-th and $(-m)$-th terms in the sum, we can express the
leading order of the sum as

$\sum_{m=1}^{K}\frac{1}{{\rm i}}\left(\log\frac{L(\mu_{j}-\mu_{j-m})+{\rm
i}\ell(\mu_{j}^{2}+\delta)}{L(\mu_{j}-\mu_{j-m})-{\rm
i}\ell(\mu_{j}^{2}+\delta)}+\log\frac{L(\mu_{j}-\mu_{j+m})+{\rm
i}\ell(\mu_{j}^{2}+\delta)}{L(\mu_{j}-\mu_{j+m})-{\rm
i}\ell(\mu_{j}^{2}+\delta)}\right),$ (B.12)

using

$\begin{split}&\frac{1}{{\rm i}}\left(\log\frac{L(\mu_{j}-\mu_{j-m})+{\rm
i}\ell(\mu_{j}^{2}+\delta)}{L(\mu_{j}-\mu_{j-m})-{\rm
i}\ell(\mu_{j}^{2}+\delta)}+\log\frac{L(\mu_{j}-\mu_{j+m})+{\rm
i}\ell(\mu_{j}^{2}+\delta)}{L(\mu_{j}-\mu_{j+m})-{\rm
i}\ell(\mu_{j}^{2}+\delta)}\right)\\\ &=\frac{1}{{\rm
i}}\log\frac{b^{2}m^{2}-[{\rm
i}\ell(\mu_{j}^{2}+\delta)-\frac{c_{3}m^{2}}{2L}]^{2}}{c_{2}^{2}m^{2}-[{\rm
i}\ell(\mu_{j}^{2}+\delta)+\frac{c_{3}m^{2}}{2L}]^{2}}\\\
&=\frac{2c_{3}\ell(\mu_{j}^{2}+\delta)}{c_{2}^{2}L}\left(1-\frac{1}{\frac{c_{2}^{2}}{\ell^{2}(\mu_{j}^{2}+\delta)^{2}}+1}\right)+\mathcal{O}\left(\frac{1}{L^{2}}\right).\end{split}$
(B.13)

The first part can be combined with the sum for $|m|>K$, since

$\frac{2\ell(\mu_{j}\mu_{j-m}+\delta)}{L(\mu_{j}-\mu_{j-m})}+\frac{2\ell(\mu_{j}\mu_{j+m}+\delta)}{L(\mu_{j}-\mu_{j+m})}\simeq\frac{2c_{3}(\mu_{j}^{2}+\delta)}{c_{2}^{2}L}.$
(B.14)

Taking the limit $K\to\infty$ (beware that $K/L\to 0$), for the second part we
have

$-\sum_{m=1}^{\infty}\frac{2c_{3}\ell(\mu_{j}^{2}+\delta)}{c_{2}^{2}L\left[\frac{c_{2}^{2}}{\ell^{2}(\mu_{j}^{2}+\delta)^{2}}+1\right]}=\frac{c_{3}\ell(\mu_{j}^{2}+\delta)}{c_{2}^{2}L}\left[1-\frac{\pi\ell(\mu_{j}^{2}+\delta)}{c_{2}}\coth\left(\frac{\pi\ell(\mu_{j}^{2}+\delta)}{c_{3}}\right)\right].$
(B.15)

Substituting back in the values in Eq. (B.10), we will obtain the finite-size
correction in Eq. (3.18).

In addition, the finite-size correction in terms of $\zeta$ variable takes the
form

$\frac{\pi\rho^{\prime}(\zeta)\ell^{2}(1+\delta\zeta^{2})^{2}}{L}\coth\left[\pi\ell(1+\delta\zeta^{2})\rho(\zeta)\right]+\mathcal{O}\left(\frac{1}{L^{2}}\right).$
(B.16)

##  Appendix C Useful formulae for elliptic functions

We collect several useful functions and formulae used in the derivations in
Section 5.2.

We begin by defining the elliptic integral of the first kind

$K({\rm k}^{2})=\int_{0}^{1}\frac{{\rm d}x}{\sqrt{(1-x^{2})(1-{\rm
k}^{2}x^{2})}}.$ (C.1)

The Jacobi elliptic function $\mathrm{sn}(x,{\rm k}^{2})$ is defined as the
inverse of the elliptic integral of the first kind,

${\rm w}=\mathrm{sn}(x,{\rm k}^{2}),\quad x=\int_{0}^{\rm w}\frac{{\rm
d}z}{\sqrt{(1-z^{2})(1-{\rm k}^{2}z^{2})}},$ (C.2)

and, without ambiguity, we can put $\mathrm{sn}(x,{\rm
k}^{2})=:\mathrm{sn}(x)$. Other types of Jacobi elliptic functions can be
defined in a similar way,

${\rm w}=\mathrm{cn}(x,{\rm k}^{2}),\quad x=\int_{\rm w}^{1}\frac{{\rm
d}z}{\sqrt{(1-z^{2})(1-{\rm k}^{2}+{\rm k}^{2}z^{2})}},$ (C.3)

and

${\rm w}=\mathrm{dn}(x,{\rm k}^{2}),\quad x=\int_{\rm w}^{1}\frac{{\rm
d}z}{\sqrt{(1-z^{2})(z^{2}+{\rm k}^{2}-1)}},$ (C.4)

such that

${\rm sn}^{2}x+\mathrm{cn}^{2}x=1,\quad{\rm k}^{2}{\rm
sn}^{2}x+\mathrm{dn}^{2}x=1.$ (C.5)

When shifting the argument by one quarter of the period of $K({\rm k^{2}})$,
we have

$\mathrm{cn}\left(x+K({\rm k^{2}})\right)=-\sqrt{1-{\rm
k}^{2}}\frac{\mathrm{sn}(x)}{\mathrm{dn}(x)},\quad\mathrm{cn}\left(x+K({\rm
k^{2}})\right)=\sqrt{1-{\rm k}^{2}}\frac{1}{\mathrm{dn}(x)}.$ (C.6)

In addition, we also make use of theta functions to express the spin field.
The most important one here is

$\vartheta_{3}(z,\tau)=\sum_{n=-\infty}^{+\infty}e^{{\rm i}\pi\tau n^{2}+2{\rm
i}zn}.$ (C.7)

For a more detailed exposition and other properties of elliptic functions we
refer the reader to Refs. [89, 90].

##  Appendix D Numerical tests

###  D.1 Gaudin norm

We present the data for several numerical checks. Firstly, we computed the
Gaudin norm of a one-cut solution without a condensate. Secondly, we include a
condensate, and consider two regimes: (i) isotropic interaction with
$\delta=0$ and (a) anisotropic regimes with $\delta>0$. Case (i) without a
condensate has been studied in Ref. [31], and we use it as a benchmark. Case
(a) is more interesting, as it enables a non-trivial quantitative confirmation
of our proposal for determining the location of a condensate or additional
contours, cf. Sec. 6.2.1.

We next present our numerical results for the Gaudin norm computed on the one-
cut solution with mode number $n=1$ and filling fraction $\nu=0.1$, for both
cases (i) and (a). For this choice of parameters, there is no condensates
involved. A linear fit on the finite-size numerical data yields

$\log\mathcal{N}-\frac{1}{2}\log
L=0.00714654(1)\,L+0.068763(9),\quad\delta=0,$ (D.1)

and

$\log\mathcal{N}-\frac{1}{2}\log L=0.0083405(3)\,L+0.083756(1),\quad\delta=1.$
(D.2)

Figure 12: Logarithm of the Gaudin norm, shown for the one-cut solution with
$n=1$ and $\nu=0.1$. Green circles (black crosses) show numerical results for
$\delta=0$ ($\delta=1$), respectively. Linear fits are indicated by dashed
lines. Figure 13: Logarithm of the Gaudin norm, shown for the one-cut solution
with $n=1$, $\nu=\frac{1}{3}$. Green (black) circles show numerical results
for $\delta=0$ ($\delta=1$). Linear fits are indicated by dashed lines.

Comparing these results to those obtained from the functional integral
approach, cf. Eq. (7.2) (denoted by $C_{1}$ in the table below), we have

|  $C_{1}$ numerical |  $C_{1}$ functional
---|---|---
$\delta=0$ |  0.00714654(1) |  0.007156(1)
$\delta=1$ |  0.0083405(3) |  0.008383(8)

We can see that the functional approach (only requiring the knowledge of the
density contours) yields very accurate results in both the isotropic (i) and
anisotropic (a) case.

Next up, we analyse the cases (i) and (a) with an extra condensate, computing
the Gaudin norm for a one-cut solution with mode number $n=1$ and filling
fraction $\nu=\frac{1}{3}$. In the functional integral approach, this amounts
to compute the integral in Eq. (7.2) along a contour $\mathcal{C}$ comprising
of three parts,

$\mathcal{C}=\mathcal{C}_{1}+\mathcal{C}_{2}+\mathcal{C}_{\rm cond}.$ (D.3)

Here $\mathcal{C}_{1}$ pertains to the original contour with density
$\rho_{1}(\zeta)$ in Eq. (6.7) (green dashed line in panel (b) in Fig. 7),
whereas contour $\mathcal{C}_{2}$ has density $\rho_{2}(\zeta)$, depicted in
Eq. (6.7) by yellow dashed line in panel (b) of Fig. 7. Finally,
$\mathcal{C}_{\rm cond}$ is the straight condensate contour with density
$\rho_{cond}=\frac{i}{1+\delta\zeta^{2}}$.

This time, a linear extrapolation of the numerical finite-size data yields

$\log\mathcal{N}-\frac{1}{2}\log L=0.091273(6)\,L+1.92813(4),\quad\delta=0,$
(D.4)

and

$\log\mathcal{N}-\frac{1}{2}\log L=0.082597(1)L+2.09809(5),\quad\delta=1,$
(D.5)

while comparing to the results of the functional integral approach, see Eq.
(7.2) ($C_{1}$ in the table below), we obtain

|  $C_{1}$ numerical |  $C_{1}$ functional
---|---|---
$\delta=0$ |  0.091273(6) |  0.091121(9)
$\delta=1$ |  0.082597(1) |  0.081761(2)

In spite of an extra condensate, the functional integral method yields very
accurate results. Even more importantly, this check provides a robust
confirmation for the additional condensate contour(s). We note that any
different contour, e.g. the usual arc-shaped contour without a condensate,
produces an appreciable numerical mismatch.

###  D.2 Slavnov overlap

(a)

(b)

Figure 14: (a) Integration contours for computing the overlap coefficient
between a one-cut state and the vacuum descendant. Square-root branch cut of
$\mathfrak{p}(\zeta)$ is indicated by blue dashed line, whereas all the
additional branch cuts due to the dilogarithm function are marked by green
dashed line. The integration contour, marked in red (orange) on the upper
(bottom) Riemann sheets, escape to infinity on the bottom sheet. (b) Logarithm
of the overlap coefficient between the one-cut solution with $n=1$ and
$\nu=0.1$ and a vacuum descendant state. Green circles (black crosses) show
numerical data for interaction parameter $\delta=0$ ($\delta=1$),
respectively, with dashed lines corresponding to linear fits.

In this section, we present a numerical check of an overlap formula. Here we
compute the overlap between a semi-classical Bethe eigenstate with a single
cut and a vacuum descendant state, which is a “domain-wall state” of the form
$|\downarrow\cdots\downarrow\uparrow\cdots\uparrow\rangle$. Again, we perform
computations for both the isotropic case (i) at $\delta=0$ and for the
anisotropic interaction (a) by setting the anisotropy parameter to $\delta=1$.

The (unnormalised) overlap can be obtained from the general Algebraic Bethe
ansatz determinant formula due to Slavnov with help of L’Hôpital rule
(presented previously in e.g. [91]),

$\mathcal{V}=\langle\phi|\\{\vartheta\\}\rangle=\prod_{l=1}^{M}\sin\left(\vartheta_{l}+i\frac{\eta}{2}\right)^{M}\frac{\det
H}{\prod_{j<k}\sin(\vartheta_{j}-\vartheta_{k})},$ (D.6)
$H_{ab}=\cot\left(\vartheta_{a}-i\frac{\eta}{2}\right)^{b}-\cot\left(\vartheta_{a}+i\frac{\eta}{2}\right)^{b},$
(D.7)

taking the “domain-wall state”
$|\phi\rangle=|\downarrow\cdots\downarrow\uparrow\cdots\uparrow\rangle$ with
$M$ down-turned spins and lattice size $L$. In the isotropic limit, we obtain

$\mathcal{V}=\langle\phi|\\{\lambda\\}\rangle=\prod_{l=1}^{M}\left(\lambda_{l}+\frac{i}{2}\right)^{M}\frac{\det
H}{\prod_{j<k}(\lambda_{j}-\lambda_{k})},$ (D.8)
$H_{ab}=\left(\frac{1}{\lambda_{a}-i/2}\right)^{b}-\left(\frac{1}{\lambda_{a}+i/2}\right)^{b}.$
(D.9)

In the following we shall ignore the phase and consider only the absolute
value of the overlap. We expect, similarly as previously for the Gaudin norm,
the following behavior at large $L$,

$\log\left|\mathcal{V}\right|=C_{2}L+\frac{1}{2}\log L+\mathcal{O}(1),\quad
C_{2}\in\mathcal{O}(1).$ (D.10)

The results of computations are shown in Fig. 14. By numerically fitting the
finite-size data, we obtained

$\log|\mathcal{V}|-\frac{1}{2}\log L=-0.144278(5)L-0.272812(7),\quad\delta=0,$
(D.11)

and

$\log|\mathcal{V}|-\frac{1}{2}\log L=-0.143827(7)L-0.254592(6),\quad\delta=1.$
(D.12)

Before we can repeat the computation using Eq. (7.6), we need to find the
“quasi-momentum” corresponding to the vacuum descendant $|\phi\rangle$. With
aid of Algebraic Bethe ansatz, a vacuum descendant state (with no
inhomogeneities) is given by [91]

$|\phi\rangle=|\underbrace{\downarrow\cdots\downarrow}_{M}\uparrow\cdots\uparrow\rangle=\lim_{\xi_{j}\to
0}\prod_{j=1}^{M}B(\xi_{j})|0\rangle,\quad|0\rangle=|\uparrow\cdots\uparrow\rangle,$
(D.13)

where $|0\rangle$ is the ferromagnetic Bethe vacuum
$|\uparrow\cdots\uparrow\rangle$, and $B(\lambda)$ is the magnon excitation
operator corresponding to the upper off-diagonal element of the quantum
monodromy matrix.

The density of the “off-shell Bethe roots” can then be expressed as

$\rho_{\phi}(\zeta)=\nu_{1}\delta(\zeta),\quad\xi_{1},\cdots\xi_{M}\to 0.$
(D.14)

With no loss of generality, we set subsequently the classical period to
$\ell=1$. The quasi-momentum associated to $|\phi\rangle$ then reads

$\mathfrak{p}_{\phi}(\zeta)=\ell\int_{\mathcal{C}_{\phi}}d\zeta^{\prime}\frac{\rho(\zeta^{\prime})(1+\delta\zeta\zeta^{\prime})}{\zeta-\zeta^{\prime}}-\frac{\ell}{2\zeta}=\frac{\nu_{1}-1/2}{\zeta}.$
(D.15)

Now we are ready to employ the functional integral formula (7.6). The
appropriate choice of contours is shown in panel (a) in Fig. 14 131313The
rationale behind this choice is to avoid the branch cuts of the dilogarithm
function. More details on this can be found in [33].. Comparing the two
computations of coefficient $C_{2}$, we find

|  $C_{2}$ numerical |  $C_{2}$ functional
---|---|---
$\delta=0$ |  -0.144278(5) |  -0.144485(3)
$\delta=1$ |  -0.143827(7) |  -0.142267(2)

Once again the computation using formula (7.6) works quite well, both in the
isotropic and anisotropic cases.

##  Appendix E Numerical recipes

###  E.1 Numerical solution to Bethe equations

We outline how to numerically solve for the Bethe roots to equations (3.3)
(for finite but possibly large system length $L$) for a specific class of
quantum eigenstates that in the thermodynamic limit become one-cut classical
solution. To this end, we employ the algorithm described in Section 7 of Ref.
[13] for the rational Bethe equations (i.e. for isotropic interaction,
$\delta=0$).

We use this method in combination with another method, given in Appendix C of
Ref. [71], where the solution to the isotropic chain is used as the initial
condition for the Newton-Raphson iteration during which the anisotropy
parameter gets gradually increased. However, while this procedure works quite
well for the simplest case of mode number $n=1$, we could not achieve good
convergence for mode numbers $n>2$ and consequently could not perform any
benchmark on classical solutions with two or more cuts.

###  E.2 Determining branch points from filling fractions and mode numbers

We describe a numerical procedure to determine the branch points from a given
set of moduli, that is the mode numbers and filling fractions. The method is
completely general and applies to solutions with an arbitrary number of cuts.
For simplicity however, we demonstrate it below on the class of two-cut
solutions.

There are four branch points $\\{\mu_{1}$, $\bar{\mu}_{1}$, $\mu_{2}$,
$\bar{\mu}_{2}\\}$ that appear in complex-conjugate pairs. Thus, there are in
total four real parameters (real and imaginary components of each branch
point). The finite-gap solution is parametrised by equivalently four
parameters, i.e. mode numbers and filling fractions of both branch cuts,
relating to the previous four parameters in a nonlinear manner. In addition,
the solution must be periodic with the period $\ell$, which adds an additional
constraint.

We would like to remark that, unlike in the one-cut case, the determination is
highly nonlinear, related to elliptic functions and integrals. Hence, there is
not a simple analytic closed-form formula available in this case. Instead, we
are going to use the following numerical procedure:

  * •

We first fix the real part of branch points of the first cut to
$a\equiv\mathrm{Re}\,\mu_{1}=1$.

  * •

We next scan a range of values $\mathrm{Im}\mu_{1}\in(b_{1},b_{2})$ and
$\mathrm{Re}\,\mu_{2}\in(c_{1},c_{2})$, and find the value of
$\mathrm{Im}\,\mu_{2}$ that yields the required mode numbers $n_{1}$ and
$n_{2}$. More specifically, for any $\mathrm{Im}\,\mu_{2}$, we compute
$\frac{d\mathfrak{p}}{\ell}$ by demanding the $\mathcal{A}$-cycle to vanish.
Since the classical period reads

$\ell=\frac{2\pi n_{1}}{\int_{\mathcal{B}_{1}}{\rm d}\mathfrak{p}/\ell},$
(E.1)

we can numerically determine $\mathrm{Im}\,\mu_{2}$ from the requirement

$\int_{\mathcal{B}_{2}}{\rm d}\mathfrak{p}=2\pi n_{2}.$ (E.2)

  * •

Having done the above, we can readily compute following quantities,

$\tilde{\nu}_{1}=\oint_{\mathcal{A}_{1}}\frac{{\rm
d}\mathfrak{p}}{\mu},\quad\tilde{\nu}_{2}=\oint_{\mathcal{A}_{2}}\frac{{\rm
d}\mathfrak{p}}{\mu},$ (E.3)

which moreover depend on $\mathrm{Im}\,\mu_{1}$ and $\mathrm{Re}\,\mu_{2}$.

  * •

By requiring that $\tilde{\nu}_{1}=\nu_{1}$, we obtain a “curve” in the plane
spanned by $\mathrm{Im}\mu_{1}$ and $\mathrm{Re}\,\mu_{2}$. By finally
requesting also that $\tilde{\nu}_{1}=\nu_{2}$, we are left with a single
point, say $(b,c)$. The last point, call it $d$, corresponds to
$\mathrm{Im}\,\mu_{2}$.

  * •

We have thus determined to complex branch points $\mu_{1}=1+bi$,
$\mu_{2}=c+di$ (alongside their complex conjugates) which yields the
prescribed classical period $\ell$, mode number $n_{1}$, $n_{2}$ and filling
fraction $\nu_{1}$, $\nu_{2}$ of a general two-solution.

In making a comparison with the Bethe root distributions, we normally prefer
to set the classical period to $\ell=1$. In this case we simply divide the
above branch points by $\ell$, see Eq. (E.1), that is

$\mu_{1,n}=\frac{1}{\ell}+\frac{b}{\ell}i,\quad\mu_{2,n}=\frac{c}{\ell}+\frac{d}{\ell}i,$
(E.4)

or equivalently in terms of spectral parameter $\zeta=1/\mu$,

$\zeta_{1,n}=\frac{1}{\bar{\mu}_{1,n}},\quad\zeta_{2,n}=\frac{1}{\bar{\mu}_{2,n}}.$
(E.5)

## References

  * [1] L. D’Alessio, Y. Kafri, A. Polkovnikov and M. Rigol,  _From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics_ ,  Adv. Phys. 65(3), 239 (2016),  10.1080/00018732.2016.1198134.
  * [2] C. Gogolin and J. Eisert,  _Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems_ ,  Rep. Prog. Phys. 79(5), 056001 (2016),  10.1088/0034-4885/79/5/056001.
  * [3] F. H. L. Essler and M. Fagotti,  _Quench dynamics and relaxation in isolated integrable quantum spin chains_ ,  J. Stat. Mech. 2016(6), 064002 (2016),  10.1088/1742-5468/2016/06/064002.
  * [4] A. De Luca and G. Mussardo,  _Equilibration properties of classical integrable field theories_ ,  J. Stat. Mech. 2016(6), 064011 (2016),  10.1088/1742-5468/2016/06/064011.
  * [5] R. Vasseur and J. E. Moore,  _Nonequilibrium quantum dynamics and transport: from integrability to many-body localization_ ,  J. Stat. Mech. 2016(6), 064010 (2016),  10.1088/1742-5468/2016/06/064010.
  * [6] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg and M. Žnidarič,  _Finite-temperature transport in one-dimensional quantum lattice models_ ,  arXiv preprint arXiv:2003.03334 (2020).
  * [7] V. Alba and P. Calabrese,  _Entanglement and thermodynamics after a quantum quench in integrable systems_ ,  Proc. Natl. Acad. Sci. U.S.A. 114(30), 7947 (2017),  10.1073/pnas.1703516114.
  * [8] O. A. Castro-Alvaredo, B. Doyon and T. Yoshimura,  _Emergent Hydrodynamics in Integrable Quantum Systems Out of Equilibrium_ ,  Phys. Rev. X 6, 041065 (2016),  10.1103/PhysRevX.6.041065.
  * [9] B. Bertini, M. Collura, J. De Nardis and M. Fagotti,  _Transport in Out-of-Equilibrium $XXZ$ Chains: Exact Profiles of Charges and Currents_,  Phys. Rev. Lett. 117, 207201 (2016),  10.1103/PhysRevLett.117.207201.
  * [10] J. A. Minahan and K. Zarembo,  _The Bethe-ansatz for $\mathcal{N}=4$ super Yang-Mills_,  J. High Energy Phys. 2003(03), 013 (2003),  10.1088/1126-6708/2003/03/013.
  * [11] M. Kruczenski,  _Spin Chains and String Theory_ ,  Phys. Rev. Lett. 93, 161602 (2004),  10.1103/PhysRevLett.93.161602.
  * [12] V. A. Kazakov, A. Marshakov, J. A. Minahan and K. Zarembo,  _Classical/quantum integrability in AdS/CFT_ ,  J. High Energy Phys. 2004(05), 024 (2004),  10.1088/1126-6708/2004/05/024.
  * [13] T. Bargheer, N. Beisert and N. Gromov,  _Quantum stability for the Heisenberg ferromagnet_ ,  New J. Phys. 10(10), 103023 (2008),  10.1088/1367-2630/10/10/103023.
  * [14] L. D. Faddeev and L. A. Takhtajan,  _Hamiltonian Methods in the Theory of Solitons_ ,  Springer Berlin Heidelberg,  10.1007/978-3-540-69969-9 (1987).
  * [15] G. Misguich, K. Mallick and P. L. Krapivsky,  _Dynamics of the spin- $\frac{1}{2}$ Heisenberg chain initialized in a domain-wall state_,  Phys. Rev. B 96, 195151 (2017),  10.1103/PhysRevB.96.195151.
  * [16] O. Gamayun, Y. Miao and E. Ilievski,  _Domain-wall dynamics in the Landau-Lifshitz magnet and the classical-quantum correspondence for spin transport_ ,  Phys. Rev. B 99, 140301 (2019),  10.1103/PhysRevB.99.140301.
  * [17] G. Misguich, N. Pavloff and V. Pasquier,  _Domain wall problem in the quantum XXZ chain and semiclassical behavior close to the isotropic point_ ,  SciPost Phys. 7, 25 (2019),  10.21468/SciPostPhys.7.2.025.
  * [18] M. Žnidarič,  _Spin Transport in a One-Dimensional Anisotropic Heisenberg Model_ ,  Phys. Rev. Lett. 106, 220601 (2011),  10.1103/PhysRevLett.106.220601.
  * [19] M. Ljubotina, M. Žnidarič and T. Prosen,  _Spin diffusion from an inhomogeneous quench in an integrable system_ ,  Nat. Comm. 8(1) (2017),  10.1038/ncomms16117.
  * [20] E. Ilievski, J. De Nardis, M. Medenjak and T. Prosen,  _Superdiffusion in One-Dimensional Quantum Lattice Models_ ,  Phys. Rev. Lett. 121, 230602 (2018),  10.1103/PhysRevLett.121.230602.
  * [21] S. Gopalakrishnan and R. Vasseur,  _Kinetic Theory of Spin Diffusion and Superdiffusion in $XXZ$ Spin Chains_,  Phys. Rev. Lett. 122, 127202 (2019),  10.1103/PhysRevLett.122.127202.
  * [22] J. De Nardis, M. Medenjak, C. Karrasch and E. Ilievski,  _Anomalous Spin Diffusion in One-Dimensional Antiferromagnets_ ,  Phys. Rev. Lett. 123, 186601 (2019),  10.1103/PhysRevLett.123.186601.
  * [23] J. De Nardis, S. Gopalakrishnan, E. Ilievski and R. Vasseur,  _Superdiffusion from Emergent Classical Solitons in Quantum Spin Chains_ ,  Phys. Rev. Lett. 125, 070601 (2020),  10.1103/PhysRevLett.125.070601.
  * [24] T. Prosen and B. Žunkovič,  _Macroscopic Diffusive Transport in a Microscopically Integrable Hamiltonian System_ ,  Phys. Rev. Lett. 111, 040602 (2013),  10.1103/PhysRevLett.111.040602.
  * [25] E. Ilievski and J. De Nardis,  _Microscopic Origin of Ideal Conductivity in Integrable Quantum Models_ ,  Phys. Rev. Lett. 119, 020602 (2017),  10.1103/PhysRevLett.119.020602.
  * [26] J. De Nardis, D. Bernard and B. Doyon,  _Diffusion in generalized hydrodynamics and quasiparticle scattering_ ,  SciPost Phys. 6, 49 (2019),  10.21468/SciPostPhys.6.4.049.
  * [27] S. Gopalakrishnan, D. A. Huse, V. Khemani and R. Vasseur,  _Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems_ ,  Phys. Rev. B 98, 220303 (2018),  10.1103/PhysRevB.98.220303.
  * [28] B. Sutherland,  _A Brief History of the Quantum Soliton with New Results on the Quantization of the Toda Lattice_ ,  Rocky Mt. J. Math. 8(1/2), 413 (1978).
  * [29] B. A. Dubrovin,  _Theta functions and non-linear equations_ ,  Russ. Math. Surv. 36(2), 11 (1981),  10.1070/rm1981v036n02abeh002596.
  * [30] E. D. Belokolos, A. I. Bobenko, V. Z. Enol’skii, A. R. Its and M. V. B,  _Algebro-geometric approach to nonlinear integrable equations_ ,  Springer (1994).
  * [31] N. Gromov, A. Sever and P. Vieira,  _Tailoring three-point functions and integrability III. Classical tunneling_ ,  J. High Energy Phys. 2012(7) (2012),  10.1007/jhep07(2012)044.
  * [32] I. Kostov,  _Classical Limit of the Three-Point Function of $N=4$ Supersymmetric Yang-Mills Theory from Integrability_,  Phys. Rev. Lett. 108, 261604 (2012),  10.1103/PhysRevLett.108.261604.
  * [33] I. Kostov,  _Three-point function of semiclassical states at weak coupling_ ,  J. Phys. A 45(49), 494018 (2012),  10.1088/1751-8113/45/49/494018.
  * [34] E. Bettelheim and I. Kostov,  _Semi-classical analysis of the inner product of Bethe states_ ,  J. Phys. A 47(24), 245401 (2014),  10.1088/1751-8113/47/24/245401.
  * [35] N. Gromov,  _Introduction to the Spectrum of $N=4$ SYM and the Quantum Spectral Curve_ (2017),  1708.03648.
  * [36] O. Babelon, D. Bernard and M. Talon,  _Introduction to Classical Integrable Systems_ ,  Cambridge Monographs on Mathematical Physics. Cambridge University Press,  10.1017/CBO9780511535024 (2003).
  * [37] H. A. Bethe,  _Zur Theorie der Metalle. i. Eigenwerte und Eigenfunktionen der linearen Atomkette_ ,  Zeit. für Physik 71, 205 (1931),  10.1007/BF01341708.
  * [38] V. E. Korepin, N. M. Bogoliubov and A. G. Izergin,  _Quantum Inverse Scattering Method and Correlation Functions_ ,  Cambridge University Press,  10.1017/cbo9780511628832 (1993).
  * [39] M. Gaudin and J.-S. Caux,  _The Bethe Wavefunction_ ,  Cambridge University Press,  10.1017/cbo9781107053885 (2009).
  * [40] R. Baxter,  _Exactly solved models in statistical mechanics_ ,  Elsevier (2016).
  * [41] B. Sutherland,  _Low-Lying Eigenstates of the One-Dimensional Heisenberg Ferromagnet for any Magnetization and Momentum_ ,  Phys. Rev. Lett. 74, 816 (1995),  10.1103/PhysRevLett.74.816.
  * [42] A. Dhar and B. Sriram Shastry,  _Bloch Walls and Macroscopic String States in Bethe’s Solution of the Heisenberg Ferromagnetic Linear Chain_ ,  Phys. Rev. Lett. 85, 2813 (2000),  10.1103/PhysRevLett.85.2813.
  * [43] Y. Miao,  _An exact solution to asymptotic Bethe equation_ ,  J. Phys. A 54(19), 194001 (2020),  10.1088/1751-8121/abf1ad.
  * [44] B. S. Shastry and A. Dhar,  _Solution of a generalized Stieltjes problem_ ,  J. Phys. A 34(31), 6197 (2001),  10.1088/0305-4470/34/31/313.
  * [45] T. Prosen and E. Ilievski,  _Families of quasilocal conservation laws and quantum spin transport_ ,  Phys. Rev. Lett. 111, 057203 (2013),  10.1103/PhysRevLett.111.057203.
  * [46] T. Prosen,  _Quasilocal conservation laws in XXZ spin-1/2 chains: Open, periodic and twisted boundary conditions_ ,  Nucl. Phys. B. 886, 1177 (2014),  10.1016/j.nuclphysb.2014.07.024.
  * [47] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler and T. Prosen,  _Complete generalized gibbs ensembles in an interacting theory_ ,  Phys. Rev. Lett. 115, 157201 (2015),  10.1103/PhysRevLett.115.157201.
  * [48] E. Ilievski, M. Medenjak, T. Prosen and L. Zadnik,  _Quasilocal charges in integrable lattice systems_ ,  J. Stat. Mech. 2016(6), 064008 (2016),  10.1088/1742-5468/2016/06/064008.
  * [49] B. A. Dubrovin,  _Inverse problem for periodic finite zoned potentials in the theory of scattering_ ,  Funct. Anal. its Appl. 9(1), 61 (1975),  10.1007/bf01078183.
  * [50] B. A. Dubrovin, V. B. Matveev and S. P. Novikov,  _Non-linear Equations of Korteweg-De Vries Type, Finite-zone Linear Operators, And Abelian Varieties_ ,  Russ. Math. Surv. 31(1), 59 (1976),  10.1070/rm1976v031n01abeh001446.
  * [51] E. Date, M. Jimbo, M. Kashiwara and T. Miwa,  _Landau-Lifshitz equation: solitons, quasi-periodic solutions and infinite-dimensional Lie algebras_ ,  J. Phys. A 16(2), 221 (1983),  10.1088/0305-4470/16/2/006.
  * [52] R. F. Bikbaev, A. I. Bobenko and A. R. Its,  _Landau-Lifshitz equation, uniaxial anisotropy case: Theory of exact solutions_ ,  Theor. Math. Phys. 178(2), 143 (2014),  10.1007/s11232-014-0135-4.
  * [53] E. K. Sklyanin,  _Separation of Variables: New Trends_ ,  Prog. Theor. Phys. Supp. 118, 35–60 (1995),  10.1143/PTPS.118.35.
  * [54] A. R. Its and V. B. Matveev,  _Schrödinger operators with finite-gap spectrum and N-soliton solutions of the Korteweg-de Vries equation_ ,  Theor. Math. Phys. 23(1), 343 (1975),  10.1007/bf01038218.
  * [55] N. Beisert, A. A. Tseytlin and K. Zarembo,  _Matching quantum strings to quantum spins: One-loop vs. finite-size corrections_ ,  Nucl. Phys. B. 715(1-2), 190 (2005),  10.1016/j.nuclphysb.2005.03.030.
  * [56] N. Beisert, M. Staudacher, J. A. Minahan and K. Zarembo,  _Stringing spins and spinning strings_ ,  J. High Energy Phys. 2003(09), 1 (2003),  10.1088/1126-6708/2003/09/010.
  * [57] M. G. Forest and J.-E. Lee,  _Geometry and Modulation Theory for the Periodic Nonlinear Schrodinger Equation_ ,  In C. Dafermos, J. L. Ericksen, D. Kinderlehrer and M. Slemrod, eds., _Oscillation Theory, Computation, and Methods of Compensated Compactness_ , pp. 35–69. Springer New York, New York, NY,  ISBN 978-1-4613-8689-6 (1986).
  * [58] E. Brézin and V. Kazakov,  _Universality of Correlations of Levels with Discrete Statistics_ ,  Commun. Math. Phys. 214(1), 233 (2000),  10.1007/s002200000259.
  * [59] D. J. Gross and E. Witten,  _Possible third-order phase transition in the large- $N$ lattice gauge theory_,  Phys. Rev. D 21, 446 (1980),  10.1103/PhysRevD.21.446.
  * [60] S. R. Wadia,  _$N=\infty$ phase transition in a class of exactly soluble model lattice gauge theories_,  Phys. Lett. B 93(4), 403 (1980),  10.1016/0370-2693(80)90353-6.
  * [61] M. R. Douglas and V. A. Kazakov,  _Large N phase transition in continuum QCD2_ ,  Phys. Lett. B 319(1-3), 219 (1993),  10.1016/0370-2693(93)90806-s.
  * [62] P. Zinn-Justin,  _Six-vertex model with domain wall boundary conditions and one-matrix model_ ,  Phys. Rev. E 62, 3411 (2000),  10.1103/PhysRevE.62.3411.
  * [63] F. Colomo and A. G. Pronko,  _Third-order phase transition in random tilings_ ,  Phys. Rev. E 88, 042125 (2013),  10.1103/PhysRevE.88.042125.
  * [64] D. Flassig, A. Franca and A. Pritzel,  _Large- $N$ ground state of the Lieb-Liniger model and Yang-Mills theory on a two-sphere_,  Phys. Rev. A 93, 013627 (2016),  10.1103/PhysRevA.93.013627.
  * [65] L. Piroli and P. Calabrese,  _Local correlations in the attractive one-dimensional Bose gas: From Bethe ansatz to the Gross-Pitaevskii equation_ ,  Phys. Rev. A 94, 053620 (2016),  10.1103/PhysRevA.94.053620.
  * [66] R. Dijkgraaf, D. Orlando and S. Reffert,  _Quantum crystals and spin chains_ ,  Nucl. Phys. B. 811(3), 463 (2009),  10.1016/j.nuclphysb.2008.11.027.
  * [67] R. Dijkgraaf, D. Orlando and S. Reffert,  _Relating field theories via stochastic quantization_ ,  Nucl. Phys. B. 824(3), 365 (2010),  10.1016/j.nuclphysb.2009.07.018.
  * [68] R. Hagemans and J.-S. Caux,  _Deformed strings in the Heisenberg model_ ,  J. Phys. A 40(49), 14605 (2007),  10.1088/1751-8113/40/49/001.
  * [69] C. Marboe and D. Volin,  _Fast analytic solver of rational Bethe equations_ ,  J. Phys. A 50(20), 204002 (2017),  10.1088/1751-8121/aa6b88.
  * [70] N. A. Slavnov,  _Calculation of scalar products of wave functions and form factors in the framework of the alcebraic Bethe ansatz_ ,  Theor. Math. Phys. 79(2), 502 (1989),  10.1007/bf01016531.
  * [71] Y. Jiang and J. Brunekreef,  _On the semi-classical limit of scalar products of the XXZ spin chain_ ,  J. High Energy Phys. 2017(3) (2017),  10.1007/jhep03(2017)012.
  * [72] C. Babenko,  _Semiclassical expansion of the Bethe scalar products in the XXZ spin chain_ (2017),  1703.01123.
  * [73] L. D. Faddeev and R. M. Kashaev,  _Quantum Dilogarithm_ ,  Mod. Phys. Lett. A 09(05), 427 (1994),  10.1142/s0217732394000447.
  * [74] R. M. Kashaev and T. Nakanishi,  _Classical and Quantum Dilogarithm Identities_ ,  SIGMA (2011),  10.3842/sigma.2011.102.
  * [75] F. A. Smirnov,  _Form Factors in Completely Integrable Models of Quantum Field Theory_ ,  World Scientific,  10.1142/1115 (1992).
  * [76] M. Jimbo and T. Miwa,  _Algebraic Analysis of Solvable Lattice Models_ ,  American Mathematical Society,  10.1090/cbms/085 (1994).
  * [77] Hagemans, Rob and Caux, Jean-Sébastien and Maillet, Jean-Michel,  _How to Calculate Correlation Functions of Heisenberg Chains_ ,  In _AIP Conference Proceedings_. AIP,  10.1063/1.2222273 (2006).
  * [78] A. Klauser, J. Mossel, J.-S. Caux and J. van den Brink,  _Spin-Exchange Dynamical Structure Factor of the $S=1/2$ Heisenberg Chain_,  Phys. Rev. Lett. 106, 157205 (2011),  10.1103/PhysRevLett.106.157205.
  * [79] F. A. Smirnov,  _Quasiclassical study of form-factors in finite volume_ (1998),  hep-th/9802132.
  * [80] H. Flaschka, M. G. Forest and D. W. McLaughlin,  _Multiphase averaging and the inverse spectral solution of the Korteweg—de Vries equation_ ,  Commun. Pure Appl. Math. 33(6), 739 (1980),  10.1002/cpa.3160330605.
  * [81] P. Calabrese, F. H. L. Essler and M. Fagotti,  _Quantum quench in the transverse-field ising chain_ ,  Phys. Rev. Lett. 106, 227203 (2011),  10.1103/PhysRevLett.106.227203.
  * [82] J. De Nardis, B. Wouters, M. Brockmann and J.-S. Caux,  _Solution for an interaction quench in the Lieb-Liniger Bose gas_ ,  Phys. Rev. A 89, 033601 (2014),  10.1103/PhysRevA.89.033601.
  * [83] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler and T. Prosen,  _Complete generalized gibbs ensembles in an interacting theory_ ,  Phys. Rev. Lett. 115, 157201 (2015),  10.1103/PhysRevLett.115.157201.
  * [84] J.-S. Caux and F. H. L. Essler,  _Time Evolution of Local Observables After Quenching to an Integrable Model_ ,  Phys. Rev. Lett. 110, 257203 (2013),  10.1103/PhysRevLett.110.257203.
  * [85] J.-S. Caux,  _The Quench Action_ ,  J. Stat. Mech. 2016(6), 064006 (2016),  10.1088/1742-5468/2016/06/064006.
  * [86] A. Cavaglià, N. Gromov and F. Levkovich-Maslyuk,  _Separation of variables and scalar products at any rank_ ,  J. High Energy Phys. 2019(9) (2019),  10.1007/jhep09(2019)052.
  * [87] N. Gromov, F. Levkovich-Maslyuk, P. Ryan and D. Volin,  _Dual separated variables and scalar products_ ,  Phys. Lett. B 806, 135494 (2020),  10.1016/j.physletb.2020.135494.
  * [88] R. Hernández, E. López, Á. Periáñez and G. Sierra,  _Finite size effects in ferromagnetic spin chains and quantum corrections to classical strings_ ,  J. High Energy Phys. 2005(06), 011 (2005),  10.1088/1126-6708/2005/06/011.
  * [89] N. I. Akhiezer,  _Elements of the Theory of Elliptic Functions_ ,  American Mathematical Society (1990).
  * [90] H. McKean and V. Moll,  _Elliptic Curves: Function Theory, Geometry, Arithmetic_ ,  Cambridge University Press,  10.1017/CBO9781139174879 (1997).
  * [91] J. Mossel and J.-S. Caux,  _Relaxation dynamics in the gapped XXZ spin-1/2 chain_ ,  New J. Phys. 12(5), 055028 (2010),  10.1088/1367-2630/12/5/055028.


soft open fences

# Schottky’s theorem in $\mathbb{C}^{n}$

P.V.Dovbush

###### Abstract.

The aim of this note is to give a proof of the Schottky theorem in general
domains in $\mathbb{C}^{n}$. The proof is short and works for the cases $n=1$
and $n>1$ at the same time.

###### Key words and phrases:

Marty’s criterion, Schottky’s theorem, normal families, holomorphic functions
of several complex variables

###### 1991 Mathematics Subject Classification:

32A19

Let $\Omega$ be a bounded domain in $\mathbb{C}^{n}$. A family
$\mathcal{F}\subset H(\Omega)$ is normal in $\Omega$ if every sequence of
functions $\\{f_{j}\\}\subseteq F$ contains either a subsequence which
converges to a limit function $f$ uniformly on each compact subset of
$\Omega$, or a subsequence which converges uniformly to $\infty$ on each
compact subset.

The notion of a normal family has played an important role in the development
of complex function theory. We refer to [2] and the references therein for
various results on normal families and applications in case $n=1$. A series of
important theorems surround the Montel criterion for normality, thereby
constituting the so-called Montel cycle. The famous Schottky theorem is one of
them.

Scottky’s Theorem can be given a completely elementary proof, based on nothing
more than Montel’s Theorem and general properties of normal families, which
works for the cases $n=1$ and $n>1$ at the same time.

######  Theorem 0.1.

[Schottky Theorem]. Suppose $f(z)$ is holomorphic in
$\Omega\subset\mathbb{C}^{n}$, $n\geq 1$, which range omits $0$ and $1$ there.
If $a\in\Omega$ and $|f(a)|<M$, then for every $\varepsilon>0$ sufficiently
small we have

(0.1) $|f(z)|<C\textrm{ for all }z\in\Omega,\ \ \delta(z)\geq\varepsilon,$

where $\delta(z)$ denotes the Euclidean distance from $z\in\Omega$ to
$\partial\Omega$. Here $C$ is a positive constant depending only on $f(a)$ and
$\varepsilon$.

The proof of Theorem 0.1 is based on the following lemma.

######  Lemma 0.2.

Let $\mathcal{F}$ be a normal family of holomorphic functions in a domain
$\Omega\subset\mathbb{C}^{n}$, $n\geq 1$. Let $M$ be a positive number and
$a\in\Omega$. Assume that for each function $f\in\mathcal{F}$, we have

$\left|f(a)\right|<M.$

Then the family $\mathcal{F}$ is uniformly bounded on each compact subset of
$\Omega$.

###### Proof.

Let $K$ be a compact set of points belonging to $\Omega$. Assume that
$\mathcal{F}$ is not uniformly bounded on $K$. Then to each positive integer
$k$, corresponds a function $f_{k}\in\mathcal{F}$ such that

$\max_{z\in K}\left|f_{k}(z)\right|>k.$

The sequence $\\{f_{k}\\}$ can not converge on compact $\\{a\\}\cup K$ to a
holomorphic function or to $\infty$ in $\Omega$. This contradicts the
hypothesis that $\mathcal{F}$ is a normal family. ∎

###### Proof of Theorem 0.1.

Let $\mathbf{G}$ denote the family of all functions holomorphic on $\Omega$
which omit the values $0$ and $1$. According to Montel’s theorem [7, Theorem
1.4] the family $\mathbf{G}$ is normal. By Lemma 0.2, the family $\mathbf{G}$
is uniformly bounded on
$\overline{\Omega}_{\varepsilon}=\left\\{z\in\Omega,\,\delta(z)\geq\varepsilon\right\\}$,
hence there is a positive constant $C$ such that (0.1) holds for each function
$f(z)$ of $\mathcal{F}$. This constant $C$ has then the required property. ∎

The method of proof in Theorem 0.1 does not give any information about what
the constant $C$ is.

######  Remark 0.3.

In the case $n=1$ the Shottky theorem has numerous proofs. We refer to [2] for
an exposition (the history, methods and references) of the theorem. Schottky’s
original theorem [1] did not give an explicit bound for $f$. Let $K(f(0),r)$
denote the best possible bound in Theorem 0.1. Various authors have dealt with
the problem of giving an explicit estimate for this bound see Jenkins [3];
Hempel [4] gave some bounds whose constants are in some sense the best
possible.

It has been proved by Murali Rao [8, Theorem 7.17, p. 118] that Shottky’s
Theorem holds with

$C=\tan(\arctan M+(c(1-r^{2}))^{-1}),$

where $c$ is some constant. The proof is based on the work of Minda and
Schober [9]. We shall give an alternative argument based on the classical
Montel Theorem, and Marty’s criterion for normality.

Recall that to every $a$ in the unit disc
$\Delta:=\left\\{z\in\mathbb{C}:\left|z\right|<1\right\\}$ corresponds an
automorphism $\varphi_{a}$ of the disc that interchanges $a$ and 0, namely
$\varphi_{a}(w):=(w-a)/(1-\overline{a}w)$. Note for later reference that
$\varphi_{a}(0)=-a$ and $\varphi_{a}^{\prime}(0)=1-\left|a\right|^{2}$. Define
$Aut(\Delta):=\left\\{\varphi_{a},a\in\Delta\right\\}$.

The following theorem is the main result of this note.

######  Theorem 0.4.

(see [8, Theorem 7.17, p. 118]) Suppose $f(z)$ is holomorphic in the unit disc
$\Delta$ whose range omits $0$ and $1$. If $|f(0)|<M$, then for every
$r\in(0,1)$ we have

(0.2) $|f(z)|<\tan(\arctan M+2rL/(1-r^{2}))\textrm{ for all
}\left|z\right|<r,r\in\interval[openright]{0}{1}.$

###### Proof.

Let $\mathbf{G}$ denote the family of all functions holomorphic on the unit
disc $\Delta$ which omit the values $0$ and $1$. By Montel’s Theorem [5, p.
218] the family $\mathbf{G}\circ Aut(\Delta)$ is normal. Marty’s Theorem [5,
p. 216] yields a non-negative constant $L$ such that for any
$f\circ\varphi_{a}\in\mathbf{G}\circ Aut(\Delta)$

$\frac{\left|f^{\prime}(\varphi_{a}(0))\right|}{1+\left|f(\varphi_{a}(0))\right|^{2}}\leq
L.$

The constant $L$ in the inequality above does not depend on
$f\circ\varphi_{a}$. The exact value $L$ is not known.

Using the chain rule for differentiation and replacing $-a$ by $z$ we obtain

$\frac{\left|f^{\prime}(z)\right|(1-\left|z\right|^{2})}{1+|f(z)|^{2}}\leq
L\textrm{ for all }z\in\Delta.$

Since $f$ never takes the value $0$, the function $t\to|f(tz)|$ is
continuously differentiable on $\interval[open]{0}{1}$ for any fixed
$z\in\Delta[r]=\left\\{\left|z\right|<r\right\\}$ and

$\frac{d}{dt}\Big{(}\arctan(|f(tz)|)\Big{)}\leq\frac{2\left|f^{\prime}(tz)\right|\left|z\right|}{1+|f(tz)|^{2}}\leq\frac{2rL}{1-r^{2}}.$

Integrating this from $t=0$ to $1$, we deduce

$\Big{|}\arctan(|f(z)|)-\arctan(|f(0)|)\Big{|}\leq\int_{0}^{1}\Big{|}\frac{d}{dt}(\arctan(|f(tz)|))\Big{|}dt\leq\frac{2rL}{1-r^{2}}$

and hence

$\arctan(|f(z)|)\leq\arctan(|f(0)|)+\frac{2rL}{1-r^{2}}<\arctan
M+\frac{2rL}{1-r^{2}}$

which is the theorem with $C=\tan(\arctan M+2rL/(1-r^{2}))$. This completes
the proof of Theorem 0.4. ∎

######  Remark 0.5.

The same proof works for $n>1$.

This type of argument can also be used to prove:

######  Theorem 0.6.

[Picard Little Theorem]. At most one complex number is absent from the range
of a non-constant entire function ($=$ holomorphic in $\mathbb{C}$).

###### Proof.

Assume, to the contrary, that an entire function $f$ omits two distinct values
from its range. Composing with a linear fractional transformation, we may also
assume that the omitted values are $0$ and $1$.

The function $f(2^{k}\lambda)$, where $k$ is positive integer, is holomorphic
in $\Delta$, $f(2^{k}\lambda)$ takes the same values in $\Delta$ as does $f$
in $\left|z\right|<2^{k}$, and $f(2^{k}\lambda)$ omits the values $0$ and $1$.

The family $\left\\{f(2^{k}\varphi_{a}),\varphi_{a}\in
Aut(\Delta),k=1,2\ldots\right\\}$, being contained in the family of functions
avoiding $0$ and $1$, is normal in $\Delta$ by Montel’s theorem [5, p. 218].
By Marty’s Theorem [5, p. 216], there exists $L>0$ such that for any
$f(2^{k}\varphi_{a})$

$\frac{\left|f^{\prime}(2^{k}\varphi_{a}(0))\right|}{1+\left|f(2^{k}\varphi_{a}(0))\right|^{2}}\leq
L.$

After some manipulations, using the chain rule for differentiation, the
equalityes $\varphi_{a}(0)=-a$ and
$\varphi^{\prime}_{a}(0)=1-\left|a\right|^{2}$, we find

$\frac{2^{k}\left|f^{\prime}(-2^{k}a)\right|(1-\left|a\right|^{2})}{1+\left|f(-2^{k}a)\right|^{2}}\leq
L.$

Replacing in this inequality $-2^{k}a$ by $z$ we get

$\frac{2^{k}\left|f^{\prime}(z)\right|(1-\left|z/2^{k}\right|^{2})}{1+\left|f(z)\right|^{2}}\leq
L.$

Thus

$\frac{\left|f^{\prime}(z)\right|}{1+|f(z)|^{2}}\leq\frac{2^{k}L}{2^{2k}-\left|z\right|^{2}}.$

Keeping $z$ fixed, let $k$ tend to infinity to deduce that $f^{\prime}(z)=0$
for every point $z$ in the complex plane. But this clearly implies that $f(z)$
is a constant, not the non-constant entire function. The obtained
contradiction proves the theorem. ∎

Immediate (and easy) corollary of Picar’s Little Theorem is Liouville’s
Theorem:

######  Corollary 0.7.

Every bounded entire function is constant.

Theorem 0.6 is a remarkable generalization of Liouville’s theorem. It is
simple to find entire functions whose range is the entire $\mathbb{C}$;
nonconstant polynomials, for instance. The exponential function is an example
of an entire function whose range omits only one value, namely zero. But there
does not exist a nonconstant entire function whose range omits two different
values.

The fundamental theorem of algebra is a consequence of Liouville’s Theorem
(for the proof see, for example, [8, Corollary 4, p. 98]).

An alternate program for obtaining the above theorems (and not only) in
elementary fashion was indicated by Zalcman [6, p. 817].

######  Remark 0.8.

Many of the theorems on the theory of normal families for holomorphic
functions of one variable may be transferred without essential changes to the
case of functions of several complex variables.

Slice functions is a useful tool that allow us to apply facts from the
function theory of unit disc $\Delta$ to questions in the unit ball
$B\subset\mathbb{C}^{n}$. Suppose $f(z)$ is holomorphic in the unit ball whose
range omits $0$ and $1$. Fix a point $z\in B$ with nonzero norm
$r=\left|z\right|$ and consider the holomorphic function of one complex
variable that sends a number $\lambda$ in the unit disk to the image
$f(\lambda z/r)$. This slice function omits the values $0$ and $1$, takes the
value $f(0)$ when $\lambda=0$, and takes the value $f(z)$ when $\lambda=r$.
Then $\left|f(z)\right|<C$, where $C$ is the constant in Theorem 0.4
corresponding to $f(0)$ and $r$. Thus Schottky’s theorem holds in the unit
ball with the same constant as in dimension one.

## References

  * [1] F. Schottky. “On Picard’s theorem and Borel’s inequality”. German. In: Berl. Ber. 1904 (1904), pp. 1244–1263.
  * [2] Robert B. Burckel. Classical analysis in the complex plane. Cornerstones. Springer, New York, $[2021]$ , pp. xxix+1123. isbn: 978-1-0716-1963-6; 978-1-0716-1965-0. doi: 10.1007/978-1-0716-1965-0.
  * [3] J. A. Jenkins. “On explicit bounds in Schottky’s theorem”. In: Canadian J. Math. 7 (1955), pp. 76–82. issn: 0008–414X. doi: 10.4153/CJM-1955-010-4.
  * [4] Joachim A. Hempel. “Precise bounds in the theorems of Schottky and Picard”. In: J. London Math. Soc. (2) 21.2 (1980), pp. 279–286. issn: 0024–6107. doi: 10.1112/jlms/s2-21.2.279.
  * [5] Lawrence Zalcman. “Normal families: new perspectives”. In: Bull. Amer. Math. Soc. (N.S.) 35.3 (1998), pp. 215–230. issn: 0273-0979. doi: 10.1090/S0273-0979-98-00755-1.
  * [6] Lawrence Zalcman. “A heuristic principle in complex function theory”. In: Amer. Math. Monthly 82.8 (1975), pp. 813–817. issn: 0002-9890. doi: 10.2307/2319796.
  * [7] P. V. Dovbush. “On normal families in $C^{n}$”. In: Complex Var. Elliptic Equ. 67.1 (2022), pp. 1–8. issn: 1747–6933. doi: 10.1080/17476933.2020.1797703.
  * [8] Murali Rao et al. Complex analysis. (2nd Edition). World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2015, pp. x+414. isbn: 978–981–4579–599. doi: 10.1142/9062.
  * [9] David Minda and Glenn Schober. “Another elementary approach to the theorems of Landau, Montel, Picard and Schottky”. In: Complex Variables Theory Appl. 2.2 (1983), pp. 157–164. issn: 0278–1077. doi: 10.1080/17476938308814039.


#  Universal emergence of local Zipf’s law

Davide Cugini  Dipartimento di Fisica, Università di Pavia, via Bassi 6,
27100, Pavia, Italy  André Timpanaro  Universidade Federal do ABC, 09210-580
Santo André, Brazil  Giacomo Livan  Dipartimento di Fisica, Università di
Pavia, via Bassi 6, 27100, Pavia, Italy  Department of Computer Science,
University College London, 66-72 Gower Street, WC1A 6EA London, UK  Giacomo
Guarnieri  Dipartimento di Fisica, Università di Pavia, via Bassi 6, 27100,
Pavia, Italy

###### Abstract

A plethora of natural and socio-economic phenomena share a striking
statistical regularity, that is the magnitude of elements decreases with a
power law as a function of their position in a ranking of magnitude. Such
regularity is commonly known as Zipf’s law, and plenty of problem-specific
explanations for its emergence have been provided in different fields. Yet, an
explanation for Zipf’s law ubiquity is currently lacking. In this paper, we
demonstrate from first principles that Zipf’s behavior naturally emerges as a
local approximation to the order statistics generated by any ranking process.
We validate our results against several relevant examples.

Figure 1: Panel (a) represents a parent distribution used to extract $N=30$
i.i.d. samples, depicted as dots. The same samples are ordered and reported in
panel (b), together with the order statistics expectation value depicted as a
dashed line. The filled area in panel (b) represents the region within the
expected standard deviation from the dashed line.

Introduction. –  Many natural and social phenomena, once ranked according to a
certain feature $x$, e.g. their magnitude or frequency, are found to follow an
empirical power law $x(r)\sim A/(B+r)^{\alpha}$ (with $r$ being the rank).
This is known as Zipf’s law (or, equivalently, Pareto distribution) [1, 2, 3],
and is observed rather ubiquitously, e.g. in magnitudes of earthquakes [4],
cosmic structures [5], cities populations [6], financial returns [7], words
usage in a language [8], choices of chess openings [9], visualizations of web
pages [10], neural activity [11, 12], the number of citations of scientific
papers [13] and genes expression [14], to cite but a few.
This universality has spurred an intense research aimed at understanding the
underlying mechanisms that could lead to the emergence of Zipf’s law in the
organization of such extremely diverse complex phenomena. Apart from the
above-mentioned scenarios where specific models were formulated to explain the
observed occurrence on a case-by-case basis, there have recently been efforts
to derive Zipf’s law from less fine-tuned conditions, including hidden
variable models [15], stochastic processes [16], dynamical constraint [17],
least-effort principle [18], and criticality in multivariate data [19, 20].
Despite all these advancements, a single general derivation that would account
for the ubiquitous emergence of Zipf’s law is still missing. This is precisely
what we provide in this work. Our result is achieved by making use of notions
from the order statistics, a field of fundamental importance in understanding
the distributional properties of sample data [21, 22, 23].

First of all, we derive a new general expression for all the cumulants of the
order statistics obtained from data drawn from a given probability
distribution and then, building on this result, we demonstrate that any such
order statistics is always _at least locally_ described by a Zipf’s law. The
implication of our main result is then a rigorous proof of the observed
ubiquity of Zipf’s law emergence in Nature. We finally showcase our findings
on three paradigmatic examples, two of which representing important benchmarks
in the field (i.e. Miller’s typing monkey [24] and the Barabási-Albert model
[25]) and finally arguably the most ubiquitous statistical model, i.e. the one
generated by a Gaussian distribution, showing also in this case how the
associated order statistics locally gives rise to a Zipf’s law.

Order Duality Relationships. – We start by considering $N$ independently and
identically distributed (i.i.d.) draws from a probability distribution $p(x)$,
henceforth called parent p.d.f., and rank them in non-increasing order, i.e.,
$X_{r}\geq X_{r+1}$, $\forall\ r<N$. The resulting sorted set is usually
referred to as the ‘order statistics’ [23]. At this point, the goal can be
properly framed as a so-called ’backward problem’ [26], namely the
identification of the parent p.d.f. giving rise to a given Zipf order
statistics. A first crucial observation to make is that the sample size $N$
can be very different and considering a given $X_{r}$ can have a completely
different meaning depending on such size. For example, given a sample size
$N=9$, $X_{5}$ would correspond to the median value, while for a larger sample
size $N=100$ it would belong to the top $5\%$ elements, thus far above the
median. In order to circumvent this, we introduce a new variable $\lambda=r/N$
which represents the relative ranking, and denote $X(\lambda)\equiv
X_{r=N\lambda}$. In the limit $N\to\infty$, $\lambda$ becomes a continuous
variable in the interval $\left(0,1\right]$. Our first main result is a closed
analytic correspondence between all the statistical cumulants of $X(\lambda)$
and the cumulative distribution associated with the parent p.d.f.
$F(x)=\int_{-\infty}^{x}dx^{\prime}p(x^{\prime})$; see Supplementary Material,
Section A. We stress that $p(x)$ can be any analytic p.d.f. . These results,
which incorporate the ‘backward problem’, are called ‘Order duality
relationships’ (ODRs). Here in the main text we explicitly report the
expression for the first two cumulants, i.e. average and variance, which read

$\displaystyle\langle X(\lambda)\rangle$
$\displaystyle=F^{-1}\left[1-\frac{N}{N+1}\lambda\right]+\mathcal{O}\left[1/N\right],$
(1) $\displaystyle\mathrm{Var}\left[X(\lambda)\right]$
$\displaystyle=\frac{N\lambda\left(N[1-\lambda]+1\right)}{p^{2}\left(\langle
X(\lambda)\rangle\right)(N+2)(N+1)^{2}}+\mathcal{O}\left[1/N^{2}\right].$

It can be noted that, for $N\to+\infty$, the average value of the order
statistics tends to a constant, while the variance decreases as $\sim 1/N$.

An immediate consequence of the ODRs is that, for any parent p.d.f. of the
form $p(x)\propto x^{-\beta}$ ($\beta\neq 1$), the resulting order statistics
follows a Zipf’s law with $\alpha=(\beta-1)^{-1}$, in agreement with the
result obtained in [26]. It is worth remarking, however, that at the
exceptional value $\beta=1$, the order statistics manifests an exponential
falloff. These considerations already provide new insights on several
benchmarking models, e.g. the Barabási-Albert algorithm [25], as we will
detail in the following Section.

It is worth mentioning that the generality of ODRs may seem in contrast with
the occasional empirically observed deviations from Zip’s law at
$\lambda\approx 0$ [26]. However, Eq. (1) do not only agree but provide a
recovery of Zipf’s law also in this regime. The key observation is that the
variance of $X(\lambda\approx 0)$ in these instances can diverge as
$\lambda/p^{2}(\langle X(\lambda)\rangle)$ for $\lambda\to 0$. Therefore such
deviations become compatible with the ideal expected Zipf’s behavior within
the predicted statistical error.

Local Zipf emergence. – Notwithstanding, Eq. (1) alone still do not account
for the apparent ubiquity of Zipf’s law. To this aim, a first analysis of
empirical evidence suggests that such behavior is expected to emerge locally
when the very top fraction of relative ranking events are considered, which
correspond to $\lambda\approx 0$. For example, in the case of earthquakes one
typically focuses on the detection of very intense magnitude events while
neglecting the huge amount of tiny ground vibrations (which could hardly be
systematically detected anyway). Another example would be to consider only the
top $100$ players in a given sport, e.g. soccer or chess.

Our main result is not only to rigorously prove this statement but also that
the same conclusion applies in the vicinity of any generic $\lambda_{0}$, i.e.
for $\lambda\approx\lambda_{0}$. In particular, we prove that, for any fixed
$\lambda_{0}$, it is possible to find a value of
$\alpha\equiv\alpha(\lambda_{0})$ such that

$\langle
X(\lambda)\rangle=\frac{A}{\left(B+\lambda\right)^{\alpha}}+\mathcal{O}\left[(\lambda-\lambda_{0})^{3}\right]\,,$
(2)

where $A$ and $B$ are constants (we provide their exact but rather cumbersome
expressions in the Supplemental Material, Section B). This result demonstrates
that the order statistics of any analytic parent p.d.f. is locally
indistinguishable, up to third order in the relative ranking, from a Zipf’s
distribution (see Supplemental Material, Section B, for details).

Applications. –  In the following, we showcase our results on three
paradigmatic examples, chosen in light of their interdisciplinarity and wide
applicability in several contexts. The first two, in particular, highlight the
close relationship between the tail exponents of a power-law parent p.d.f. and
the corresponding Zipf’s law. The third example, instead, shows the
aforementioned emergence of local Zipf behaviour for data drawn from an
associated Gaussian parent p.d.f. .

1\. Miller’s Typing monkey. We begin by benchmarking our formalism to re-
derive one of the classical results associated with Zipf’s law, known as
Miller’s typing monkey [24]. Let $m$ denote the number of letters of a given
alphabet and $q_{s}$ be the probability of hitting the space-bar, thus
terminating a word. Each letter is assumed to have the same probability
$(1-q_{s})/m$ to be typed and we remark that any possible string of letters
will be considered to be a valid word. In this gedankenexperiment, one
imagines to have a ’dictionary’, consisting of all the existing words, a
number of which are typed by the monkey. One then records the frequency
(relative number of times) each word is typed. At the end, one sorts the typed
words from the most frequently typed to the least one, i.e. build the order
statistics. Our goal now is to apply our results in order to predict the
latter from Eq. (1). In this thought experiment, the ranking feature $x$ is
the words’ frequency and thus the parent p.d.f. $p(x)$ expresses the
distribution of the frequency values that have been recorded. We can obtain an
expression for $p(x)$ by observing that the frequency $x$ of a typed word is
strictly connected to the word length $y$ by

$x(y)=q_{s}[1-q_{s}]^{y}m^{-y}=q_{s}\mathrm{exp}\left[-\mathrm{ln}\left(\frac{m}{1-q_{s}}\right)y\right]\
.$ (3)

Since the probability for a word in the dictionary to have exactly $y$ letters
is $P(y)\propto m^{y}$, a straightforward change of variable leads to the
desired result

$p(x)=P(y)\frac{dy}{dx}\propto x^{-(1+b/a)}\ ,$ (4)

where $a\equiv\mathrm{ln}\left[m(1-q_{s})^{-1}\right]$ and
$b\equiv\mathrm{ln}\,m$. Since the parent distribution is a power-law
$p(x)\propto x^{-\beta}$ with exponent

$\beta=1+\frac{b}{a},$ (5)

the theory we developed guarantees that the associated expectation value for
the order statistics also approaches a Zipf’s law $\langle
X(\lambda)\rangle\propto\lambda^{-\alpha}$ with

$\alpha=(\beta-1)^{-1}=\frac{a}{b}.$ (6)

This result agrees with and extends the result originally obtained by Miller
[24], which only described the average behavior of groups of words with equal
length.

Figure 2: Ranking of nodes of a network generated through Barabási-Albert
algorithm, where nodes are ranked basing on their number of edges $X$. In the
plot $X$ has been normalized with the total number of edges in the network.
Numerical results are reported for $N=\\{10,20,50,100\\}$. The coloured shaded
regions represent the associated standard deviations, obtained from Eq. (1),
while the dotted line represents the theoretical prediction in the
$N\to\infty$ limit.

2\. Barabási-Albert model Several real world networks turn out to be scale-
free [27, 28]. These include, e.g., the internet [29], cellular structures
[30] and scientific coauthorship networks [31], among many others. This is
equivalent to the statement that a randomly picked node of the network has a
probability

$p(x)\propto x^{-\beta},\quad\beta>0$ (7)

to have a number of edges $x$. The Barabási-Albert (BA) model provides an
algorithm, based on the concepts of growth and preferential attachment [25,
28], that allows to construct a network with this feature. It is known that,
independently from initial conditions, BA networks exhibit $\beta=3$ as the
number of nodes grows very large. More complex models have been proposed in
order to produce networks with $\beta\neq 3$ [32]. We can now identify $p(x)$
as the parent distribution that assigns to each node its number of edges.
Using the ODRs Eq. (1), we can predict the behaviour $\langle
X(\lambda)\rangle\propto\lambda^{-\frac{1}{\beta-1}}$ for the corresponding
order statistics. In particular, for the BA model we obtain $\langle
X(\lambda)\rangle\propto\lambda^{-1/2}$. We test such a result numerically by
using the BA algorithm to independently generate $100$ networks with a total
number of nodes $N$, which are then sorted based on their number of edges.
This allows us to calculate the mean and the standard deviation of the order
statistics. We report the results for $N=\\{10,20,50,100\\}$ in Fig. 2, where
we correctly observe the convergence to the analytical result obtained from
Eq. 1 in the $N\to\infty$ limit, with a standard deviation that decays as
$\sim 1/\sqrt{N}$.

3\. Normal distribution. We now turn our attention to showcasing the notion of
local Zipf behaviour. Given its ubiquity — due to the Central Limit Theorem —
the order statistics generated by a Gaussian distribution is of paramount
importance, so in the following we shall focus on a Gaussian parent p.d.f.:

$p(x)=\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2}\ .$ (8)

Using the ODRs together with Eq. (2) one obtains that the order statistics
produced by ordering $N$ i.i.d. random variables drawn from $p(x)$ manifests a
local Zipf behaviour almost everywhere. The only exception is in fact the
local behaviour around $\langle X(\lambda)\rangle=\pm 1$, i.e., at one
standard deviation from the mean, that is exponential (see Supplemental
Material, Section C).

In Fig. 3 we report the order statistics’ expectation value for the top $100$
draws in samples of $N=10^{6}$ draws from the above Gaussian p.d.f. (solid
line), i.e., the order statistics in the top $0.01\%$ of the ranking. Together
with that, we show the corresponding order statistics obtained from Eq. (2)
with two different values of $\lambda_{0}$, namely $\lambda_{0}^{L}=10^{-5}$
and $\lambda_{0}^{R}=9\cdot 10^{-5}$ (highlighted by vertical lines). As it
can be seen, in both cases the empirical order statistics is very well
approximated by the local Zipf behaviour predicted by Eq. (1) and Eq. (2). Let
us remark that ‘zooming’ on the upper percentiles of a ranking — just like we
did in this example — corresponds to most practical situations where one is
dealing with very large samples, e.g., ranking the $\sim$$10^{9}$ people on
the planet based on their wealth.

Figure 3: Order statistics associated with a Gaussian parent distribution.
Numerical results are reported together with shades representing their
statistical error. The plot represents the top $10^{2}$ ranked draws out of
$N=10^{6}$ total samples. The theoretically predicted local power-law
behaviour around $\lambda_{0}^{L}=10^{-5}$ and $\lambda_{0}^{R}=9\cdot
10^{-5}$ is reported as dashed lines together with the statistical
fluctuations as obtained from Eq. (1).

Conclusions. –  In our work we have obtained a novel analytical correspondence
between a generic probability distribution and its associated order
statistics. We then made use of this result, which we called Order Duality
Relationships (ODRs), to rigorously show that any analytic distribution gives
rise to a Zipf’s law at least locally, with an error that is suppressed as
$N^{-3}$, with $N$ being the sample size. This provides a demonstration of the
universality of Zipf’s behavior in Nature, which we then proceeded to showcase
in three paradigmatic examples.

The first model we considered is the Miller’s typing monkey problem, that
allowed us to compare our results with another analytical benchmark. The
second is the Barabási-Albert model, for which we managed to predict the order
statistics of the numerically generated samples. Finally, we moved our
attention to the case of Gaussian parent p.d.f., that is of extreme interest
for most real-word applications. In this case we further verified that the
local power-like approximation of the parent p.d.f. leads to a satisfactory
description of the order statistics, whenever a small fraction of the full
ranking is taken into account.

## References

  * Newman [2005] M. E. Newman, Contemporary physics 46, 323 (2005).
  * Saichev  _et al._ [2009] A. I. Saichev, Y. Malevergne, and D. Sornette,  _Theory of Zipf’s law and beyond_ , Vol. 632 (Springer Science & Business Media, 2009).
  * Clauset  _et al._ [2009] A. Clauset, C. R. Shalizi, and M. E. Newman, SIAM review 51, 661 (2009).
  * Sornette  _et al._ [1996] D. Sornette, L. Knopoff, Y. Kagan, and C. Vanneste, Journal of Geophysical Research: Solid Earth 101, 13883 (1996).
  * De Marzo  _et al._ [2021a] G. De Marzo, F. S. Labini, and L. Pietronero, Astronomy & Astrophysics 651, A114 (2021a).
  * Gabaix [1999] X. Gabaix, American Economic Review 89, 129 (1999).
  * Gabaix [2009] X. Gabaix, Annu. Rev. Econ. 1, 255 (2009).
  * Piantadosi [2014] S. T. Piantadosi, Psychonomic bulletin & review 21, 1112 (2014).
  * Blasius and Tönjes [2009] B. Blasius and R. Tönjes, Physical Review Letters 103 (2009), 10.1103/physrevlett.103.218701.
  * Adamic and Huberman [1999] L. Adamic and B. A. Huberman, (1999).
  * Mora and Bialek [2011] T. Mora and W. Bialek, Journal of Statistical Physics 144, 268 (2011).
  * Tyrcha  _et al._ [2013] J. Tyrcha, Y. Roudi, M. Marsili, and J. Hertz, Journal of Statistical Mechanics: Theory and Experiment 2013, P03005 (2013).
  * Price [1965] D. J. D. S. Price, Science 149, 510 (1965).
  * Furusawa and Kaneko [2003] C. Furusawa and K. Kaneko, Physical review letters 90, 088102 (2003).
  * Aitchison  _et al._ [2016] L. Aitchison, N. Corradi, and P. Latham, PLoS Computational Biology 12 (2016), 10.1371/journal.pcbi.1005110.
  * James  _et al._ [2018] C. James, S. Azaele, A. Maritan, and F. Simini, Physical Review E (2018), 10.1103/PhysRevE.98.032408.
  * De Marzo  _et al._ [2021b] G. De Marzo, A. Gabrielli, A. Zaccaria, and L. Pietronero, Physical Review Research 3, 013084 (2021b).
  * Wang [2020] Q. Wang, Chaos, Solitons and Fractals (2020), 10.1016/j.chaos.2021.111489.
  * Schwab  _et al._ [2014] D. J. Schwab, I. Nemenman, and P. Mehta, Physical Review Letters 113 (2014), 10.1103/physrevlett.113.068102.
  * Marsili and Roudi [2022] M. Marsili and Y. Roudi, Physics Reports 963, 1 (2022).
  * Rényi [1953] A. Rényi, Acta Math. Acad. Sci. Hung 4, 48 (1953).
  * Arnold  _et al._ [2008] B. C. Arnold, N. Balakrishnan, and H. N. Nagaraja,  _A First Course in Order Statistics_ (Society for Industrial and Applied Mathematics, 2008).
  * David and Nagaraja [2004] H. A. David and H. N. Nagaraja,  _Order statistics_ (John Wiley & Sons, 2004).
  * Miller [1957] G. A. Miller, The American journal of psychology 70, 311 (1957).
  * Barabási and Albert [1999] A.-L. Barabási and R. Albert, science 286, 509 (1999).
  * Cristelli  _et al._ [2012] M. Cristelli, M. Batty, and L. Pietronero, Scientific reports 2, 812 (2012).
  * Caldarelli [2007] G. Caldarelli,  _Scale-free networks: complex webs in nature and technology_ (Oxford University Press, 2007).
  * Albert and Barabási [2002] R. Albert and A.-L. Barabási, Reviews of modern physics 74, 47 (2002).
  * Pastor-Satorras and Vespignani [2001] R. Pastor-Satorras and A. Vespignani, Physical review letters 86, 3200 (2001).
  * Albert [2005] R. Albert, Journal of cell science 118, 4947 (2005).
  * Newman [2001] M. E. Newman, Proceedings of the national academy of sciences 98, 404 (2001).
  * Krapivsky  _et al._ [2000] P. L. Krapivsky, S. Redner, and F. Leyvraz, Phys. Rev. Lett. 85, 4629 (2000).

Supplemental Material

In this Supplemental Material we provide the proofs and all additional details
necessary to obtain the results presented in the main text. In particular,
Section A is devoted to the analytical derivation od the Order Duality
Relationships, viz. Eq. (1) of the main text; in Section B we demonstrate the
local indistinguishability of a generic order statistics from a Zipf’s law,
viz. Eq. (2) of the main text; Section C finally contains additional details
and considerations about the Zipf’s behavior associated to a Gaussian parent
p.d.f. .

##  Appendix A Section A: Order Duality Relations

In this section we prove the Order Duality Relations in Eq. (1). At this
scope, we first focus on the simple scenario of a uniformly distributed parent
p.d.f., whose results are then extended to more generic parent distributions.
Let $g(u)$ be the uniform distribution defined by

$g(u)=\begin{cases}1,&\mathrm{if}\;u\in[0,1]\\\ 0,&\mathrm{otherwise}\\\
\end{cases}$ (S1)

and let $\\{u_{r}\\}_{r\in\\{1,2,...,N\\}}$ be a set of $N$ independent random
variables distributed according to $g(u)$. We denote with $\\{U_{r}\\}$ the
order statistic obtained by sorting the $\\{u_{r}\\}$ in non-increasing order.
The p.d.f. for the r-th order statistic $U_{r}$ is

$g_{r}(u)=\frac{N!}{(N-r)!(r-1)!}[1-u]^{r-1}u^{N-r}\,.$ (S2)

The expectation value for the k-th moment can be obtained by straightforward
calculation and is given by

$\langle
U_{r}^{k}\rangle=\frac{N!}{(N+k)!}\frac{(N+k-r)!}{(N-r)!}\quad\quad\forall
k\in\mathds{N}.$ (S3)

We hereby introduce the change of variable

$\lambda\equiv\frac{r}{N}\in\left\\{\frac{1}{N},\frac{2}{N},..1\right\\}$ (S4)

and use the notation $U(\lambda)\equiv U_{\lambda N}$, in terms of which

$\langle
U^{k}(\lambda)\rangle=\frac{N!}{(N+k)!}\frac{(N[1-\lambda]+k)!}{(N[1-\lambda])!}\,.$
(S5)

This change of variable allows us to study the $N\to\infty$ limit while
keeping the ratio $\lambda=r/N$ fixed. This has as a direct consequences

$\langle U(\lambda)\rangle=1-\frac{N}{N+1}\lambda\,,$ (S6)
$\left\langle\left(U(\lambda)-\left\langle
U(\lambda)\right\rangle\right)^{2}\right\rangle=\frac{\lambda
N(N[1-\lambda]+1)}{(N+1)^{2}(N+2)}$ (S7)

and

$\left\langle\left(U(\lambda)-\left\langle
U(\lambda)\right\rangle\right)^{p}\right\rangle=\mathrm{o}\left(\frac{1}{N}\right)\quad\forall
p\in\mathbb{N},\,p>2,$ (S8)

that will be extensively used in what follows. In particular eq. (S7) implies
that the variance of $U(\lambda)$ is suppressed as $1/N$ in the limit of
$N\to\infty$. We now move our attention to the more generic framework of
independent random variables $\\{x_{r}\\}_{r\in\\{1,...,N\\}}$ whose
associated analytic p.d.f. is $p(x)$. Let
$F(x)=\int_{-\infty}^{x}dx^{\prime}\,p(x^{\prime})$ be the correspondent
cumulative density function (c.d.f.) that can be inverted since $p(x)$ is
analytic. Then, the random variables $u_{r}\equiv F(x_{r})$ are uniformly
distributed between 0 and 1 and therefore the results obtained above hold. The
independent variables

$X_{r}\equiv F^{-1}\left(U_{r}\right)$ (S9)

obtained from the order statistics $\\{U_{r}\\}$ are themselves order
statistics, i.e.

$X_{r}\geq X_{r+1}.$ (S10)

Using the same notation introduced for the uniform parent scenario, we write
$X(\lambda)=X_{\lambda N}$. Starting from eq. (S9) and Taylor-expanding around
$\langle U(\lambda)\rangle$ we obtain

$\langle
X^{k}(\lambda)\rangle=\left\langle\left[F^{-1}\left(U(\lambda)\right)\right]^{k}\right\rangle=\left\langle\left(\Sigma_{j=0}^{\infty}\frac{1}{j!}\frac{d^{j}F^{-1}}{dU^{j}}\Big{|}_{\langle
U(\lambda)\rangle}\cdot\left(U(\lambda)-\langle
U(\lambda)\rangle\right)^{j}\right)^{k}\right\rangle$ (S11)

We can know use eq. (S8) to obtain

$\boxed{\langle X(\lambda)\rangle=F^{-1}\left[\langle
U(\lambda)\rangle\right]+\mathrm{o}\left(1\right)}\,,$ (S12)
$\boxed{\left\langle\left(X(\lambda)-\left\langle
X(\lambda)\right\rangle\right)^{2}\right\rangle=\frac{\left\langle\left(U(\lambda)-\left\langle
U(\lambda)\right\rangle\right)^{2}\right\rangle}{p^{2}\left(\left\langle
X(\lambda)\right\rangle\right)}+\mathrm{o}\left(\frac{1}{N}\right)}$ (S13)

that, together with eq. (S6) and eq. (S7), give the Order Duality Relations.

##  Appendix B Section B: Local Zipf behaviour

We hereby start from Eq. (S12) to show that any order statistics can be
locally described by a Zipf Law. In particular, once fixed
$\lambda_{0}\in[0,1]$, one gets

$\frac{N}{N+1}\left(\lambda_{0}-\lambda\right)=\int_{\langle
X(\lambda_{0})\rangle}^{\langle X(\lambda)\rangle}p(x)dx.$ (S14)

In what follows we assume that both $\langle X(\lambda)\rangle$ and $\langle
X(\lambda_{0})\rangle$ are positive, analogous results can be obtained for
negative values after the redefinition $x\mapsto-x$. We consider the change of
variables

$\begin{cases}x=e^{\eta}\\\ p(x)=e^{z(\eta)}\end{cases}$ (S15)

From the the assumption of $p(x)$ to be analytic it follows that also
$z(\eta)$ is analytic. As a consequence $z(\eta)$ is smooth and its Taylor
expansion around any point converges in an interval. Notice that the parent
distribution $p(x)$ is a power-law iff the derivatives $z^{(i)}(H_{0})$ vanish
$\forall i\geq 2$. Indeed, let $p(x)\propto x^{-\beta}$, then

$\frac{p(x)}{p(x_{1})}=\left(\frac{x}{x_{1}}\right)^{-\beta}$ (S16)

After the change of variables defined in Eq. (S15), the previous equation
becomes

$z(\eta)=z(\eta_{1})-\beta(\eta-\eta_{1}),\quad\quad x_{1}\equiv
e^{\eta_{1}},$ (S17)

that is a linear behaviour and therefore $z^{(i)}(H_{0})=0,\;\forall i\geq 2$.
We now define $H(\lambda)=\mathrm{ln}(\langle X(\lambda)\rangle)$ and
$H_{0}=\mathrm{ln}(\langle X(\lambda_{0})\rangle)$, obtaining

$\begin{split}\frac{N}{N+1}\left(\lambda_{0}-\lambda\right)&=\int_{H_{0}}^{H(\lambda)}d\eta\,\mathrm{exp}\left[z(\eta)+\eta\right]=\\\
&=\mathrm{exp}\left[z(H_{0})+H_{0}\right]\int_{H_{0}}^{H(\lambda)}d\eta\,\mathrm{exp}\left[(z^{(1)}(H_{0})+1)(\eta-
H_{0})+\sum_{i=2}^{\infty}\frac{z^{(i)}(H_{0})}{i!}\,[\eta-
H_{0}]^{i}\right].\end{split}$ (S18)

The idea is now to Taylor-expand both sides around $\lambda=\lambda_{0}$ and
compare them term by term. From the first order we obtain

$\frac{dH}{d\lambda}\Big{|}_{\lambda_{0}}=-\frac{N}{N+1}\mathrm{exp}\left[-\left(z(H_{0})+H_{0}\right)\right]\,.\,$
(S19)

Moving to the second order we get

$\frac{d^{2}H}{d\lambda^{2}}\Big{|}_{\lambda_{0}}=-(z^{(1)}(H_{0})+1)\left(\frac{N}{N+1}\right)^{2}\mathrm{exp}\left[-2\left(z(H_{0})+H_{0}\right)\right]\,.\,$
(S20)

Finally, we repeat the calculation for the third order, obtaining

$\frac{d^{3}H}{d\lambda^{3}}\Big{|}_{\lambda_{0}}=-\left[2(z^{(1)}(H_{0})+1)^{2}-z^{(2)}(H_{0})\right]\left(\frac{N}{N+1}\right)^{3}\mathrm{exp}\left[-3\left(z(H_{0})+H_{0}\right)\right]\,.\,$
(S21)

Observe that if we approximate $p(x)$ with a power law around $\langle
X(\lambda_{0})\rangle$, i.e. we set $z^{(i)}(H_{0})=0$ for all $i\geq 2$, the
result obtained for $H(\lambda)$ from Eq. (S18) doesn’t change, up to third
order corrections. Since $\langle
X(\lambda)\rangle=\mathrm{exp}\left[H(\lambda)\right]$, the same property
holds for $\langle X(\lambda)\rangle$. Therefore, one can approximate the
p.d.f. as

$\begin{split}p(x)&\approx\mathrm{exp}\left[z\left(H_{0}\right)+z^{(1)}\left(H_{0}\right)\left(\eta-
H_{0}\right)\right]\\\ &=p\left(\langle
X(\lambda_{0})\rangle\right)\left[\frac{x}{\langle
X(\lambda_{0})\rangle}\right]^{z^{(1)}\left(H_{0}\right)}\end{split}$ (S22)

and use it to solve Eq. (S18), obtaining an expression for $\langle
X(\lambda)\rangle$ that is known to be correct up to third order terms.
Explicitly, one obtains

$\langle X(\lambda)\rangle=\left[\langle
X(\lambda_{0})\rangle^{z^{(1)}(H_{0})+1}-(z^{(1)}(H_{0})+1)(\lambda-\lambda_{0})\right]^{1/(z^{(1)}(H_{0})+1)}+\mathcal{O}\left[\left(\lambda-\lambda_{0}\right)^{3}\right]$
(S23)

that can be cast as in Eq. (2) $\langle
X(\lambda)\rangle=/\left(B+\lambda\right)^{\alpha}+\mathcal{O}\left[(\lambda-\lambda_{0})^{3}\right]\,$
upon the following identifications:

$\displaystyle\alpha$
$\displaystyle\equiv\alpha(\lambda_{0})=-1/(z^{(1)}(H_{0})+1)\,,$ (S24)
$\displaystyle A$ $\displaystyle\equiv\alpha^{\alpha}\,,$ (S25) $\displaystyle
B$ $\displaystyle\equiv\alpha\langle
X(\lambda_{0})\rangle^{-1/\alpha}-\lambda_{0}\,.$ (S26)

##  Appendix C Section C: Gaussian parent p.d.f.

Here we focus on the order statistics associated to a Gaussian parent p.d.f.

$p(x)=\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2}.$ (S27)

In order to perform the change of variables in Eq. (S15) we restrict our study
to either $x\geq 0$ or $x<0$. In any case we define

$\begin{cases}\absolutevalue{x}=e^{\eta}\\\ p(x)=e^{z}(\eta),\end{cases}$
(S28)

from where

$z(\eta)=-\frac{1}{2}\left[e^{2\eta}+\mathrm{ln}2\pi\right].$ (S29)

The local approximation of $p(x)$ with a power law is equivalent to consider

$z(\eta)\approx z(H_{0})+z^{(1)}(H_{0})(\eta-H_{0}).$ (S30)

We proved that the associated order statistics will locally follow a Zipf’s
law whenever $z^{(1)}(H_{0})\neq-1$, i.e. iff

$-e^{2H_{0}}\neq-1\implies H_{0}\neq 0\implies\absolutevalue{\left\langle
X(\lambda_{0})\right\rangle}\neq 1.$ (S31)

In other words, the only two points of the order statistics with an
exponential local behaviour, instead of a Zipf’s one, are those were
$\left\langle X(\lambda_{0})\right\rangle=\pm 1$, i.e. those corresponding to
one standard deviation from the mean value of the Gaussian distribution.