The Newcomb–Benford law: Scale invariance and a simple Markov process based on it

In the late 19th century, an astronomer and mathematician visits his institution's library and consults a table of logarithms to perform certain astronomical calculations. As on previous occasions, he is struck by the fact that the first pages (those corresponding to numbers that start at 1) are much more worn than the last ones (corresponding to numbers that start at 9). Intrigued, this time he decides not to overlook the matter. He closes his eyes to concentrate, sketches a few calculations on a piece of paper, and finally smiles. He has found the answer and it turns out to be enormously simple and elegant.

A little over half a century later, a physicist and electrical engineer, who is unaware of his predecessor's discovery, observes the same curious phenomenon on the pages of logarithm tables and arrives at the same conclusion. Both have understood that, in a long list of records ${rn}$ obtained from measurements or observations, the fraction p_d of records beginning with the significant digit $d=1,2,…,9$ is not $pd=1/9$⁠, as one might naively expect, but rather follows a logarithmic law. More specifically,

The numeric values of p_d are shown in the second column of Table I. We see that the records that start with 1, 2, or 3 account for around 60% of the total, while the other six digits must settle for the remaining 40%.

Our 19th century character is Simon Newcomb (Fig. 1) and he published his discovery in a modest two-page note.¹ The second character is Frank Benford (Fig. 2) and he wrote a 22-page article² in which, in addition to mathematically justifying Eq. (1), he showed its validity in the analysis of more than 20 000 first digits taken from sources as varied as river areas, populations of American cities, physical constants, atomic and molecular weights, specific heats, numbers taken from newspapers or the Reader's Digest, postal addresses, …, and the series $n−1, n$⁠, n², or $n!$⁠, among others, with $n∈[[1,100]]$⁠.

With such overwhelming evidence, it is not surprising that Eq. (1) is usually known as Benford's law (or first-digit law), even though it was found nearly sixty years earlier by Newcomb. This is one more manifestation of the so-called Stigler's law, according to which no scientific discovery is named after the person who discovered it in the first place. In fact, as Stigler himself points out,³ the law that bears his name was actually spelled out in a similar way twenty-three years earlier by the American sociologist Robert K. Merton. In order not to fall completely into Stigler's law, many authors refer to Eq. (1) as Newcomb–Benford's law and that is the convention (by means of the acronym NBL) that we will follow in this article.

While the literature on the NBL in specialized journals is vast,⁴ and several books on the topic are available,^5,6 the number of papers in general or science education journals is much scarcer.^7–9 In this paper, apart from revisiting the NBL and illustrating it with a few examples, we construct a Markov-chain model inspired by the invariance property of the NBL under the operation of a change of scale.

The remainder of the paper is organized as follows. The connection between the NBL (and some of its generalizations) with the property of scale invariance is formulated in Sec. II. This is followed by a few examples in Sec. III. The most original part is contained in Sec. IV, where our Markov-chain model is proposed and solved. Finally, the paper is closed in Sec. V with some concluding remarks. For the interested reader,  Appendices A–C contain the most technical and mathematical parts of Secs. II and IV.

Oftentimes, when one first speaks to a friend, relative, or even a colleague about the NBL, their first reaction is usually of skepticism. Why is the first digit not evenly distributed among the nine possible values? A simple argument shows that, if a robust distribution law exists, it cannot be the uniform distribution whatsoever.

Imagine a long list of river lengths, mountain heights, and country surfaces, for example. It is possible that the lengths of the rivers are in km, the heights of the mountains in m, and the surfaces of countries in km², but they could also be expressed in miles, feet, or acres, respectively. Will the distribution p_d depend on whether we use some units or others, or even if we mix them? It seems logical that the answer should be negative; that is, the p_d distribution should be (statistically) independent of the units chosen; in other words, it is expected that the p_d distribution is invariant under a change of scale. The uniform distribution $pd=19$ obviously does not follow that property of invariance. Suppose we start from a dataset in which all the values of the first digit are equally represented. If we multiply all the records in the dataset by 2, we can see that those records that started before with 1, 2, 3, and 4 then start with either 2 or 3, either 4 or 5, either 6 or 7, and either 8 or 9, respectively. On the other hand, all those that started with 5, 6, 7, 8, or 9 start now with 1. All those possibilities are schematically shown in Fig. 3. Therefore, if $pd=19$ initially, then $p1=59$ and $p2+p3=p4+p5=p6+p7=p8+p9=19$ after doubling all the records, thus destroying the initial uniformity.

Thus, the most identifying hallmark of the NBL is that it must be applied to records that have units or, as Newcomb himself writes,¹ “As natural numbers occur in nature, they are to be considered as the ratios of quantities.”

While relatively formal mathematical proofs of the NBL can be found in the literature,^10–12 in  Appendix A, we present a sketch of a simpler, physicist-style derivation of the law by imposing invariance under a change of scale.¹³ 

Equation (1) can be generalized beyond the first digit. The probability that the m first digits of a record match the ordered string $(d1,d2,…,dm)$⁠, where $d1∈{1,2,…,9}$ and $di∈{0,1,2,…,9}$ if $i≥2$⁠, is given by (see  Appendix A)

As an example, the probability that the first three digits of a record form precisely the string $(3,1,4)$ is $p3,1,4=log101+1/314=0.001 38$⁠.

Once we have $pd1,d2,…,dm$⁠, we can calculate the probability $pd(m)$ that the mth digit is d, regardless of the values of the preceding $m−1$ digits, just by summing for all possible values of those preceding $m−1$ digits,

In Table I, the law for the first digit, p_d, is accompanied by the laws for the second, third, and fourth digits, obtained from Eqs. (2) and (3). As the digit moves to the right, the probability distribution becomes less and less disparate.

In Fig. 3 we saw that, when multiplying a dataset ${rn}$ by 2, part of the records that started with $d=1,2,3,4$⁠, specifically those that start with the first two digits $(d,0), (d,1), (d,2), (d,3)$⁠, or $(d,4)$⁠, will start with 2d, while the rest will start with $2d+1$⁠. Let us call α_d the fraction of records that, initially starting with $d=1,2,3,4$⁠, start with 2d by doubling all the records. Thus,

If the dataset fulfills the NBL, then one has

We will use these values later in Sec. IV.

The applications and verifications of the NBL are numerous and cover topics as varied and prosaic as the study of the genome,¹⁴ atomic,¹⁵ nuclear,^8,16 and particle¹⁷ physics, astrophysics,¹⁸ quantum correlations,¹⁹ toxic emissions,²⁰ biophysics,²¹ medicine,²² dynamical systems,²³ distinction of chaos from noise,²⁴ statistical physics,²⁵ scientific citations,²⁶ tax audits,^5,27 electoral²⁸ or scientific²⁹ frauds, gross domestic product,³⁰ stock market,^9,31 inflation data,³² climate change,³³ world wide web,³⁴ internet traffic,³⁵ social networks,³⁶ textbook exercises,³⁷ image processing,³⁸ religious activities,³⁹ dates of birth,⁴⁰ hydrology and geology,⁴¹ fragmentation processes,⁴² the first letters of words,⁴³ or even COVID-19.⁴⁴ Other examples can be seen in the link of Ref. 45. In this section, we will present some additional examples.

Let us start with one of the situations that Benford himself studied in his classic paper:² city populations. Using data from the Spanish National Institute of Statistics,⁴⁶ we have considered the 2019 population of the 165 municipalities in the province of Badajoz (plus the total population of the province of Badajoz), of the 223 municipalities of the province of Cáceres (plus the total population of the province of Cáceres), and the total population of the 388 municipalities of the Spanish region of Extremadura, which encompasses the provinces of Badajoz and Cáceres (plus the total populations of the provinces of Badajoz and Cáceres). We have also considered the population (according to the 2016 census) of the $8 110$ Spanish municipalities. For all these lists of records, we have analyzed the frequency of those starting with $d=1,2,…9$ and the results are compared in Fig. 4. There is a good general agreement between the population data and the NBL predictions. This is interesting, since it is not obvious that the distribution of the number of inhabitants of municipalities should be invariant under a change of scale.

Let us now turn to two examples from astronomy. In the first one, we take the distance from Earth (in light-years and in parsecs) to the 300 brightest stars.⁴⁷ In the second case, the data are the daily number of sunspots from 1818 to the present.⁴⁸ As seen in Fig. 5, the distances between our planet and the 300 brightest stars agree reasonably well with the NBL, despite the fact that the list is not excessively long (only 300 data points) and that there are “local” deviations (for example, $p6<p7$ in the two choices of units). This general agreement was to be expected, since the distribution of digits in distances (which are expressed in units) is a clear example of invariance under a change of scale. However, in the case of the daily number of sunspots, quantitative (although not qualitative) differences are observed with the NBL, especially in the d = 1, 3, 4, and 5 cases. It should be noted that, although the series is rather long (more than $59 000$ records, after excluding days without data or with zero spots), each record only has one, two, or three digits (the maximum number of sunspots was 528 and corresponded to August 26, 1870), thus producing a certain bias to records starting with 1. Moreover, the daily number of sunspots cannot be expressed in units and therefore may not be scale-invariant.

Finally, we have analyzed the prices of $1 016$ items from a chain of fashion retailers⁴⁹ and of $1 373$ products from a chain of hypermarkets.⁵⁰ The results are shown in Fig. 6. In this case, the discrepancies with the NBL are more pronounced. Although the highest frequencies occur for d = 1 and d = 2, the observed values of p_d do not decrease monotonically with increasing d. In the case of the fashion retailers, we have $p4>p3$ and $p9>p8>p6>p7$⁠; in the prices of the chain of hypermarkets, $p8>p9>p7>p6$⁠. In principle, one might think that, since they can be expressed in euro, pound, dollar, peso, ruble, yen, …, prices should satisfy the property of invariance under a change of scale inherent to the NBL. However, commercial and artificial pricing strategies must be superimposed on this invariance, which generates relevant deviations with respect to the NBL.

As already said, NBL, Eq. (1), is invariant under a change of scale; that is, if we start from a dataset ${rn}$ that fulfills the NBL and multiply all the records by a constant λ (other than a fractional power of 10), the resulting dataset ${λrn}$ still fulfills the NBL.

It is tempting to conjecture that the NBL should be an attractor of this process. This would mean that, if we started from an initial dataset ${rn(0)}$ that does not fulfill the NBL and generated new sets ${rn(t)}={λtrn(0)}$ at times $t=1,2,…$ by multiplying each successive dataset by λ, the first-digit distribution ${pd(t)}$ of the generated sets would converge toward the NBL. However, this is not the case. As illustrated in Fig. 7 for the case λ = 2 and an initial dataset of 10⁴ random numbers, the time-dependent distribution ${pd(t)}$ exhibits a quasiperiodic oscillatory pattern around the NBL distribution without any apparent amplitude attenuation. In fact, since $210=1024≃103$⁠, the distribution nearly returns to the uniform initial one at times $t=10,20,30,…$⁠. This behavior is reminiscent of the Poincaré recurrence time in mechanical systems and the associated Zermelo paradox about irreversibility.⁵¹ 

In the transformation ${rn(t)}→{rn(t+1)=2rn(t)}$⁠, the first-digit distribution changes, according to Fig. 3, as

At this point, we construct a simplified dynamical model in which the four unknown and time-dependent parameters α_d in Eqs. (6) are replaced by fixed constants. In matrix notation,

where $p(t)=(p1(t),p2(t),…,p9(t))†$ is a column vector (⁠$†$ denotes transposition) and

is a fixed square matrix. Equation (7) fits the canonical description of Markov chains,⁵² where $W$ is the so-called transition matrix, and ${αd}$ correspond to transition probabilities. In this way, we may forget about the meaning of ${pd(t)}$ as the first-digit distribution of the dataset ${rn(t)}$ and look at the numerals $1,2,…,9$ as labels of nine distinct internal states of a certain physical system, which follow each other in the prescribed order sketched by Fig. 3.

Note that the matrix $W$ is singular, that is, not invertible. This implies the irreversible character of the transition ${pd(t)}→{pd(t+1)}$⁠. The stationary solution $p∗$ to Eq. (7) satisfies the condition $p∗=W·p∗$⁠. Such a stationary solution will be an attractor of our Markov process if $limt→∞p(t)=p∗$ for any initial condition $p(0)$⁠. This property, along with the exact solution of Eq. (7), is proved in  Appendix B. If we choose for α_d the values given by Eqs. (5), the stationary solution coincides exactly with the NBL, as could be expected. This will be the choice made in the rest of this section.

To illustrate the irreversible evolution of $p(t)$ toward $p∗$⁠, we are going to consider two different initial conditions. First, we start from a uniform distribution, that is, $pd(0)=19$⁠. The result is shown in Fig. 8, where we see that the evolution is oscillatory (see  Appendix B for an explanation) and the oscillations are rapidly damped to attain the stationary solution. As a second example, we take an NBL inverted initial distribution, that is, $pd(0)=p10−d∗$⁠, so that, initially, state 9 is the most probable and state 1 is the least probable. In this second example, as shown in Fig. 9, the initial oscillations are of greater amplitude but, as before, the stationary distribution is practically reached after a few iterations. Comparison between Figs. 7 and 8 shows that the main difference between our Markov model and the non-Markovian evolution of ${pd(t)}$ in a real dataset is that the oscillation amplitudes are attenuated in the model and not in the original framework.

It seems convenient to characterize the evolution of the set of probabilities ${pd(t)}$ obeying the Markov process [Eq. (7)] toward the attractor ${pd∗}$ by means of a single parameter that, in addition, evolves monotonically, thus representing the irreversibility of evolution. It is expected that these properties are verified by the Kullback–Leibler divergence,⁵³ which in our case is defined as

This quantity represents the opposite of the Shannon entropy⁵⁴ of ${pd(t)}$⁠, except that it is measured with respect to the stationary distribution ${pd∗}$⁠. In the context of our Markov-chain model, $DKL$ is related to information theory. On the other hand, the mathematical structure of $DKL$ can be extended to physical systems, thus providing a statistical-mechanical basis to the thermodynamic entropy,^54,55 as exemplified below.

Figure 10 shows the evolution of $DKL(t)$ for the same initial conditions as in Figs. 8 and 9. Both cases confirm the desired monotonic evolution of $DKL(t)$⁠. Also, the asymptotic decay $DKL(t)→0$ occurs essentially exponentially with a rate independent of the initial probability distribution. This is proved in  Appendix C, as well as the monotonicity property,

with the equality not holding for two successive times unless $pd(t)=pd∗$⁠, in which case $DKL=0$⁠.

Thus, we can say that the NBL in our Markov model plays a role analogous to the equilibrium state in thermodynamics and statistical mechanics. In this analogy, the “entropy” of the out-of-equilibrium system would be (except for a constant) $S=−DKL$⁠, so that S increases irreversibly in the evolution toward equilibrium.

This analogy is put forward in Table II, where we compare a system described by the Markov chain [Eq. (7)] to a classical dilute gas as an example of a well-known physical system. In both cases, a statistical description is constructed by introducing the adequate probability distribution: the velocity distribution function (gas) and the probability of the internal state d (Markov chain); while $f(v→,t)$ is continuous in both $v→$ and t, $pd(t)$ is discrete in d and t. The evolution equation for the probability distribution function is the Boltzmann equation (gas, under the molecular chaos ansatz⁵⁶) and Eq. (7) (Markov chain). The equilibrium state is characterized by the Maxwell–Boltzmann distribution⁵⁷ $fMB(v→)$ (gas) and the NBL distribution $pd∗$ (Markov chain). In both classes of systems, the nonequilibrium entropy functional S(t) can be defined, within a constant, as the opposite of the Kullback–Leibler divergence^53,58 of the equilibrium distribution from the time-dependent one. Boltzmann's celebrated H-theorem^56,59 (gas) and the result presented in Eq. (10) (Markov chain) show that the entropy S(t) never decreases, irreversibly evolving toward its equilibrium value.

One of the main goals of this article was to show that, contrary to what might be initially thought, the first significant digit d of a dataset extracted from measurements or observations of the real world is not evenly distributed among the nine possible values, but typically the frequency is higher for d = 1 and decreases as d increases. The NBL, Eq. (1), gives a mathematical expression to this empirical fact, although it does not always need to be rigorously verified due to statistical fluctuations (as happens with populations in Fig. 4 and with distances in Fig. 5), limitations of the sample (as in the sunspot case in Fig. 5), or an artificial bias (as in the case of prices of articles in Fig. 6). It is to be expected that, except for unavoidable statistical fluctuations, the law is fulfilled in datasets in which quantities are measured in units, so that the distribution of the first digit is independent of the units chosen due to its invariance under a change of scale. More generally, the NBL is satisfied when the mantissa of the logarithms of the considered data is uniformly distributed. That makes lists not directly related to physical quantities, such as Fibonacci numbers or powers of 2, also satisfy the NBL.

Gaining inspiration from the scale invariance property of the NBL, we have constructed a Markov-chain model for a nine-state system that evolves in time according to the scheme depicted in Fig. 3. The initial-value model can be exactly solved, the solution showing an irreversible relaxation toward the NBL probability distribution. Moreover, we have proved that the associated (relative) entropy monotonically increases, in analogy with the second law of thermodynamics in physical systems.

Until the 1970s (which is when scientific pocket calculators began to be used), physicists used tables of logarithms (or their application in slide rules) for small everyday scientific calculations. Those calculations are nowadays performed on pocket calculators, cellular phones, or personal computers with a wide variety of existing mathematical programs. Since the data that are manipulated in physics are extracted from “real” situations, such as experiments, models, physical constants, …, we can conclude, as a tribute to Newcomb and Benford and their logarithmic tables, that keyboard button 1 will be the one that presents the greatest wear and tear, while that of 9 will be the least used.

A.S. acknowledges financial support from Grant No. FIS2016-76359-P/AEI/10.13039/501100011033 and the Junta de Extremadura (Spain) through Grant No. GR18079, both partially financed by Fondo Europeo de Desarrollo Regional funds.

Consider a long list of records ${rn}$ that, without loss of generality for the matter at hand, we will assume positive. Each record can be written in the form $rn=xn×10kn$⁠, where k_n is an integer and $xn∈[1,10)$ is the significand. Obviously, it is the distribution of the significand that is relevant for the NBL. The significand x_n is directly related to the mantissa μ_n of the decimal logarithm of r_n, that is, $log10rn=kn+μn$⁠, where $μn=log10xn∈[0,1)$⁠.

Let $Px(x)dx$ be the probability that the significand lies between x and $x+dx$⁠, so that the normalization condition is $∫110dx Px(x)=1$⁠. The probability that the first significant digit of any record is d is then given by the integral,

More generally, if we consider an ordered string $(d1,d2,…,dm)$ made up of the first m digits, where $d1∈{1,2,…,9}$ and $di∈{0,1,2,…,9}$ if $i≥2$⁠, it is obvious that the records whose m first digits match the string $(d1,d2,…,dm)$ will be those whose significand x is greater than or equal to $d1+d2×10−1+⋯+dm×10−(m−1)$ and less than $d1+d2×10−1+⋯+(dm+1)×10−(m−1)$⁠. Consequently,

If the distribution $Px(x)$ is actually invariant under a change of scale, that means that $Px(λx)=f(λ)Px(x)$ with arbitrary λ. Taking into account the normalization condition in the form $∫λ10λd(λx) Px(λx)=1$⁠, it follows that necessarily $f(λ)=λ−1$⁠, that is, $Px(λx)=λ−1Px(x)$⁠. Differentiating both sides of the equation with respect to λ and then taking λ = 1, we easily obtain $xPx′(x)=−Px(x)$⁠, which, according to Euler's theorem, implies that $Px(x)$ is a homogeneous function of degree −1, that is, $Px(x)∝x−1$⁠. The constant of proportionality is obtained from the normalization condition, finally yielding

This is the unique distribution of significands that is invariant under a change of scale. From Eq. (A3), and applying Eqs. (A1) and (A2), it is straightforward to obtain NBL, Eq. (1), and its generalization, Eq. (2). As a consistency test of Eq. (2), it is easy to check that

It is interesting to note that the inverse law for the significand, Eq. (A3), implies a uniform law for the mantissa (and vice versa). Let $Pμ(μ)dμ$ be the probability that the mantissa lies between μ and $μ+dμ$⁠. Since $Pμ(μ)dμ=Px(x)dx$ and $dμ=(x−1/ ln 10)dx$⁠, Eq. (A3) gives us $Pμ(μ)=1$⁠. In Newcomb's words,¹ “The law of probability of the occurrence of numbers is such that all mantissæ of their logarithms are equally probable.” An immediate consequence is that, if μ is a random variable uniformly distributed between 0 and 1, then the random variable $x=10μ$ fulfills the NBL. This is in fact an easy way to generate a list of random records matching the NBL.

There are deterministic series that also satisfy the NBL. Suppose the series ${rn=aαn+bβn,n=1,2,…}$⁠, where a, b, α, and β are real numbers with $a>0, α>β≥0$⁠, and $log10α=irrational$⁠.¹² In that case, $limn→∞ log10rn=n log10α+log10a$ has a uniformly distributed mantissa, so ${rn}$ satisfies the NBL. That includes, for example, the series ${2n}, {3n}$⁠, and ${Fn}$⁠, where $Fn=15[φn−(−φ−1)n]$ are the Fibonacci numbers, $φ≡12(1+5)$ being the golden ratio. Similarly, the series ${n!}$ also verifies the law.¹² 

Another important property is that if a list ${rn}$ fulfills the NBL, so does the list ${rna}$⁠, a being a real number. Indeed, if $log10rn=kn+μn$⁠, the mantissa μ_n being uniformly distributed, then the mantissa of $log10rna=a(kn+μn)$ is also evenly distributed. This is directly related to the fact that the NBL is not only invariant under a change of scale but also under base change,¹¹ as would be expected, given the artificial character of the decimal base. To see it, let us assume a base b and build the list ${rna}$⁠, with $a=log10b$⁠, from a list ${rn}$ that fulfills the NBL. In that case, $rna=yn×bkn$⁠, where $yn=xna∈[1,b)$ is the significand of $rna$ in the base b. The probability distribution $Py(y)$ is related to the distribution $Px(x)$ through the equation $Py(y)dy=Px(x)dx$⁠, so that Eq. (A3) leads to

Therefore, NBL in an arbitrary base b takes the form,

Note first that Eqs. (6) verify the normalization condition $∑d=19pd(t+1)=∑d=19pd(t)=1$⁠. Therefore, only eight of the probabilities ${pd,d=1,2,…,9}$ are linearly independent, so we can eliminate one of them. If we choose $p9=1−∑d=18pd$⁠, Eq. (6a) gives us $p1(t+1)=1−p1(t)−p2(t)−p3(t)−p4(t)$⁠. Although it is not strictly necessary, it is mathematically convenient to split the column vector $p(t)$ into the subsets $pI(t)≡(p1(t),p2(t),p3(t),p4(t))†$⁠, $pII(t)≡(p5(t),p6(t),p7(t),p8(t))†$⁠, and $p9(t)$⁠. As a consequence, the matrix equation (7) can be decomposed into

where $q=(1,0,0,0)†$ and

In this way, one deals with two independent 4 × 4 matrices instead of the 9 × 9 transition matrix introduced in Eq. (8). Moreover, only the equation for the projected vector $pI$ in Eq. (B1) needs to be solved. Once the solution is obtained, the solution for the complementary projected vector $pII$ is straightforward.

By recursion, it is easy to check that the solution to the initial-value problem posed by Eq. (B1) is

The stationary solution will be an attractor if $limt→∞pI(t)=pI∗$ and $limt→∞pII(t)=pII∗$ for any initial condition $pI(0)$⁠. According to Eqs. (B3), this implies $limt→∞At=0$⁠.

To check the above condition, let us obtain the four eigenvalues ${ai,i=0,1,2,3}$ of $A$⁠. It is easy to see that the characteristic equation is $a(α1α2+a+a2+a3)=0$⁠. Therefore, the eigenvalues are $a0=0$ and

Consequently

where

From Eqs. (B5), it can be verified that $|a1|<|a2,3|<1$ for $0<α1α2<1$⁠, so that $limt→∞ait=0$ for i = 1, 2, 3. Therefore, $limt→∞Dt=0$ and hence $limt→∞At=0$⁠. This proves the character of the stationary distribution $p∗$ as an attractor of the dynamics. Moreover, since both the real eigenvalue (a₁) and the real part of the complex eigenvalues (⁠$a2,3$⁠) are negative, the evolution of each $pd(t)$ to $pd∗$ is oscillatory, as can be observed in Figs. 8 and 9. Note that, in the analysis above, the eigenvalues 0 (double), $1−α3$⁠, and α₄ of the matrix $B$ are not needed.

If we choose $αd=12$⁠, then Eqs. (B4) provide the stationary solution $p1∗=413≃0.308, p2=p3=213≃0.154, p4=p5=p6=p7=113≃0.077, p8=p9=126≃0.038$⁠. These values are not too different from those of true NBL, as can be seen from Table I. On the other hand, the most natural choice is provided by Eqs. (5), in which case the stationary solution coincides exactly with the NBL.

The Kullback–Leibler divergence is defined by Eq. (9). In order to analyze its asymptotic decay, let us consider times that are long enough so that the deviations $δpd(t)≡pd(t)−pd∗$ can be considered small. In this regime, we can expand Eq. (9) in a power series and retain the dominant terms. The result is

On the other hand, for times long enough, $|a1|t≪|a2,3|t$ (note that $|a1|=0.4261$ and $|a2,3|=0.8692$⁠), so that, according to Eqs. (B3), $δpd(t)∼|a2,3|t$⁠. Thus,

This asymptotic behavior is represented in Fig. 10.

Finally, let us prove Eq. (10). According to Eq. (9),

Evolution equations (6) and the stationarity condition $p∗=W·p∗$ can be rewritten as

Therefore,

Given the stationary values ${pd∗,d=5,…,9}$⁠, the difference $ΔDKL(t+1)$ is a function of the five probabilities ${pd(t),d=5,…,9}$⁠. To find the maximum value of $ΔDKL(t+1)$⁠, we take the derivative