A Gaussian-Distributed Quantum Random Number Generator Using Vacuum Shot Noise

Abstract

Among all the methods of extracting randomness, quantum random number generators are promising for their genuine randomness. However, existing quantum random number generator schemes aim at generating sequences with a uniform distribution, which may not meet the requirements of specific applications such as a continuous-variable quantum key distribution system. In this paper, we demonstrate a practical quantum random number generation scheme directly generating Gaussian distributed random sequences based on measuring vacuum shot noise. Particularly, the impact of the sampling device in the practical system is analyzed. Furthermore, a related post-processing method, which maintains the fine distribution and autocorrelation properties of raw data, is exploited to extend the precision of generated Gaussian distributed random numbers to over 20 bits, making the sequences possible to be utilized by the following system with requiring high precision numbers. Finally, the results of normality and randomness tests prove that the generated sequences satisfy Gaussian distribution and can pass the randomness testing well.

1. Introduction

Random numbers are of extreme importance for a great range of applications from scientific to engineering fields, including statistical sampling, numerical simulation, lottery and cryptography. A typical example is the quantum key distribution (QKD), in which true random numbers are essential to guarantee its unconditional security [1,2,3,4]. The algorithm-based classical pseudo-random number generators have been widely applied for their simple implementation and extremely high generation rate [5]. However, the inherent determinacy of pseudo-random number generator makes it substantially deterministic and predictable, which leads to the failure in satisfying theoretically requirements of secure communication systems. Aside from the algorithmic method, extracting randomness from objective physical processes is feasible. An outstanding alternative is a quantum random number generator (QRNG), which exploits the intrinsic random nature of quantum mechanics [6,7], acts as a promising method in generating truly random numbers.

Practical (conventional) QRNG schemes could generate relatively large amount of quantum random numbers with high generation rate via utilizing easily accessible commercial devices. Schemes of various quantum random sources have been demonstrated, including discrete ones by measuring photon path [8,9,10], photon arrival time [11,12,13,14,15,16], photon number distribution [17,18], as well as continuous types of phase noise of lasers [19,20,21,22,23,24,25] and Raman scattering [21], intensity fluctuation of amplified spontaneous noise (ASE) [26,27,28,29], quadrature fluctuation of vacuum shot noise (VSN) [30,31,32,33,34], and other potential quantum state [35,36].

Practical schemes are proposed under assumptions only if the QRNG model works well, indicating the system is fully trusted. This rigorous condition can hardly be fulfilled, and device-independent (DI) protocols are proposed for closing the loophole. DI QRNG verify the randomness physically, taking the violation of Bell’s inequality [37,38] as a judgment [39,40]. Later, two branches are researched for alternative proposes, namely randomness extraction [41,42,43] and randomness amplification [44,45]. While DI protocols sacrifice too much on feasibility, a third choice which compromises between practical scheme and DI protocol is proposed. These semi-DI protocols merely make a reasonable assumption on critical devices [46,47,48,49,50,51], pursuing practical security instead of unconditional security.

The application fields of Gaussian RNG are diverse, in which the most significant application is a simulation, ranging from Monte Carlo method to simulation of communication channels and noises, biology, psychology, and so on. Specific to quantum information, Gaussian RNG provides Gaussian distributed random numbers for the modulation of coherent states in continuous-variable QKD systems [4,52,53]. However, all the previous QRNG schemes provide uniformly distributed random numbers. Despite the universality that uniformly distributed random numbers could be converted to any distributions mathematically, the conversion process itself somehow costs much time and resources. An even higher potential risk is the process is that approximate in principle [54], which may lead to the defects of performance in applications. In fact, most of the continuous-variable quantum random sources, owing to the central limit theorem, feature Gaussian distributed signals in the time domain, including vacuum shot noise and phase noise of laser. Hence, it is possible for hardware-based schemes, naturally including QRNG, to utilize the Gaussian distribution profile and directly generate random numbers as required.

In this paper, a practical scheme directly generating Gaussian distributed quantum random numbers is proposed. Here “directly” means there are no conversion steps from the uniform distribution to the Gaussian distribution, however, the scheme is not post-processing free. Firstly, we point out the inherent difference in entropy estimation for Gaussian distribution versus uniform distribution. Practical issues of sampling devices are discussed for entropy estimation and system optimization. Secondly, a novel post-processing method is proposed, which takes a step further from the recursive method in classical Gaussian distributed RNG [55]. It is designed to remove the impacts of classical noise in the system, along with fulfilling the precision and auto-correlation requirements from applications. Finally, an experimental setup is demonstrated to show the feasibility of this scheme, using vacuum fluctuation of the quantum state as a quantum random source, and the implementation has passed tests both on normality and randomness.

The structure of this article is described as follows. In Section 2, firstly we discuss the difference in entropy estimation between Gaussian and uniform distribution, followed by the analysis on the impacts of practical sampling device to the system, namely sampling range and sampling resolution. In Section 3, a novel post-processing method is proposed to overcome the disadvantage of low precision in sampling, and substantially eliminate the impacts of electronic noise. In Section 4, an experimental setup is demonstrated, as well as the optimization and post-processing operation on a practical system. Finally, the test results for both normality and randomness are shown.

2. Analysis of Gaussian Distribution QRNG Scheme

2.1. Gaussian Random Source and Entropy Estimation

2.1.1. Vacuum Fluctuation

In principle, most quantum random sources with Gaussian distributed signals in the time domain can be applied in our scheme. Particularly, for the following excellent features, we choose vacuum fluctuation of the quantum state as the random source. Firstly, vacuum shot noise is caused by vacuum fluctuation, thus the randomness of the pure state is secured. Secondly, it is a Gaussian state, which means the measurement of either position or momentum quadrature $\widehat{x}$ or $\widehat{p}$ in a pure state will always follow a Gaussian distribution. Finally, it is identical, which means additional vacuum fluctuation introduced by devices, such as the beam splitter, will not affect the randomness of the quantum source.

The Wigner function of vacuum fluctuation is as follows:

$$W_0(x,p)=\frac{1}{\pi }exp(-x^2-p^2).$$

(1)

As a quasi-probability function, one can repeatedly measure either $\widehat{x}$ or $\widehat{p}$ quadrature, given fixed phase difference $\theta$ between the vacuum and LO signals. Taken $\widehat{x}$ quadrature as an example, the probability density distribution (PDF) of detected signal should be:

$${\left|{\psi }_0\left(x\right)\right|}^2={\int }_{-\infty }^{+\infty }{W}_0(x,p)dp=\frac{1}{\sqrt{\pi }}exp(-x^2),$$

(2)

which is perfectly Gaussian distributed, with mean value $\mu =0$ and variance ${\sigma }^2=1/2$ centering at the origin in phase space.

Security is always an important issue to a cryptographic system, including quantum random number generator, compared to its classical counterpart. While there definitely exist some risks of leaking information to an adversary in randomness extraction, modeling of vacuum fluctuation also takes advantage of its property, of which it could never be tampered even by the most powerful adversary under the limitation of physical laws. Hence, unlike traditional applications of classical RNG, where the noises are usually treated as introduced by the system itself, we could regard any noise in the system as introduced by the eavesdropper (Eve) in the QRNG system, in attempt to reach a lower bound in entropy estimation.

For homodyne detection, signals of two balanced arms are subtracted to supress the common mode noise, while the amplification factor is decided by the system:

$$\begin{array}{ccc}\hfill \Delta n& =& {\widehat{n}}_2-{\widehat{n}}_1=\widehat{a}{\widehat{a}}_{LO}^{†}+{\widehat{a}}_{LO}{\widehat{a}}^{†}\hfill \\ \hfill V_{samp}& \propto & A\Delta n\hfill \end{array}$$

(3)

$$\begin{array}{cc}=& A\left|{\alpha }_{LO}\right|(\widehat{x}cos\theta +\widehat{p}sin\theta )\hfill \end{array}$$

(4)

where $\widehat{a},{\widehat{a}}^{†}$ are annihilation and creation operators, and $\widehat{n}=\widehat{a}{\widehat{a}}^{†}$ is photon number operator. $V_{samp}$ is the signal at sampling device (after subtraction), A is the amplification factor of the system excluding LO signal, $\left|{\alpha }_{LO}\right|$ indicates the X quadrature of LO signal, and $\theta$ is the phase difference between vacuum and LO signal.

2.1.2. Entropy Estimation

As a conventional scheme, entropy estimation should be done before randomness extraction. The most significant difference between uniform and Gaussian distribution, from the perspective of information theory, is that the information entropy $H\left(X\right)$ should have different maximal value under different constraints.

In order to eavesdrop most information under classical scenario, Eve’s best strategy is figuring out ${max}_{x\in X}{p}_i$, the highest probability of a single bin in a random variable X, which directly related to the minimal entropy (min-entropy):

$$H\left(X\right)=-{log}_d\underset{x\in X}{max}{p}_i$$

(5)

where d is the base of logarithmic function that defines whether the signal is binary, decimal and so on. For binary information, we often define $d=2$. However, if the signal precision n is more than one bit, it could be also treated as $d=2^n$.

According to Equation (5), uniform distribution possesses the highest min-entropy with no constraints. Noticing that for the continuous case, the classical entropy $H\left(X\right)$ always goes to infinity for ideal sampling device with infinite sampling range and precision. Therefore, we assume the total amount of information a single signal carries is 1, thus normalize the maximal value of information entropy rate. This method is also adopted for Gaussian distribution entropy estimation in the following analysis.

Meanwhile, in applications under certain constraints, namely the quadratic quantity of energy (or power) of signal is fixed, one is expecting a different PDF. This conclusion is naturally derived from the property that among all distributions with the same variance ${\sigma }_m^2$, Gaussian distribution possess the highest information entropy $H\left(X\right)$:

$$H\left(X\right)={\int }_{-\infty }^{+\infty }-p\left(x\right){log}_{d}p\left(x\right)dx,$$

(6)

where $p\left(x\right)$ is PDF of ideal Gaussian distribution.

This property indicates that, if the variance of continuous noise signal ${\sigma }_m^2$ is observable and steady, Gaussian distribution, instead of uniform distribution, could achieve a higher entropy. Fortunately, the variance of total noise is indeed measurable in a QRNG scheme, and perfectly matches the assumption. Therefore when adopting a Gaussian distributed random source, and the output is supposed to be Gaussian distributed, the information quantity acquired is significantly reduced during the conversion phase of uniform distributed RNG schemes.

To achieve a Gaussian distributed QRNG scheme, we should adopt the goodness of fit (GoF) test essentially, to verify whether the PDF of our samples are sufficiently close to Gaussian distribution. For Gaussian distribution, there are several specific methods, namely Kolmogorov-Smirnov test and Anderson-Darling test, which will be introduced detailed in Appendix A.

2.2. Impact of Sampling Device

In the entropy estimation phase, a similar idea of “worst-case scenario” to its uniform counterpart could be adopted. Alice loses some entropy due to the sampling device, while Eve may acquire original information from ideal Gaussian distribution. As a continuous distribution, either infinite sampling range or precision is not practical, and will cause the entropy to be infinite, hence we should set conditions considering the performance of the practical device.

Classical Gaussian RNG often set $\pm 10\sigma$ as the bounds in high multiple-sigma test [56], it seems reasonable to follow this assumption. Meanwhile, sampling precision can hardly exceed 20 bits for current commercial analog-to-digital converter (ADC). Practical issues of range, precision and depth will be discussed in detail.

Despite the classical noise, our scheme is still a trusted device scheme, where the extractable randomness of the scheme is described as:

$$R_{dis}=I(A:B),$$

(7)

where $R_{dis}$ refers to the generation rate of a QRNG, with “dis” indicates the discretized samples which may lose some information, and $I(A:B)$ is the mutual information between the authorized users (Alice, Bob) in a cryptographic system. In QRNG scheme, specifically speaking, Alice is the random source and Bob is the randomness extractor. Apparently, QRNG could be (and in most occasions is) local, while Alice plays both roles of the sender (random source) and receiver (randomness extractor), thus $I(A:B)$ is actually determined by the entropy $H\left(A\right)$ of measured classical data.

Particularly, in the vacuum fluctuation scheme we demonstrated below, the variance of total noise ${\sigma }_m^2$ and the variance of classical noise ${\sigma }_c^2$ can be observed by separately turn on/off the LO signal. Since we believe the quantum noise Q and classical noise E are independent from each other, the min-entropy of quantum noise is a conditional entropy, with classical noise part E is given by [33]:

$$\begin{array}{cc}\hfill H_{min}\left(M\right|E)& =-{log}_2\left[\underset{e\in E}{max}\underset{m\in M}{max}{P}_{M|E}\left(m\right|e)\right]\hfill \\ & =-{log}_{2}max\left(c_1,c_2\right)\hfill \end{array}$$

(8)

where $M,E$ are the random variable of total measured noise and electric noise, $m,e$ mean the specific measured value. $R,\delta$ means the sampling range and sampling resolution respectively. $c_1=\frac{1}{2}\left[\mathrm{erf}\left(\frac{e_{max}-R+3\delta /2}{2{\sigma }_q}\right)+1\right],c_2=\mathrm{erf}\left(\frac{\delta }{2\sqrt{2}{\sigma }_q}\right)$ refer to the two possible values that could be the maximal $p_i$ in Equation (5), and $e_{max}$ is the maximal possible electric noise.

Things differ a little under the Gaussian scenario comparing to the analysis in Ref. [33]. In the uniform scenario, the optimization for the system is setting $c_1=c_2$ to achieve maximal value of min-entropy. However, if we adopt $c_1=c_2$ in the Gaussian scheme, the raw data will definitely fail the GoF test. Therefore, we have to analyze the impact of the sampling device under the restriction of the GoF test, where there always exists $c_1<c_2$.

2.2.1. Sampling Range

The sampling device will change the instantaneous voltages beneath (above) the lower (upper) threshold into $V_{min}$ (or $V_{max}$). Parameter k is the ratio between sampling range R and the deviation of signal $\sigma$. Finite sampling range will truncate probability distribution $P(x\ge \left|k\sigma \right|)$ outside the range $\pm k\sigma$, another consequence is a significant defect at the tails, causing the PDF non-Gaussian.

We define a parameter called normalized min-entropy in our analysis. Supposing an ideal Gaussian distributed random variable, and the information carried by the variable is described as $H_{ideal-min}$ before normalized to 1. When taking practical sampling device into consideration, the distribution is changed and entropy is estimated by Equation (8), however, it should be monitored by the GoF test. In the following analysis of sampling range and resolution, the utmost assumption is that the signal should satisfy normality, meanwhile, this is also the assumption of the post-processing method below. Therefore, the min-entropy of distribution from a practical system should initially pass the GoF test, before it can be normalized according to the ideal case, which could be described as:

$$H_{norm-min}=H_{min}\left(M\right|E)/H_{ideal-min}$$

(9)

Figure 1 shows the relationship between sampling range and entropy $H_{norm-min}$. Cases $R\le \pm 3.5\sigma$ are discarded for all precisions, due to these cases feature defected PDF and frequently fail the GoF test (with default significance level at $\alpha =0.01$). However, lower sampling precision n with too large a sampling range will also fail the GoF test (as the curve $n=12$ stops at $R=\pm 4.6\sigma$), since the discretization effect is notably increased for lower precision cases.

• If k is too small, $V_{min}$ ($V_{max}$) will occur too often, making the random variable more predictable, and reducing entropy $H_{dis}\left(X\right)$. Furthermore, the worse profile of Gaussian distribution has a higher possibility to fail the GoF test, which does not match our requirement in post-processing and applications;
• If k is too large, most signals will locate in a small range of sample bins, making the most significant bits (MSB) of samples more predictable, and also reducing entropy $H_{dis}\left(X\right)$. On the other hand, many sampling bins are unoccupied, wasting the ability of devices and substantially reduce the sampling precision.

For cases that successfully pass the GoF test, the normalized min-entropy decreases as the sampling range increases. The curves do have a period of rapid increase under no constraint assumption within the range of $R\le \pm 3\sigma$, however, these cases are rejected by the GoF test. From the view of variance, as long as the raw data pass GoF test, the higher k value is, the lower normalized min-entropy is, which shows great significance on the matching of signals and range of sampling device.

2.2.2. Sampling Resolution

Finite sampling resolution $\delta$ will result in information loss of probability distribution inside the minimal discrete sampling interval, i.e., resolution: $\left[x_i-1/2\delta ,x_i+1/2\delta \right](i\in [0,2^n-1])$. Intuitively, entropy grows monotonously as the precision n increases. If n is too small, too many detailed information is lost, and we can hardly extract random numbers after entropy estimation.

Figure 2 shows the relationship between precision and entropy $H_{norm-min}$. For the same reason discussed in sampling range analysis, cases $R\le \pm 3.5\sigma$ are discarded. Despite precision below $n=12$ will frequently fail the GoF test due to a strong discretization effect, we estimate the entropy to show the trend of entropy curve.

2.2.3. Sampling Depth

Sampling depth (maximal samples in a single buffer) mainly affects the practical system on the GoF test. As the test statistic shows, the AD test is distribution-free, but sample data size L related. An identical distribution with sample length 10 times larger would lead to approximately linearly increased test statistic, while the critical value remains the same. This is due to it possessing a larger sample space, any violation on the PDF becomes more significant to be detected by the GoF test. Therefore, a larger buffer should have a better PDF for raw data to pass the GoF test.

According to the three factors discussed above, we consider $R=\pm 4\sigma$ to be the optimal sampling range for noise-free cases, while the precision and depth should be as high as possible, which is not so crucial in uniform occasions. It is highly recommended that if one wants to achieve a Gaussian distribution QRNG, sampling precision should be at least 12 bits for decent performances. Although the sampling range of the practical device is often fixed, one can adjust the amplification factor $A\left|{\alpha }_{LO}\right|$ to adapt the range, aiming at achieving better performance. However, noise introduced by the system with a variance of ${\sigma }_c^2$ often alters the PDF and optimization condition. If the noise introduced by the system is not crucial enough to change the PDF, the following post-processing method could significantly reduce its influence.

3. Post-Processing

Post-processing is an essential part in QRNG scheme. It is adopted to remove the impacts of classical noise in the system as well as the imperfections caused by finite sampling. Most of the post-processing methods can also improve the probability distribution of the raw data.

The Toeplitz matrix hashing method [23,57,58] is widely acknowledged as the most effective method in QRNG post-processing. However, the whole method aims at uniform distribution generation [59], hence does not meet our requirement. Here we propose another post-processing method originating from recursive method [55] adopted in classical Gaussian distribution RNG schemes.

The recursive method takes the essence of Gaussian distribution that is, the summation of any amount of Gaussian distributed variables is still Gaussian distributed:

$$Y=\sum _i{k}_i{X}_i,$$

(10)

while the original Gaussian variables $X_i$ satisfy $X_i=N({\mu }_i,{\sigma }_i^2)$, the output Y should satisfy $Y=N({\sum }_i{k}_i{\mu }_i,{\sum }_i{k}_i^2{\sigma }_i^2)$.

Traditional central limit theorem (CLT) of non-Gaussian cases is only valid for a large amount of independent identical distributed (i.i.d.) variables. On the contrary, we notice that the recursive method takes merely four elements as Equation (11) shows. By adopting the recursive method, one could avoid the risk that raw data of distinguishable non-Gaussian variables are converted to identical Gaussian distribution.

The original transfer matrix $T_{\mathrm{rec}}$ is derived from following operations:

$$\begin{array}{cc}\hfill a_1^{\prime }& =c-a_3,\hfill \\ \hfill a_2^{\prime }& =c-a_2,\hfill \\ \hfill a_3^{\prime }& =a_1-c,\hfill \\ \hfill a_4^{\prime }& =a_4-c,\hfill \end{array}$$

(11)

where $c=\frac{1}{2}{\sum }_{i=1}^4{a}_i$. Thus we can denote the relationship between input and output vectors $A_i,A_i^{\prime }$, as well as the operating matrix $T_{rec}$ ($\frac{1}{2}$ is normalization coefficient):

$$\begin{array}{ccc}\hfill A_i^{\prime }& =& T_{\mathrm{rec}}{A}_i,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill T_{rec}& =& \frac{1}{2}\left[\begin{array}{cccc}1& 1& -1& 1\\ 1& -1& 1& 1\\ 1& -1& -1& -1\\ -1& -1& -1& 1\end{array}\right]\hfill \end{array}$$

(13)

The output of the recursive method possesses the perfect auto-correlation property. However, it cannot extend the precision of a single number. To make full utilization of all significant bits from different raw data in precision extension, adding random numbers from i.i.d. Gaussian distributed variables with different weight is an effective method.

The post-processing method includes two steps. Firstly, we utilize the m-MSB (most significant method) as pre-processing. When entropy estimation phase introduced in Section 2 is done, the value m utilized in m-MSB processing is:

$$m=\lfloor H_{min}\left(M\right|E)\rfloor .$$

(14)

Then we should adopt an operation that could achieve precision extension based on the matrix in Equation (13). Noticing that since the raw data has passed GoF test, the condition in Equation (10) is satisfied.

Assuming X is the original variable from ADC, and we could divide X into groups of Gaussian distributed variables $X_i$, before taking operation as Equation (10) shows. As an example where we divide the raw data into $l=4$ groups, consecutive four random numbers $x_{4i-3},x_{4i-2},x_{4i-1},x_{4i}$ will form a vector $A_i$ before operating by the matrix. Particularly in Equation (10), suppose $k_i=2^{-i}$, then every adjacent raw data in vector $A_i$ shifts only 1 more bit, thus the summation has a precision of $n=m+l-1$ bit, while $m,n$ are the precision of variables $X_i$ and Y respectively, and l is the number of groups.

Combining the analyses above, we modify $T_{rec}$, adding different weights in the matrix similar to the original method:

$$\begin{array}{cc}\hfill a_1^{″}& =1/2a_1+1/4a_2-1/8a_3+1/16a_4,\hfill \\ \hfill a_2^{″}& =1/16a_1-1/2a_2+1/4a_3+1/8a_4,\hfill \\ \hfill a_3^{″}& =1/8a_1-1/16a_2-1/2a_3-1/4a_4,\hfill \\ \hfill a_4^{″}& =-1/4a_1-1/8a_2-1/16a_3+1/2a_4,\hfill \end{array}$$

(15)

Thus we can denote ($k_{\mathrm{NC}}$ is normalized coefficient):

$$\begin{array}{ccc}\hfill A_i^{″}& =& S_{rec}{A}_i,\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill S_{rec}& =& k_{\mathrm{NC}}\left[\begin{array}{cccc}1/2& 1/4& -1/8& 1/16\\ 1/16& -1/2& 1/4& 1/8\\ 1/8& -1/16& -1/2& -1/4\\ -1/4& -1/8& -1/16& 1/2\end{array}\right]\hfill \end{array}$$

(17)

Noticing that, the structure of $S_{rec}$ is very much similar to the original structure of $T_{rec}$. Both of them share two rows/columns with three positive and one negative element, and others with three negative and one positive. This type of structure is convenient for expansion to a $8\times 8$ or even larger size of $T_{rec}$ [55]. For $S_{rec}$, the expansion method is similar, as long as obeying the rules discussed below.

A crucial difference between the original recursive method and our modified method is that, since we introduce different (absolute) value in the operating matrix, the auto-correlation coefficient will not remain flat. Therefore, we can only extract one number from $A_i^{″}$ (of $n=m+l-1$ bit precision, where m is the precision of $A_i$, n is the precision of $A_i^{″}$ as the final output, and l is the size of $S_{rec}$), while in the original case all numbers of $A_i^{\prime }$ (of m bit precision) could be extracted.

The recursive method post-processing operation can be designed, hence it is definitely more flexible than Equation (10). Utilizing matrix for a precision extension instead of simply adding i.i.d. Gaussian distributed variables have several merits:

• Elements in the matrix, which are the weights in Equation (10), is not fixed, as long as they obey fundamental rules. For $4\times 4$ matrix, each row/column should have 3 (1) positive and 1 (3) negative elements, and the position should not be the same; the absolute value of each row and column should not be the same either. Thus there is a group of $S_{rec}$ with hundreds of possible matrices;
• The size of the matrix can be designed, which indicates how many raw numbers will be used to generate a final number. We take the $4\times 4$ matrix as the simplest example for a demonstration. However, when the precision after m-MSB pre-processing is inadequate, and a larger matrix should be made. For instance, in the following section of implementation, we generate 12-bit Gaussian distribution numbers from 5-bit pre-processed data, by utilizing an $8\times 8$ matrix. If the matrix size is larger, it has a potential for even higher precision, such as five-bit pre-processed data with a $16\times 16$ matrix will generate 20-bit Gaussian distribution random numbers for high multiple-sigma applications.
• The values of matrix elements can also be designed, which indicate shifted bits of the pre-processed data. In the discussion above, weights of adjacent numbers always follow the power of $1/2$, which means that adjacent numbers in $A_i$ should shift one bit in the summation operation. However, if we change $1/2$ to $1/4$, it means that adjacent numbers in $A_i$ should shift two bits. Remember that according to Equation (17), a normalized coefficient $k_{\mathrm{NC}}$ should be carefully calculated to match the designation, making sure that the input and output share the same variance.

Due to these merits above, one can design his/her own $S_{rec}$ matrix for alternative experimental setup and application requirements. Furthermore, these properties leave huge space for further introduction of pre-generated random seed. It is possible to prepare several operating matrices and, based on the random seed that generated before or even feedback from real-time QRNG scheme, alter the post-processing operation in real-time.

Table 1. The relative entropy $H\left(p\right(x\left)\right|q\left(x\right))$ between unknown distribution $p\left(x\right)$ (with normalized variance) and reference $q\left(x\right)$ after post-processing. All data unit is ${10}^{-5}$. $H_{\mathrm{rel}}=0$ means the unknown distribution is identical with $q\left(x\right)$. We assume $p\left(x\right)$ a standard Gaussian distribution with minor classical noise of Quantum-to-Classical Noise Ratio (QCNR) ranges between 3–20 dB. In order to highlight the smooth effect on profile, data is designed for small size with $n_{\mathrm{tot}}=10M$. The residual relative entropy after post-processing is possibly due to the finite size effect of this calculation method.

	Normal		t-Dist.		Uniform		Rayleigh
QCNR(dB)	Before	After	Before	After	Before	After	Before	After
3	1.2225	1.1653	1.2966	1.3338	64.036	4.6896	179.40	1.4993
6	1.2582	1.2320	1.4416	1.4348	9.0556	1.4507	39.991	1.2510
10	1.1920	1.2031	1.2478	1.2741	1.4064	1.3917	4.5185	1.1799
20	1.2717	1.2455	1.2132	1.2510	1.1764	1.1996	1.2150	1.1964

shows the relative entropy $H\left(p\right(x\left)\right|q\left(x\right))$ between $p\left(x\right)$ and $q\left(x\right)$. $p\left(x\right)$ is quasi-Gaussian distributed, mixing ideal Gaussian distribution with several types of classical noise of small variance. $q\left(x\right)$ is the reference of standard Gaussian distribution. It is clear that the post-processing method dramatically reduces the impact of noise for low Quantum-to-Classical Noise Ratio (QCNR) cases, especially for those noises which a not Gaussian, regardless of the profile of raw data. However, noise is still distinguishable from a standard Gaussian distribution.

4. Implementation and Results

4.1. Experimental Setup

We experimentally demonstrate our scheme and the setup is described as follows (as shown in Figure 3). The local oscillator (LO) is 1550 nm distributed feedback laser (NKT Basic E15, linewidth 100 Hz) with adjustable output power up to 15 mW, connecting to an external variable optical attenuator (VOA) precisely setting the amplification factor of the LO signal. Vacuum shot noise, physically provided by blocking one input port of a 50:50 beam splitter (BS), interferes with the LO light. The signals are sent to a well-tuned homemade AC coupling homodyne detector (measurement bandwidth limited to 100 MHz by low-pass filter) to measure the noise. Following circuits including an analog-to-digital converter (ADC, ADS5400, sampling frequency 200 MHz, sampling precision 12 bits and input voltage range 1.5 V peak-to-peak), a field-programmable gate array (FPGA, KC705 evaluation board) that realizes randomness extraction and data precision adjustment. The power spectral density function of total noise and classical noise is shown in Figure 4.

To obtain better performance, the power of LO light was examined by setting different LO power with fixed steps. When LO light was off, the vacuum fluctuation can be ignored, and classical noise contributes to the output with variance of ${\sigma }_c^2$, which is quite steady. The variance of total noise ${\sigma }_m^2$ increases as LO power (after VOA) getting stronger. The linear region ends when LO power increases at around 9.5 mW, and finally saturates at around 13 mW.

Our system features a high QCNR to obtain more potential information from the signal. By setting LO power slightly less than saturation at around 12 mW (6 mW for each branch of the balanced detector), we have acquired 12-bit raw data after ADC, and calculate the variance of signals when LO light is on/off, representing the total and classical noise respectively. Noticing that, all the units mentioned here are sampling bins, and according to the ability of our sampling device, one sampling bin roughly equals 0.366 mV.

Firstly, due to the fixed sampling range, $V_{\mathrm{range}}=2^{12}$ and peak-to-peak value is $V_{\mathrm{p}-\mathrm{p}}=200$, around 3-bit MSB is discarded. Since the variance of total noise is ${\sigma }_m^2=1200.7$, classical noise is ${\sigma }_c^2=82.5$, variance of quantum noise can be calculated: ${\sigma }_q^2=1118.2$, thus the maximal QCNR is defined by:

$$\gamma =\mathrm{QCNR}=\frac{{\sigma }_m^2-{\sigma }_c^2}{{\sigma }_c^2}=13.55,$$

(18)

with $\mathrm{QCNR}=13.55\left(11.3\mathrm{dB}\right)$, the classical noise after normalization is $\epsilon =1/\mathrm{QCNR}=0.074$.

As QCNR indicates, classical noise can only fluctuate in a small range of voltage. The MSB part of the residual sample is more likely to be affected by quantum noise, while the LSB part is affected by both quantum and classical noise, which is opposite to uniform distributed occasions.

We adopted entropy estimation initially in the post-processing phase. Our ADC has a sampling range of 1.5 V peak-to-peak and sampling precision of 12 bits, thus the quantization error is ${(\delta /12)}^2=9.3132\times {10}^{-10}{V}^2$. While the LO is turned off, measured voltage variance is ${\sigma }_c^2=1.11\times {10}^{-5}{V}^2=82.5{\delta }^2$, and the total measured voltage variance is ${\sigma }_m^2=1.61\times {10}^{-4}{V}^2=1200.7{\delta }^2$. Since the requirement of passing GoF test, the system always works under safety condition $c_1<c_2$ in Equation (8), hence the min-entropy is determined by the middle of the distribution, i.e., $H_{min}=-{log}_2\left(\mathrm{erf}\left(\delta /2\sqrt{2}{\sigma }_q\right)\right)=6.39$ bits. Therefore, the rest 12-3-6 = 3-bit LSB is doubtful for security aspect, and its influence should be eliminated by post-processing. To make our scheme more reserved, we keep five bits per signal from the highest non-zero MSB as the pre-processed data in precision extension. Hence, the output has a precision of 5 + (8 −1) = 12 bit per signal, while the generation rate by number is 1/8 of the original sampling rate, i.e., 25 M samples per second.

We compare our scheme in generation rate with traditional method of uniform distribution QRNG plus inverse CDF conversion post-processing. Under the condition of same implementation settings, namely sampling rate $f_s$, sampling precision n and min-entropy (extractable quantity of randomness) $H\left(x\right)$, traditional method can generate $f_{s}n$ raw data, and around $f_{s}n·H\left(x\right)$ final data in uniform distribution with estimated entropy $H\left(x\right)$, thus the generation rate of k-bit Gaussian distributed number is $f_{s}n·H\left(x\right)/k$. On the other hand, our scheme provide $f_s$ raw data, and around $f_s/4$ or $f_s/8$ final data of Gaussian distributed numbers. Considering the practical condition $n=12$, $H\left(x\right)=0.6\sim 0.8$ and $k=12\sim 32$, $n·H\left(x\right)/k$ and $1/4$ are approximately at the same order of magnitude. Hence, the generation rate of two schemes are leveled, but our scheme has avoided the enormous time cost to calculate the accurate Gaussian distributed value, or space cost to store the huge library of inverse CDF conversion in post-processing [55].

4.2. Test Results

Normality Tests

Initially, the random sequences after post-processing should pass the normality test. The fitting result is shown in Figure 5. Random sequences also pass several goodness of fit tests, the test result is shown in

Table 2. The result of commonly used matlab functions for normality test. Random sequences are normalized to $N(0,1)$ for the convenience of matlab tests, so there are no terms of variance test.

Function	Mean	AD Test	JB Test	t-Test
Calculated result	−$3.6066\times {10}^{-4}$	p = 0.4788	p = 0.3678	p = 0.2023
Confidence Interval	[$-0.0036$, 0.0036]	NULL	NULL	NULL
Hypothesis value	$H=0$	$H=0$	$H=0$	$H=0$
Status	Pass	Pass	Pass	Pass

.

We calculate the $3\sigma$ threshold of bias $e\left(n\right)$ and auto-correlation $a_k\left(n\right)$ under Gaussian distribution. $3\sigma$ criterion is a rough threshold indicating the bias $e\left(n\right)$ and auto-correlation coefficient $a_k\left(n\right)$ of a finite sample from ideal random sequence, should only exceed the reference by a probability of $1-\mathrm{erf}\left(\frac{3}{\sqrt{2}}\right)=0.3\%$.

The $3\sigma$ criterion originates from the central limit theorem (CLT). The traditional description of CLT indicates that, the summation $S_n$ of a large amount of i.i.d. variables $\left\{X_i\right\}$ should always have asymtotic behavior to Gaussian distribution:

$$Z_n=\frac{S_n-E\left(S_n\right)}{\sqrt{D\left(S_n\right)}}\to N(0,1),$$

(19)

where $S_n=\frac{1}{n}{\sum }_i{X}_i$.

As long as we can derive the mean value $\mu$ and variance ${\sigma }^2$ of certain test statistic, the $3\sigma$ threshold is determined. These two statistics can be described as:

$$\begin{array}{cc}\hfill e\left(n\right)& =\frac{1}{N}\sum _{i=1}^N(s_i-\overline{s}),\hfill \\ \hfill a_k\left(n\right)& =\frac{{\sum }_{i=1}^N(s_i-\overline{s})(s_{(i+k)\mathrm{mod}\mathrm{N}}-\overline{s})}{{\sum }_{i=1}^N{(s_i-\overline{s})}^2},\hfill \end{array}$$

(20)

while for the Gaussian distribution, there exists: $\overline{s}=\mu =0,\frac{1}{N}{\sum }_{i=1}^N{(s_i-\overline{s})}^2={\sigma }^2=1$. Hence the simplification of Equation (20) is:

$$\begin{array}{cc}\hfill e\left(n\right)& =\frac{1}{N}\sum _{i=1}^N{s}_i,\hfill \\ \hfill a_k\left(n\right)& =\frac{{\sum }_{i=1}^N{s}_i{s}_{(i+k)\mathrm{mod}\mathrm{N}}}{{\sum }_{i=1}^N{s}_i^2},\hfill \end{array}$$

(21)

One can easily derive that for Gaussian distribution, bias follows distribution $e\left(n\right)\sim N(\mu ,\frac{1}{n}{\sigma }^2)$, while auto-correlation follows distribution $a_k\left(n\right)\sim N(0,\frac{1}{n})$ (and free from k), both of which are normal distribution. We utilize the threshold of $3\sigma$ criterion to test our Gaussian distribution QRNG, and the result of auto-correlation coefficient is shown in Figure 6.

In addition, since no test suites for Gaussian distributed random numbers are proposed, we converted some random sequences into uniform distribution for randomness test. The conversion is done by CDF method discussed in Appendix B, and result is shown in

Table 3. The result of NIST-STS test after CDF conversion of Gaussian distributed random sequences. For the test environment of significance level $\alpha =0.01$ and block number of $n=400$, The p-value should be over $p=0.01$ threshold (uniformity version over $0.0001$) and proportion should be within the range $[0.9750,1]$.

Test Name	p-Value	Proportion	Status
Frequency	0.811993	394	Success
Block Frequency	0.719747	396	Success
Cumulative Sums	0.785103(KS)	395.5(avg)	Success
Runs	0.270275	396	Success
Longest Run	0.788728	397	Success
Rank	0.375313	396	Success
FFT	0.272297	395	Success
Non-overlapping	0.647530(KS)	394(avg)	Success
Overlapping	0.830808	396	Success
Universal	0.451234	393	Success
Approx. Entropy	0.739918	397	Success
Excursions	0.726852(KS)	392(avg)	Success
Excursions Var.	0.670396(KS)	395(avg)	Success
Serial	0.589359(KS)	392.5(avg)	Success
Complexity	0.124115	392	Success

.

5. Conclusions

We proposed a QRNG scheme generating random numbers with a Gaussian distribution based on vacuum fluctuation of a quantum state, a theoretically proved Gaussian distributed random source. We analyzed the impacts of practical issues in the QRNG system, including sampling range, resolution and depth of the sampling device, along with the optimization method. A novel and flexible post-processing method is proposed, inspired from the classical RNG scheme, to extend the precision of a single number to 12, or even over 20 bits, where the property of Gaussian distributed PDF and the auto-correlation coefficient is maintained at the cost of generation rate. The generated random sequence simultaneously pass normality test focusing on distribution, as well as widely acknowledged NIST-STS test suite of randomness (after converted to uniformly distributed sequences). We experimentally demonstrated the scheme based on vacuum shot noise with conventional devices at a generation rate 25M of sample per second.

Our scheme takes advantage of the Gaussian distributed profile of quantum random sources. Impacts of practical issues, strictly monitored by the GoF test, could not essentially alter the profile of Gaussian distribution, and consequently eliminated by the designed fast post-processing method. We provide a novel method generating Gaussian distribution random numbers effectively.

We have to admit that, despite other QRNG schemes would face consumption in uniform-Gaussian conversion procedure, the generation rate in our system is questionably inadequate for a practical continuous-variable QKD system [52,53]. However, we demonstrate the feasibility of such kind of QRNG, and two factors limiting the generation rate, both of which have huge space to improve. Firstly, the frequency in our system is quite low, due to the limitation of detector bandwidth in our system. By using a balanced detector and ADC with higher bandwidth, the generation rate can be further improved by at least one order of magnitude. Secondly, despite the amplification, the amplified vacuum fluctuation is still too small compared with the sampling range, thus fails to make full use of the ADC.

Security is another issue that is extremely significant to the QRNG system. Despite that, we have estimated the min-entropy in the trusted device scenario and operate accordingly, it is not totally clear whether the MSB method in the post-processing phase eliminates the classical noise substantially. The security issue of Gaussian-distributed schemes needs further discussion, and we are keen on tracing related works.