# Normal forms, universal scaling functions, and extending the validity of the
RG
Contribution to the book commemorating Michael Fisher’s contributions to
condensed-matter physics and critical phenomena.

James P. Sethna  LASSP, Cornell University, Ithaca, NY 14853, USA  David
Hathcock  IBM T. J. Watson Research Center, Yorktown Heights, NY 10598  Jaron
Kent-Dobias  INFN, Unità di Roma 1, 00185 Rome, Italy  Archishman Raju  Simons
Centre for the Study of Living Machines, National Centre for Biological
Sciences, Tata Institute of Fundamental Research, Bengaluru 560065, India

###### Abstract

Our community has a deep and sophisticated understanding of phase transitions
and their universal scaling functions. We outline and advocate an ambitious
program to use this understanding as an anchor for describing the surrounding
phases. We explain how to use normal form theory to write universal scaling
functions in systems where the renormalization-group flows cannot be
linearized. We use the 2d Ising model to demonstrate how to calculate high-
precision implementations of universal scaling functions, and how to extend
them into a complete description of the surrounding phases. We discuss
prospects and challenges involved into extending these early successes to the
many other systems where the RG has successfully described emergent scale
invariance, making them invaluable tools for engineers, biologists, and social
scientists studying complex systems.

##  1 Introduction

Half a century ago, Michael Fisher, together with Wilson and Kadanoff,
introduced the renormalization group to analyze systems with emergent, fractal
scale invariance. For five decades, physicists have applied these techniques
to equilibrium phase transitions, avalanche models, glasses and disordered
systems, the onset of chaos, plastic flow in crystals, surface morphologies,
etc. But these tools have not made a substantial impact on engineering and
biology. We believe it is our duty to make these tools accessible to the
broader science community.

$\bullet$ We need to provide them with tools that allow them to describe not
only the critical point, but properties of systems that exhibit incipient
scale-invariant fluctuations yet are far from the critical region. These
demand that we understand corrections to scaling, which become more important
farther from the critical point. In section 3, we build on Michael’s work on
analytic corrections to scaling with Aharony [1, 2] and his work on the
complex analytic properties of the Onsager solution [3] to extend our
understanding of the 2d Ising critical point to a full description of the
ferromagnetic and paramagnetic phases [4].

$\bullet$ We need not only universal power laws, but a complete description of
all behavior of the material. These demand convenient access to accurate
universal scaling functions that govern behavior involving more than two
quantities at a time. Michael taught us about these powerful tools in his
pedagogical reviews and lectures. So, magnetization as a function of
temperature goes as $M\sim t^{\beta}$, but magnetization as a function of
temperature and external field goes as $M\sim
t^{\beta}\mathcal{M}(h/t^{\beta\delta})$, where $\mathcal{M}(X)$ is an
universal function, in principle predicted from the renormalization group. In
section 2, we discuss the correct way [5] of writing universal scaling
functions for systems where logarithms and exponentials are found in addition
to power laws (e.g., in the upper and lower critical dimensions). In section
4, we present a user-friendly solution [6] for $\mathcal{M}(X)$ in the 2d
Ising model, with a systematic expansion that captures the correct singular
behaviors with seven-digit accuracy.

$\bullet$ We need rapidly converging methods that can connect our deep
understanding of the singularities in phases and in universal scaling
functions to quantitative predictions which make the best use of limited
information away from the critical points. In modern numerical analysis, one
can integrate or approximate analytic functions of one variable with
exponential accuracy [7, secs. 4.6,5.8.1] so long as one knows the
singularities at the endpoints. Chebyshev, Gauss, and Romberg methods have
superseded Simpson’s rule and its generalizations for approximating and
integrating analytic functions, and can be adapted to capture known
singularities. We know that all properties are analytic in phases, and our
mission for five decades has been to understand the singularities between
phases. Can we explain the phases using the critical points? The exponential
convergence in sections 3 and 4 demonstrate that the normal form theory can do
so.

Focusing on the 2d Ising model allows us to show proof of principle, but can
we aspire to similar progress for other, less well studied critical points? We
suggest this as a key task for the next stage of research in critical
phenomena and emergent scale invariance. In section 5, we discuss progress by
those using NPFRG (non-perturbative functional renormalization-group) methods
[8], which have found broad application not only in equilibrium thermodynamic
systems, but in avalanche models, quantum systems, turbulence, …They
explicitly coarse-grain and renormalize a system systems in fixed dimension,
implicitly calculating the universal scaling functions and the system-specific
behavior far from criticality. They have mostly been used to extract high-
precision critical exponents, amplitude ratios, and proofs of success. If we
can organize these calculations into user-friendly universal scaling
functions, we may provide experimentalists, simulators, and theorists in a
variety of fields with tools to describe matter beyond systems tuned to a
single point on the phase diagram.

In the final section 6, we embark into deep issues in the renormalization
group and the prospects and challenges they provide to our mission to make the
theory of critical phenomena an organizing principle for much of science.

##  2 Normal form theory, RG flows, and logs

The renormalization group takes an enormous leap of abstraction – studying
emergent scale invariance as a flow under coarse-graining in the space of all
possible systems. This reduces the problem to the study of fixed points for
differential equations in an infinite-dimensional space (e.g., free energy
$f$, temperature $t=(T-Tc)$, field $h$, other parameters $u$, …). Near the
critical point for the Ising model, coarse-graining and rescaling the free
energy $f$ by a factor $e^{\ell}$ is described by the differential equations

$\displaystyle df/d\ell$
$\displaystyle=df+At^{2}+\mathrm{other~{}nonlinear\dots}$ (1) $\displaystyle
dt/d\ell$ $\displaystyle=\lambda_{t}t+\mathrm{nonlinear\dots}$ $\displaystyle
dh/d\ell$ $\displaystyle=\lambda_{h}h+\mathrm{nonlinear\dots}$ $\displaystyle
du/d\ell$ $\displaystyle=\lambda_{u}u+\mathrm{nonlinear\dots}$

The first step in most treatments of these RG flows is to linearize the flow
near the fixed point. Positive eigenvalues $\lambda_{t}$ and $\lambda_{h}$
correspond to relevant operators like field $h$ and temperature $t$; negative
eigenvalues $\lambda_{u}$ correspond to ‘irrelevant’ perturbations like $u$
that provide singular corrections to scaling that are subdominant near the
critical point. Using this linearization, this treatment then argues for
universal power laws for things like the correlation length $\xi(t)\sim
t^{\nu}=t^{1/\lambda_{t}}$ at $h=0$, and universal power laws times universal
functions

$f(t,h,u)\sim
t^{d\nu}\mathcal{F}(h/t^{\beta\delta},ut^{\omega})=t^{d/\lambda_{t}}\mathcal{F}(h/t^{\lambda_{h}/\lambda_{t}},ut^{-\lambda_{u}/\lambda_{t}}),~{}\mathrm{and}$
(2) $M(t,h,u)\sim
t^{\beta}\mathcal{M}(h/t^{\beta\delta},ut^{\omega})=t^{(d-\lambda_{h})/\lambda_{t}}\mathcal{M}(h/t^{\lambda_{h}/\lambda_{t}},ut^{-\lambda_{u}/\lambda_{t}})$
(3)

for things like the field-dependent free energy and magnetization, that
describe relations between more than two quantities.

This linearization is not useful in many cases (e.g., the 1d, 2d, and 4d Ising
models and all models in their upper and lower critical dimensions). How to
systematically formulate universal scaling functions in these cases has
hitherto been mysterious. In this section, we describe the use of normal form
methods from dynamical systems theory by Raju et al. [[5]] to understand when
this linearization is possible, and how to modify the invariant arguments to
the universal scaling functions functions when it is not.

Wegner and co-workers [9, 10] in the early days justified this linearization
by changing variables to nonlinear scaling fields which transform linearly
under the renormalization group. Cardy [11] denotes these new variables
$u_{t}$ and $u_{h}$, but we shall use tildes $\widetilde{t}$ and
$\widetilde{h}$. By choosing a suitable Taylor expansion of the change of
variables,

$\displaystyle t(\widetilde{t},\widetilde{h},\widetilde{u}\dots)$
$\displaystyle=\widetilde{t}+a_{tu}\widetilde{t}\widetilde{u}+a_{th^{2}}\widetilde{t}\widetilde{h}^{2}+\dots$
(4) $\displaystyle h(\widetilde{t},\widetilde{h},\widetilde{u}\dots)$
$\displaystyle=\widetilde{h}+b_{hu}\widetilde{h}\widetilde{u}+b_{ht}\widetilde{h}\widetilde{t}+\dots$
$\displaystyle\dots{}$

the equations simplify to $d\widetilde{f}/d\ell=d\widetilde{f}$,
$d\widetilde{t}/d\ell=\lambda_{t}\widetilde{t}$,
$d\widetilde{u}/d\ell=\lambda_{u}\widetilde{u}$, etc. Aharony and Fisher [1,
2] ten years later noted that this change of variables leads to what we call
analytic corrections to scaling, again subdominant near the critical point.
The analytic corrections to scaling have power laws that involve integers and
combinations of existing critical exponents $\beta$, $\nu$, $\delta$, and
include an analytic background to the free energy; these corrections can be
written in terms of derivatives of the universal scaling function. The
aforementioned singular corrections to scaling introduce new critical
exponents, and become independent variables in the universal scaling
functions.

Dynamical systems theory [12, 5] tells us that linearizing the flow can only
be done for what are called hyperbolic fixed points. The change of variables
$f,t,h,u\to\widetilde{f},\widetilde{t},\widetilde{h},\widetilde{u}$ is
calculated one polynomial order at a time, and its radius of convergence can
be a subtle mathematical issue [13]. In section 3, we shall take on the
ambitious task of attempting to describe the entire surrounding paramagnetic
and ferromagnetic phases for the 2d Ising model. There we shall see that
convergence is tricky, but a good choice of coordinates can yield a radius of
convergence that appears to converge precisely in the range from zero to
infinite temperature.

Even for the two-dimensional Ising model, the RG fixed point cannot be
linearized. The specific heat has a logarithmic singularity: often deemed
$\alpha=0(\log)$, but incompatible with a linearized flow. This arises because
in 2d no polynomial change can remove the $At^{2}$ term in the flow of $f$ in
Eq. 1 (although rescaling the relative magnitudes of $f$ and $t$ can change
the value of $A$ to one). Somewhat messy algebra can confirm that this is due
to a integer resonance between the two linear eigenvalues
$\lambda_{f}=d=2=2/\nu=2*\lambda_{t}$ for the 2d Ising model. The simplest
normal form is thus $d\widetilde{f}/d\ell=2\widetilde{f}-\widetilde{t}^{2}$,
$d\widetilde{t}/d\ell=\widetilde{t}$, etc. This results in a singularity in
the free energy of the form $t^{2}\log(t^{2})$ which will play an important
role in sections 3 and 4.

For the 4d Ising model, the leading irrelevant operator $u$ becomes marginal,
with a zero eigenvalue $\lambda_{u}$. There it has always been clear that one
cannot linearize the RG flow. In the dynamical systems nomenclature, the flow
undergoes a transcritical bifurcation in $d=4$ (together with the same
resonance $d=2/\nu$ seen in 2d Ising above). Our analysis [5] shows that the
normal form for the RG flows is111 The nonlinear term proportional to
$\widetilde{u}\widetilde{h}$ in the equation for $d{\widetilde{h}}/d{\ell}$
needed in the general normal form is known to be zero for the 4d Ising model.
Raju [14] has shown this is because there is a redundant variable proportional
to the magnetization cubed in the renormalization group (see sec. 6).

$d{\widetilde{f}}/d{\ell}=4\widetilde{f}-\widetilde{t}^{2},~{}~{}~{}d{\widetilde{t}}/d{\ell}=2\widetilde{t}-A\widetilde{u}\widetilde{t},~{}~{}~{}d{\widetilde{h}}/d{\ell}=3\widetilde{h},~{}~{}~{}d{\widetilde{u}}/d{\ell}=-\widetilde{u}^{2}+D\widetilde{u}^{3}.$
(5)

The universal scaling of the free energy in a system of length $L$ does not
take the usual scaling form
$f(t,h,u,L)=L^{-d}\mathcal{F}(X,Y,Z)+f_{a}(t,h,u,L)$ with
$X=\widetilde{t}L^{2}$, $Y=\widetilde{h}L^{3}$, and
$Z=\widetilde{u}/L^{\omega/\nu}$. First, we have quite unusual scaling
variables

$X=\widetilde{t}L^{2}\left(W(yL^{1/D})/(1/(D\widetilde{u})-1)\right)^{-A},Y=\widetilde{h}L^{3},~{}~{}\mathrm{and}~{}~{}Z=yL^{1/D}$
(6)

Second, the free energy has a more complex form

$f(t,h,u,L)=L^{-d}f(X,Y,Z)-W(Z)^{-A}\left(\frac{W(Z)^{-A}}{1-A}-\frac{1}{A}\right)+f_{a}(t,h,u,L).$
(7)

Here $W(x)$ is the Lambert $W$ or product log function, $\widetilde{u}$ is the
marginal quartic term in the Landau free energy, $y$ is a messy known function
of $\widetilde{u}$, $A$ and $D$ are the amplitudes of nonlinear terms in the
renormalization group flow (Eq. 5), and $f_{a}$ is a non-singular, analytic
function near the critical point.

Eq. 6 captures the complete, correct singularity for the 4d Ising model; the
traditional log and log-log corrections in the specific heat and the
susceptibility arise from expansions of the Lambert $W$ function for large
arguments. At bifurcations like $d=4$ and resonances as at $d=2$, normal form
theory dictates the nature of the singularity at the critical point.

It also tells us that the free energy in the surrounding phase can be found by
changing variables. So, for example, we know that the liquid-gas critical
point at high pressures and temperatures is in the 3d Ising universality
class. Hence the liquid and gas phase properties should be given by

$f(T,P)=\widetilde{t}(T,P)^{3\nu}\mathcal{F}_{\mathrm{3dIsing}}(\widetilde{h}(T,P)/\widetilde{t}(T,P)^{\beta\delta},\widetilde{u}(T,P)t^{\omega},\dots)+f_{a}(T,P),$
(8)

where $f_{a}$, $\widetilde{t}$, $\widetilde{h}$, and $\widetilde{u}$ are
analytic in temperature and pressure near the critical point. As usual, from
the free energy (and a corresponding scaling form for the correlation
function) one can derive all of the equilibrium and linear-response properties
of the resulting phases. We shall implement a change of coordinates like this
for the 2d Ising model in section 3.

Figure 1: Avalanche sizes for the 2d $T=0$ random-field Ising model (from
[[15]]). (a) Avalanches for disorders $w=0.5$, $1.0$, and $5.0$; each color is
a separate avalanche. (b) Scaling collapse of the area-weighted avalanche size
distribution. Note the factor of ten range in disorder $w$, and the factor of
2000 range in typical avalanche size $\Sigma$. Here the invariant scaling
combination $\Sigma(w)=\Sigma_{s}(B+1/w)^{-Bd_{f}+C}\exp(d_{f}/w)$ is not a
ratio of power laws, but is derived directly from the nonlinear normal form of
the renormalization-group equations in the lower critical dimension (as in Eq.
5). The area-weighted avalanche size distribution is thus
$sA(s|w)=(s/\Sigma(w))^{x}\mathcal{A}((s/\Sigma(w))^{y})$, with $\mathcal{A}$
a universal scaling function.

We highlight the use of our normal-form methods [5] to solve a twenty-year-old
puzzle – the unusual scaling of avalanches in the non-equilibrium 2d random-
field Ising model [15] (see Fig. 1). Over 25 years ago, we used the random-
field Ising model to study avalanches. We could understand the scaling in 3,
4, 5, 7, and 9 dimensions [16], but two dimensions made no sense. Changing the
arguments from ratios of power laws to the invariant functions determined from
normal form theory was one part of the puzzle.222Using random lattices to
suppress faceting was the other obstacle.

The distribution of avalanche sizes in the random-field Ising model is cut off
at size $\Sigma(w)$ depending on the disorder $w$. In three dimensions and
higher it takes the traditional form $\Sigma(w)=w^{-d_{f}\nu}$. In two
dimensions, normal form theory predicts that
$\Sigma=(B+1/w)^{-Bd_{f}+C}\exp(d_{f}/w)$, with $d_{f}$ a universal critical
exponent, and $B$ and $C$ being universal constants associated with
irremovable nonlinear terms in the renormalization-group flow. Using this,
Fig. 1 does an excellent job explaining the behavior.

##  3 Changing coordinates to describe phases: Matching Onsager

As foreshadowed in section 2, normal form theory tells us how to transform
universal scaling forms to describe the entire phases surrounding a critical
point. This section summarizes our recent results [4], which implements this
for the 2d Ising model in zero field, where Onsager’s exact solution [17] for
the free energy enables quantitative validation of the results.

Given the asymptotic scaling form for the free energy (or any other quantity)
in normal form coordinates,
$\widetilde{f}(\widetilde{t},\widetilde{h},\widetilde{u})$, the free energy as
a function of physical temperature, field, and irrelevant perturbations is
$f(t,h,u)=\widetilde{f}(\widetilde{t}(t,h,u),\widetilde{h}(t,h,u),\widetilde{u}(t,h,u))+f_{a}(t,h,u)$.
The functions $\widetilde{t}(t,h,u)$ and $\widetilde{h}(t,h,u)$, are given by
precisely the analytic change of variables that inverts Eq. 4, mapping the
nonlinear scaling variables back to their physical counterparts. Additionally,
we must also add terms $f_{a}(t,h,u)$ accounting for the analytic background
of the free energy. Note that the change of coordinates are linear at lowest
order: for example, $\widetilde{t}\sim t$ and $\widetilde{h}\sim h$.
Therefore, $f(t,h,u)\sim\widetilde{f}(t,h,u)$ near the critical point; in
other words, the normal form free energy $\widetilde{f}$ is the asymptotic
scaling form near the critical point, usually computed using RG and related
techniques.

We have already heard that the 2d Ising model in zero field has a nonlinear
normal form due to the resonance between the free energy and temperature,
$d\widetilde{f}/d\ell=2\widetilde{f}-\widetilde{t}^{2}$,
$d\widetilde{t}/d\ell=\widetilde{t}$, and that these flow equations give rise
to a logarithmic singularity,
$\widetilde{f}=-\widetilde{t}^{2}\log(\widetilde{t}^{2})$. Following the
procedure outlined above, the free energy as a function of temperature is
simply

$f(t)=-\widetilde{t}(t)^{2}\log(\widetilde{t}(t)^{2})+f_{a}(t).$ (9)

This result is consistent with Onsager’s exact solution, which is known to
have the form $a(t)\log t^{2}+b(t)$ for some analytic functions $a$ and $b$.
Comparing these expressions for the free energy, we need to find the
coordinate change $\widetilde{t}(t)$ and analytic background $f_{a}(t)$ by
solving

$a(t)=-\widetilde{t}(t)^{2},\quad\quad
b(t)=\widetilde{t}(t)^{2}\log(\widetilde{t}(t)^{2}/t^{2})+f_{a}(t).$ (10)

Note that because $\widetilde{t}(t)$ is linear to lowest order in t, the term
$\log\left(\widetilde{t}(t)^{2}/t^{2}\right)$ is indeed analytic.

We recently [4] computed the free energy Eq. 9, by perturbatively expanding
Onsager’s exact free energy around the critical point. A key question is the
radius of convergence of the coordinate transformation $\widetilde{t}(t)$ and
$f_{a}(t)$ to the normal form. Unlike Taylor expansions about analytic points,
the radius of convergence of this normal-form analytic expansion about a
singular point is not simply the distance in the complex plane to the next-
nearest singularity. One might hope that physics would govern the convergence
– perhaps the distance to zero temperature, infinite temperature, or the
nearest other phase transition. Indeed, in each of the expansions discussed
below, the critical point is closest to zero temperature, with this distance
determining the radius of convergence.

Our investigations showed the importance of the choice of variable used to
parameterize the distance to the critical point (Fig. 2(a)). Expanding
$\widetilde{t}$ and $f_{a}$ in temperature $t=(T-T_{c})$, for example,
converges in an estimated range $-T_{c}<t<T_{c}$ (or $0<T<2T_{c}$). This is
not the only natural choice, however. The low temperature expansions for the
Ising model are expressed in powers of $X=\exp(-2/T)$. $X$ is also a natural
variable in that the zeros of the 2d Ising partition function in the complex
$X$ plane form a circle passing through $X_{c}$, as Fisher explained [3].
Using the Onsager solution to expand $\widetilde{t}$ and $f_{a}$ in terms of
$x=X-X_{c}$ yields a radius of convergence that extends all the way from zero
temperature to $X=2X_{c}$ (corresponding to $T\approx 4.7\,T_{c}$), but fails
to describe higher temperatures.

Here we used special properties of the 2d Ising model to identify a new
variable

$V=\frac{5-3\sqrt{2}+X}{1+\sqrt{2}+X}\quad\quad v=V-V_{c},$ (11)

which allows our estimated radius of convergence to cover the full physical
temperature range, from zero to infinity (Fig. 2(a)). Our coordinate $v$,
unlike $t$ or $x$, respects the self-dual symmetry of the 2d Ising model.
Fisher’s circle of zeros in $X$ (and thus $x$) breaks this self-dual symmetry.
The linear fractional transformation in Eq. 11 precisely unwraps this circle
of zeros into a straight line, extending the self-dual symmetry to the complex
plane. The circle of zeros becomes the branch cut of the logarithm in the
scaling function for the free energy (Eq. 9).

Figure 2: Capturing the entire phase with analytic corrections, 2d Ising (from
[[4]]). (a) Radius of convergence depends on the expansion variable. Expanding
Onsager’s exact free energy (black line) in $t$, $x$, and $v$ (colored solid
lines) leads to increasing radii of convergence, with the convergence of the
$v$-expansion estimated to cover the entire physical temperature range. Each
expansion shown is $20^{\text{th}}$ order. Determining analytic corrections in
$v$ by matching to the low-temperature expansion does even better, accurately
reproducing the free energy across all temperatures even at low orders (red
dashed line, $6^{\text{th}}$ order). (b) Exponential convergence. Adding
analytic corrections to the universal scaling function at $T_{c}$ by fitting
to the low-temperature expansion results in exponential convergence to the
Onsager solution as we add more terms.

Fig. 2a compares the expansions $f(t)$, $f(x(t))$, and $f(v(t))$ to the exact
free energy. In accordance with the discussion above, we see improved
convergence for the $x$ and $v$ expansions. The expansion in $v$ numerically
appears to have a radius of convergence that extends to the entire physical
range, zero to infinite temperature. Physics here does determine the range of
convergence of normal form theory.

The above results leave two questions: can we improve our approximation for
the free energy further, such that it converges, even at low orders, across
all temperatures? Furthermore, can we approximate the coordinate change and
analytic background without knowledge of the exact solution or the nonlinear
terms in the RG flows? To resolve these challenges, we again expand the free
energy in $v$ using the form given in Eq. 9, but fix the expansion
coefficients of $\widetilde{t}$ and $f_{a}$ by matching to low-temperature
expansions of the free energy, instead of expanding the exact solution.
Importantly, this approach requires minimal knowledge about the critical point
(only the asymptotic scaling form), with most of the information coming from
deep within the low-temperature phase. Because the matching guarantees the
correct low temperature behavior and has the correct log-singularity at the
critical point built into the expansion, we see uniform convergence across all
temperatures. For example, by sixth order, the approximation differs from the
true free energy by at most $0.5\%$ (red dashed line in Fig. 2a) and we see
exponential convergence as we add additional terms (Fig. 2b).

Our approach for extending critical scaling forms to the neighboring phases
should naturally generalize to unsolved statistical physics models and
experimental systems. For example, low-temperature expansions are relatively
easy to compute in all dimensions and for a variety of systems and can be used
to approximate the analytic corrections to established critical scaling forms.
Candidates for this method include the 3d Ising model, where critical
exponents are known to high precision from conformal bootstrap and the 2d
Ising model in a field, whose scaling function is computed to high precision
in the section 4 below. Finally, it is credible that the liquid phase, long-
known as being challenging because there is no ‘small parameter’, could be
described as a perturbation of the liquid-gas critical point, and described as
the Ising critical point plus analytic corrections (determined, for example,
by matching to virial expansions).

##  4 2d Ising critical point in a field

Figure 3:  Comparison between parametric approximations to the scaling
functions with and without the appropriate singularities (from [[6]]). Scaling
functions for the magnetization and susceptibility plotted as functions of the
scaling invariant $h|t|^{-\beta\delta}$. The errors in the blue curves (Kent-
Dobias) are estimated to be roughly $10^{-7}$. Caselle’s earlier work [18]
that we build upon has significant errors near $h=0$ (red point); these are
due to the essential singularity at $h=0$ and $t<0$ that we address. These
discrepancies grow larger with higher derivatives of the free energy, as shown
in the right panel of Fig. 5.

Ever since Onsager solved the zero-field 2d Ising model, people have searched
for a high-precision approximate solution for the 3d Ising model, and
occasionally also for the 2d Ising model in a field. The most successful of
these attempts make use of parametric coordinates to interpolate behavior in
temperature and field around the critical point in such a way that the
singularity is naturally incorporated. [19, 18] These coordinates, first
introduced many decades ago, [20] are polar-like, with a radius-like component
$R$ that controls the proximity to the critical point, and an angle-like
coordinate $\theta$ that rotates around the temperature–field plane. When
taken to parameterize the temperature– _magnetization_ plane, such coordinates
can be periodic in $\theta$ like true polar coordinates, something once
studied by Michael and collaborators to describe the coexistence region. [21]
However, when parameterizing the temperature– _field_ plane, there is an
inevitable cut resulting from the discontinuity between the up and down states
at low temperatures, and $\theta$ is not a periodic coordinate but ends at the
abrupt transition at some $\theta_{0}$.

The power of adopting these coordinates is that, as a function of $\theta$,
the critical singularity does not appear, and one can reasonably well-
approximate the scaling function in this coordinate as a simple analytic
function. Using this principle, Caselle _et al._ did a remarkably careful job
of creating a function that satisfied lots of known properties and measured
values (in yellow of Fig. 3). Their method was doomed to slow convergence at
small external fields, though, because while the singularity of the critical
point was indeed removed by the coordinate change, other more subtle
singularities in the free energy remain, including a key essential singularity
as one crosses the $h=0$ line for $t<0$. Fig. 3 shows a comparison between the
approximate form of Caselle _et al._ and that of our work including this
singularity.

Figure 4: Contour in the complex $\theta$ plane to generate a free energy
scaling form with the correct singularities (from [[6]]).

In order to rigorously realize the idea of a free energy scaling function
broken into singular and analytic pieces, two singularities in the scaling
function must be accounted for in addition to the critical one. One is ancient
and well known to the experts: the Yang–Lee edge singularity, for $T>T_{c}$ in
the complex plane.[22, 23] The other has only recently been realized to be
relevant, but has a nice physical picture. As one crosses $h=0$, the
equilibrium magnetization jumps. But if you cross just a little bit, the up-
spin magnetized state is metastable, lasting for a good while until a bubble
of down spins forms and grows to flip the system. (Think of supersaturated
vapor – 101% humidity – and the nucleation of raindrops.) Here the surface
tension cost between the up and down regions for small droplets is bigger than
the bulk free energy gain of aligning the interior with the external field.
The energy barrier $B$ to reach large sizes diverges diverges as $h$ vanishes,
$B/k_{B}T=c/h$ for some constant $c$. Just like in quantum mechanics, where
the lifetime of a state gives the energy an imaginary part, the free energy
becomes complex for $h<0$ with an essential singularity $\mathrm{Im}{f}\sim
e^{(c/h)}$. One can use a Kramers–Kronig transform to see that the real part
also has an essential singularity, influencing $f$ for $h>0$ as well.

Figure 5:  Convergence of the universal scaling function (from [[6]]). Left:
The difference between the $n$th order approximation for the free energy
scaling function in a field and the 6th order approximation as a function of
the scaling invariant $th^{-1/\beta\delta}$. The 5th and 6th order
approximations differ by at most $2\times 10^{-6}$. Higher derivatives behave
similarly. Right: The Taylor series coefficients $\mathcal{F}_{-}^{(m)}(0)$ of
the free energy scaling function at the abrupt transition line (the location
of the essential singularity) as a function of derivative $m$. By
incorporating the essential singularity into the scaling function explicitly,
the approximate forms have zero radius of convergence, matching numeric
measurements of these coefficients.

We can incorporate these two singularities into the universal scaling function
by first generating a ‘simplest’ functional form that has the correct
singularities, and then changing variables
$\theta\to\widetilde{\theta}=g(\theta)$ and adding an overall analytic
function $\mathcal{F}_{a}(\theta)=G(\theta)$, in a precise analogy with the
analytic corrections we introduce to match the universal scaling function to
describe the surrounding phases (sec. 3).

The ‘simplest’ form is taken from writing the most singular part of the
imaginary part of the free energy, and then taking advantage of a
Kramers–Kronig like relation to find the corresponding real part. The
requisite contour integral in the complex-$\theta$ plane is shown in Fig. 4.

Though complicated, this procedure pays dividends. By incorporating these
singularities and matching the power series terms in $g(\theta)$ and
$G(\theta)$ (analogous to $\widetilde{t}(t)$ and $f_{a}(t)$ in section 3), we
are able to achieve exponential convergence of the scaling function to exactly
known values at $t=0$, and also achieve exponential convergence of
derivatives. The left-hand side of Fig. 5 shows this convergence in the free
energy itself as a function of the scaling invariant $th^{-\beta\delta}$ by
subtracting our best converged approximation $\mathcal{F}^{[6]}$ from the
lower-order approximations $\mathcal{F}^{[n]}(\theta)$ for
$n\in\\{2,3,4,5\\}$. These provide evidence that our seven-digit convergence
at $t=0$ extends to the whole scaling function. On the right of this plot, we
can see the origin of this good behavior: the series expansion for the free
energy around the abrupt transition point has zero radius of convergence, but
this is captured naturally by our approximate scaling form.

##  5 Interpolating scaling functions between dimensions

Figure 6: Interpolating between dimensions. A schematic showing the dimensions
and spin components where traditional power-law scaling will break down. This
leads to logarithms and exponentials replacing power laws, and rich but
complicated invariant scaling combinations as arguments for the universal
scaling functions [5]. On the other hand, these complexities also promise to
inform and accelerate convergence of the universal scaling functions for
dimensions in between.

Figure 6 illustrates the traditional thermodynamic critical points. Over five
decades, these have been extensively explored using $\epsilon$-expansions,
$1/N$ expansions, cluster expansions, etc. The universal critical exponents
are known essentially exactly from conformal bootstrap methods.333In the same
spirit, Kent-Dobias’ solution [6] for the scaling function in 2d in section 4
is also essentially exact: both are well-defined algorithms that generate an
exponentially converging approximation in a form useful for applications.
After all, this is what we call the ‘exact’ solution for $\sin(x)$. Can we do
the same for the universal scaling functions as functions of dimension $d$ and
the number of spin components $N$? And for (say) the random-field Ising model
[24], or turbulence [25]?

Interpolation between dimensions could have several benefits. First, much is
known analytically in two and four dimensions about the universal scaling
functions that is only known numerically in three. Second, there are important
features in scaling functions and their corrections to scaling that are a
clear foreshadowing of properties in other dimensions. So the leading
correction to scaling in 3d is the echo of the marginal variable of the
Wilson-Fisher transcritical bifurcation in 4d. And, in a more dramatic
example, the universal scaling function for the avalanche size distribution in
the 3d random-field Ising model has a striking feature (it grows by a factor
of ten before cutting off exponentially [16]), which is clearly related to the
unusual scaling of the avalanche sizes in two dimensions [15] (Fig. 1).

In ordinary differential equations, normal form theory provides not only the
behavior at the transition, but also an unfolding of the behavior near the
bifurcation (up to a smooth coordinate transformation). Here we may expect to
use this, not just to describe the phases near the critical points, but
(following the lead of Wilson and Fisher) to describe the universal scaling
functions as they evolve between dimensions.

How would we interpolate scaling functions, when even the arguments of the
scaling functions vary with dimension (as in the $W$-functions we needed in
section 2 for the 4d Ising model)? Consider taking the 4d flow equations 5
about the mean-field fixed-point, and keeping the unremovable nonlinear terms
to describe nonlinear flows about the mean-field fixed point in $d$
dimensions:

$d{\widetilde{f}}/d{\ell}=d\widetilde{f}-\widetilde{t}^{2},~{}~{}~{}d{\widetilde{t}}/d{\ell}=2\widetilde{t}-A\widetilde{u}\widetilde{t},~{}~{}~{}d{\widetilde{u}}/d{\ell}=(4-d)\widetilde{u}-\widetilde{u}^{2}+D\widetilde{u}^{3},~{}~{}~{}$
(12)

Things to note. 1) The normal-form procedure of removing the nonlinear terms
allows us to set the nonlinear terms in any way we wish. 2) We know the
Wilson-Fisher fixed point for Eqs. 12 at
$(\widetilde{t}^{*},\widetilde{u}^{*})$ must have a linearization whose
eigenvalues $1/\nu$ and $\omega$ that, as a function of dimension, matches the
conformal bootstrap values (or the exact values in 2d). This constrains the
two nonlinear terms $A(d)$ and $D(d)$. Serendipitously, the
$\widetilde{t}^{2}$ resonance term in the 4d free energy flow is also needed
to get the logarithmic specific heat in 2d. A high-precision universal scaling
function for the free energy using the nonlinear scaling variables
$\widetilde{f}$, $\widetilde{t}$ and $\widetilde{u}$ in Eq. 12 would describe
the universal crossover scaling from mean field to short-range magnetization
in all dimensions between two and four. One could then linearize that scaling
function about the 3d Wilson-Fisher fixed point to extract the traditional
scaling functions.

How would one calculate such scaling functions? Here serendipity strikes in
recent advances in non-perturbative functional renormalization-group (NPFRG)
calculations [8], whose critical exponents [26] are almost competitive with
those of conformal bootstrap, but which in the process also compute a
functional form for the coarse-grained free energy as a function of
magnetization.

Matthew Tissier and colleagues at the Jussieu campus of the Sorbonne and the
first author have been exploring how to calculate high-precision universal
scaling functions by combining the strengths of the normal-form theory with
the systematically improvable scaling functions of NPFRG. In initial work, we
have found that the NPFRG for the 4d Ising model should indeed have a scaling
function whose arguments are of the exotic form in Eq. 6. We are now learning
to solve the partial differential equations for $\widetilde{t}(t,u)$
$\widetilde{u}(t,u)$ and $\widetilde{f}(t,u)$ in 3d to extract the universal
scaling functions from the simplest of NPFRG models. We hope then to use the
technology in section 4 to demonstrate how to tabulate universal scaling
functions in a form convenient for experimentalists and simulators.

The NPFRG method has had striking success in high-precision calculations of
equilibrium thermodynamic systems, disordered thermodynamic and avalanche
models, fully developed turbulence, quantum many-body systems, and in QCD and
electroweak models [8]. Our hope and vision is to inspire our colleagues to
gather and tabulate their results into universal scaling functions that can
describe the behaviors both near and far from points of emergent scale
invariance, in a way accessible to engineers, biologists, and social
scientists studying complex systems.

##  6 Prospects and challenges for future work

We have summarized what we believe to be promising indications that we can
indeed generate usable, high-precision universal scaling functions in a
systematic way, and can extend them systematically into high-precision
descriptions of the surrounding phases. If this continues to work, one
imagines our community will be able to provide quantitative theories of
liquids, materials plasticity and failure, fluctuations in biological systems,
and potentially turbulence and glass behavior. In all these cases, we must
mesh our understanding of the universal, emergent critical behavior and our
system-dependent knowledge of properties away from criticality to describe the
entire phases. Will this work in practice?

To show that our method works to high precision, we chose to use a problem
where we knew a great deal about the answer: the 2d Ising model. The reader
should legitimately be skeptical that this is the easiest case. In particular,
the 2d Ising mode is special in that it (1) is self-dual, (2) is analytically
solvable in zero field, and (3) has no singular corrections to scaling. Let us
discuss causes for concern that our methods will become general tools, and
reasons for optimism.

(1) The 2d Ising model is self-dual. In section 3, we found exponential
convergence of the Ising free energy $f(t)$ in zero field using a power series
of the Onsager solution about the critical point. To do so, we expanded in a
variable $v(t)$ that was fine-tuned to the self-dual and complex analytic
properties of the 2d Ising model. We then used the low-temperature expansion
to avoid any use of Onsager’s solution. Will we need to find the ‘right’
variable in 3d to proceed? It turns out that matching to properties outside
the radius of convergence stabilizes the expansion. In preliminary work,
matching both high and low temperature series using the standard variable
$x=X-X_{c}=e^{-2/T}-e^{-2/T_{c}}$ in 2d (ignoring duality), we found
exponential convergence in the entire temperature region [4].

(2) The 2d Ising model in zero field has an exact solution, and also has an
exact solution at $T_{c}$ as a function of field.444It may seem in retrospect
clear that the essential singularity at $h=0$ would frustrate a search for an
exact solution in a field below $T_{c}$. But will the external field prevent
us from finding an exponentially converging solution using analytic
corrections to scaling $\widetilde{t}(t,h)$, $\widetilde{h}(t,h)$, and
$f_{a}(t,h)$, by matching the 2d universal scaling function from sec. 4 to the
low $T$/high $h$ and high $T$/low $h$ cluster expansions, as we did for zero
field in sec. 3?

We have an exact solution for the magnetization – the first derivative with
respect to field. This has allowed us to determine $\widetilde{t}$ and
$\widetilde{h}$ to linear order in $h$. The susceptibility $\chi(T)$ at zero
field is not known exactly, but a remarkable amount is known [27]: there are
formal expansions strongly indicating subdominant logarithmic corrections to
scaling, and a possible essential boundary in the complex plane at the Fisher
circle of zeros in the complex temperature plane [28, 29, 30, 27]. Log
corrections to scaling are known to arise from resonances between relevant and
irrelevant eigenvalues (sec. 2); the observed logs could arise from irrelevant
variables whose resonances contribute only in $O(h^{2})$ to the free energy.
It appears that we can use these results [27] not only to determine
$\widetilde{t}$ and $\widetilde{h}$ to quadratic order in $h$, but also to
extract information about irrelevant variables and their singular corrections
to scaling in 2d (see below).

The universal scaling function in 3d presumably should succumb to the same
kind of systematic approximation that we used in sec. 4 for 2d. All the
ingredients appear to be available.555And, of course, high-precision NPFRG
calculations of the critical exponents [26] implicitly have calculated the
universal scaling function. The Yang-Lee singularity in 3d [31] has yielded to
a high-precision NPFRG calculation. The essential singularity calculation
depends only on the value of the surface tension near the critical point. In
2d, we used estimates of derivatives of the universal scaling function near
$h=0$ for $T>T_{c}$ and $T<T_{c}$. The latter (shown in Fig. 5b) are
challenging to calculate using NPFRG methods (perhaps because of the essential
singularity), but the former are now available [32]. Indeed, these authors
used their estimates and the traditional Schofield coordinates to estimate the
scaling function; adding the proper droplet singularity and Yang-Lee
singularities could be straightforward and hopefully will yield convincing
exponential convergence.

Can we expect to extend a high-precision 3d Ising scaling funciton to describe
the in-field behavior in the entire surrounding phases? Using analytic
corrections to extend the critical singularity could break down at the
roughening transition.666While one might argue that the behavior of an
interface between spin up and spin down regions should not affect the bulk
free energy, the (incredibly subtle) essential singularity as one crosses
$h=0$ (sec. 4) depends on the surface energy of the critical droplet, which
will itself have an essential singularity as a function of temperature as it
develops facets at the roughening temperature. However, this subtle transition
should in principle also impede the use of Dlog Padé methods using the low-
temperature expansion of the magnetization to describe the critical-point
behavior, which works well in practice. In any case, the roughening transition
does not arise in isotropic systems, so these concerns will not interfere with
finally extracting a quantitative theory of liquids by adiabatic continuation
of our understanding of the 3d Ising critical point.

(3) The Onsager solution exhibits no singular corrections to scaling. First,
this statement is more subtle than it seems. Onsager’s solution has a pure
logarithmic singularity, but conceivably there could be irrelevant eigenvalues
with integer exponents that would lead to analytic terms in the free energy
and magnetization. Also, there are irrelevant anisotropies in the correlation
functions and the susceptibilities due to the square lattice, indeed with
integer exponents. Even more subtle, these anisotropies are singular
corrections to scaling in Wilson’s RG formulation, but are clearly present at
the fixed point of the real-space RG. Any coarse-graining procedure
maintaining a square lattice will have a short-range square anisotropy.

Experts will remember that some of the irrelevant operators in the RG are
redundant, and these redundant operators often involve flows between different
fixed points in the same universality class.777We have recently looked at this
question [33] in the context of the period doubling onset of chaos in iterated
maps $f(x)$. About half of the irrelevant eigenvalues are negative powers of
the relevant critical exponent $\alpha$, and are redundant: one can set the
amplitudes of their correction to scaling to zero using a coordinate change
$x=\phi(y)$, so instead of $\alpha f(f(x/\alpha))$ one considers the
renormalization group $\alpha\phi(f(f(\phi^{-1}(x)/\alpha))$. We called these
gauge irrelevant, as they depend on the choice of coordinate. A similar gauge
invariance arises for the quasiperiodic onset of chaos [34, 35] that we
studied when the first author was a postdoc. The idea of redundant variables
had been first explored by Wegner [36], who pointed out that a RG
transformation would always have certain operators that are redundant; they
correspond to infinitesimal redefinitions of the field, and do not contribute
to the singular part of the free energy.888How does one think of these
redundant operators in the discrete theory? Murthy and Shankar had found a
clever way to think about redundant operators for the discrete Ising model
[37]. They essentially used a version of Kadanoff’s RG projection operator
without reducing the degrees of freedom. They found evidence that these
redundant operators had observable consequences in simulations of the 2d Ising
model [38]. Evidence for redundant operators had already been found in the 3d
Ising model [39].

The first author’s first graduate student, Mohit Randeria, worked for Michael
Fisher before shifting to my group. The two of them wrote a comment on a paper
by Swendsen [40] claimed that one could modify the RG transformation to move
the fixed point anywhere on the critical surface. Randeria and Fisher asserted
that this could not in general be correct, because the amplitudes of the
singular corrections to scaling are zero at the fixed points. To quote them,

> “One fixed point may be mapped into another by a change of redundant
> variables and two RGs, say $\mathcal{R}$ and $\widehat{\mathcal{R}}$, may
> produce formally different fixed points by this mechanism. Nevertheless, a
> general point on a critical manifold cannot be transformed into a fixed
> point nor vice versa.”

In the context of the 2d Ising model, there seem to be no singular corrections
in the exact solution. Barma and Fisher had predicted that there should be
logarithmic corrections to scaling coming from irrelevant variables [41]. They
had also found evidence for an irrelevant operator with exponent $-4/3$ [42].
Conformal field theory predicts a large number of irrelevant variables which
come from descendant operators but does not predict the $-4/3$ exponent. It
has been conjectured in the past that all of these descendent operators are
redundant [43]. This may explain why they do not lead to any logarithmic
corrections to scaling (which otherwise are generically expected [14]). We
would consider all of these descendant operators as constituting ‘gauge’
corrections to scaling, which can be removed by an appropriate coordinate
choice, whereas the $-4/3$ exponent would be a genuine singular correction to
scaling.

In sec. 3 we succeeded in solving for the free energy of the square-lattice
Ising model at all temperatures by applying a normal-form change of variables
to the universal scaling function for the critical point. Is that success due
to the lack of (non-redundant) singular corrections to scaling for this
particular universality class?999It is believed that the powers of logarithms
in the 2d Ising susceptibility vanish in isotropic systems. Since all
irrelevant operators in 2d have rational eigenvalues, and the redundant ones
cannot generate resonant logarithms, one suspects that all the isotropic
corrections to scaling are redundant. It is clear that in 3d we shall need to
incorporate the (genuine) singular corrections to scaling in order to
accurately describe the behavior away from the fixed point. Indeed, these
contributions are subdominant near the fixed point, meaning conversely that
they grow faster as one leaves the fixed point than do the relevant operators.
Thus we will need a universal scaling function that depends not only on
$\widetilde{t}$ and $\widetilde{h}$, but on potentially an infinite family of
irrelevant variables $\widetilde{u}$. Can we nonetheless hope for exponential
accuracy?

The NPFRG methods mentioned in sec. 5, as it happens, precisely implement a
flow with many irrelevant eigendirections at the fixed point. Applying our
normal form analysis, we have been able to linearize this flow in the simplest
case, allowing us to extract a functional form for the leading correction
subdominant as $\widetilde{t}^{\omega}$. It seems likely that one could
generate exponential convergence for all properties of the phases by combining
increasingly sophisticated NPFRG calculations, adding more irrelevant
operators to our scaling functions, and matching to higher orders of low and
high temperature expansions (or virial expansions, $1/N$ expansions …).

Finally, there are intriguing hints of deep connections between normal form
theory, redundant operators, and the construction of universal scaling
functions.

We noted in sec. 4 that the procedure we used to generate the universal
scaling function recapitulated exactly the same steps as the procedure we used
in sec. 3 to generate the Onsager solution from the critical point. In both
cases, we knew the singularities of the functions involved, and fit parameters
in two functions (an analytic change of variables and an additive analytic
background) to match known properties measured separately. If this systematic
procedure can be generally used as a best practice, rapid progress could be
made.

There is also a striking analogy between redundant operators (removing some
singular corrections to scaling by reparameterizing the space of predictions)
and our normal form theory (removing analytic corrections to scaling by
reparameterizing the control parameters). Can we simultaneously remove both
with a joint transformation? How are they related?

In thermodynamics, whether you write the free energies as a function of the
external field and temperature or as a function of magnetization and
temperature is considered a matter of choice. A Legendre transformation from
control parameters $(f,t,h)$ to $(u,t,m)$ swaps the control parameter $h$ with
the prediction $m$: does it swap normal-form transformations to redundant
ones? (Indeed, the NPFRG methods coarse-grain the fixed-magnetization Gibbs
free energy $u(t,m)$, rather than Wilson’s free energy at fixed field.) When
we use our normal form theory to remove nonlinear terms in the RG flows by
changing $h$ to $\widetilde{h}(t,h)$, we change the predicted magnetization.
Is that change a redundant one?

Hankey and Stanley [44] showed that if you Legendre transform the generalized
homogeneous functions we use for hyperbolic fixed points, you get another
homogeneous function. Later, others, studying finite-size scaling of the Ising
model [45, 46] in the microcanonical ensemble, predicted that it would be
possible to speak of the entropy as a scaling function of the energy instead
of the free energy as a function of the temperature.

Our preliminary work on this question ran into a difficulty, generalizing
these ideas to the nonlinear normal forms necessary in the upper and lower
critical dimensions. We have found that it is possible to Legendre transform
RG flows at the linear level, but nonlinear flows generically lead to non-
analyticities in the RG flow [14, 47].101010As noted in the footnote on page
1, even in four dimensions the RG flow of the magnetization is linearizable.
Presumably this means that implementing an RG using flows in $(f,t,h)$ and
$(u,t,m)$ will be possible even in four dimensions, while we find $(f,t,h)$ to
$(S,E,h)$ leads to non-analytic RG flows. Whether this clarifies or confuses
the correspondence between redundant and normal-form changes of variables
needs to be explored.

As a long-time colleague of Michael Fisher, the lead author hopes that Fisher
would have enjoyed this deep plunge into the complex history of the field and
our ambitions for the future. He also hopes that Fisher’s friends, colleagues,
and collaborators who are contributing to this book, and of course our
colleagues of the future who are reading it, will find this chapter more
illuminating than obscure, and more useful than misleading or misguided.

## Acknowledgments

JPS would like to thank Matthieu Tissier for his insights and patience
explaining the NPFRG methods, and for continuing collaborations in this work,
Colin Clement for careful early work on the 2d Ising model, and Jacques Perk
for helpful correspondence. JPS benefited by funding from NSF DMR-1719490 and
CNRS (the French National Center for Scientific Research). DH was partially
supported by an NSF Graduate Research Fellowship Grant No. DGE-2139899. JK-D
acknowledges Simons Foundation Grant No. 454943. AR acknowledges support from
the Department of Atomic Energy, India under project No. RTI4006, and the
Simons Foundation Grant No. 287975.

## References

  * [1] A. Aharony and M. E. Fisher, University in analytic corrections to scaling for planar Ising models, _Phys. Rev. Lett._ 45, 679–682 (Sep, 1980).  doi: 10.1103/PhysRevLett.45.679.  URL https://link.aps.org/doi/10.1103/PhysRevLett.45.679.
  * [2] A. Aharony and M. E. Fisher, Nonlinear scaling fields and corrections to scaling near criticality, _Phys. Rev. B_. 27, 4394–4400 (Apr, 1983).  doi: 10.1103/PhysRevB.27.4394.  URL https://link.aps.org/doi/10.1103/PhysRevB.27.4394.
  * [3] M. E. Fisher.  The nature of critical points.  In _Boulder Lectures in Theoretical Physics: Statistical Physics, Weak Interactions, Field Theory_ , vol. 7c, pp. 1–159. University of Colorado Press, (1964).
  * [4] D. Hathcock and J. P. Sethna, Phase transitions beyond criticality: normal form and analytic corrections for the 2d ising model, _Manuscript in preparation_.
  * [5] A. Raju, C. B. Clement, L. X. Hayden, J. P. Kent-Dobias, D. B. Liarte, D. Z. Rocklin, and J. P. Sethna, Normal form for renormalization groups, _Phys. Rev. X_. 9, 021014 (Apr, 2019).  doi: 10.1103/PhysRevX.9.021014.  URL https://link.aps.org/doi/10.1103/PhysRevX.9.021014.
  * [6] J. Kent-Dobias and J. P. Sethna.  Smooth and global Ising universal scaling functions.  https://arxiv.org/abs/1707.03791, (2021).
  * [7] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, _Numerical Recipes 3rd Edition: The Art of Scientific Computing_. (Cambridge University Press, 2007).
  * [8] N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. Pawlowski, M. Tissier, and N. Wschebor, The nonperturbative functional renormalization group and its applications, _Physics Reports_. 910, 1–114, (2021).  ISSN 0370-1573.  doi: https://doi.org/10.1016/j.physrep.2021.01.001.  URL https://www.sciencedirect.com/science/article/pii/S0370157321000156.  The nonperturbative functional renormalization group and its applications.
  * [9] F. J. Wegner, Corrections to scaling laws, _Physical Review B_. 5, 4529, (1972).
  * [10] F. J. Wegner and E. K. Riedel, Logarithmic corrections to the molecular-field behavior of critical and tricritical systems, _Physical Review B_. 7, 248, (1973).
  * [11] J. Cardy, _Scaling and renormalization in statistical physics_. vol. 5, (Cambridge university press, 1996).
  * [12] S. Wiggins, _Introduction to applied nonlinear dynamical systems and chaos_. vol. 2, (Springer Science & Business Media, 2003).
  * [13] J. Guckenheimer.  (private communication).
  * [14] A. Raju.  _Aspects of the renormalization group_.  PhD thesis, Cornell University, (201).  https://sethna.lassp.cornell.edu/pubPDF/RajuPhD.pdf.
  * [15] L. X. Hayden, A. Raju, and J. P. Sethna, Unusual scaling for two-dimensional avalanches: Curing the faceting and scaling in the lower critical dimension, _Phys. Rev. Research_. 1, 033060 (Oct, 2019).  doi: 10.1103/PhysRevResearch.1.033060.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.033060.
  * [16] O. Perkovic, K. A. Dahmen, and J. P. Sethna, Disorder-induced critical phenomena in hysteresis: Numerical scaling in three and higher dimensions, _Physical Review B_. 59, 6106–6119, (1999).
  * [17] L. Onsager, Crystal statistics. i. a two-dimensional model with an order-disorder transition, _Phys. Rev._ 65, 117–149 (Feb, 1944).  doi: 10.1103/PhysRev.65.117.  URL https://link.aps.org/doi/10.1103/PhysRev.65.117.
  * [18] M. Caselle, M. Hasenbusch, A. Pelissetto, and E. Vicari, The critical equation of state of the two-dimensional Ising model, _Journal of Physics A: Mathematical and General_. 34(14), 2923–2948 (3, 2001).  doi: 10.1088/0305-4470/34/14/302.  URL https://doi.org/10.1088%2F0305-4470%2F34%2F14%2F302.
  * [19] R. Guida and J. Zinn-Justin, 3d Ising model: the scaling equation of state, _Nuclear Physics B_. 489(3), 626–652, (1997).
  * [20] P. Schofield, Parametric representation of the equation of state near a critical point, _Physical Review Letters_. 22(12), 606–608 (3, 1969).  doi: 10.1103/physrevlett.22.606.  URL https://doi.org/10.1103%2Fphysrevlett.22.606.
  * [21] M. E. Fisher, S.-y. Zinn, and P. J. Upton, Trigonometric models for scaling behavior near criticality, _Physical Review B_. 59(22), 14533–14545 (6, 1999).  doi: 10.1103/physrevb.59.14533.  URL https://doi.org/10.1103%2Fphysrevb.59.14533.
  * [22] C. N. Yang and T. D. Lee, Statistical theory of equations of state and phase transitions. i: Theory of condensation, _Physical Review_. 87(3), 404–409 (8, 1952).  doi: 10.1103/physrev.87.404.  URL https://doi.org/10.1103%2Fphysrev.87.404.
  * [23] T. D. Lee and C. N. Yang, Statistical theory of equations of state and phase transitions. ii: Lattice gas and ising model, _Physical Review_. 87(3), 410–419 (8, 1952).  doi: 10.1103/physrev.87.410.  URL https://doi.org/10.1103%2Fphysrev.87.410.
  * [24] I. Balog, G. Tarjus, and M. Tissier, Dimensional reduction breakdown and correction to scaling in the random-field Ising model, _Phys. Rev. E_. 102, 062154 (Dec, 2020).  doi: 10.1103/PhysRevE.102.062154.  URL https://link.aps.org/doi/10.1103/PhysRevE.102.062154.
  * [25] A. Gorbunova, G. Balarac, L. Canet, G. Eyink, and V. Rossetto, Spatio-temporal correlations in three-dimensional homogeneous and isotropic turbulence, _Physics of Fluids_. 33(4), 045114, (2021).  doi: 10.1063/5.0046677.  URL https://doi.org/10.1063/5.0046677.
  * [26] G. De Polsi, I. Balog, M. Tissier, and N. Wschebor, Precision calculation of critical exponents in the O(N) universality classes with the nonperturbative renormalization group, _Phys. Rev. E_. 101, 042113 (Apr, 2020).  doi: 10.1103/PhysRevE.101.042113.  URL https://link.aps.org/doi/10.1103/PhysRevE.101.042113.
  * [27] Y. Chan, A. J. Guttmann, B. Nickel, and J. Perk, The ising susceptibility scaling function, _Journal of Statistical Physics_. 145, 549–590, (2011).
  * [28] B. Nickel, On the singularity structure of the 2d ising model susceptibility, _Journal of Physics A: Mathematical and General_. 32(21), 3889–3906 (1, 1999).  doi: 10.1088/0305-4470/32/21/303.  URL https://doi.org/10.1088%2F0305-4470%2F32%2F21%2F303.
  * [29] W. Orrick, B. Nickel, A. Guttmann, and J. H. Perk, The susceptibility of the square lattice ising model: new developments, _J. Stat. Phys._ 102, 795–841, (2001).
  * [30] M. Assis, J. L. Jacobsen, I. Jensen, J.-M. Maillard, and B. M. McCoy, Analyticity of the ising susceptibility, _Journal of Physics A: Mathematical and Theoretical_. 50(36), 365203 (August, 2017).  doi: 10.1088/1751-8121/aa8128.  URL https://doi.org/10.1088%2F1751-8121%2Faa8128.
  * [31] A. Connelly, G. Johnson, F. Rennecke, and V. V. Skokov, Universal location of the Yang-Lee edge singularity in O(N) theories, _Phys. Rev. Lett._ 125, 191602 (Nov, 2020).  doi: 10.1103/PhysRevLett.125.191602.  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.191602.
  * [32] G. D. Polsi, G. Hernández-Chifflet, and N. Wschebor, Precision calculation of universal amplitude ratios in o$(n)$ universality classes: Derivative expansion results at order $\mathcal{o}(\partial^{4})$, _Physical Review E_. 104(6), 064101 (12, 2021).  doi: 10.1103/physreve.104.064101.  URL https://doi.org/10.1103%2Fphysreve.104.064101.
  * [33] A. Raju and J. P. Sethna.  Reexamining the renormalization group: Period doubling onset of chaos.  https://arxiv.org/abs/1807.09517, (2018).
  * [34] D. Rand, S. Ostlund, J. Sethna, and E. D. Siggia, Universal transition from quasiperiodicity to chaos in dissipative systems, _Phys. Rev. Lett._ 49, 132–135 (Jul, 1982).  doi: 10.1103/PhysRevLett.49.132.  URL http://link.aps.org/doi/10.1103/PhysRevLett.49.132.
  * [35] S. Ostlund, D. Rand, J. P. Sethna, and E. D. Siggia, Universal properties of the transition from quasi-periodicity to chaos in dissipative systems, _Physica D_. 8, 303–342, (1983).
  * [36] F. Wegner, Some invariance properties of the renormalization group, _Journal of Physics C: Solid State Physics_. 7(12), 2098, (1974).
  * [37] G. Murthy and R. Shankar, Redundant operators for ising spins, _Physical Review B_. 32(9), 5851, (1985).
  * [38] R. Shankar and R. Gupta, Clear evidence of redundant operators in monte carlo studies of the ising model, _Physical Review B_. 32(9), 6084, (1985).
  * [39] G. Pawley, R. Swendsen, D. Wallace, and K. Wilson, Monte carlo renormalization-group calculations of critical behavior in the simple-cubic ising model, _Physical Review B_. 29(7), 4030, (1984).
  * [40] R. H. Swendsen, Optimization of real-space renormalization-group transformations, _Physical review letters_. 52(26), 2321, (1984).
  * [41] M. Barma and M. E. Fisher, Corrections to scaling and crossover in two-dimensional ising and scalar-spin systems, _Physical Review Letters_. 53(20), 1935, (1984).
  * [42] M. Barma and M. E. Fisher, Two-dimensional ising-like systems: Corrections to scaling in the klauder and double-gaussian models, _Physical Review B_. 31(9), 5954, (1985).
  * [43] H. W. Blöte and M. P. den Nijs, Corrections to scaling at two-dimensional ising transitions, _Physical Review B_. 37(4), 1766, (1988).
  * [44] A. Hankey and H. E. Stanley, Systematic application of generalized homogeneous functions to static scaling, dynamic scaling, and universality, _Physical Review B_. 6(9), 3515, (1972).
  * [45] R. C. Desai, D. W. Heermann, and K. Binder, Finite-size scaling in a microcanonical ensemble, _Journal of statistical physics_. 53, 795–823, (1988).
  * [46] M. Kastner, M. Promberger, and A. Hüller, Microcanonical finite-size scaling, _Journal of Statistical Physics_. 99, 1251–1264, (2000).
  * [47] C. B. Clement.  _Inference and Modeling in Physics: Image Reconstruction, Spin Glasses, and the Renormalization Group_.  PhD thesis, Cornell University, (2019).  https://sethna.lassp.cornell.edu/pubPDF/ClementPhD.pdf.


# How many digits are needed?

I. W. Herbst, J. Møller, A. M. Svane

###### Abstract

Let $X_{1},X_{2},...$ be the digits in the base-$q$ expansion of a random
variable $X$ defined on $[0,1)$ where $q\geq 2$ is an integer. For
$n=1,2,...$, we study the probability distribution $P_{n}$ of the (scaled)
remainder $\sum_{k=n+1}^{\infty}X_{q}q^{n-k}$: If $X$ has an absolutely
continuous CDF then $P_{n}$ converges in the total variation metric to
Lebesgue measure on the unit interval; under certain smoothness conditions we
establish exponentially fast convergence of $P_{n}$ and its PDF $f_{n}$; and
we give examples of these results. The results are extended to the case of a
multivariate random variable defined on a unit cube.

Keywords:  asymptotic distribution; exponential convergence rate; extended
Newcomb-Benford law; multivariate digit expansion; remainder of a digit
expansion; total variation distance; uniform distribution

2020 Mathematics Subject Classification:  60F25; 62E17; 37A50

##  1 Introduction

Let $X$ be a random variable so that $0\leq X<1$, and for $x\in\mathbb{R}$,
let $F(x)=\mathrm{P}(X\leq x)$ be the cumulative distribution function (CDF)
of $X$. For a given integer $q\geq 2$, we consider the base-$q$ transformation
$T_{q}:[0,1)\mapsto[0,1)$ given by

$T_{q}(x)=xq-\lfloor xq\rfloor$ (1)

where $\lfloor\cdot\rfloor$ is the floor function (so $\lfloor xq\rfloor$ is
the integer part of $xq$). For $n=1,2,...$, let
$T_{q}^{n}=T_{q}\circ\cdots\circ T_{q}$ denote the composition of $T_{q}$ with
itself $n$ times and define

$X_{n}=\lfloor T_{q}^{n-1}(X)q\rfloor$ (2)

where $T_{q}^{0}(X)=X$. Then

$X=\sum_{n=1}^{\infty}X_{n}q^{-n}$ (3)

is the base-$q$ expansion of $X$ with digits $X_{1},X_{2},...$. Note that $X$
is in a one-to-one correspondence to the first $n$ digits $(X_{1},...,X_{n})$
together with $T_{q}^{n}(X)=\sum_{k=n+1}^{\infty}X_{k}q^{n-k}$, which is the
remainder multiplied by $q^{n}$. Let $\mu$ denote Lebesgue measure on $[0,1)$,
$P_{n}$ the probability distribution of $T_{q}^{n}(X)$ and $F_{n}$ its CDF, so
$X$ follows $P_{0}$ and has CDF $F_{0}=F$. The following facts are well-known
(see [3]):

  1. (a)

$P_{0}=P_{1}$ (i.e., invariance in distribution under $T_{q}$) is equivalent
to stationarity of the process $X_{1},X_{2},...$.

  2. (b)

$P_{0}=P_{1}$ and $F$ is absolutely continuous if and only if $P_{0}=\mu$.

  3. (c)

$P_{0}=\mu$ if and only if $X_{1},X_{2},...$ are independent and uniformly
distributed on $\\{0,1,...,q-1\\}$.

Items (a)-(c) together with the fact that $T_{q}$ is ergodic with respect to
$\mu$ are used in metric number theory (see [7] and the references therein) to
establish properties such as ‘for Lebesgue almost all numbers between 0 and 1,
the relative frequency of any finite combination of digits of a given length
$n$ and which occurs among the first $m>n$ digits converges to $q^{-n}$ as
$m\rightarrow\infty$’ (which is basically the definition of a normal number in
base-$q$, cf. [2]). To the best of our knowledge, less (or perhaps no)
attention has been paid to the asymptotic behaviour of the (scaled) remainder
$T_{q}^{n}(X)$ as $n\rightarrow\infty$. This paper fills this gap.

By ergodicity (see e.g. [4]), since the event
$A=\\{x\in[0,1)\,|\,\lim_{n\rightarrow\infty}F_{n}(x)=x\\}$ is
$T_{q}$-invariant, either $\mu(A)=0$ or $\mu(A)=1$, but we shall prove even
more when $F$ is absolutely continuous: In Section 2 we show that
$\lim_{n\rightarrow\infty}d_{\mathrm{TV}}(P_{n},\mu)=0$ where
$d_{\mathrm{TV}}$ is the total variation metric (as given later in (7)).
Because of the convergence of $P_{n}$ to $\mu$, in an experiment, if a
realization of $X$ is observed and the first $n$ digits are kept, and if (so
far) the only model assumption is absolutely continuity of $F$, then the
remainder rescaled by $q^{n}$ is approximately uniformly distributed when $n$
is large. Thus it may be wise to increase $n$ if it is found experimentally
that the remainder is not distributed approximately uniformly. Furthermore, in
Section 2 we study the convergence rate of $d_{\mathrm{TV}}(P_{n},\mu)$ and
other related properties, and we illustrate our results in connection to
various specific choices of $F$, including the case where $F$ follows from the
extended Newcomb-Benford law (Example 1). Finally, in Section 3 we generalize
our results to the situation where $X$ is extended to a multivariate random
variable with values in the $k$-dimensional unit cube $[0,1)^{k}$ and each of
the $k$-coordinates of $X$ is transformed by $T_{q}$.

We plan in a future paper to study the asymptotic behaviour of the remainder
in other expansions, including a certain base-$\beta$ expansion of a random
variable, namely when $q$ is replaced by $\beta=(1+\sqrt{5})/2$ (the golden
ratio) in all places above.

##  2 The asymptotic distribution of the remainder in the base-$q$ expansion

Let again the situation be as in (1)–(3). The following lemma is true in
general (i.e., without assuming $F$ is absolutely continuous). As in [3], we
define a base-$q$ fraction in $[0,1)$ to be a number of the form
$\sum_{k=1}^{n}j_{k}q^{-k}$ with
$(j_{1},...,j_{n})\in\mathbb{\\{}0,1,...,q-1\\}^{n}$ and $n\in\mathbb{N}$.

######  Lemma 2.1.

If $F$ has no jump at any base-$q$ fraction in $[0,1)$ then for every
$x\in[0,1]$,

$F_{n}(x)=\sum_{j=0}^{q^{n}-1}F(q^{-n}(j+x))-F(q^{-n}j).$ (4)

###### Proof.

Clearly, (4) holds for $x=1$, so let $0\leq x<1$. For
$j_{1},...,j_{n}\in\\{0,1,...,q-1\\}$ and $j=\sum_{i=1}^{n}j_{i}q^{n-i}$, the
event that $X_{1}=j_{1},...,X_{n}=j_{n}$, and $T_{q}^{n}(X)\leq x$ is the same
as the event that $q^{-n}j\leq X<q^{-n}(j+1)$ and $X\leq q^{-n}(j+x)$. Hence,
since $0\leq x<1$,

$F_{n}(x)=\sum_{j=0}^{q^{n}-1}\mathrm{P}(q^{-n}j\leq X\leq q^{-n}(j+x))$

whereby (4) follows since $F(x)$ has no jumps at the base-$q$ fractions. ∎

The property that $F$ has no jump at any base-$q$ fraction is of course
satisfied when $F$ is continuous.

For the remainder of this section we assume that $X$ has a probability density
function (PDF) $f$ concentrated on $(0,1)$, meaning that $F$ is absolutely
continuous with $F(x)=\int_{-\infty}^{x}f(t)\,\mathrm{d}t$ for all
$x\in\mathbb{R}$. Then, by (4), $F_{n}$ is absolutely continuous with PDF

$f_{n}(x)=q^{-n}\sum_{j=0}^{q^{n}-1}f(q^{-n}(j+x))$ (5)

for $0<x<1$.

We start by considering a special case of $f$ where convergence of $P_{n}$ is
obtained within a finite number of steps.

######  Proposition 2.2.

Suppose, for some integer $m\geq 1$, $f$ is constant on each of the open
intervals $(jq^{-m},(j+1)q^{-m})$, $j=0,1,...,q^{m}-1$. Then $P_{n}=\mu$ for
$n=m,m+1,...$.

###### Proof.

By (5) and the assertion in the proposition, $f_{m}$ is constant on $(0,1)$,
so $P_{m}=\mu$ and hence by invariance of $\mu$ under $T_{q}$, $P_{n}=\mu$ for
$n=m,m+1,...$. ∎

We need some notation for the following main results. For a real, measurable
function $g$ defined on $(0,1)$, denote its $L_{1}$\- and supremum-norm by
$\|g\|_{1}=\int_{0}^{1}|g(t)|\,\mathrm{d}t$ and
$\|g\|_{\infty}=\sup_{x\in(0,1)}|g(x)|$, respectively, and denote the
corresponding $L_{1}$-space by $L_{1}(0,1)=\\{g\,|\,\|g\|_{1}<\infty\\}$. Let
$\bar{L}_{1}(0,1)=\\{g\,|\int_{0}^{1}g(t)dt=1,\allowbreak\|g\|_{1}<\infty\\}$
be the subset of functions with finite $L_{1}$-norm and integral over $[0,1]$
equal one, and $\bar{L}_{1}^{\prime}(0,1)\subset\bar{L}_{1}(0,1)$ its subset
of differentiable functions $g$ such that $\|g^{\prime}\|_{\infty}<\infty$.
For $g\in\bar{L}_{1}^{\prime}(0,1)$, $n\in\mathbb{N}$, $j=0,1,...,q^{n}-1$,
and $0<x<1$, define $g_{n,j}^{\prime}(x)=g^{\prime}(x)$ if
$q^{-n}j<x<q^{-n}(j+1)$ and $g_{n,j}^{\prime}(x)=0$ otherwise, and define

$g_{n}(x)=q^{-n}\sum_{j=0}^{q^{n}-1}g(q^{-n}(j+x)).$ (6)

Henceforth, we also think of $f$ as an element of $\bar{L}_{1}(0,1)$. Denote
$\mathcal{B}$ the class of Borel subsets of $[0,1)$. The total variation
distance between two probability measures $\nu_{1}$ and $\nu_{2}$ defined on
$\mathcal{B}$ and with PDFs $g_{1}$ and $g_{2}$, respectively, is given by

$d_{\mathrm{TV}}(\nu_{1},\nu_{2})=\sup_{A\in\mathcal{B}}|\nu_{1}(A)-\nu_{2}(A)|=\frac{1}{2}\|g_{1}-g_{2}\|_{1},$
(7)

see e.g. Lemma 2.1 in [8].

######  Theorem 2.3.

If $f\in\bar{L}_{1}(0,1)$ and $g\in\bar{L}_{1}^{\prime}(0,1)$ then

$d_{\mathrm{TV}}(P_{n},\mu)\leq\frac{1}{2}\|f-g\|_{1}+\frac{1}{6}q^{-2n}\sum_{j=0}^{q^{n}-1}\|g_{n,j}^{\prime}\|_{\infty}\leq\frac{1}{2}\|f-g\|_{1}+\frac{1}{6}q^{-n}\|g^{\prime}\|_{\infty}.$
(8)

In particular,

$\lim_{n\rightarrow\infty}d_{\mathrm{TV}}(P_{n},\mu)=0$ (9)

and we have the following sharper convergence results. If
$f\in\bar{L}_{1}^{\prime}(0,1)$ then $P_{n}$ converges exponentially fast:

$d_{\mathrm{TV}}(P_{n},\mu)\leq\frac{1}{6}q^{-2n}\sum_{j=0}^{q^{n}-1}\|f_{n,j}^{\prime}\|_{\infty}\leq\frac{1}{6}q^{-n}\|f^{\prime}\|_{\infty}.$
(10)

If $\|f\|_{\infty}<\infty$ and $f$ is continuous except for finitely many
points, then

$|f_{n}(x)-1|\to 0\quad\mbox{uniformly for $x\in(0,1)$}.$ (11)

If $f$ is twice differentiable then we have the following improvement of (10):

$d_{\mathrm{TV}}(P_{n},\mu)=\frac{1}{8}q^{-2n}\left|\sum_{j=0}^{q^{n}-1}f^{\prime}(\xi_{nj})\right|+O(q^{-2n})\leq\frac{1}{8}q^{-n}\|f^{\prime}\|_{\infty}+O(q^{-2n})$
(12)

where $\xi_{n,j}\in(q^{-n}j,q^{-n}(j+1))$ is arbitrary.

Before proving this theorem we need the following lemma.

######  Lemma 2.4.

Let $f\in\bar{L}_{1}(0,1)$, $g\in\bar{L}_{1}^{\prime}(0,1)$. For every
$x\in(0,1)$,

$|g_{n}(x)-1|\leq
q^{-2n}\left(x^{2}-x+\frac{1}{2}\right)\sum_{j=0}^{q^{n}-1}\|g^{\prime}_{n,j}\|_{\infty}\leq\frac{1}{2}q^{-n}\|g^{\prime}\|_{\infty},$
(13)

and

$\int_{0}^{1}|f_{n}(x)-g_{n}(x)|\,\mathrm{d}x\leq\|f-g\|_{1}.$ (14)

If $g$ is twice differentiable on $(0,1)$ with
$\|g^{\prime\prime}\|_{\infty}<\infty$ then for every $x\in(0,1)$,

$g_{n}(x)-1=q^{-2n}\left(x-\frac{1}{2}\right)\sum_{j=0}^{q^{n}-1}g^{\prime}(\xi_{n,j})+O(q^{-2n}),$
(15)

where each $\xi_{n,j}\in(q^{-n}j,q^{-n}(j+1))$ is arbitrary.

######  Remark 1.

Of course, (13) and (15) hold with $g_{n}$ replaced by $f_{n}$ if $f$ is
differentiable respectively twice differentiable with
$\|f^{\prime\prime}\|_{\infty}<\infty.$ For Example 2 below it is useful to
realize that in (15), $q^{-n}\sum_{j=0}^{q^{n}-1}g^{\prime}(\xi_{n,j})$ is a
Riemann sum for the integral $\int_{0}^{1}g^{\prime}(t)\mathrm{d}t$.

###### Proof.

Let $x\in(0,1)$. From (6) we have

$\displaystyle g_{n}(x)-1$
$\displaystyle=\sum_{j=0}^{q^{n}-1}\int_{q^{-n}j}^{q^{-n(j+1)}}[g(q^{-n}(j+x))-g(t)]\,\mathrm{d}t$
$\displaystyle=\sum_{j=0}^{q^{n}-1}\int_{0}^{1}q^{-n}[g(q^{-n}(j+x))-g(q^{-n}(j+y))]\,\mathrm{d}y.$
(16)

If $g$ is differentiable on $(0,1)$ with $\|g^{\prime}\|_{\infty}<\infty$, we
get by the mean value theorem,

$|g(q^{-n}(j+x))-g(q^{-n}(j+y))|\leq\|g^{\prime}_{n,j}\|_{\infty}q^{-n}|x-y|,$

which yields the bound

$|g_{n}(x)-1|\leq
q^{-2n}\int_{0}^{1}|x-y|\,\mathrm{d}y\sum_{j=0}^{q^{n}-1}\|g^{\prime}_{n,j}\|_{\infty}=q^{-2n}\left(x^{2}-x+\frac{1}{2}\right)\sum_{j=0}^{q^{n}-1}\|g^{\prime}_{n,j}\|_{\infty}.$
(17)

Thereby (13) follows. Moreover,

$\displaystyle\int_{0}^{1}|f_{n}(x)-g_{n}(x)|\,\mathrm{d}x$
$\displaystyle\leq\sum_{j=0}^{q^{n}-1}\int_{0}^{1}q^{-n}|f(q^{-n}(j+x))-g(q^{-n}(j+x))|\,\mathrm{d}x$
$\displaystyle=\|f-g\|_{1}$ (18)

whereby (14) follows. If $g$ is twice differentiable on $(0,1)$ with
$\|g^{\prime\prime}\|_{\infty}<\infty$, the mean value theorem gives

$g(q^{-n}(j+x))-g(q^{-n}(j+y))=g^{\prime}(\xi_{x,y})q^{-n}(x-y)=g^{\prime}(\xi_{n,j})q^{-n}(x-y)+O(q^{-2n}),$

where $\xi_{x,y}\in(q^{-n}j,q^{-n}(j+1))$ depends on $x$ and $y$ and
$\xi_{n,j}\in(q^{-n}j,q^{-n}(j+1))$ is arbitrary. The second equality was
obtained by applying the mean value theorem to
$g^{\prime}(\xi_{x,y})-g^{\prime}(\xi_{n,j})$. Inserting this into (16) yields

$g_{n}(x)-1=q^{-2n}\sum_{j=0}^{q^{n}-1}g^{\prime}(\xi_{n,j})\int_{0}^{1}(x-y)\,\mathrm{d}y+O(q^{-2n}),$

which reduces to (15). ∎

We are now ready for the proof of Theorem 2.3.

###### Proof.

We have

$\displaystyle d_{\mathrm{TV}}(P_{n},\mu)$
$\displaystyle=\frac{1}{2}\int_{0}^{1}|f_{n}(x)-1|\,\mathrm{d}x$
$\displaystyle\leq\frac{1}{2}\int_{0}^{1}|f_{n}(x)-g_{n}(x)|\,\mathrm{d}x+\frac{1}{2}\int_{0}^{1}|g_{n}(x)-1|\,\mathrm{d}x$
$\displaystyle\leq\frac{1}{2}\|f-g\|_{1}+\frac{1}{6}q^{-2n}\sum_{j=0}^{q^{n}-1}\|g^{\prime}_{nj}\|_{\infty}$
(19)
$\displaystyle\leq\frac{1}{2}\|f-g\|_{1}+\frac{1}{6}q^{-n}\|g^{\prime}\|_{\infty}$

where we get the equality from (7) and the second inequality from (13), (14),
and since $\int_{0}^{1}\left(x^{2}-x+1/2\right)\,\mathrm{d}x=1/3$. Thereby (8)
is verified. Taking $n\to\infty$ in (8) and using that
$\bar{L}^{\prime}_{1}(0,1)$ is dense in $\bar{L}_{1}(0,1)$, we get (9).
Equation (10) follows from (8) by setting $g=f$.

For the proof of (11) we suppose $f$ is continuous except at
$x_{1},\ldots,x_{m}\in(0,1)$ and set $x_{0}=0$ and $x_{m+1}=1$. Let $\delta>0$
and

$\displaystyle I_{n}$
$\displaystyle=\\{j\in\\{0,1,\ldots,q^{n}-1\\}\mid\exists
i\in\\{0,1,\ldots,m+1\\}:|q^{-n}j-x_{i}|<\delta\\},$ $\displaystyle J_{n}$
$\displaystyle=\\{0,1,\ldots,q^{n}-1\\}\backslash I_{n}.$

By (16),

$\displaystyle f_{n}(x)-1$ $\displaystyle=\sum_{j\in
I_{n}}\int_{q^{-n}j}^{q^{-n}(j+1)}(f(q^{-n}(j+x))-f(t))\,\mathrm{d}t$
$\displaystyle+\sum_{j\in
J_{n}}\int_{q^{-n}j}^{q^{-n}(j+1)}(f(q^{-n}(j+x))-f(t))\,\mathrm{d}t.$ (20)

Given $\varepsilon>0$, we choose $\delta$ so that
$\delta<\varepsilon/(6(m+2)\|f\|_{\infty})$. Then, since the cardinality of
$I_{n}$ is at most $(m+2)(2q^{n}\delta+1)$, the first sum in (20) is bounded
by

$\left(\frac{2q^{n}\varepsilon}{6\|f\|_{\infty}}+m+2\right)q^{-n}\|f\|_{\infty}=\frac{\varepsilon}{3}+(m+2)q^{-n}\|f\|_{\infty}<\frac{\varepsilon}{2}$

for $n$ sufficiently large. Moreover, for $n$ large enough, the second sum in
(20) is bounded by $\varepsilon/2$ since $f$ is uniformly continuous on
$(0,1)\backslash\bigcup_{i=0}^{m+1}(x_{i}-\delta/2,x_{i}+\delta/2)$, which is
a closed set. Thus, for large enough $n$, $|f_{n}(x)-1|<\varepsilon$ which
gives (11) since $\varepsilon>0$ is arbitrary.

To prove (12) we use (15) with $g$ replaced by $f$. Then

$\int_{A}(f_{n}(t)-1)\,\mathrm{d}t=q^{-2n}\int_{A}\left(t-\frac{1}{2}\right)\,\mathrm{d}t\sum_{j=0}^{q^{n}-1}f^{\prime}(\xi_{nj})+O(q^{-2n}).$

We have

$\sup_{A\in\mathcal{B}}\bigg{|}\int_{A}\left(t-\frac{1}{2}\right)\,\mathrm{d}t\bigg{|}=\frac{1}{2}\sup_{A\in\mathcal{B}}\bigg{|}\int_{A}(2t-1)\,\mathrm{d}t\bigg{|}=\frac{1}{4}\int_{0}^{1}(2t-1)\,\mathrm{d}t=\frac{1}{8}$

where the second identity follows from (7). This gives (12). ∎

######  Example 1.

Any number $y\not=0$ can uniquely be written as $y=sq^{k}(y_{0}+y_{f})$ where
$s=s(y)\in\\{\pm 1\\}$ is the sign of $y$, $k=k(y)\in\mathbb{Z}$ determines
the decimal point of $y$ in base-$q$, $y_{0}=y_{0}(y)\in\\{1,...,q-1\\}$ is
the leading digit of $y$ in base-$q$, and $y_{f}=y_{f}(y)\in[0,1)$ is the
fractional part of $y$ in base-$q$; $y_{0}+y_{f}$ is the so-called significand
of $y$ in base-$q$. Correspondingly, consider any real-valued random variable
$Y\not=0$ (or just $\mathrm{P}(Y=0)=0$), so (almost surely)
$Y=Sq^{K}(X_{0}+X)$ where $S=s(Y)$, $K=k(Y)$, $X_{0}=y_{0}(Y)$, and
$X=y_{f}(Y)$ are random variables. Let $X_{1},X_{2},...$ be the digits of $X$
in the base-$q$ expansion, cf. (3). We call $X_{0},X_{1},X_{2},...$ the
significant digits of $Y$ in base-$q$. By definition $Y$ satisfies the
extended Newcomb-Benford law if

$\mathrm{P}(X_{0}=x_{0},...,X_{n}=x_{n})=\log_{q}\left(1+1\bigg{/}\sum_{j=0}^{n}q^{n-j}x_{j}\right)$
(21)

for $n=0,1,...$ and any $x_{0}\in\\{1,...,q-1\\}$ and
$x_{j}\in\\{0,1,...,q-1\\}$ with $1\leq j\leq n$. Equivalently, the log-
significand of $Y$ in base-$q$, $\log_{q}(X_{0}+X)$, is uniformly distributed
on $[0,1)$ (Theorem 4.2 in [1]). Then $X$ has a CDF and a PDF given by

$F(x)=\sum_{j=1}^{q-1}\log_{q}(1+x)-\log_{q}j,\quad
f(x)=\sum_{j=1}^{q-1}\frac{1}{\ln q}\frac{1}{1+x},$ (22)

for $0\leq x\leq 1$.

The extended Newcomb-Benford law applies to a wide variety of real datasets,
see [5, 1] and the references therein. The law is equivalent to appealing
scale invariance properties: Equation (21) is equivalent to that $Y$ has scale
invariant significant digits (Theorem 5.3 in [1]) or just that there exists
some $d\in\\{1,...,\allowbreak q-1\\}$ such that $\mathrm{P}(y_{0}(aY)=d)$
does not depend on $a>0$ (Theorem 5.8 in [1]). Remarkably, for any positive
random variable $Z$ which is independent of $Y$, if the extended Newcomb-
Benford law is satisfied by $Y$, it is also satisfied by $YZ$ (Theorem 8.12 in
[1]).

For the remainder of this example, suppose (21) is satisfied. Then, for
$n=1,2,...$ and $x_{n}\in\\{0,1,...,q-1\\}$,

$\mathrm{P}(X_{n}=x_{n})=\log_{q}\left(\prod_{j=1}^{q-1}\prod_{i=1}^{q^{n-1}}\left(1+\frac{1}{jq^{n}+(i-1)q+x_{n}}\right)\right)$

where $\mathrm{P}(X_{n}=x_{n})$ is a decreasing function of $x_{n}$. It is
well-known that the distribution of $X_{n}$ quickly tends to the uniform
distribution on $\\{0,1,...,q-1\\}$ as $n$ grows: The fast convergence is in
accordance with Theorem 2.3 stating that $T_{q}^{n}(X)$ converges to a uniform
distribution on $[0,1)$ and hence the first digit $X_{n}$ of $T_{q}^{n}(X)$ is
approximately uniformly distributed on $\\{0,1,...,q-1\\}$ when $n$ is large.
The left part of Figure 1 shows plots of
$\mathrm{P}(X_{n}=0)-\mathrm{P}(X_{n}=q-1)$ versus $n$ for $q=2,3,5,10$
indicating fast convergence to uniformity and that the convergence speed
increases with $q$. The right part of Figure 1 illustrates the stronger
statement in (11) that the PDF $f_{n}$ of $T_{q}^{n}(X)$ converges uniformly
to the uniform PDF.

Figure 1: The left plot shows $\mathrm{P}(X_{n}=0)-\mathrm{P}(X_{n}=q-1)$ as a
function of $n$ for various values of $q$ when $f$ is as in (22). The plot on
the right shows $f_{n}$ when $q=2$ and $n=0,\ldots,5$.

To further illustrate the fast convergence, we drew a sample of 1000
observations with CDF (22) and made a $\chi^{2}$ goodness-of-fit test for
uniformity of $X_{n}$. Considering a significance level of 0.05, the rejection
rate for 10.000 repetitions is shown in Table 1. Such a $\chi^{2}$ test can
also be used as a test for uniformity of the remainder $T_{q}^{n-1}(X)$. A
more refined test can be performed by basing the goodness-of-fit test on
$2^{k}$ combinations of the first $k$ digits $(X_{n},\ldots,X_{n+k-1})$. The
result is shown in Table 1 for $k=1,2,3$. When $n=1$ we always rejected the
hypothesis that $(X_{n},\ldots,X_{n+k-1})$ is uniformly distributed, when
$n=2$ the rejection rate decreases as $k$ grows and it is 0.067 for $k=3$, and
when $n\geq 3$ the rejection rates are close to 0.05 as expected if the
hypothesis is true. When we instead tried with a sample of 100 observations,
even when $n=1$ the test had almost no power for $k=1,2,3$.

$n$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
$k=1$ | 1.000 | 0.094 | 0.050 | 0.054 | 0.052 | 0.051 | 0.054 | 0.055
$k=2$ | 1.000 | 0.081 | 0.047 | 0.050 | 0.051 | 0.053 | 0.052 | 0.047
$k=3$ | 1.000 | 0.067 | 0.050 | 0.049 | 0.049 | 0.050 | 0.052 | 0.052
Table 1: Rejection rate for a $\chi^{2}$ goodness-of-fit test for uniformity
of $(X_{n},\ldots,X_{n+k-1})$.

######  Example 2.

To illustrate how the convergence rate in Theorem 2.3 depends on the
smoothness of $f$, let $f(t)=\alpha t^{\alpha-1}$ be a beta-density with shape
parameters $\alpha>0$ and $1$. Then, $f\in\bar{L}_{1}^{\prime}(0,1)$ if and
only if $\alpha=1$ or $\alpha\geq 2$. Of course, $P_{n}$ and $\mu$ agree if
$\alpha=1$. For $q=2$, Figure 2 shows plots of $d_{\mathrm{TV}}(P_{n},\mu)$
and $\ln(d_{\mathrm{TV}}(P_{n},\mu))$ versus $n$ when
$\alpha=0.1,0.5,1,1.5,5,10$ as well as a plot of
$\ln(\frac{1}{8}q^{-2n}\sum_{j=0}^{q^{n}-1}\|f^{\prime}_{n,j}\|_{\infty})-\ln(d_{\mathrm{TV}}(P_{n},\mu))$
(cf. (12)) versus $n$ when $\alpha=2,5,10$. For the calculation of
$d_{\mathrm{TV}}(P_{n},\mu)$ observe that $f_{n}^{\prime}$ is $<0$ if
$\alpha<1$ and $>0$ if $\alpha>1$, so $f_{n}(x_{0})=1$ for some unique
$x_{0}\in(0,1)$, and hence since $F_{n}(0)-0=F_{n}(1)-1=0$,

$d_{\mathrm{TV}}(P_{n},\mu)=\frac{1}{2}\|f_{n}-1\|_{1}=\frac{1}{2}\bigg{|}\int_{0}^{x_{0}}(f_{n}(t)-1)\mathrm{d}t\bigg{|}+\frac{1}{2}\bigg{|}\int_{x_{0}}^{1}(f_{n}(t)-1)\mathrm{d}t\bigg{|}=|F_{n}(x_{0})-x_{0}|.$

We used the Newton-Raphson procedure to find $x_{0}$ (the procedure always
converges).

The first plot in Figure 2 shows that for all values of $\alpha$,
$d_{\mathrm{TV}}(P_{n},\mu)$ goes to zero, as guaranteed by Theorem 2.3. The
second plot indicates that for $\alpha>1$, $d_{\mathrm{TV}}(P_{n},\mu)$ decays
exponentially at a rate independent of $\alpha$, while for $\alpha<1$, the
decay is also exponential, but with a slower rate. The graphs in the third
plot seem to approach zero, indicating that for $\alpha\geq 2$, the rate of
decay is indeed as given by (12), which holds since $f^{\prime\prime}$ is
bounded. In the middle plot, the decay rate also seems to be $q^{-n}$ for
$\alpha=1.5$, though this is not guaranteed by Theorem 2.3. To see why the
rate $q^{-n}$ also holds for $1<\alpha<2$, we argue as follows. In (16), (17),
and (19), we may refine to the cases $j=0$ and $j>0$ (observing that
$\|f^{\prime}_{n,j}\|_{\infty}<\infty$ when $j>0$) to obtain the following
modification of (10),

$d_{\mathrm{TV}}(P_{n},\mu)\leq
q^{-n}\left(\frac{1}{2}\|f_{n,0}\|_{\infty}+\frac{1}{6}q^{-n}\sum_{j=1}^{q^{n}-1}\|f^{\prime}_{n,j}\|_{\infty}\right).$

Furthermore, since $|f^{\prime}|$ is decreasing for $\alpha<2$,
$\sum_{j=1}^{q^{n}-1}\|f_{n,j}^{\prime}\|_{\infty}q^{-n}$ is a lower Riemann
sum for the improper Riemann integral
$\int_{0}^{1}|f^{\prime}(t)|\,\mathrm{d}t$, which exists and is finite when
$1<\alpha<2$. Consequently, for every $x\in(0,1)$,

$d_{\mathrm{TV}}(P_{n},\mu)\leq
q^{-n}\left(\frac{1}{2}\|f_{n,0}\|_{\infty}+\frac{1}{6}\|f^{\prime}\|_{1}\right).$

Figure 2: The first two figures show $d_{\mathrm{TV}}(P_{n},\mu)$ and
$\log(d_{\mathrm{TV}}(P_{n},\mu))$, respectively, as a function of $n$ for
$q=2$ and various values of $\alpha$. The last plot shows the difference
between
$\log(\frac{1}{8}q^{-2n}\sum_{j=0}^{q^{n}-1}\|f^{\prime}_{n,j}\|_{\infty})$
and $\log(d_{\mathrm{TV}}(P_{n},\mu))$ for three values of $\alpha\geq 2$.

As in Example 1, we tested for uniformity of the remainder $T_{q}^{n-1}(X)$ by
a $\chi^{2}$ goodness-of-fit test for uniform distribution of the $k=3$ first
digits $(X_{n},X_{n+1},X_{n+2})$ again using 10.000 replications of samples of
1000 observations from a beta-distribution with $\alpha=0.1,0.5,1.5,2$. Table
2 shows that for $\alpha=0.1$, uniformity is rejected in all samples for all
$n$ indicating that the distribution of the remainder remains far from uniform
even for $n=8$. For $\alpha=0.5$, the rejection rate reaches the 0.05 level
for $n=5$, while for $\alpha=1.5$, this happens already for $n=2$ and for
$\alpha=5$ it happens around $n=3$ or $n=4$. For $\alpha>1$ close to 1, the
results are comparable to those for the Benford law in Example 1, while for
large $\alpha$ and $\alpha<1$, the rejection rate is higher indicating slower
convergence.

$n$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
$\alpha=0.1$ | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
$\alpha=0.5$ | 1.000 | 1.000 | 0.416 | 0.078 | 0.052 | 0.051 | 0.048 | 0.049
$\alpha=1.5$ | 1.000 | 0.123 | 0.049 | 0.048 | 0.051 | 0.052 | 0.047 | 0.051
$\alpha=5$ | 1.000 | 0.932 | 0.059 | 0.049 | 0.048 | 0.048 | 0.049 | 0.050
Table 2: Rejection rate for a $\chi^{2}$ goodness-of-fit test for uniformity
of $(X_{n},\ldots,X_{n+2})$ in a beta-distribution for various values of
$\alpha$.

######  Remark 2.

In conclusion, Examples 1 and 2 demonstrate that the answer to the title of
our paper (‘How many digits are needed?’) of course depend much on $q$ (in
Example 1, the higher $q$ is, the less digits are needed) and on how much $f$
deviates from the uniform PDF on $[0,1)$ (in Example 2, the more skew $f$ is,
the more digits are needed). Moreover, as this deviation increases or the
sample size decreases, the $\chi^{2}$ goodness-of-fit test as used in the
examples becomes less powerful; alternative tests are discussed in [6].

##  3 The multivariate case

Theorem 2.3 extends as follows. For a given positive integer $k$, let now
$X=(X_{1},...,X_{k})$ be a $k$-dimensional random variable with values in the
unit cube $[0,1)^{k}$ so that its CDF $F(x_{1},...,x_{n})=\mathrm{P}(X_{1}\leq
x_{1},...,X_{k}\leq x_{k})$ is absolutely continuous, and denote its
multivariate PDF by $f$. Extend the function $T_{q}$ to be a function
$T_{q}:[0,1)^{k}\mapsto[0,1)^{k}$ so that
$T_{q}(x_{1},...,x_{k})=(T_{q}(x_{1}),...,T_{q}(x_{k}))$. For $n=1,2,...$,
denote the multivariate CDF of $T_{q}^{n}(X)$ by $F_{n}$. For a real Lebesgue
integrable function $g$ defined on $(0,1)^{k}$, let
$\|g\|_{1}=\int_{(0,1)^{k}}|g(t)|\,\mathrm{d}t$ and let $L_{1}((0,1)^{k})$ be
the set of such functions $g$ (i.e., $\|g\|_{1}<\infty$). For a real
$k$-dimensional function $g=(g_{1},...,g_{k})$ defined on $(0,1)^{k}$, let
$\|g\|_{\infty}=\sup_{x\in(0,1)^{k}}\sqrt{g_{1}(x)^{2}+...+g_{k}(x)^{2}}$.
Define the set $\bar{L}_{1}((0,1)^{k})=\\{g\in
L^{1}((0,1)^{k})\,|\,\int_{(0,1)^{k}}g(t)\,\mathrm{d}t=1\\}$ and
$\bar{L}_{1}^{\prime}((0,1)^{k})$ as its subset of differentiable functions
$g$ with gradient

$\nabla g(x_{1},\ldots,x_{k})=\left(\frac{\partial g}{\partial
x_{1}}(x_{1},\ldots,x_{k}),\ldots,\frac{\partial g}{\partial
x_{k}}(x_{1},\ldots,x_{k})\right)\quad$

such that $\|\nabla g\|_{\infty}<\infty$. Thus, for $k=1$,
$\bar{L}_{1}^{\prime}((0,1)^{k})=\bar{L}_{1}^{\prime}(0,1)$ as used in Theorem
2.3. For $g\in\bar{L}_{1}^{\prime}((0,1)^{k})$, $n\in\mathbb{N}$,
$j:=(j_{1},...,j_{k})\in\\{0,1,...,q-1\\}^{k}$, and
$x=(x_{1},...,x_{k})\in(0,1)^{k}$, define $\nabla g_{n,j}(x)=\nabla g(x)$ if
$q^{-n}j_{i}<x_{i}<q^{-n}(j_{i}+1)$ for $i=1,...,k$ and $\nabla g_{n,j}(x)=0$
otherwise, and define

$F_{g}(x)=\int_{0}^{x_{1}}\cdots\int_{0}^{x_{k}}g(t_{1},...,t_{k})\,\mathrm{d}t_{1}\cdots\,\mathrm{d}t_{k}.$

For notational convenience, we can consider $F$ and $F_{n}$ to be functions
defined on $(0,1)^{k}$, so $F=F_{f}$. Let $e=(1,\cdots,1)$, that is $x$ with
each component equal to $1$, and as a short hand notation write
$\sum_{j=0}^{(q^{n}-1)e}...$ for
$\sum_{j_{1}=0}^{q^{n}-1}\cdots\sum_{j_{k}=0}^{q^{n}-1}...$, and for a real
function $g$ defined on $(0,1)^{k}$, $n=1,2,...$, and $x\in(0,1)^{k}$, let

$g_{n}(x)=q^{-nk}\sum_{j=0}^{(q^{n}-1)e}g(q^{-n}(j+x)).$

Then, as in (5), we see that $F_{n}$ is absolutely continuous with PDF
$f_{n}$. Finally, let $P_{n}$ be the probability distribution with CDF
$F_{n}$, $\mu$ Lebesgue measure on $[0,1)^{k}$, and
$d_{\mathrm{TV}}(P_{n},\mu)$ the total variation distance between these
measure (where (7) extends to the multivariate case with obvious
modifications).

######  Theorem 3.1.

If $g\in\bar{L}_{1}^{\prime}((0,1)^{k})$ then

$\displaystyle d_{\mathrm{TV}}(P_{n},\mu)$
$\displaystyle\leq\frac{1}{2}\|f-g\|_{1}+\frac{1}{2}{\sqrt{\frac{k}{3}}}q^{-n(k+1)}\sum_{j=0}^{(q^{n}-1)e}\|\nabla
g_{n,j}\|_{\infty}$ (23)
$\displaystyle\leq\frac{1}{2}\|f-g\|_{1}+\frac{1}{2}{\sqrt{\frac{k}{3}}}q^{-n}\|\nabla
g\|_{\infty}.$ (24)

In particular,

$\lim_{n\rightarrow\infty}d_{\mathrm{TV}}(P_{n},\mu)=0.$

Furthermore, if $f\in\bar{L}_{1}^{\prime}((0,1)^{k})$ then $P_{n}$ converges
exponentially fast:

$d_{\mathrm{TV}}(P_{n},\mu)\leq\frac{1}{2}\sqrt{\frac{k}{3}}q^{-n(k+1)}\sum_{j=0}^{(q^{n}-1)e}\|\nabla
f_{n,j}\|_{\infty}\leq\frac{1}{2}{\sqrt{\frac{k}{3}}}q^{-n}\|\nabla
f\|_{\infty}.$

Finally, if $\|f\|_{\infty}<\infty$ and $f$ is continuous except for finitely
many points, then

$|f_{n}(x)-1|\to 0\quad\mbox{uniformly for $x\in(0,1)^{k}$}.$ (25)

###### Proof.

Let $x\in(0,1)^{k}$ and $g\in\bar{L}_{1}^{\prime}((0,1)^{k})$. As in (16),

$\displaystyle g_{n}(x)-1$
$\displaystyle=\sum_{j=0}^{(q^{n}-1)e}\int_{q^{-n}j_{1}}^{q^{-n(j_{1}+1)}}\cdots\int_{q^{-n}j_{k}}^{q^{-n(j_{k}+1)}}[g(q^{-n}(j+x))-g(t)]\,\mathrm{d}t$
$\displaystyle=\sum_{j=0}^{(q^{n}-1)e}\int_{(0,1)^{k}}q^{-nk}[g(q^{-n}(j+x))-g(q^{-n}(j+y))]\,\mathrm{d}y.$

By the mean value theorem,

$|g(q^{-n}(j+x))-g(q^{-n}(j+y))|\leq\|\nabla g_{n,j}\|_{\infty}q^{-n}\|x-y\|$

where $\|\cdot\|$ is usual Euclidean distance. We estimate

$\int_{(0,1)^{k}}\|x-y\|\,\mathrm{d}y\leq\left(\int_{(0,1)^{k}}\|x-y\|^{2}\,\mathrm{d}y\right)^{1/2}\leq\frac{\sqrt{k}}{3}$

which yields the bound

$|g_{n}(x)-1|\leq q^{-n(k+1)}\frac{\sqrt{k}}{3}\sum_{j=0}^{(q^{n}-1)e}\|\nabla
g_{n,j}\|_{\infty}\leq q^{-n}\frac{\sqrt{k}}{3}\|\nabla g\|_{\infty}.$

As in (18) we have

$\int_{0}^{1}\left|f_{n}(x)-g_{n}(x)\right|\,\mathrm{d}x\leq\|f-g\|_{1}.$

Combining the last two estimates gives

$\displaystyle 2d_{\mathrm{TV}}(P_{n},\mu)$
$\displaystyle\leq\|f-g\|_{1}+\sqrt{\frac{k}{3}}q^{-n(k+1)}\sum_{j=0}^{(q^{n}-1)e}\|\nabla
g_{n,j}\|_{\infty}$
$\displaystyle\leq\|f-g\|_{1}+\sqrt{\frac{k}{3}}q^{-n}\|\nabla g\|_{\infty}$

whereby (23) follows.

To show (25), we assume for simplicity that $f$ only has one discontinuity at
$x_{0}=(x_{0,1},\ldots,x_{0,k})\in(0,1)^{k}$. The case of more than one
discontinuity can be treated as in the proof of Theorem 2.3. Let $\delta>0$
and define

$\displaystyle I_{n}$
$\displaystyle=\\{(j_{1},\ldots,j_{k})\in\\{0,1,\ldots,q^{n}-1\\}^{k}\mid\exists
i:j_{i}<\lfloor q^{n}\delta\rfloor\vee j_{i}>\lceil q^{n}(1-\delta)\rceil\\},$
$\displaystyle J_{n}$
$\displaystyle=\\{(j_{1},\ldots,j_{k})\in\\{0,1,\ldots,q^{n}-1\\}^{k}\mid\max_{i}|j_{i}-q^{n}x_{0,i}|<\lfloor
q^{n}\delta\rfloor\\},$ $\displaystyle K_{n}$
$\displaystyle=\\{0,\ldots,q^{n}-1\\}^{k}\backslash(I_{n}\cup J_{n}).$

Then

$\displaystyle f_{n}(x)-1$ $\displaystyle=\sum_{j\in
I_{n}}\int_{(0,1)^{k}}q^{-nk}[f(q^{-n}(j+x))-f(q^{-n}(j+y))]\,\mathrm{d}y$
$\displaystyle+\sum_{j\in
J_{n}}\int_{(0,1)^{k}}q^{-nk}[f(q^{-n}(j+x))-f(q^{-n}(j+y))]\,\mathrm{d}y$
$\displaystyle+\sum_{j\in
K_{n}}\int_{(0,1)^{k}}q^{-nk}[f(q^{-n}(j+x))-f(q^{-n}(j+y))]\,\mathrm{d}y.$
(26)

If $\varepsilon>0$ is given, we can choose $\delta$ such that each term in
(26) is less than $\varepsilon/3$. This follows as in the proof of Theorem 2.3
by noting that the cardinality of $I_{n}$ is at most $\delta q^{n(k-1)}$, the
cardinality of $J_{n}$ is at most $2\delta q^{n}$, and $f$ is bounded on
$(0,1)^{k}$ and uniformly continuous on the closed set
$[\delta/2,1-\delta/2]^{k}\backslash C(x_{0})$ where $C(x_{0})$ denotes the
cube of sidelength $\delta$ centered at $x_{0}$. ∎

## Acknowledgements

Supported by The Danish Council for Independent Research — Natural Sciences,
grant DFF – 10.46540/2032-00005B.

## References

  * [1] Berger, A. and Hill, T.P. (2015). An Introduction to Benford’s Law. Princeton University Press, Princeton, NJ.
  * [2] Borel, E. (1909). Les probabilités dénombrables et leurs applications arithmétiques. Rendiconti del Circolo Matematico di Palermo 27, 247–271.
  * [3] Cornean, H., Herbst, I., Møller, J., Støttrup, B.B., Studsgaard, K.S. (2022). Characterization of random variables with stationary digits. J. Appl. Prob. 59, 931–947.
  * [4] Dajani, K., Kalle, C. (2021). A First Course in Ergodic Theory. Chapman and Hall/CRC.
  * [5] Hill, T.P. (1998). The first digit phenomenon. American Scientist. 86, 358–363.
  * [6] Morrow, J. (2014). Benford’s Law, families of distributions and a test basis. CEP Discussion Paper No 1291, Centre for Economic Performance, London School of Economics and Political Science.
  * [7] Schweiger, F. (1995). Ergodic Theory of Fibered Systems and Metric Number Theory. Clarendon Press, Oxford.
  * [8] Tysbakov, A. (2009). Introduction to Nonparametric Estimation. Springer Series in Statistics, Springer - Verlag, New York.