⇦Chinese Physics Letters, 2019, Vol. 36, No. 7, Article code 070301 Experimental Point-to-Multipoint Plug-and-Play Measurement-Device-Independent Quantum Key Distribution Network * Guang-Zhao Tang (唐光召)**, Shi-Hai Sun (孙仕海)**, Chun-Yan Li (李春燕) Affiliations College of Science, National University of Defense Technology, Changsha 410073 Received 21 January 2019, online 20 June 2019 *Supported by the National Natural Science Foundation of China under Grant Nos 11674397 and 61671455.
**Corresponding author. Email: tang130124@163.com; shsun1983@126.com Citation Text: Tang G Z, Sun S H and Li C Y 2019 Chin. Phys. Lett. 36 070301    Abstract Measurement-device-independent quantum key distribution (MDI-QKD) offers a practical way to realize a star-type quantum network. Previous experiments on MDI-QKD networks can only support the point-to-point communication. We experimentally demonstrate a plug-and-play MDI-QKD network which can support the point-to-multipoint communication among three users. Benefiting from the plug-and-play MDI-QKD architecture, the whole network is automatically stabilized in spectrum, polarization, arrival time, and phase reference. The users only need the encoding devices, which means that the hardware requirements are greatly reduced. Our experiment shows that it is feasible to establish a point-to-multipoint MDI-QKD network. DOI:10.1088/0256-307X/36/7/070301 PACS:03.67.Dd, 03.67.Hk, 89.70.Cf © 2019 Chinese Physics Society Article Text The fundamental laws of quantum mechanics ensure the information-theoretical security of quantum key distribution (QKD) for remote users.[1,2] Until now, various QKD systems and networks have been demonstrated by tremendous experiments.[3–6] Most of the QKD networks[3–5] apply the trusted relays to form the point-to-point links between a quantum transmitter and a quantum receiver. To reduce the hardware requirement and to increase the scalability, the point-to-multipoint architecture[6] is suggested to share those expensive resources in QKD systems, especially single-photon detectors (SPDs). A critical shortcoming of these QKD networks is that the relays have to be assumed to be trustworthy. As soon as any one of the relays is dishonest, the security of the whole network will have failed. Moreover, a QKD network with trusted relays still suffers from various loopholes originating from the gap between the theoretical models and the practical devices. Among those loopholes, the detection-related loopholes are considered as the main source of attacks.[7–9] The measurement-device-independent QKD (MDI-QKD)[10–20] can not only eliminate all security loopholes on detection, but also provide a perfect star-type network architecture. In MDI-QKD, the users send the independently prepared encoding states to an untrusted measurement site (Charlie). Charlie measures the signals received and reveals the Bell states obtained. Then, the users broadcast the intensity and basis settings to generate the final secure key. The architecture of a traditional MDI-QKD network is shown in Fig. 1(a). An experiment[20] has been designed to show the feasibility of the traditional MDI-QKD network. It randomly switches any two users to the untrusted measurement site every two hours. Therefore, it still belongs to the point-to-point links. On the other hand, it requires complicate stabilization systems to calibrate the spectrum, polarization, arrival time, and phase reference of signals. It may become difficult to implement when more users are added. The plug-and-play MDI-QKD network,[16] which is shown in Fig. 1(b), contains the potential to establish a self-stabilized MDI-QKD network, which can support the point-to-multipoint communication.
cpl-36-7-070301-fig1.png
Fig. 1. MDI-QKD network. (a) The architecture of traditional MDI-QKD network. (b) The architecture of plug-and-play MDI-QKD network. Enc: encoder; Cir: circulator; FM: Faraday mirror.
In this study, we make a proof-of-principle demonstration of a self-stabilized plug-and-play MDI-QKD network which can support the point-to-multipoint communication. By encoding the signals launched by SynL1 and SynL3 subsequently, David can communicate with Alice and Bob simultaneously. The signal laser and SPDs are in the charge of Charlie. Alice, Bob, and David only need the encoding devices. Benefiting from the plug-and-play MDI-QKD architecture, the whole network is self-stabilized in spectrum, polarization, arrival time, and phase reference. Our network simplifies the stabilization systems and presents a practical way to realize a point-to-multipoint star-type MDI-QKD network. We make a proof-of-principle experiment to demonstrate the feasibility of the plug-and-play MDI-QKD network. The experimental setup is shown in Fig. 2(a). Charlie controls a signal laser (1550 nm), three synchronization lasers (SynL,1310 nm), and four SPDs. The whole system is working at the repetition rate of 1 MHz. The homemade signal laser is internally modulated into a pulse train with a width of 2 ns (FWHM). We use an asymmetric Mach–Zehnder interferometer (AMZI) to separate the pulses of the signal laser into two time bins with 25 ns time delay. Alice, Bob, and David only have the modulation devices including PMs and AMs for encoding. For the $X$ basis, the key bit is encoded into the relative phase, 0 or $\pi$, by PM1. For the $Z$ basis, the key bit is encoded into the time bin, 0 or 1, by AM1. PM2 is used for the phase randomization. We use a sawtooth wave with a repetition rate of 15 kHz to make the global phase of each optical pulse randomize in the range of $[0, 2\pi]$.
cpl-36-7-070301-fig2.png
Fig. 2. (a) Experimental setup of the plug-and-play MDI-QKD network. Charlie controls all the lasers, including a signal laser (1550 nm) and three synchronization lasers (1310 nm). The optical switch unit is used to assign the signals. When the pulses of the SynLs are reflected back by the Faraday mirror (1310 nm), a homemade photoelectric detector (PD) is utilized to detect them and generate the driving clock for the signal laser. Alice, Bob, and David encode the signals and send them back to Charlie. Charlie uses four single-photon detectors (SPDs) to realize the Bell-state measurement. ConSys: control system; Drsys: driving system; SynL: synchronization laser; Attn: attenuator; PC: polarization controller; AMZI: asymmetric Mach-Zehnder interferometer; OSU: optical switch unit; WDM: wavelength division multiplexer; PM: phase modulator; AM: amplitude modulator; BS: beam splitter; PBS: polarizing beam splitter; CIR: circulator; FR: Faraday rotator; SPDs: single-photon detectors. (b) Schematic of the driving system. PD: photoelectric detector. (c) Schematic of the optical switch unit. (d) Schematic of the optical switch.[21]
Our network can support communication for any two users. For example, if Alice wants to communicate with Bob, we can use SynL1 and SynL3 to generate the signals of Bob and Alice, which has been demonstrated in experiment.[18] Here we put emphasis on the point-to-multipoint communication among three users, namely, David communicates with Alice and Bob simultaneously. David needs to encode the signals launched by SynL1 and SynL3 independently. Alice and Bob need to encode the signals launched by SynL2. We use the setup shown in Fig. 2(b) to generate the driving clock of the signal laser and the optical switch unit. The synchronization pulses of SynL1 and SynL3 are sent from Charlie to Alice and Bob correspondingly. After being reflected back by a Faraday mirror (1310 nm), they are detected by the photoelectric detectors (PD). The output of PD4 is used to drive the signal laser (1550 nm) to generate the signal pulses of David. Similarly, the synchronization pulses of SynL2 are sent from Charlie to David, reflected back by the Faraday mirror (1310 nm), and detected by the PDs. The output of PD4 is used to drive the signal laser (1550 nm) to generate the signal pulses of Alice and Bob. Figure 2(c) shows the diagram of the optical switch unit (OSU), which is used to assign the signals to the correct users. We use the setup shown in Fig. 2(d) to serve as the optical switch (OS).[21] By modulating the relative phase shift ($\phi$) between the two subsequent passages of the light (s-path and l-path) through the phase modulator, we can realize an arbitrary optical splitting ratio between output1 and output2. The extinction ratio can exceed 20 dB with an outstanding feature of stabilization.[21] We modulate $\phi$ to 0 (or $\pi$), the optical pulses output at output1 (or output2). In our network, all users share the signal laser and the AMZI, which means that there is naturally no mismatch in spectral, in pulse waveform, and in phase-reference-frame at all. As to the polarization mode, the plug-and-play architecture can automatically compensate for the birefringence effects. The temporal mode difference between David and Alice, David and Bob can be expressed as $$\begin{alignat}{1} \Delta t^{\rm AD}=\,&(t_{1310}^{\rm C\rightleftarrows D}+t_{1550}^{\rm C\rightleftarrows A})-(t_{1310}^{\rm C\rightleftarrows A}+t_{1550}^{\rm C\rightleftarrows D}) \\ =\,&\Delta t_{0}^{\rm AD}+(1/v_{1550}-1/v_{1310})\Delta L^{\rm AD},~~ \tag {1} \end{alignat} $$ $$\begin{alignat}{1} \Delta t^{\rm BD}=\,&(t_{1310}^{\rm C\rightleftarrows D}+t_{1550}^{\rm C\rightleftarrows B})-(t_{1310}^{\rm C\rightleftarrows B}+t_{1550}^{\rm C\rightleftarrows D}) \\ =\,&\Delta t_{0}^{\rm BD}+(1/v_{1550}-1/v_{1310})\Delta L^{\rm BD},~~ \tag {2} \end{alignat} $$ where $\Delta t_{0}^{\rm AD(BD)}=(1/v_{1550}-1/v_{1310})(L_{0}^{\rm C\rightleftarrows D}-L_{0}^{\rm C\rightleftarrows A(B)})$, and $\Delta L^{\rm A(B)D}=\Delta L^{\rm C\rightleftarrows D}-\Delta L^{\rm C\rightleftarrows A(B)}$, $L^{\rm C\rightarrow A(B,D)}$ represents the fiber length between Charlie and Alice (Bob, David), and $\Delta L=\alpha_{\rm T}L^{0}\Delta T$ with $\alpha_{\rm T}=5.4\times10^{-7}/^{\circ}\!$C being the thermal expansion coefficient of fiber, and $\Delta T$ representing the change of temperature. The second part of Eqs. (1) and (2) are negligible.[18] Therefore, the arrival time differences of signals between David and Alice, David and Bob are constants, which can be compensated by adjusting the time delay between the SynLs with high-resolution (10 ps) delay chips. Our network is automatically stabilized in spectrum, polarization, arrival time, and phase reference, while the previous experiment[20] needs the polarization stabilization system, the phase stabilization system, and the wavelength calibration system to realize the high quality Hong–Ou–Mandel (HOM) interference. Thus the stabilization systems are greatly simplified in our experiment. At the measurement site, four commercial InGaAs SPDs (two id201 of id Quantique; two single-photon detectors of QuantumCTek) are used to record the coincidences. For SPDs of id201, we chose the efficiency to be $10\%$ with a gate width of 2.5 ns. The dead time is $10\,µ$s with a dark count rate of $6\times10^{-6}$ per gate. As to the SPDs of QuantumCTek, we adjust the parameters (efficiency, gate width, and dead time) to make the performance in accordance with that of the SPDs of id201. We modulate the signal pulses of Alice and David (Bob and David) into three different intensities, namely, signal state intensity ($\mu=0.3$ for Alice and David, $\mu=0.4$ for Bob and David), decoy state intensity ($\nu=0.1$) and vacuum state intensity ($\omega=0.01$). The experimental gains are listed in Tables 1 and 3, and the quantum bit error rates (QBER) are listed in Tables 2 and 4. The QBERs of $Z$-basis are due to the extinction ratio of AM1 and the background counts (Rayleigh backscattering), while in $X$-basis, the vacuum and multiphoton components of weak coherent states cause accidental coincidences which introduce a typical error rate of $25\%$.
Table 1. Experimental values of gains $Q_{I_{\rm A}I_{\rm D}}^{Z(X)}$($10^{-5}$). Here $I_{\rm A}$ and $I_{\rm D}$ are the optical intensities of Alice and David.
$Z$-basis $X$-basis
$I_{\rm A}$
$I_{\rm D}$ $\mu$ $\nu$ $\omega$ $\mu$ $\nu$ $\omega$
$\mu$ 0.892 0.2971 0.04314 1.372 0.9287 0.7347
$\nu$ 0.2894 0.0993 0.0143 0.6725 0.2706 0.1068
$\omega$ 0.0413 0.0134 0.0015 0.325 0.0535 0.00455
Table 2. Experimental values of QBERs. Here $I_{\rm A}$ and $I_{\rm D}$ are the optical intensities of Alice and David.
$Z$-basis $X$-basis
$I_{\rm A}$
$I_{\rm D}$ $\mu$ $\nu$ $\omega$ $\mu$ $\nu$ $\omega$
$\mu$ 0.0138 0.0310 0.1730 0.267 0.3626 0.4740
$\nu$ 0.0233 0.0329 0.1556 0.2905 0.2954 0.4128
$\omega$ 0.1014 0.1011 0.1222 0.4439 0.3885 0.3878
Table 3. Experimental values of gains $Q_{I_{\rm B}I_{\rm D}}^{Z(X)}$ ($10^{-5}$). Here $I_{\rm B}$ and $I_{\rm D}$ are the optical intensities of Bob and David.
$Z$-basis $X$-basis
$I_{\rm B}$
$I_{\rm D}$ $\mu$ $\nu$ $\omega$ $\mu$ $\nu$ $\omega$
$\mu$ 1.242 0.3605 0.0572 1.925 1.037 1.139
$\nu$ 0.3742 0.1314 0.0195 0.7739 0.2592 0.1315
$\omega$ 0.0527 0.0171 0.0023 0.428 0.0757 0.0090
Table 4. Experimental values of QBERs. Here $I_{\rm B}$ and $I_{\rm D}$ are the optical intensities of Bob and David.
$Z$-basis $X$-basis
$I_{\rm B}$
$I_{\rm D}$ $\mu$ $\nu$ $\omega$ $\mu$ $\nu$ $\omega$
$\mu$ 0.0151 0.0384 0.1628 0.277 0.339 0.462
$\nu$ 0.0287 0.0450 0.1569 0.294 0.286 0.388
$\omega$ 0.1412 0.1422 0.1957 0.460 0.382 0.420
In the asymptotic case, the secure key rate is given by[10,22] $$ R\geq q\{Q_{\mu\mu,11}^{Z,L}[1-H(e_{11}^{X,U})]-Q_{\mu\mu}^{Z}fH(E_{\mu\mu}^{Z})\},~~ \tag {3} $$ where $q$ is the proportion of pulses when both Alice and Bob send out signal states in the $Z$ basis, $Q_{\mu\mu}^{Z}$ and $E_{\mu\mu}^{Z}$ are the overall gain and error rate when Alice and Bob send the signal states in the $Z$ basis, and they can be directly measured from the experiment: $Q_{\mu\mu,11}^{Z,L}=\mu^{2}e^{-2\mu}Y_{11}^{Z,L}$, where $Y_{11}^{Z,L}$ is a lower bound of the yield of single photon states in the $Z$ basis, $e_{11}^{X,U}$ is an upper bound of the QBER of the single photon states in the $X$ basis, $f=1.16$ is the error correction efficiency, and $H(e)=-e\log_{2}e-(1-e)\log_{2}(1-e)$ is the binary Shannon entropy function. Here $Y_{11}^{Z,L}$ and $e_{11}^{X,U}$ can be estimated using an analytical method with two decoy states according to Ref.  [22], and $Y_{11}^{Z,L} $ is given by $$\begin{align} Y_{11}^{Z,L}=\,&[(\mu_{a}^{2}-\omega_{a}^{2})(\mu_{b}-\omega_{b})Q_{M1}^{Z} -(\nu_{a}^{2}-\omega_{a}^{2})(\nu_{b}\\ &-\omega_{b})Q_{M2}^{Z}]/[(\mu_{a}-\omega_{a})(\mu_{b}-\omega_{b}) (\nu_{a}\\ &-\omega_{a})(\nu_{b}-\omega_{b})(\mu_{a}-\nu_{a})],~~ \tag {4} \end{align} $$ where $Q_{M1}^{Z}=Q_{\nu_{a}\nu_{b}}^{Z}e^{(\nu_{a}+\nu_{b})} +Q_{\omega_{a}\omega_{b}}^{Z}e^{(\omega_{a}+\omega_{b})} -Q_{\nu_{a}\omega_{b}}^{Z}e^{(\nu_{a}+\omega_{b})} -Q_{\omega_{a}\nu_{b}}^{Z}e^{(\omega_{a}+\nu_{b})}$, and $Q_{M2}^{Z}=Q_{\mu_{a}\mu_{b}}^{Z}e^{(\mu_{a}+\mu_{b})} +Q_{\omega_{a}\omega_{b}}^{Z}e^{(\omega_{a}+\omega_{b})} -Q_{\mu_{a}\omega_{b}}^{Z}e^{(\mu_{a}+\omega_{b})} -Q_{\omega_{a}\mu_{b}}^{Z}e^{(\omega_{a}+\mu_{b})}$. We can obtain $e_{11}^{X,U}$ as follows: $$\begin{align} e_{11}^{X,U}=\,&\frac{1}{(\nu_{a}-\omega_{a})(\nu_{b}-\omega_{b})Y_{11}^{X,L}} \\ &\cdot [Q_{\nu_{a}\nu_{b}}^{X}E_{\nu_{a}\nu_{b}}^{X}e^{(\nu_{a}+\nu_{b})} +Q_{\omega_{a}\omega_{b}}^{X}E_{\omega_{a}\omega_{b}}^{X}e^{(\omega_{a}+\omega_{b})} \\ &-Q_{\nu_{a}\omega_{b}}^{X}E_{\nu_{a}\omega_{b}}^{X}e^{(\nu_{a}+\omega_{b})} -Q_{\omega_{a}\nu_{b}}^{X}E_{\omega_{a}\nu_{b}}^{X}e^{(\omega_{a}+\nu_{b})}],~~ \tag {5} \end{align} $$ where $Y_{11}^{X,L}$ can be achieved with a similar method for $Y_{11}^{Z,L}$. The parameters we estimated are listed in Tables 5 and 6.
Table 5. Parameters estimated in the process of secure key rate estimation: $Q_{M1(2),A}^{Z(X)}$ ($10^{-5}$), and $Y_{11,A}^{Z(X),L}$ ($10^{-4}$).
$Q_{M1}^{Z(X)}$ $Q_{M2}^{Z(X)}$ $Y_{11}^{Z(X),L}$
$Z$-basis 0.0919 1.512 0.734
$X$-basis 0.1563 1.060 2.275
Table 6. Parameters estimated in the process of secure key rate estimation: $Q_{M1(2),B}^{Z(X)}$ ($10^{-5}$), and $Y_{11,B}^{Z(X),L}$ ($10^{-4}$).
$Q_{M1}^{Z(X)}$ $Q_{M2}^{Z(X)}$ $Y_{11}^{Z(X),L}$
$Z$-basis 0.1219 1.690 1.650
$X$-basis 0.0945 1.932 1.136
About $1.2\times10^{11}$ signal pulses are sent out in our experiment. We chose $q=\frac{1}{36}$ for the three optical intensity states which are prepared with the same probability. Therefore, we obtain that $Y_{11,AD}^{Z,L}=7.34\times10^{-5}$ and $e_{11,AD}^{X,U}=14.65\%$, $Y_{11,BD}^{Z,L}=1.65\times10^{-4}$ and $e_{11,BD}^{X,U}=5.54\%$. The secure key rate for Alice and David can reach $1\times10^{-8}$ bits per pulse, and the secure key rate for Bob and David is $1.83\times10^{-7}$ bits per pulse. We remark that, according to the security analysis in Ref.  [21] the energy and arrival time of signal pulses should be monitored precisely to acquire the certain information about the photon-number distribution and the timing mode in plug-and-play MDI-QKD network. In our proof-of-principle demonstration, however, we cut down the energy of signal pulses to reduce the Rayleigh backscattering (RBS). The intensity detector cannot monitor the energy of signal pulses successfully with such a low pulse energy level. However, it does not affect the feasibility of our plug-and-play MDI-QKD network, since we can adopt the scheme of pulse trains to reduce the RBS.[23] The reduction of the RBS is demonstrated with an arbitrary waveform generator (Tektronix AWG7082C) worked at 1 MHz repetition rate in gated mode. Another waveform generator (KEYSIGHT 33250A) is used to offer the gate signals (5 kHz with a width of 5$\,µ$s). The background counts can reach $10^{-6}$ per gate. On the other hand, the performance can be significantly improved by optimizing the parameters and applying the high-speed SPDs. It has been demonstrated that the transmission distance and secure key rate are substantially improved using the superconducting nanowire single photon detectors (SNSPD) with an efficiency over $40\%$.[14] In conclusion, we have made a proof-of-principle demonstration of a self-stabilized plug-and-play MDI-QKD network over an asymmetric channel setting. In our network, David can communicate with Alice and Bob simultaneously by encoding the signals launched by SynL1 and SynL3 subsequently. Our experiment presents a practical way to realize a star-type MDI-QKD network which greatly reduces the hardware requirements for users. References Quantum cryptography based on Bell’s theoremThe SECOQC quantum key distribution network in ViennaMetropolitan all-pass and inter-city quantum communication networkField test of quantum key distribution in the Tokyo QKD NetworkA quantum access networkExperimental demonstration of phase-remapping attack in a practical quantum key distribution systemHacking commercial quantum cryptography systems by tailored bright illuminationPassive Faraday-mirror attack in a practical two-way quantum-key-distribution systemMeasurement-Device-Independent Quantum Key DistributionReal-World Two-Photon Interference and Proof-of-Principle Quantum Key Distribution Immune to Detector AttacksExperimental Measurement-Device-Independent Quantum Key DistributionExperimental Demonstration of Polarization Encoding Measurement-Device-Independent Quantum Key DistributionMeasurement-Device-Independent Quantum Key Distribution over 200 kmPhase-Reference-Free Experiment of Measurement-Device-Independent Quantum Key DistributionMeasurement-device-independent quantum communication with an untrusted sourceQuantum key distribution without detector vulnerabilities using optically seeded lasersExperimental asymmetric plug-and-play measurement-device-independent quantum key distributionTime-Bin Phase-Encoding Measurement-Device-Independent Quantum Key Distribution with Four Single-Photon DetectorsMeasurement-Device-Independent Quantum Key Distribution over Untrustful Metropolitan NetworkProof-of-concept of real-world quantum key distribution with quantum framesPractical aspects of measurement-device-independent quantum key distributionQuantum key distribution based on phase encoding in long-distance communication fiber
[1]Bennett C H and Brassard G 1984 Int. Conf. Comput. Syst. Signal Process (Bangalore India) (New York: IEEE) p 175
[2] Ekert A K 1991 Phys. Rev. Lett. 67 661
[3] Peev M et al 2009 New J. Phys. 11 075001
[4] Chen T Y et al 2010 Opt. Express 18 27217
[5] Sasaki M et al 2011 Opt. Express 19 10387
[6] Fröhlich B et al 2013 Nature 501 69
[7] Xu F et al 2010 New J. Phys. 12 113026
[8] Lydersen L et al 2010 Nat. Photon. 4 686
[9] Sun S H et al 2011 Phys. Rev. A 83 062331
[10] Lo H K et al 2012 Phys. Rev. Lett. 108 130503
[11] Rubenok A et al 2013 Phys. Rev. Lett. 111 130501
[12] Liu Y et al 2013 Phys. Rev. Lett. 111 130502
[13] Tang Z Y et al 2014 Phys. Rev. Lett. 112 190503
[14] Tang Y L et al 2014 Phys. Rev. Lett. 113 190501
[15] Wang C, Song X T, Yin Z Q, Wang S, Chen W, Zhang C M, Guo G C and Han Z F 2015 Phys. Rev. Lett. 115 160502
[16] Xu F 2015 Phys. Rev. A 92 012333
[17] Comandar L C, Lucamarini M, Fröhlich B, Dynes J F et al 2016 Nat. Photon. 10 312
[18] Tang G Z, Sun S H, Xu F, Chen H, Li C Y and Liang L M 2016 Phys. Rev. A 94 032326
[19] Tang G Z, Sun S H, Chen H, Li C Y and Liang L M 2016 Chin. Phys. Lett. 33 120301
[20] Tang Y L, Yin H L, Zhao Q, Liu H et al 2016 Phys. Rev. X 6 011024
[21] Lucio-Martinez I, Chan P, Mo X, Hosier S and Tittel W 2009 New J. Phys. 11 095001
[22] Xu F, Curty M, Qi B and Lo H K 2013 New J. Phys. 15 113007
[23] Sun S H, Ma H Q, Han J J, Liang L M and Li C Z 2010 Opt. Lett. 35 1203
⇦Chinese Physics Letters, 2019, Vol. 36, No. 7, Article code 070201 A Proof of First Digit Law from Laplace Transform * Mingshu Cong (丛明舒)1, Bo-Qiang Ma (马伯强)1,2** Affiliations 1School of Physics and State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871 2Center for High Energy Physics, Peking University, Beijing 100871 Received 22 March 2019, online 20 June 2019 *Supported by the National Natural Science Foundation of China under Grant No 11475006.
**Corresponding author. Email: mabq@pku.edu.cn Citation Text: Cong M S and Ma B Q 2019 Chin. Phys. Lett. 36 070201    Abstract The first digit law, also known as Benford's law or the significant digit law, is an empirical phenomenon that the leading digit of numbers from real world sources favors small ones in a form $\log(1+{1}/{d})$, where $d=1, 2,\ldots, 9$. Such a law has been elusive for over 100 years because it has been obscure whether this law is due to the logical consequence of the number system or some mysterious mechanism of nature. We provide a simple and elegant proof of this law from the application of the Laplace transform, which is an important tool of mathematical methods in physics. It is revealed that the first digit law originates from the basic property of the number system, thus it should be attributed as a basic mathematical knowledge for wide applications. DOI:10.1088/0256-307X/36/7/070201 PACS:02.30.Uu, 02.50.-r, 02.50.Cw © 2019 Chinese Physics Society Article Text The first digit law, which is also called the significant digit law or Benford's law, was first noticed by Newcomb in 1881,[1] and then re-discovered independently by Benford in 1938.[2] It is an empirical observation that the first digits of natural numbers are preferentially small rather than a uniform distribution as might be expected. More accurately, the probability that a number begins with digit $d$, where $d=1,2,\ldots,9$, can be expressed as $$ P_d=\log\Big(1+\frac{1}{d}\Big),~~d=1, 2,\ldots, 9,~~ \tag {1} $$ as shown in Fig. 1.
cpl-36-7-070201-fig1.png
Fig. 1. Benford's law of the first digit distribution, from which we can see that the probability of finding numbers with leading digit 1 is larger than that with $2,\ldots, 9$, respectively.
Empirically, the populations of countries, the areas of lakes, the lengths of rivers, the Arabic numbers on the front page of a newspaper,[2] physical constants,[3] the stock market indices,[4] file sizes in a personal computer,[5] etc., all conform to the peculiar law well. Benford's law has been verified to hold true for a vast number of examples in various domains, such as economics,[4] social science,[6] environmental science,[7] biology,[8] geology,[9] astronomy,[10] statistical physics,[11,12] nuclear physics,[13–17] particle physics,[18] and some dynamical systems.[19–21] Also, there have been many explorations on applications of the law in various fields, mainly to detect data and judge their reasonableness, such as in distinguishing and ascertaining fraud in taxing and accounting[22–24] and falsified data in scientific experiments.[25] Benford's law has several elegant properties. It is scale-invariant,[26,27] which means that the law does not depend on any particular choice of units. This law is also base-invariant,[28–30] which means that it is independent of the base $b$ with a general form $$ P_d=\log_{b}\Big(1+\frac{1}{d}\Big),~d=1, 2,\ldots, {b-1},~{\rm for}~b \ge 2.~~ \tag {2} $$ The law is also found to be power-invariant,[18] i.e., any power ($\neq 0$) on numbers in the data set does not change the first digit distribution. Although there have been many studies on Benford's law,[31] the underlying reason for the success of this law has remained elusive for more than 100 years. It was unclear whether Benford's Law was due to some unknown mechanism of nature or whether it is merely a logical consequence of the human number system. However, the situation has changed due to the appearance of a general derivation of Benford's law from the application of the Laplace transform,[32] where a strict version of Benford's law is derived as composed of a Benford term and an err term. The Benford term explains the prevalence of Benford's law and the err term leads to derivations from the law with four categories of number sets. It is the purpose of this work to provide a simpler and more elegant version of the derivation of Benford's law compared to Ref.  [32]. Through this derivation, it is easier to understand the rationality of Benford's law. We reveal that the first digit law can be derived as the main term from the Laplace transform. This explains why Benford's law is so successful for many number sets. We perform similar analyses on the regularities of the second digit and $i$th-significant digit distributions, and extend the law to a more general rule for the first several digit distribution. We also estimate the error term and point out conditions for the validity of this law. For simplification, we constrain ourselves to the decimal system first. Let $F(x)$ be an arbitrary normalized probability density function defined on positive real number set $\mathbb{R}^+$. Here we use the capital letter $F$ instead of the lowercase one, as opposed to the convention. Of course, in the real case, the variable $x$ may be negative or bounded, but this is not harmful to our derivation. When $x$ can be negative, we can use the probability density function of its absolute value, keeping results unchanged. In the decimal system, the probability $P_d$ of finding a number whose first digit is $d$ is the sum of the probability that it is contained in the interval $[d \cdot 10^n, (d+1) \cdot 10^n)$ for all integer $n$, therefore $P_d$ can be expressed as $$ P_d=\sum_{n=-\infty}^{\infty} \int_{d \cdot 10^n}^{(d+1) \cdot 10^n} F(x) dx,~~ \tag {3} $$ which can also be rewritten as $$ P_d=\int_{0}^{\infty} F(x)g_d(x) dx,~~ \tag {4} $$ with $g_d(x)$ being a new density function whose meaning will be clear in the following. Here the lowercase letter is used, due to conventions for Laplace transform in the following. Adopting the Heaviside step function $$ \eta(x)= \left\{\begin{matrix} 1, & {\rm if~} x \ge 0, \\ 0, & {\rm if~} x < 0, \end{matrix} \right.,~~ \tag {5} $$ we can write $g_d(x)$ as $$ g_d(x)=\sum_{n=-\infty}^{\infty}[\eta (x-d \cdot 10^n)- \eta (x-(d+1) \cdot 10^n)].~~ \tag {6} $$ Based on the above discussion, we can understand to some extent why numbers prefer smaller first digits. Naively one might think that the nine digits in the decimal system play the same roles, but they define different density $g_d(x)$ as shown above, thus behave differently in the decimal system. For better illustration, we draw the images of $g_1(x)$ and $g_2(x)$ in the interval $[1,30)$, as shown in Fig. 2, from which we notice that the two density functions have totally different shapes. Neither of them can simply be a translation or an expansion of the other. All the above derivations are rigorous. In fact, using Eq. (3) or (4), we can calculate $P_d$ for any given $F(x)$ numerically. Usually, it does not strictly fit in with Eq. (1). In this sense, Benford's law is not a rigorous 'law' with strong predictive power. However, using the technique of Laplace transform, we show in the following that Benford's law is a rather good approximation for those well-behaved probability density functions.
cpl-36-7-070201-fig2.png
Fig. 2. Images of $g_1(x)$ (upper) and $g_2(x)$ (lower), from which we notice that the gap indicated by the red line in $g_2(x)$ is wider than that in $g_1(x)$. This shows that the distribution of $g_1(x)$ is denser than $g_2(x)$ in the whole number range.
We now prove that if a probability density function has an inverse Laplace transform, it satisfies Benford's law well. Recalling the complex inversion formula,[33] if $F(x)$, extended to the complex plane, satisfies: (1) $F(x)$ is analytic on $\mathbb{C}$ except for a finite number of isolated singularities, (2) $F(x)$ is analytic on the half plane $\{x|{\rm Re}z>0\}$, (3) there are positive constants $M$, $R$, and $\beta$ such that $|F(x)|\le M/|x|^{\beta}$ whenever $|z|\geq R$, $F(x)$ will have an inverse Laplace transform. We can say a probability density function to be 'well-behaved' if it satisfies these three conditions and its inverse Laplace transform is smooth enough, i.e., without violent oscillation. Exponential functions, some fractional functions, and a handful of other common functions are all well-behaved. Thus the derivation in the following has wide application. In what follows, we assume that $F(x)$ is well-behaved. Let $f(t)$ be the inverse Laplace transform of $F(x)$, and $G(t)$ be the Laplace transform of $g(x)$, i.e., $$\begin{align} F(x)=\,&\int_{0}^{\infty} f(t)e^{-tx} dt,~~ \tag {7} \end{align} $$ $$\begin{align} G(t)=\,&\int_{0}^{\infty} g(x)e^{-tx} dx.~~ \tag {8} \end{align} $$ Laplace transforms have the following property $$\begin{align} \int_{0}^{\infty} F(x)g(x) dx =\,& \int_{0}^{\infty}dx g(x) \int_{0}^{\infty} f(t)e^{-tx} dt \\ =\,&\int_{0}^{\infty}dt f(t) \int_{0}^{\infty} g(x)e^{-tx} dx \\ =\,& \int_{0}^{\infty} f(t)G(t) dt,~~ \tag {9} \end{align} $$ which means that Laplace transform can act on either the function $f$ or $g$ with the above integral result keeping unchanged. To derive the left-hand side of the above equation, we would like to calculate the right-hand side instead. Because it is comparably convenient to calculate the Laplace transform of $g_d(x)$, $$\begin{alignat}{1} G_d(t)=\,&\int_{0}^{\infty} g_d(x)e^{-tx} dx \\ =\,&\sum_{n=-\infty}^{\infty} \int_{d \cdot 10^n}^{(d+1) \cdot 10^n}e^{-tx} dx \\ =\,&\frac{1}{t}\sum_{n=-\infty}^{\infty}(e^{-td \cdot 10^n}-e^{-t(d+1) \cdot 10^n}),~~ \tag {10} \end{alignat} $$ which can be treated as a function of two variables $d$ and $t$. Although $d$ is defined on the decimal digit set ${1,2,\ldots,9}$, it can be extended to the whole real axis. Therefore, $G_d(t)$ is a continuous function of $d$ as well as $t$. A technique to evaluate $G_d(t)$ is to calculate its partial derivative with respect to $d$ approximately, and then integrate the partial derivative to derive the result $$\begin{align} \frac{\partial G_d(t)}{\partial d} =\,&\sum_{n=-\infty}^{\infty}(-10^n e^{-td \cdot 10^n}+10^n e^{-t(d+1) \cdot 10^n}) \\ \simeq\,& \int_{ -\infty}^{\infty}(- 10^x e^{-td \cdot 10^x}+ 10^x e^{-t(d+1) \cdot 10^x}) dx \\ =\,& \frac{1}{\rm {ln}10}\int_{0}^{\infty}(-e^{-tdy}+e^{-t(d+1)y}) dy \\ =\,& \frac{1}{\rm{ln}10}\Big(-\frac{1}{td}+\frac{1}{t(d+1)}\Big).~~ \tag {11} \end{align} $$ There is one and only one approximation, i.e., we adopt an integration to replace a summation. Because $G_d(t)\to 0$ when $d\to \infty$, Eq. (11) can be integrated to yield $$ G_d(t)\simeq \frac{1}{t} \log_{10}\Big(1+\frac{1}{d}\Big).~~ \tag {12} $$ Then using Eq. (9), we obtain $$\begin{align} P_d =\,& \int_{0}^{\infty} F(x)g_d(x) dx \\ =\,&\int_{0}^{\infty} G_d(t)f(t) dt \\ \simeq & \int_{0}^{\infty} \frac{f(t)}{t} \log_{10}\Big(1+\frac{1}{d}\Big) dt \\ =\,& \log_{10}\Big(1+\frac{1}{d}\Big) \int_{0}^{\infty} \frac{f(t)}{t} dt \\ =\,& \log_{10}\Big(1+\frac{1}{d}\Big),~~ \tag {13} \end{align} $$ where we have used the following normalization condition of $f(t)$, $$\begin{align} 1 =\,& \int_{0}^{\infty} F(x) dx \\ =\,& \int_{0}^{\infty} dx\int_{0}^{\infty} f(t)e^{-tx} dt \\ =\,& \int_{0}^{\infty} dt f(t)\int_{0}^{\infty} e^{-tx} dx \\ =\,& \int_{0}^{\infty} \frac{f(t)}{t} dt.~~ \tag {14} \end{align} $$ Equation (13) is exactly the first digit law for the decimal system. Thus we show that well-behaved functions satisfy Benford's law approximately. A more rigorous derivation without the approximately equal signs in Eqs. (11)-(13) can be found in Ref.  [32]. Compared to Ref.  [32], the method provided above accords with our intuition better. In fact, unnecessarily complicated treatments are introduced to guarantee the strictness of the proof in Ref.  [32]. For example, a logarithmic scale is adopted after Laplace transform, merely to derive Eq. (12) of Ref.  [32], which corresponds to Eq. (13) in this study. Equation (13), though approximately holds, is set up on the original linear scale, thus manifests itself as a property of the direct Laplace transform, instead of the logarithmic Laplace transform which bears less intuitive physical meanings. In this study, no logarithmic transform is required to derive Benford's law. According to derivations so far, we can already explain the rationality of Benford's law through a clear chain of logic, as follows: (1) The integral of the product of $F(x)$ and $g(x)$ equals the integral of the product of the inverse Laplace transform of $F(x)$ and the Laplace transform of $g(x)$, i.e., Eq. (9). (2) The Laplace transform of $g(x)$ approximately equals the Benford term divided by $t$, i.e., Eq. (12). (3) The normalization condition of $F(x)$ guarantees that the integral of the inverse Laplace transform of $F(x)$ divided by $t$ equals $1$, i.e., Eq. (14). (4) Therefore, the integral of the product of $F(x)$ and $g(x)$ approximately equals the Benford term, i.e., Eq. (13). Such a chain of logic is not apparent in Ref.  [32]. The second significant digit law was also given by Newcomb.[1] In the decimal system, it is $$\begin{align} &P({\rm 2nd~digit~being~}d)\\ =\,&\sum_{k=1}^9\log_{10}(1+(10k+d)^{-1}),\\ &~d=0,1,\ldots,9.~~ \tag {15} \end{align} $$ Hill derived a general $i$th-significant digit law:[30] letting $D_i$ ($D_1,D_2,\ldots$) denote the $i$th-significant digit (with base 10) of a number (e.g., $D_1(0.0314)=3$, $D_2(0.0314)=1$, $D_3(0.0314)=4$), then for all positive integers $k$ and all $d_j \in {0,1,\ldots,9} $, $j=1,2,\ldots, k$, one has $$\begin{align} &P(D_1=d_1,\ldots,D_k=d_k)\\ &=\log_{10}[1+(\sum_{i=1}^k d_i \cdot10^{k-i})^{-1}].~~ \tag {16} \end{align} $$ We propose here a general form of digit law, and show that both the second significant digit law and the general $i$th-significant digit law are only corollaries of this general form. We calculate $P_{b,d,l,k}$, which is the probability that the integer composed of the first $k$ digits (base $b$) of an arbitrary number (e.g., for the number $0.0314$ and $k=2$, this integer is 31) is between $d$ and $d+l$ ($b^{k-1} \le d < d+l < b^k$). Correspondingly we introduce the density function $g_{b,d,l,k}(x)$ as $$ g_{b,d,l,k}(x)=\sum_{n=-\infty}^{\infty}[\eta (x-d \cdot b^n)-\eta (x-(d+l) \cdot b^n)],~~ \tag {17} $$ where the right-hand side is independent of $k$ (while $k$ puts restrictions on $d$ and $l$). Thus we can omit the subscript $k$ in the following. A similar technique gives the Laplace transform of $g_{b,d,l}(x)$, $$ G_{b,d,l}(t)\simeq \frac{1}{t} \log_{b}\Big(1+\frac{l}{d}\Big).~~ \tag {18} $$ Thus we arrive at the general significant digit law $$ P_{b,d,l,k}= \log_{b}\Big(1+\frac{l}{d}\Big).~~ \tag {19} $$ We find that Benford's law (2) corresponds to a special case of this general form for $k=1$ and $l=1$, whereas Hill's general $i$th-significant law (Eq. (16)) corresponds to the case for $b=10$, $d=\sum_{i=1}^k d_i \cdot 10^{k-i}$ and $l= 1$. Newcomb's second significant digit law can be considered as a corollary of Hill's law according to the addition principle in probability theory. We now calculate the error brought by our replacement of the summation to the integration in Eq. (11). Since Eq. (4) is always an accurate expression, the total error is $$ {\it \Delta}_{{\rm total}}=\int_{0}^{\infty} F(x)g_{b,d,l}(x) dx-\log_{b}\Big(1+\frac{l}{d}\Big).~~ \tag {20} $$ If we define $$ {\it \Delta}_{b,d,l}(t)=tG_{b,d,l}(t)-\log_{b}\Big(1+\frac{l}{d}\Big),~~ \tag {21} $$ the total error can be written as $$\begin{align} {\it \Delta}_{{\rm total}} =\,& \int_{0}^{\infty}\frac{ f(t)}{t}[t G_{b,d,l}(t) -\log_{b}\Big(1+\frac{l}{d}\Big)]dt \\ =\,& \int_{0}^{\infty}\frac{ f(t)}{t}{\it \Delta}_{b,d,l}(t)dt.~~ \tag {22} \end{align} $$ Checking the definitions of the two terms of ${\it \Delta}_{b,d,l}$, we can find that the variables of both of them can be multiplied by $b$ and the results remain unchanged, i.e., ${\it \Delta}_{b,d,l}$ is scale invariant. Hence $$ {\it \Delta}_{b,d,l}(bt)={\it \Delta}_{b,d,l}(t).~~ \tag {23} $$ For clarity, we define $$\begin{align} t=\,&e^s,~~ \tag {24} \end{align} $$ $$\begin{align} \widetilde{{\it \Delta}}_{b,d,l}(s)=\,&{\it \Delta}_{b,d,l}(e^s),~~ \tag {25} \end{align} $$ $$\begin{align} \widetilde{f}(s)=\,&f(e^s).~~ \tag {26} \end{align} $$ The corresponding normalization condition is $$\begin{align} \int_{-\infty}^{+\infty}\widetilde{f}(s) ds=1,~~ \tag {27} \end{align} $$ and the property Eq. (23) becomes $$ \widetilde{{\it \Delta}}_{b,d,l}(s+\ln {b})=\widetilde{{\it \Delta}}_{b,d,l}(s).~~ \tag {28} $$ Clearly, $\widetilde{{\it \Delta}}_{b,d,l}$ is a function of period $\ln b$. Furthermore, according to the result for exponential distribution in Ref.  [34] (Corollary 2, $\widetilde{f}(s)$ here is exactly $h_1(x)$ in Ref.  [34], $|\widetilde{{\it \Delta}}_{10,d,1}(s)|$ is $|h_1(x)-\log_{10}(1+\frac{1}{d} )|$ in the equation of Corollary 2), a rather good estimation can be made when $b=10$ and $d= l=1$, $$ 0.029 < \max|\widetilde{{\it \Delta}}_{10,1,1}(s) | < 0.03.~~ \tag {29} $$ We notice that the total error can be expressed as $$\begin{align} {\it \Delta}_{{\rm total}}=\int_{-\infty}^{+\infty}\widetilde{{\it \Delta}}_{b,d,l}(s)\widetilde{f}(s)ds,~~ \tag {30} \end{align} $$ where $\widetilde{f}(s)$ is dependent on $F(x)$ ultimately. In most cases, the correlation between $\widetilde{f}(s)$ and $\widetilde{{\it \Delta}}_{b,d,l}(s)$ is small, thus is the total error. Therefore, Benford's law can be a rather good approximation. However, if $\widetilde{f}(s)$ is close to a periodic function with the exact period $\ln b$, or $\widetilde{f}(s)$ changes signs very fast between positive and negative numbers (this may happen when $F(x)$ is artificially chosen, as the case of telephone numbers in a given region), the small $\widetilde{{\it \Delta}}_{b,d,l}(s)$ is counted and accumulated for many times, therefore the correlation becomes larger. Similar problems also exist for some special types of probability density functions, whose inverse Laplace transforms oscillate violently between positive and negative numbers, e.g., uniform distributions or normal distributions with small variances. Number sets drawn from such distributions, e.g., heights or ages of people, though being natural, still violate Benford's law. By arguing this, we point out that although the above derivation seems quite general, it cannot be universally true. More rigorous discussions about the err term with general applications to four types of number sets can be found in Ref.  [32]. A special case is when the integral of $\widetilde f(s)$ is not only convergent to 1, but also absolutely convergent to a positive real number $M$, then $$\begin{align} {\it \Delta}_{{\rm total}} \le\,&\int_{-\infty}^{+\infty}|\widetilde{{\it \Delta}}_{10,1,1}(s)||\widetilde{f}(s)|ds \\ \le\,& \int_{-\infty}^{+\infty}0.03|\widetilde{f}(s)|ds \\ =\,& 0.03M.~~ \tag {31} \end{align} $$ If $f(s)$ is a positive or negative definite function, it is absolutely integrable. Such an $F(x)$ is called the completely monotonic function in mathematics. This means that $M$ is 1, thus ${\it \Delta}_{{\rm total}}$ is not greater than 0.03. Consequently Benford's law is a good estimation. For example, when $F(x)=\frac{2}{\sqrt{x}}e^{-\sqrt{x}}$, $f(t)=\frac{2}{\sqrt{\pi t}}e^{-\frac{1}{4t}} >0$, and when $F(x)=\frac{4}{(x+1)^5}$, $f(x)=\frac{e^{-t}}{6}>0$. We can assert that in these cases, the total errors are less than 0.03. In fact, numerical results are 0.0005 and 0.009. This verifies our estimation. As a rule of thumb, distributions with monotonic decreasing and relatively smooth probability density functions often conform to Benford's law well.[32] Inverse Laplace transforms of such probability density functions generally change signs only for finite times, thus being absolutely convergent. To understand this, one can view inverse Laplace transform as decomposing the original probability density function into a series of exponential functions, among which some are positive and others negative. If a monotonic decreasing probability density function is relatively smooth, i.e., without a sharp change of probability density, it can be approached mainly by positive exponential distributions, therefore its inverse Laplace transform does not oscillate between positive and negative numbers very much. As an application of this rule of thumb, for non-monotonic decreasing distributions, we can transform them into monotonic decreasing distributions, resulting in better performance of Benford's law, e.g., for normal distributions, we can subtract the mean value from the original data set and obtain a monotonic decreasing distribution. The above calculations and derivations tell us that the significant digit behaviors demonstrate that although our nature has no preference to any specific number, it does have discrimination to digits in numbers as a logical consequence of human counting systems. Therefore, our results justify the conventional wisdom that the violation of Benford's law is a sign that a table of numbers is artificial or anomalous. The underlying reason for the uneven distribution of first digits is due to the basic property of digital system, but not some dynamic source behind nature as people suspected. This also explains why we can use Benford's law to distinguish anomalies or unnaturalness in artificial numbers. The mathematical expressions and derivations provided here are simple, elegant, and all with clear intuitive pictures. They are easily comprehensible. Therefore, this version of proof of Benford's law can also serve as an example for the application of the Laplace transform. The first digit law reveals an astonishing regularity in realistic numbers. We provide in this work a proof of this law from the Laplace transform, and point out the condition for the validity of the law. Compared to Ref.  [32], the derivation here is simple and elegant, and it directly reveals the rationality of the first digit law. From our work, the first digit law is due to the basic structure of the number system. Thus the first digit law is a general rule that applies to vast data sets in the natural world as well as in human social activities. It is no longer strange why Benford's law is so successful in various domains. Such a law should be regarded as a basic mathematical knowledge with great potential for vast applications. References Note on the Frequency of Use of the Different Digits in Natural NumbersBenford’s law and physical constants: The distribution of initial digitsOn the Peculiar Distribution of the U.S. Stock Indexes' DigitsHow do numbers begin? (The first digit law)Survival Distributions Satisfying Benford's LawBenford's Law and the screening of analytical data: the case of pollutant concentrations in ambient airThe number of cells in colonies of the cyanobacterium Microcystis aeruginosa satisfies Benford's lawBenford’s Law Applied to Hydrology Data—Results and Relevance to Other Geophysical DataEmpirical mantissa distributions of pulsarsThe significant digit law in statistical physicsFirst-digit law in nonextensive statisticsAn illustration of Benford's first digit law using alpha decay half livesBenford’s law and half-lives of unstable nucleiBenford's Law and β-Decay Half-LivesBenford’s law and cross-sections of A(n, $ \alpha$ )B reactionsBenford's Law in Nuclear Structure PhysicsFIRST DIGIT DISTRIBUTION OF HADRON FULL WIDTHDo dynamical systems follow Benford’s law?One-dimensional dynamical systems and Benford's lawMulti-dimensional dynamical systems and Benford's LawNot the First Digit! Using Benford's Law to Detect Fraudulent Scientif ic DataOn the Distribution of First Significant DigitsScale-Distortion Inequalities for Mantissas of Finite Data SetsThe Significant-Digit PhenomenonBase-Invariance Implies Benford's LawA Statistical Derivation of the Significant-Digit LawFirst digit law from Laplace transformBenford's law for exponential random variables
[1] Newcomb S 1881 Am. J. Math. 4 39
[2]Benford F 1938 Proc. Am. Philos. Soc. 78 551
[3] Burke J and Kincanon E 1991 Am. J. Phys. 59 952
[4] Ley E 1996 Am. Stat. 50 311
[5] Torres J et al 2007 Eur. J. Phys. 28 L17
[6] Leemis L M et al 2000 Am. Stat. 54 236
[7] Brown R J C 2005 Analyst 130 1280
[8] Costas E et al 2008 Aquat. Bot. 89 341
[9] Nigrini M and Miller S J 2007 Math. Geol. 39 469
[10] Shao L and Ma B Q 2010 Astropart. Phys. 33 255
[11] Shao L and Ma B Q 2010 Physica A 389 3109
[12] Shao L and Ma B Q 2010 Phys. Rev. E 82 041110
[13] Buck B et al 1993 Eur. J. Phys. 14 59
[14] Ni D and Ren Z 2008 Eur. Phys. J. A 38 251
[15] Ni D D et al 2009 Commun. Theor. Phys. 51 713
[16] Liu X J et al 2011 Eur. Phys. J. A 47 78
[17] Jiang H et al 2011 Chin. Phys. Lett. 28 032101
[18] Shao L and Ma B Q 2009 Mod. Phys. Lett. A 24 3275
[19] Tolle C R et al 2000 Chaos 10 331
[20] Berger A et al 2005 Trans. Am. Math. Soc. 357 197
[21] Berger A 2005 Discrete Contin. Dyn. Syst. A 13 219
[22]Nigrini M J 1996 J. Am. Tax. Assoc. 18 72
[23]Nigrini M J 1999 Intern. Auditor 56 21
[24]Rose A M and Rose J M 2003 J. Accountancy 196 58
[25] Diekmann A 2007 J. Appl. Stat. 34 321
[26] Pinkham R S 1961 Ann. Math. Stat. 32 1223
[27] Berger A et al 2008 J. Theor. Probab. 21 97
[28] Hill T P 1995 Am. Math. Mon. 102 322
[29] Hill T P 1995 Proc. Am. Math. Soc. 123 887
[30] Hill T P 1995 Stat. Sci. 10 354
[31]Berger A et al 2009 Benford Online Bibliography http://www.benfordonline.net/
[32] Cong M et al 2019 Phys. Lett. A 383 1836
[33]Marsden J and Hoffman M 1999 Basic Complex Anal. 3rd edn (New York: W. H. Freeman) p 471
[34] Engel H A and Leuenberger C 2003 Stat. Probab. Lett. 63 361