A simple low-latency real-time certifiable quantum random number generator

Quantum random numbers distinguish themselves from others by their intrinsic unpredictability arising from the principles of quantum mechanics. As such they are extremely useful in many scientific and real-world applications with considerable efforts going into their realizations. Most demonstrations focus on high asymptotic generation rates. For this goal, a large number of repeated trials are required to accumulate a significant store of certifiable randomness, resulting in a high latency between the initial request and the delivery of the requested random bits. Here we demonstrate low-latency real-time certifiable randomness generation from measurements on photonic time-bin states. For this, we develop methods to certify randomness taking into account adversarial imperfections in both the state preparation and the measurement apparatus. Every 0.12 s we generate a block of 8192 random bits which are certifiable against all quantum adversaries with an error bounded by 2−64. Our quantum random number generator is thus well suited for realizing a continuously-operating, high-security and high-speed quantum randomness beacon.

for certifying randomness against classical side information; however, the implementation of this framework is casedependent.

A. Construction of the classical model
To construct the classical model for a trial, we first temporarily assume that Eve's classical side information about the state is independent of Eve's partial knowledge of the input and measurement at the trial. Then we remark and explain that the classical model in the general case where Eve's classical side information about the state is classically correlated with Eve's partial knowledge of the input and measurement can be easily obtained as the convex closure of the model C constructed with the above independence assumption. The model C can be constructed by considering all the possible combinations of Eve's classical side information about the state and Eve's partial knowledge of the input and measurement, which is detailed in the following three steps: Step 1. We consider the case where at each trial the photon source emits a single photon with probability 1 and the probabilities P X and P Z of selecting the X-and Z-bases are fixed. Since each measurement has two possible outcomes +1 or −1 and the probabilities P X and P Z are fixed, the set of probability distributions of the input (basis choice, denoted by I) and output (measurement outcome, denoted by O) achievable at each trial is fully determined by the region accessible by the expectation values σ X and σ Z of the measurements σ X = σ · n X and σ Z = σ · n Z on an arbitrary qubit state. So, we need only to characterize the region accessible by σ X and σ Z , which is detailed in the next paragraph.
Given two arbitrary unit vectors n X and n Z , we have that (n X + n Z ) · (n X − n Z ) = 0. Hence, the two Pauli matrices σ + = σ · n + and σ − = σ · n − , where n + = (n X + n Z )/ √ 2 + 2n X · n Z and n − = (n X − n Z )/ √ 2 − 2n X · n Z , are orthogonal, that is, Tr(σ † + σ − ) = 0. So, there exists another Pauli matrix σ such that the set {1l 2 , σ + , σ − , σ }, where 1l 2 is the identity matrix of 2 × 2, is an orthonormal basis for the Hilbert space of a qubit. Then, in view of the fact that the purity Tr(ρ 2 ) of an arbitrary qubit state ρ cannot be larger than 1, we obtain the inequality Therefore, if the angle between the two directions n X and n Z is fixed to be θ, the region R 1,θ accessible by σ X and σ Z is described as We observe that the region R 1,θ has a shape of a filled ellipse, so it is convex. In practice the angle θ cannot be precisely characterized and can vary from trial to trial; however, by calibration we can ensure a bound on the deviation such that |θ − π/2| ≤ δ at each trial. Then, the region R 1 accessible by σ X and σ Z becomes the union of the convex regions R 1,θ accessible by all θ satisfying |θ − π/2| ≤ δ. Since the angle θ takes a value from a convex set and for each possible θ the region R 1,θ is convex, so is the region R 1 . Specifically, in the first quadrant where σ X ≥ 0 and σ Z ≥ 0, the boundary of the region R 1 is described as follows: when 0 ≤ σ X ≤ sin(δ), σ Z = 1; when sin(δ) < σ X < 1, ( σ X + σ Z ) 2 2+2 sin(δ) 2−2 sin(δ) = 1; when σ X = 1, 0 ≤ σ Z ≤ sin(δ).
By symmetries, we can also obtain the boundary of the region R 1 in the other quadrants, which is illustrated in Supplementary Figure 1.
We now can remark on the possibility of dealing with positive operator-valued measures (POVMs). Suppose that the two measurement directions n X and n Z are in the X-Z plane of the Bloch sphere. Without loss of generality, let n X = (cos(θ X ), sin(θ X )) and n Z = (sin(θ Z ), cos(θ Z )), where |θ X | ≤ δ X , |θ Z | ≤ δ Z and δ X + δ Z ≤ δ in order to satisfy the bound on the misalignment. Then, the expectation values ( σ · n X , σ · n Z ) are in the region R 1 . Because the region R 1 is convex, the expectation values ( M X , M Z ) of M X = k ω X,k σ · n X,k and M Z = l ω Z,l σ · n Z,l , where for all k (l), ω X,k ≥ 0 and |θ X,k | ≤ δ X (ω Z,l ≥ 0 and |θ Z,l | ≤ δ Z ), and k ω X,k = l ω Z,l = 1, are also in the region R 1 . In other words, if we randomly choose to perform the POVM {(1l 2 + M X )/2, (1l 2 − M X )/2} or {(1l 2 + M Z )/2, (1l 2 − M Z )/2} on a qubit, the resulting distribution of the output O conditional on the input I is in the set determined by the convex region R 1 . Therefore, in our construction of the classical model not only the assumed projective measurements but also the POVMs that can be decomposed into convex combinations of these projective measurements are allowed. We emphasize that the set of POVMs considered above does not include all the possible POVMs that are close to the ideal projective measurement σ 1 or σ 3 . Also, our current method does not allow Eve to hold purifications of POVMs, as it assumes that Eve can access classical side information but not quantum side information about the measurement performed at a trial. How to completely handle general POVMs with our method deserves further investigation in future work.
Step 2. We consider the case where at each trial the probability of emitting a single photon from the photon source is at least q 1,lb and the probabilities of selecting the X-and Z-bases are fixed. For the same reason as in Step 1, we need only to characterize the region accessible by the expectation values at the two measurement bases. Recall that in the single-photon subspace the two measurements are described by Pauli operators σ X and σ Z , while in the general case the two measurements, denoted by M X and M Z , may be difficult to be exactly specified but are always block-diagonal with respect to various photon-number subspaces. Because of the block-diagonal structure of M X and M Z , the expectation values M X and M Z can be expressed as convex combinations of those expectation values conditional on a varying number of photons emitted from the source. In Step 1, we already characterize the region R 1 accessible by the expectation values M X and M Z conditionally on the emission of a single photon at a trial. Moreover, since each measurement has two possible outcomes +1 or −1, the region R =1 accessible by M X and M Z conditionally on the emission of multiple photons or zero photon at a trial is at most as large as the square R s . = {−1 ≤ M X ≤ 1, −1 ≤ M Z ≤ 1}. For the model construction, we set R =1 = R s . That is, when the source emits multiple photons or zero photon we constrain the accessible region R =1 in a device-independent way, regardless of the state prepared or measurements performed. Given that the single-photon probability is at least q 1,lb , the region R accessible by M X and M Z is thus described as Here the equality in the second line of the above equation is due to the fact that R 1 ⊆ R s . Since both regions R 1 and R s are convex and R 1 ⊆ R s , the region R, as illustrated in Supplementary Figure 1, is convex. We remark that Eq. (4) implies the statement in the main text that the classical model C can be expressed as a convex combination of sub-models constructed conditionally on a varying number of photons emitted from the source.
Before moving to the next step, we recall that the expectation values M X and M Z are implicitly conditional on Eve's classical side information. Therefore, by setting R =1 = R s we pessimistically allow perfect correlations between Eve's classical side information and all the possible measurement outcomes obtained when the source emits multiple photons or zero photon. In this way, we choose to not certify the randomness contributed by these outcomes. As a consequence, our security analysis is robust against photon-number splitting attacks.
Step 3. We consider the general case where at each trial the probability of emitting a single photon from the photon source is at least q 1,lb and the probabilities P X and P Z of selecting the X-and Z-bases are not exactly known but bounded such that P X + P Z = 1 and the imbalance τ = (P X − P Z )/2 satisfies τ lb ≤ τ ≤ τ ub . The set of input distributions is thus given as S = {P(I) : τ lb ≤ (P(I = X) − P(I = Z))/2 ≤ τ ub and P(I = X) + P(I = Z) = 1}.
It is easy to see that S is a closed interval and so is convex. From Step 2, we know that the set of probability distributions of the output O conditional on the input I is characterized as follows: (6) The set C O|I is convex as each member of C O|I is obtained from a member of the convex set R by a linear transformation. To obtain the joint distribution of I and O, we consider an arbitrary combination of the input distribution P(I) ∈ S and the input-conditional distribution of outputs P(O|I) ∈ C O|I . Accordingly, the classical model C (that is, the set of distributions of I and O achievable at a trial) is given as C = P(I, O) : P(I) ∈ S, P(O|I) ∈ C O|I , and for each i and o, P( The model C is convex because both sets S and C O|I are convex. We emphasize that to construct the classical model C, we exploit only the fact that the measurements M X and M Z are block-diagonal with respect to various photon-number subspaces, and we do not use any information about the normalized state in each photon-number subspace. The measurement operators and the distribution of photon numbers emitted from the source, as well as the input distribution, may depend on some auxiliary degrees of freedom, such as spatial mode, which are not used for encoding information. If Eve can access these auxiliary degrees of freedom such as by a side-channel attack, Eve could be able to manipulate the measurements performed and the state prepared. In this case, we take advantage of the fact that in practice the measurements are block-diagonal with respect to various states for the auxiliary degrees of freedom manipulable by Eve (see the Methods section of the main text for detailed explanations). Therefore, as long as the adversarial imperfections are characterized by the parameters (q 1,lb , δ, τ lb , τ ub ), the classical model C constructed as above, as well as the PEFs constructed below, is valid, meaning that our security analysis is robust against such side-channel attacks.
We finally remark that in the general case where Eve's classical side information about the state is classically correlated with Eve's partial knowledge of the input and measurement at a trial, the joint distribution of the input I and output O of the trial (conditional on the classical side information E implicitly) can be written as where for each λ, ω Λ=λ ≥ 0, and λ ω Λ=λ = 1. Here, the random variable Λ characterizes the classical correlation between Eve's classical side information about the state and Eve's partial knowledge of the input and measurement, and for each possible value Λ = λ, the joint distribution P Λ=λ (I, O) is in the model C of Eq. (7). Therefore, the classical model for the trial will become the convex closure of the classical model C as constructed above. Since the model C is convex, the classical trial model in the general case is still C.

B. Construction of PEFs
We observe that the classical model C obtained above is convex, so we need only to ensure that the PEF inequality in Eq. (1) of the main text is satisfied at each extremal distribution of C according to Lemma 14 of Ref. [1]. However, the set C has an infinite number of extremal distributions, and so a PEF needs to satisfy an infinite number of linear constraints, where each extremal distribution of C imposes a linear constraint. In order to effectively construct PEFs by numerical optimization, instead we can find a convex polytope P that has a finite number of extreme points and contains the classical model C of interest. In this way, if a non-negative function F c : (i, o) → F c (i, o) satisfies all the PEF inequalities imposed by the finite number of extreme points of P, then the function F c : (i, o) → F c (i, o) will also satisfy all the PEF inequalities imposed by each distribution in C (see Lemma 14 of Ref. [1]). Therefore, the function F c : (i, o) → F c (i, o) is a PEF for the classical model C. Below we will first explain how to find such a convex polytope P and then present the explicit optimization problem for constructing PEFs.
For the classical model C constructed according to Eq. (7), each member P(I, O) of C is given by the product of a member P(I) in a closed interval S and a member P(O|I) in the convex set C O|I . Hence, in order to find a convex polytope P such that C ⊆ P, it suffices to find another convex polytope P O|I such that C O|I ⊆ P O|I . To achieve this goal, in view of the construction of C O|I given in Eq. (6) we just need to find a convex polygon R such that R ⊆ R . Since R = q 1,lb R 1 + (1 − q 1,lb )R s (see Eq. (4)) and the region R s is a square, we need only to find a convex polygon R 1 such that R 1 ⊆ R 1 . Since the boundary of the region R 1 is well characterized (see Eq. (3) and Supplementary  Figure 1), we can find such a convex polygon R 1 as detailed in the next paragraph.
We choose R 1 to be a convex polygon circumscribed around R 1 . For this, we observe that in each quadrant the boundary of R 1 is composed by two line segments and one elliptic curve. For example, in the first quadrant the two line segments are shown in the first and third lines of Eq. (3) and the curve is shown in the second line of Eq. (3), that is, σ X + σ Z = 2 + 2 sin(δ) cos(ψ) and σ X − σ Z = 2 − 2 sin(δ) sin(ψ), where ψ ∈ (−φ, φ) and φ = arcsin( (1 − sin(δ))/2). Consequently, if we choose (2K − 1) points in the elliptic curve of each quadrant, for example, in the first quadrant the points characterized by ψ = 0, ±φ/K, ±2φ/K, ..., ±(K − 1)φ/K, and draw tangent lines of the curve at these points, then the region R 1 bounded by these tangent lines and the line segments in the boundary of R 1 is a (8K)-sided convex polygon circumscribed around R 1 . Therefore, following the arguments of the previous paragraph, we can obtain a convex polytope P circumscribed around the convex set C. Moreover, with the increase of the value for K, the resulting convex polytopes R 1 and P become better outer-approximations of their respective convex regions R 1 and C. In our numerical construction of PEFs detailed below, we choose K = 64. Now we are ready to present the optimization problem for constructing PEFs. We emphasize that to ensure the validity of the subsequent security analysis using PEFs, a PEF for a trial needs to be constructed before preforming the trial. Since the classical model for each trial is the identical C, we can use the same PEF, F c : (i, o) → F c (i, o), for each trial, and construct this PEF before the experiment. With the same PEF for each trial, from Theorem 1 of the main text we know that if the classical s -smooth conditional min-entropy certified in a successful experiment with n trials (without taking into account the probability of success) is at least − log 2 (p). Therefore, the quantity in the left-hand side of Eq. (9) can be informally interpreted as the amount of classical s -smooth min-entropy in the outputs O n certifiable conditionally on the inputs I n and the classical side information E. Hence, before the experiment we need to choose a PEF such that the expected value of the quantity in the left-hand side of Eq. (9) is as large as possible. For this purpose, we temporarily assume that the experimental trials are independent and identically distributed (i.i.d.). Suppose that the actual distribution of each trial's input and output is P exp (I, O) and that the classical trial model C is contained in a convex polytope P with L extremal distributions P l (I, O), l = 1, 2, ..., L. Then, to construct a well-performing PEF with a proper power β c , before the experiment we can solve the following optimization problem: The above maximization is over both the power β c > 0 and the PEF F c : (i, o) → F c (i, o). When fixing the power β c , the objective function is a strictly concave function of the PEF F c and the constraints are linear in F c , so there is a unique maximum. This maximum and the corresponding PEF can be found by the sequential quadratic programming [1,2]. After solving the maximization over the PEFs with a fixed power β c , the maximization over β c can be solved by a local search. We would like to emphasize that as long as a function F c : satisfies the constraints in the optimization problem of Eq. (10), the function is a valid PEF with power β c for the classical trial model C. That is, the validity of the constructed PEF depends only on the validity of the classical model C. This has two implications: First, the PEF constructed by the above maximization is valid no matter what the actual distribution P exp (I, O) of the trial's input and output is. Hence, in practice it is not necessary to learn the actual distribution P exp (I, O) before the experiment, and in the optimization problem of Eq. (10) we can replace P exp (I, O) by a distribution P est (I, O) estimated based on the calibration data that is available before the experiment (see Supplementary Note 4 for details). Second, even if the trials are not i.i.d. (but are characterized by the same classical model C), the constructed PEF is still valid.
Finally we define the quantity Here the optimization is over both the power β c > 0 and the PEF F c : (i, o) → F c (i, o) for the classical trial model C.
In view of the arguments in the previous paragraph, R c can be identified as the asymptotic randomness-generation rate per trial witnessed by all possible PEFs when each trial has the same distribution P exp (I, O) and is characterized by the same model C, see Refs. [1,2] for detailed justifications. Moreover, the asymptotic rate R c is optimal, see Sect. VI D of Ref. [2] for a general proof. We can bound the asymptotic rate R c from below by constructing a convex polytope P with L extremal distributions such that P ⊇ C. Specifically, using the same numerical method as the above for constructing PEFs we can find a lower bound on R c . With the increase of the value for L, we obtain a tighter lower bound on R c . To obtain the results presented in Figure 1 of the main text, we chose L = 1024. This choice is justified by our observation that the obtained lower bound on R c would not increase further if a larger value for L is used.

Supplementary Note 2. Quantum probability estimation for our experiment
In this supplementary note, we present the details of how to implement the framework of quantum probability estimation for certifying randomness against quantum side information [3,4]. We note that for device-independent randomness generation, the excellent finite-data efficiency of quantum probability estimation has been demonstrated in Refs. [3][4][5]. Here we would like to apply quantum probability estimation for low-latency randomness generation with partially characterized quantum devices. In Supplementary Note 2 A and Supplementary Note 2 B, we construct the quantum model for each trial and the corresponding quantum estimation factors (QEFs) for the experiment of our interest, respectively. At each trial, the input distribution P(I) is not exactly known but is chosen from the convex set S of Eq. (5). We emphasize that quantum probability estimation is a general framework for certifying randomness against quantum side information; however, the implementation of this framework is case-dependent.
In this supplementary note, we identify classical systems with random variables. We denote quantum systems by upper-case letters in sans serif font (such as D, E). Throughout this supplementary note, E plays a distinguished role as the system carrying the quantum side information. We denote the identity operator on a quantum system by 1l. Quantum states, unless otherwise specified, refer to normalized quantum states.

A. Construction of the quantum model
To construct the quantum model for a trial, we first temporarily assume that Eve's quantum side information about the state is independent of Eve's partial knowledge of the input and measurement at the trial. Then we remark and explain that the quantum model in the general case where Eve's quantum side information about the state is classically correlated with Eve's partial knowledge of the input and measurement can be easily obtained as the convex closure of the model Q constructed with the above independence assumption. The model Q can be constructed by considering all the possible combinations of Eve's partial knowledge of the input and measurement (i.e., all the possible input distributions and all the possible measurements on the quantum system D of the device used in our experiment) and Eve's quantum side information about the state (i.e., all the possible quantum states ρ DE shared between the device's system D and Eve's system E before the trial). The construction of Q is detailed in the following two steps: Step 1. We consider the case where at each trial the photon source emits a single photon with probability 1. In this case, the device's system D 1 is described as a qubit, and the two measurements on the device are described as σ X = σ · n X and σ Z = σ · n Z where |n X · n Z | ≤ sin(δ). Each measurement is a projection-valued measure (PVM) on the system D 1 . Let Π D1,o|i be the input-dependent PVM element for the trial output O = o given the trial input I = i, where i = X or Z and o = +1 or −1. Suppose that before the trial, the joint state of the systems D 1 and E is ρ D1E . Given a trial input i, the state of the trial output O and Eve's system E is induced by performing the corresponding measurement on the initial state ρ D1E . That is, where Tr D1 is the partial trace over the system D 1 and 1l E is the identity operator on the system E. If we assume that the input I at each trial can be arbitrarily chosen with probability distribution P(I) in the convex set S of Eq. (5), then the joint state of the trial input I, the trial output O and Eve's system E is a classical-quantum state of the form where P(I) ∈ S. Since the initial state ρ D1E of the systems D 1 and E is arbitrary, the measurement directions satisfy |n X · n Z | ≤ sin(δ), and the input distributions P(I) is arbitrarily chosen from the convex set S, we obtain a set of classical-quantum states ρ IOE . This set is the quantum model Q 1 for each trial under the assumption that the device's system is a qubit. Note that the quantum model Q 1 in general is not convex but has a convex closureQ 1 . To construct QEFs for the trial model Q 1 , considering Property 2 of Ref. [3] we need only to characterize the extremal states in the convex setQ 1 and the corresponding set of sandwiched Rényi powers, which are detailed in the next two paragraphs respectively. It is obvious to see that the initial state ρ D1E should be a pure state in order to ensure the resulting state ρ IOE is extremal inQ 1 . Since the system D 1 is a qubit, by the Schmidt decomposition a pure state |ψ D1E can be written as Suppose that the two measurement directions n X and n Z are in the X-Z plane of the Bloch sphere, and let |0 and |1 be the two eigenstates of the measurement operator σ Z = σ · n Z with eigenvalues +1 and −1 respectively. Then, in the basis {|0 , |1 } the two measurements σ X and σ Z , as well as their corresponding PVM elements, are represented by real matrices. In general, the Schmidt basis states |0 and |1 can be not in the X-Z plane of the Bloch sphere, that is, the density matrix ρ D1E = |ψ D1E ψ| can be not a real matrix in the basis {|0 , |1 }. However, since all PVM elements Π D1,o|i with i = X, Z and o = ±1 are represented by real matrices in the basis {|0 , |1 }, the density matrix ρ D1E must be also represented by a real matrix in the same basis in order to ensure the resulting state ρ IOE is extremal. Therefore, the Schmidt basis states |0 and |1 should be in the X-Z plane of the Bloch sphere. In other words, to obtain the extremal states of the convex closureQ 1 of the quantum model Q 1 , we need only to consider pure initial states |ψ D1E with the Schmidt basis states being in the X-Z plane of the Bloch sphere. In particular, suppose that the Schmidt basis states are |0 = cos(γ)|0 + sin(γ)|1 and |1 = sin(γ)|0 − cos(γ)|1 , where γ ∈ [0, π), and let n X = (sin(θ), − cos(θ)) and n Z = (0, 1) be two unit vectors in the X-Z plane of the Bloch sphere, where θ ∈ [π/2 − δ, π/2 + δ]. Then, the extremal states of the convex closureQ 1 of the quantum model Q 1 can be expressed as where the sub-normalized states ρ E (i, o) of E conditional on I = i and O = o are represented in the Schmidt basis {|0 , |1 } as follows: Now we are ready to characterize the set RQ 1,e of sandwiched Rényi powers according to the extremal states of the convex setQ 1 . As the PVM elements Π D1,o|i for all i and o are rank-1 and the initial state ρ D1E is pure, we observe that the sub-normalized states ρ E (i, o) in Eq. (15) are rank-1. With this observation in mind and in view of the definition of sandwiched Rényi powers in Eq. (4) of the main text, direct calculation shows that where we abbreviate the sandwiched Rényi power . Then, the set RQ 1,e is given as We have two remarks on the constructions of the quantum model Q 1 and the corresponding QEFs as follows. First, in view of the convexity ofQ 1 ⊇ Q 1 and the joint convexity of sandwiched Rényi powers (Proposition 3 of Ref. [6]), the convex closure of the set RQ 1,e , which is denoted byRQ 1 ,e , satisfies the following property: Each element in the set R Q1 of sandwiched Rényi powers according to all the states in the model Q 1 is bounded from above by an element in the setRQ 1,e . We denote the above relation by R Q1 ≤RQ 1,e . Therefore, if a non-negative function satisfies the QEF inequalities imposed by all the elements of the setRQ 1,e , then this function will also satisfy the QEF inequalities imposed by all the elements of the set R Q1 . In other words, if a function is a QEF with respect to the sandwiched Rényi powers in the setRQ 1,e , then this function is a QEF with respect to the sandwiched Rényi powers in the set R Q1 according to the quantum model Q 1 . This observation applies to a general quantum model. In this work, we will construct the QEFs for the quantum model Q of our interest by characterizing the set of sandwiched Rényi powers according to the convex closure of Q. Second, suppose that the two measurement directions n X and n Z are in the X-Z plane of the Bloch sphere. Without loss of generality, let n X = (cos(θ X ), sin(θ X )) and n Z = (sin(θ Z ), cos(θ Z )), where |θ X | ≤ δ X , |θ Z | ≤ δ Z and δ X + δ Z ≤ δ satisfying the bound on the misalignment. Denote the PVM element for the trial output o given the trial input i by Π θ X D1,o|i=X or Π θ Z D1,o|i=Z , where the dependence on the measurement direction is made explicit. Then, given arbitrary |θ X | ≤ δ X and |θ Z | ≤ δ Z , an arbitrary initial state ρ D1E , and an input distribution P(I) in the convex set S of Eq. (5), the resulting classical-quantum state ρ is in the quantum model Q 1 . Now consider the two POVMs where for all k (l), ω X,k ≥ 0 and |θ X,k | ≤ δ X (ω Z,l ≥ 0 and |θ Z,l | ≤ δ Z ), and k ω X,k = l ω Z,l = 1. The classicalquantum state induced by performing the above two POVMs selected randomly at a trial can be expressed as a convex mixture where each state ρ is induced by performing the assumed projective measurements selected randomly at a trial and thus in the quantum model Q 1 . So, the classical-quantum state ρ IOE in Eq. (19) is in the convex closurē Q 1 of the quantum model Q 1 . Hence, the quantum model induced by the above POVMs is generally larger than the quantum model Q 1 induced by the assumed projective measurements, but both models have the same convex closure. In view of the first remark above and the remarks below Eq. (5) of the main text, a QEF works well for both a quantum model and its convex closure. Therefore, in our construction of the quantum model not only the assumed projective measurements but also the POVMs of the form in Eq. (18) are allowed. We emphasize that the set of POVMs considered above does not include all the possible POVMs that are close to the ideal projective measurement σ 1 or σ 3 . Also, our current method does not allow Eve to hold purifications of POVMs, as it assumes that Eve can access classical side information but not quantum side information about the measurement performed at a trial. How to completely handle general POVMs with our method deserves further investigation in future work.
Step 2. We consider the general case where at each trial the probability of emitting a single photon from the photon source is at least q 1,lb . In view of the fact that the measurements on the quantum system D of the device are blockdiagonal with respect to various photon-number subspaces, without loss of generality the states ρ DE shared between the device's system D and Eve's system E before the trial can be restricted to those with the same block-diagonal structure. That is, ρ DE has the direct-sum (denoted by ⊕) structure where D 1 is a qubit system corresponding to the case that the photon source emits a single photon, and D 2 is a quantum system of arbitrary dimension corresponding to the case that the photon source emits multiple photons or zero photon. Moreover, the sub-normalized state ρ D1E satisfies the constraint that Tr(ρ D1E ) ≥ q 1,lb in order to be consistent with the single-photon probability bound allowed. For the model construction, we conservatively assume that both the state of the subsystem D 2 and the measurements on D 2 are completely uncharacterized except that each measurement has two outcomes ±1. As the dimension of the system D 2 can be arbitrarily high, the measurements on D 2 can be restricted to projective ones. Let Π D2,o|i be the input-dependent PVM element on D 2 for the trial output O = o given the trial input I = i. Recall that the input-dependent PVM elements on D 1 are denoted by Π D1,o|i and that they are partially characterized. Suppose that before the trial, the joint state of the systems D and E is ρ DE of the form in Eq. (20). Then, given a trial input i, the state of the trial output O and Eve's system E induced by performing the corresponding measurement is where ρ Suppose that the input distribution P(I) is arbitrarily chosen from the convex set S of Eq. (5). The joint state of the trial input I, the trial output O and Eve's system E is a classical-quantum state of the form where ρ (k) OE|i , k = 1, 2, and P(I) ∈ S. To construct the quantum model, let us introduce the normalized statesρ Since the joint of the systems D and E has an arbitrary initial state ρ DE of the form in Eq. (20), all the possible states ρ (1) IOE achievable at a trial form the set Q 1 studied in Step 1. Denote the set of all possible statesρ (2) IOE by Q 2 and its convex closure byQ 2 . We observe that Q 2 is convex and so Q 2 =Q 2 , considering the following facts: 1) both the sub-normalized state ρ D2E and the measurements on the system D 2 are arbitrary; 2) the input distribution P(I) is arbitrarily chosen from the convex set S of Eq. (5). When the trace Tr(ρ D1E ) is fixed to a value t, all the possible states ρ IOE achievable at a trial form the set In the case of our interest, the trace Tr(ρ D1E ) takes a value at least q 1,lb , so the set of classical-quantum states ρ IOE achievable at a trial is the union of all the sets Q t where t ≥ q 1,lb . This union is the quantum model Q for each trial assuming that the single-photon probability is bounded from below by q 1,lb . Specifically, the quantum model Q can be expressed as Note that the set Q in general is not convex. We also remark that Eq. (24) implies the statement in the main text that the quantum model Q can be expressed as a direct sum of sub-models constructed conditionally on a varying number of photons emitted from the source. To construct QEFs for the trial model Q, we need to characterize the set R Q of sandwiched Rényi powers according to all the states in Q, which is detailed in the next three paragraphs. We first exploit the block-diagonal structure of the classical-quantum states in the quantum model Q to simplify the characterization of R Q . For all ρ IOE ∈ Q and for all trial inputs i and trial outputs o, we have that where ρ (k) for all inputs i and outputs o are simultaneously block-diagonal. Considering this fact and in view of the definition in Eq. (4) of the main text, direct calculation shows that the sandwiched Rényi powers satisfy In other words, a sandwiched Rényi power according to the normalized state ρ IOE = ρ IOE is given by the probabilistic mixture of sandwiched Rényi powers according to the normalized statesρ IOE is the corresponding state in the set Q k , k = 1, 2, any sandwiched Rényi power according to Q can be written as a probabilistic mixture of sandwiched Rényi powers according to Q k , k = 1, 2. Specifically, let the sets of sandwiched Rényi powers according to Q 1 and Q 2 be R Q1 and R Q2 respectively. Then we have considering that ω = Tr(ρ D1E ) ≥ q 1,lb . Recall that we already characterized the set R Q1 in Step 1. Precisely speaking, we characterized a convex setRQ 1 ,e such that R Q1 ≤RQ 1,e , see the paragraph below Eq. (17). Second, we characterize the set R Q2 . In particular, to construct QEFs we need only to characterize the set RQ 2,e of sandwiched Rényi powers according to the extremal states of the convex set Q 2 =Q 2 . So, we need to first characterize the extremal states ofQ 2 . Since there are only two binary-outcome measurements on the device, by the lemma of Ref. [7] or Lemma 2 of Ref. [8] there exists an orthonormal basis for the system D 2 where the PVM elements of the two measurements are simultaneously block-diagonal with blocks of size at most 2 × 2. Therefore, the system D 2 can be assumed to be a qubit in order to obtain the extremal states ofQ 2 . So we reduce to the case considered in Step 1. The only difference lies in that the two projective measurements on D 2 are conservatively assumed to be completely uncharacterized such that the angle θ between the two measurement directions on the Bloch sphere can be arbitrary. Hence, the extremal states ofQ 2 have the form of Eq. (14), where the sub-normalized states ρ E (i, o) are given in Eq. (15) with the parameters µ ∈ [0, π/2], γ ∈ [0, π) and θ ∈ [0, π]. According to these extremal states, Eve's probability of correctly guessing the output O given the input I at a trial can be as large as 1, meaning that no randomness can be certified according to the quantum model Q 2 =Q 2 . Recall that the quantum model Q 2 is constructed under the conservative assumption that both the state of the subsystem D 2 and the measurements on D 2 are completely uncharacterized. We thus like to remark that with the above assumption, we choose to not certify the randomness contributed by the outcomes obtained from measurements on the subsystem D 2 , and so our security analysis is robust against photon-number splitting attacks. Now we can characterize the set RQ 2,e . Similar to the set RQ 1 ,e in Eq. (17), the set RQ 2,e can be expressed as has the form of Eq. (16) with µ ∈ [0, π/2], γ ∈ [0, π) and θ ∈ [0, π]}.
In view of the convexity ofQ 2 and the joint convexity of sandwiched Rényi powers (Proposition 3 of Ref. [6]), the convex closure of the set RQ 2 ,e , which is denoted byRQ 2 ,e , satisfies the following property: Each element in the set R Q2 is bounded from above by an element in the setRQ 2 ,e . We denote the above relation by R Q2 ≤RQ 2,e . In the same sense, below we extend the bound relation "≤" to two arbitrary sets.
Third and now we are ready to characterize the set R Q . In view of the results in the above two paragraphs, we have R Q k ≤RQ k ,e , k = 1, 2, and so Here, for a fixed ω the set ωRQ 1,e + (1 − ω)RQ 2,e is convex as both setsRQ 1 ,e andRQ 2 ,e are convex. By Eqs. (17) and (28), we also have that RQ 1,e ⊆ RQ 2,e and soRQ 1,e ⊆RQ 2 ,e . Therefore, we get In view of Eqs. (27), (29) and (30), we thus get the relation whereRQ = q 1,lbRQ 1,e + (1 − q 1,lb )RQ 2,e is a convex set. The above relation is useful for constructing QEFs, as detailed in Supplementary Note 2 B. We emphasize that to construct the quantum model Q, we exploit only the fact that the measurements M X and M Z are block-diagonal with respect to various photon-number subspaces, and we do not use any information about the normalized state in each photon-number subspace. The measurement operators and the distribution of photon numbers emitted from the source, as well as the input distribution, may depend on some auxiliary degrees of freedom, such as spatial mode, which are not used for encoding information. If Eve can access these auxiliary degrees of freedom such as by a side-channel attack, Eve could be able to manipulate the measurements performed and the quantum state prepared. In this case, we take advantage of the fact that in practice the measurements are block-diagonal with respect to various states for the auxiliary degrees of freedom manipulable by Eve (see the Methods section of the main text for detailed explanations). Therefore, as long as the adversarial imperfections are characterized by the parameters (q 1,lb , δ, τ lb , τ ub ), the quantum model Q constructed as above, as well as the QEFs constructed below, is valid, meaning that our security analysis is robust against such side-channel attacks.
We finally remark that in the general case where Eve's quantum side information about the state is classically correlated with Eve's partial knowledge of the input and measurement at a trial, the joint state of the trial input I, the trial output O and Eve's system E can be written as where for each λ, ω Λ=λ ≥ 0, and λ ω Λ=λ = 1. Here, the random variable Λ characterizes the classical correlation between Eve's quantum side information about the state and Eve's partial knowledge of the input and measurement, and for each possible value Λ = λ, the classical-quantum state ρ Λ=λ IOE is in the model Q of Eq. (24). Therefore, the quantum model for the trial will become the convex closureQ of the model Q as characterized above. In view of the remarks below Eq. (5) of the main text, no matter whether the quantum model for a trial is Q orQ the set of QEFs constructed is the same. For this reason, without loss of generality the quantum model for a trial is characterized by the model Q of Eq. (24). But we emphasize that the construction of QEFs detailed in the next section and so our security analysis with respect to quantum side information work in the general case where Eve's quantum side information about the state is classically correlated with Eve's partial knowledge of the input and measurement at a trial.

B. Construction of QEFs
Considering the relation in Eq. (31) and the joint convexity of sandwiched Rényi powers (Proposition 3 of Ref. [6]), to construct QEFs for the quantum model Q it suffices to consider the QEF inequalities imposed by the extreme points of the convex setRQ. These extreme points are expressed as RQ ,e = q 1,lb RQ 1 ,e + (1 − q 1,lb )RQ 2,e , where RQ 1,e and RQ 2,e are in the sets RQ 1,e and RQ 2,e of Eq. (17) and Eq. (28) respectively. However, the setRQ has an infinite number of extreme points, and so a QEF needs to satisfy an infinite number of linear constraints, where each extreme point ofRQ imposes a linear constraint. In order to effectively construct QEFs by numerical optimization, instead we can find a convex polytope P that has a finite number of extreme points and contains the convex setRQ. In this way, if a non-negative function F q : (i, o) → F q (i, o) satisfies all the QEF inequalities imposed by the finite number of extreme points of P, then the function F q : (i, o) → F q (i, o) will also satisfy all the QEF inequalities imposed by each member of the convex setRQ (see Property 2 of Ref. [3]). Therefore, the function is a QEF for the quantum model Q. Below we will first explain how to find such a convex polytope P and then present the explicit optimization problem for constructing QEFs.
Since the extreme points ofRQ can be expressed as in Eq. (33), to find a convex polytope P such thatRQ ⊆ P we need only to find another two polytopes P 1 and P 2 such that RQ 1,e ⊆ P 1 and RQ 2,e ⊆ P 2 . Because the elements R(ρ E ) of the two sets RQ 1,e (Eq. (17)) and RQ 2,e (Eq. (28)) are expressed as the same function of parameters γ, θ and µ (they are different only in the range of the parameter θ), we can find the two polytopes P 1 and P 2 in the same way. Our strategy for finding P 1 and P 2 is outlined as follows: We express the vector R(ρ E ) as the entry-wise product of another two vectors, and then find out two convex polytopes which contain the accessible regions of these two vectors respectively. Consequently, all the accessible vectors R(ρ E ) are contained in a polytope whose extreme points are given by the entry-wise products of the extreme points of the two polytopes found above. The details of our strategy are presented in the next two paragraphs.

(34)
Recall that the input distribution P(I) is arbitrarily chosen from a closed interval S of Eq. (5), which is a convex polytope. So, secondly we need to find a convex polytope that contains all the accessible vectors q(ρ E ). As q(ρ E , i, o) = r(ρ E , i, o) αq with α q > 1 is a convex function of r(ρ E , i, o) and r(ρ E , i, o) is non-negative, it suffices to find a convex polytope that contains all the accessible vectors r(ρ E ) . = (r(ρ E , X, +1), r(ρ E , X, −1), r(ρ E , Z, +1), r(ρ E , Z, −1)). To see this, suppose that r(ρ E ) = g ω g r g with ω g ≥ 0 and g ω g = 1. Then q(ρ E ) ≤ g ω g q g , where the entries of q g are given by the α q -th powers of the entries of r g . As p is a vector with non-negative entries, if a non-negative function F q : (i, o) → F q (i, o) satisfies all the QEF inequalities imposed by p · q g , then the function F q : (i, o) → F q (i, o) satisfies the QEF inequality imposed by p · q(ρ E ). Therefore, the convex polytope with extreme points q g suffices for our purpose of constructing QEFs. To find a convex polytope containing all accessible r(ρ E ), we introduce the vectors v 1 = (1/2, 1/2), v 2 = (1/2, −1/2) and v µ = cos 2/αq (µ), sin 2/αq (µ) . Then, the expressions in Eq. (34) can be rewritten as Let us introduce another vector v γ,θ = cos(2γ + θ), cos(2γ) . As detailed in the next paragraph, for all values of γ, θ and µ we can find convex decompositions where ω k,γ,θ , ω m,µ ≥ 0, k ω k,γ,θ = m ω m,µ = 1, and all a k , b k , c m , d m are constants independent of γ, θ and µ. Then from Eqs. (35) and (36) we obtain that Eq. (37) tells us that all accessible vectors r(ρ E ) can be written as convex combinations of the vectors r km . Therefore, the convex polytope defined as the convex closure of the points r km contains all the accessible vectors r(ρ E ).
In this paragraph we find the convex decompositions in Eq. (36). We first consider the set of vectors v γ,θ = cos(2γ +θ), cos(2γ) . To obtain the elements R(ρ E ) in both sets RQ 1,e (Eq. (17)) and RQ 2,e (Eq. (28)), γ ∈ [0, π). So, the set of vectors v γ,θ is the same as the set of vectors v γ,θ = cos(2γ + θ/2), cos(2γ − θ/2) . Let σ X = cos(2γ + θ/2) and σ Z = cos(2γ − θ/2). Then the set of vectors v γ,θ (v γ,θ ) is described by the ellipse ( σ X + σ Z ) 2 2 + 2 cos(θ) which corresponds to Eq. (2) for constructing the classical trial model. For the elements R(ρ E ) in the set RQ 1,e of Eq. (17) (or in the set RQ 2,e of Eq. (28)), the parameter θ takes values in the interval [π/2 − δ, π/2 + δ] (or in the interval [0, π], corresponding to the case δ = π/2). Hence, the set of vectors v γ,θ is contained in the convex set R 1 specified below Eq. (2). We note that the shape of the convex set R 1 depends on the value of δ. In Supplementary Note 1 B we already found a convex polygon R 1 such that R 1 ⊆ R 1 . So, we can use the extreme points of R 1 as the vectors (a k , b k ) for the convex decompositions of the vectors v γ,θ . Specifically, for the vectors v γ,θ used for characterizing the set RQ 1 ,e where θ ∈ [π/2 − δ, π/2 + δ], we choose R 1 be a convex polygon with K 1 = 512 extreme points, while for the vectors v γ,θ used for characterizing the set RQ 2,e where θ ∈ [0, π] (that is, δ = π/2), R 1 becomes R s , a square with 4 extreme points (1, 1), (1, −1), (−1, 1) and (−1, −1), and so we choose R 1 to be R s . Secondly, we consider the set of vectors v µ = cos 2/αq (µ), sin 2/αq (µ) . To obtain the elements R(ρ E ) in both sets RQ 1,e (Eq. (17)) and RQ 2,e (Eq. (28)), µ ∈ [0, π/2]. So, the set of vectors v µ is represented by a curve S µ in the first quadrant starting at the point (1, 0) and ending at the point (0, 1). To obtain a convex polygon which contains the curve S µ , we choose (K 2 − 1) points in the curve S µ , for example, the points characterized by µ = φ/(2K 2 ), 2φ/(2K 2 ), ..., (K 2 − 1)φ/(2K 2 ), and draw tangent lines of the curve S µ at these points. The region bounded by these tangent lines and the line connecting the two points (1, 0) and (0, 1) is a convex polygon circumscribing around the curve S µ . Therefore, with the extreme points of this convex polygon, we can express all the vectors v µ by convex decompositions. In our numerical construction of QEFs detailed below, we choose K 2 = 64. We note that with the increase of the values for K 1 and K 2 , the resulting convex polygons become better outer-approximations of the regions accessible by the vectors v γ,θ and v µ respectively. Now we are ready to present the optimization problem for constructing QEFs. We emphasize that to ensure the validity of the subsequent security analysis using QEFs, a QEF for a trial needs to be constructed before preforming the trial. Since the quantum model for each trial is the identical Q, we can use the same QEF, F q : (i, o) → F q (i, o), for each trial, and construct this QEF before the experiment. With the same QEF for each trial, from Theorem 2 of the main text we know that if the quantum s -smooth conditional min-entropy certified in a successful experiment with n trials (without taking into account the probability of success) is at least − log 2 (p). Therefore, the quantity in the left-hand side of Eq. (39) can be informally interpreted as the amount of quantum s -smooth min-entropy in the outputs O n certifiable conditionally on the inputs I n and the quantum side information in E. Hence, before the experiment we need to choose a QEF such that the expected value of the quantity in the left-hand side of Eq. (39) is as large as possible. For this purpose, we temporarily assume that the experimental trials are i.i.d. Suppose that the actual distribution of each trial's input and output is P exp (I, O) and that the set R Q of sandwiched Rényi powers according to the quantum model Q is contained in a convex polytope P with L extreme points R l , l = 1, 2, ..., L, where the entries of R l are R l (i, o) for different inputs i = X or Z and outputs o = +1 or −1. We note that the set R Q and the polytope P as well as the extreme points of P all depend implicitly on the power parameter β q and the parameters (q 1,lb , δ, τ lb , τ ub ) which characterize the adversarial imperfections in the experiment. To construct a well-performing QEF with a proper power β q , before the experiment we can solve the following optimization problem: The above maximization is over both the power β q > 0 and the QEF F q : (i, o) → F q (i, o). When fixing the power β q , the objective function is a strictly concave function of the QEF F q and the constraints are linear in F q , so there is a unique maximum. This maximum and the corresponding QEF can be found by the sequential quadratic programming, the same algorithm as that for constructing PEFs. After solving the maximization over the QEFs with a fixed power β q , the maximization over β q can be solved by a local search. We would like to emphasize that as long as a function F q : (i, o) → F q (i, o) satisfies the constraints in the optimization problem of Eq. (40), the function is a valid QEF with power β q for the quantum trial model Q. That is, the validity of the constructed QEF depends only on the validity of the quantum model Q. This has two implications: First, the QEF constructed by the above maximization is valid no matter what the actual distribution P exp (I, O) of the trial's input and output is. Hence, in practice it is not necessary to learn the actual distribution P exp (I, O) before the experiment, and in the optimization problem of Eq. (40) we can replace P exp (I, O) by a distribution P est (I, O) estimated based on the calibration data that is available before the experiment (see Supplementary Note 4 for details). Second, even if the trials are not i.i.d. (but are characterized by the same quantum model Q), the constructed QEF is still valid.
Finally we define the quantity Here the optimization is over both the power β q > 0 and the QEF F q : (i, o) → F q (i, o) for the quantum trial model Q. In view of the arguments in the previous paragraph, R q can be identified as the asymptotic randomness-generation rate per trial witnessed by all possible QEFs when each trial has the same distribution P exp (I, O) and is characterized by the same model Q, see Refs. [3,4] for detailed justifications. Moreover, the asymptotic rate R q is optimal, see Sect. 6.5 of Ref. [4] for a general proof. We can bound the asymptotic rate R q from below by constructing a convex polytope P with L extreme points such that P ⊇ R Q . Specifically, using the same numerical method as the above for constructing QEFs we can find a lower bound on R q . With the increase of the value for L, we obtain a tighter lower bound on R q . To obtain the results presented in Figure 1 of the main text, we chose L = 675840. This choice is justified by our observation that the obtained lower bound on R q would not increase further if a larger value for L is used.

Supplementary Note 3. Imperfections in our experiment
The imperfections to be characterized depend on the specific implementation of the QRNG scheme and can be bounded by calibrating the photon source and measurement apparatus in real time. The lower bound q 1,lb on the single-photon probability for constructing the models and performing the corresponding security analyses can be obtained by monitoring the average photon number per trial given the photon-number distribution. The bounds τ lb and τ ub on the imbalance τ = (P X − P Z )/2 depend on how the measurement basis at a trial is selected. For the active-selection scheme, a random number from some pseudo-or physical randomness source is used and so a calibration of this randomness source is needed. For the passive-selection scheme like in our experiment, beam splitters are used and so a calibration of their splitting ratios is needed. The calibration of the misalignment angle δ also depends on the setup. For the active-selection scheme with polarization encoding, the measurement directions are set by a polarization rotator driven by some electric or magnetic field. Because of the fluctuation of the electric or magnetic field in practice, the measurement directions set by the rotator can vary from trial to trial. Therefore, we need to calibrate the electric or magnetic field applied. For the passive-selection scheme with time-bin encoding as in our experiment, the measurement directions are determined by the splitting ratios of the two beam splitters in the interferometer used as well as by the efficiency mismatch of the two detectors employed (see our characterization detailed below). Therefore, we need to calibrate the splitting ratios and the efficiency mismatch. The imperfections in the measurement apparatus of our experiment are characterized in Supplementary Note 3 A, while the imperfections in our photon source are in Supplementary Note 3 B.

A. Imperfections in measurement apparatus
We perform measurements on photonic time-bin states with an unbalanced Mach-Zehnder interferometer (MZI) as shown in Supplementary Figure 2 (1). The MZI has one input spatial mode, which can transmit two different temporal modes-the early and late time-bin modes labelled as e and l respectively, and has two output spatial modes labelled as a and b respectively. The difference in photon transit time between the two unbalanced paths of the MZI matches the separation between the two time-bin modes e and l. Therefore, depending on which path an input photon in the early (late) time-bin mode takes, the photon can arrive at the output of the MZI in the early or middle time-bin mode (in the middle or late time-bin mode). So, the output photon can be in three different temporal modes, the early, middle and late time-bin modes labelled as e, m and l respectively. Considering both the temporal and spatial modes, the MZI is described by an optical device with two input modes and six output modes as shown in Supplementary  Figure 2 (2). We remark that for detecting photons at the two output spatial modes, threshold detectors (which cannot resolve photon number) of finite efficiency are employed. Below we derive the measurement operators in the single-photon subspace and characterize their dependence on imperfections considering the following two cases: 1) The two detectors employed have the same detection efficiency, and 2) the two detectors at the output ports a and b have efficiencies η a and η b respectively. We emphasize that for constructing the classical and quantum models and performing the corresponding security analyses, our method requires a characterization of the measurement operators in the single-photon subspace but not those operators in the multiphoton subspace.  (1) is the description in the experiment. The MZI is composed of two beam splitters, BS1 and BS2, and has two input temporal modes (the early and late time-bin modes labelled as e and l) and two output spatial modes (labelled as a and b). The right panel (2) is the hypothetical description which helps to illustrate the relations between the input and output modes. Here we separate different combinations of the spatial and temporal modes for the input and output of the MZI. Note that the output can be in three different temporal modes, the early, middle and late time-bin modes (labelled as e, m and l).
Case 1. We assume that the detectors employed are trusted such that their detection efficiencies are the same and also independent of any degree of freedom of the incoming photons. As a consequence, the no-click events due to detection inefficiency can be treated as an attenuation of the incoming optical pulse before performing measurements with threshold detectors of perfect efficiency. (The effect of optical-pulse attenuations will be discussed in Supplementary Note 3 B.) So, without loss of generality we assume that each of the detectors employed has the perfect efficiency.
In the ideal case, each of the two beam splitters BS1 and BS2 in the MZI has the 50:50 splitting ratio. However, in practice it is difficult to ensure the ideal splitting ratio. Suppose that the two beam splitters BS1 and BS2 have transmission coefficients t 1 and t 2 respectively. The coefficients t 1 and t 2 may depend on the auxiliary degrees of freedom of the single photon which are not used for information encoding but manipulable by Eve. As explained in the Methods section of the main text, we take advantage of the assumption that the coherent superposition of states for an auxiliary degree of freedom manipulable by Eve does not play a role throughout the measurement process. (This assumption can be justified if in the setup for time-bin measurements there is no quantum interference between any pair of states for the auxiliary degree of freedom manipulable by Eve-which is true in practice as we think.) Therefore, each measurement operator on a single photon is block-diagonal with respect to various states for the auxiliary degrees of freedom, where each block is described by a qubit measurement. As a consequence, we need only to characterize the dependence of the qubit measurement operators on t 1 and t 2 , even if t 1 and t 2 are manipulable by Eve. For this, we utilize the relations between the creation operators of the input and output modes of a beam splitter. Specifically, if a beam splitter has a transmission coefficient t with two input modes 1 in and 2 in and two output modes 1 out and 2 out , where a photon in the mode 1 in is transmitted to the mode 1 out with probability t and is reflected to the mode 2 out with probability (1 − t), then the creation operators of the output modes are related with those of the input modes bŷ By exploiting such relations for the two beam splitters BS1 and BS2 in the MZI, we obtain the measurement operators, represented in the single-photon input-mode basis {|e , |l }, as follows: where r 1 = 1 − t 1 and r 2 = 1 − t 2 . Note that the sum of the above four operators is equal to the identity operator, meaning that only click events are observed when the detection efficiency is perfect and the input contains a single photon. Suppose that t i = 1/2 + τ i , i = 1, 2, where τ 1 and τ 2 characterize the imperfections of the two beam splitters BS1 and BS2. Then, from Eq. (43) we get the following inequalities where 1l 2 is the identity matrix of 2 × 2. Therefore, if we assume that the events e and l are due to the time-bin Z-basis measurement and that the events (m, a) and (m, b) are due to the superposition X-basis measurement on the input time-bin qubit, then the probabilities P X and P Z of selecting the X-and Z-bases at a trial satisfy that When the event e or (m, a) happens, we specify the outcome as +1; when the event l or (m, b) happens, the outcome is −1. With these outcome assignments, the (unnormalized) input-dependent measurement operators in the X-and Z-bases are where σ 1 = 0 1 1 0 and σ 3 = 1 0 0 −1 are two of the Pauli matrices. In view of the above operators and according to Property 3 of Ref. [3], for constructing PEFs and QEFs we can use the classical and quantum models induced by the (unnormalized) input-dependent measurement operators The corresponding (normalized) input-dependent measurement operators are where θ = cos −1 − 2τ 2 / 4τ 2 2 + (1 − 4τ 2 1 )(1 − 4τ 2 2 ) . So, the two measurements along the X-and Z-bases can be represented by two directions n X and n Z on the Bloch sphere, where n X = (sin(θ), 0, cos(θ)) and n Z = (0, 0, 1). Note that even if the imperfections τ 1 and τ 2 are not exactly known but can only be bounded, we can still bound the probabilities P X and P Z as well as the angle θ. Therefore, we allow Eve to manipulate the imperfections in the two beam splitters used via side-channel attacks, as long as the manipulations satisfy the assumed imperfection bounds.
In a word, the imperfection in the first beam splitter BS1 determines the imbalance, τ = (P X − P Z )/2, between the probabilities P X and P Z of selecting the X-and Z-bases, and the imperfections in both beam splitters BS1 and BS2 affect the misalignment angle, δ = max θ |θ − π/2|. Case 2. We assume that the two detectors at the output ports a and b have efficiencies η a and η b , which are independent of the temporal modes of the incoming photons but can be controlled by Eve via manipulating other degrees of freedom of the photons. For the same reason as in Case 1, we need only to characterize the dependence of the qubit measurement operators on the transmission coefficients t 1 and t 2 and on the efficiencies η a and η b , even if these parameters are manipulable by Eve. A practical detector with efficiency η can be described as a beam splitter with transmission coefficient η followed by an ideal detector with perfect efficiency [9]. In view of this fact and by exploiting the relations in Eq. (42) for the two beam splitters in the MZI as well as for the beam splitters describing the two practical detectors, we get the following measurement operators represented in the single-photon input-mode basis {|e , |l } τ lb = 1 2 By the same argument as in the first case and in view of Eq. (50), for constructing PEFs and QEFs we can use the classical and quantum models induced by the (normalized) input-dependent measurement operators where θ = cos −1 (χ 3 / χ 2 1 + χ 3 3 ). So, like in the first case, the two measurements along the X-and Z-bases can be represented by two directions n X and n Z on the Bloch sphere, where n X = (sin(θ), 0, cos(θ)) and n Z = (0, 0, 1). Note that even if the imperfections τ 1 , τ 2 and ∆η/η are not exactly known but can only be bounded, we can still bound the probabilities P X and P Z as well as the angle θ. Therefore, in this case we allow Eve to manipulate, via sidechannel attacks, the efficiency mismatch in addition to the imperfections in the two beam splitters used, as long as the manipulations satisfy the assumed imperfection bounds. Also, in contrast to the first case, here both the imbalance τ = (P X − P Z )/2 and the misalignment angle δ = max θ |θ − π/2| depend on all the considered imperfections τ 1 , τ 2 and ∆η/η.
To calibrate the imperfections in our measurement apparatus, we send a weak pulse in the early time-bin mode, the same as the pulse used for randomness generation, into the MZI. The pulse can arrive at the output of the MZI in the early (e) or middle (m) time-bin mode, and then be detected by a threshold detector at either the output port a or b. We perform the calibration in the absence of Eve's manipulation, and we assume that the parameters to be calibrated are stable over time. Specifically, we assume that the two detectors at the output ports a and b have efficiencies η a and η b , and that the two beam splitters BS1 and BS2 in the MZI have transmission coefficients t 1 and t 2 , respectively. In order to handle the cases where the parameters to be calibrated fluctuate over time or can be manipulated by Eve via side-channel attacks, we later set conservative bounds on the calibrated imperfections for performing security analysis. The probabilities of detections at each possible output mode (a combination of the temporal and spatial modes) depend on the parameters t 1 , t 2 , η a and η b , as shown in Supplementary Table I(1). Therefore, we can estimate these parameters based on experimental results. In particular, we used 360 s of experimental data for calibration. Denote the number of counts by the detector a for the early time-bin mode e of the output over 360 s of experiment time by n(a, e); the numbers of counts at the other output modes are denoted in a similar way. These results are presented in Supplementary Table I( 2). In view of the fact that the number of counts at an output mode is proportional to the corresponding detection probability, we obtain the following relations Plugging the experimental results in Supplementary Table I(2) into the above relations, we get t 1 = 0.462 and t 2 = 0.531. So, the splitting ratios between the reflection and transmission coefficients for the two beam splitters BS1 and BS2 in the MZI are 53.8:46.2 and 46.9:53.1, respectively. In addition, we can estimate the mismatch between the efficiencies η a and η b of the two detectors employed. Specifically, in view of Supplementary Table I(1) we have n(a, e) + n(a, m) n(b, e) + n(b, m) = η a η b Combining Eqs. (57) and (58) and then plugging the experimental results in Supplementary Table I(2), we get that η a /η b = 1.024 and so ∆η/η = 0.0119. This estimated efficiency-mismatch ratio is consistent with the detection efficiencies η a = (60 ± 2)% and η b = (59 ± 2)% measured before the experiment. This estimation suggests the potential presence of efficiency mismatch in our measurement setup. Now we can translate the above estimated imperfections in our measurement apparatus into the imperfection bounds required by our security-analysis method. If we consider the first case discussed above, where the efficiency mismatch is assumed to be absent, then according to Eqs. (45) and (48) we estimate that |τ | ≤ 0.038 and δ ≤ 3.565 • . If we consider the second case where the effect of efficiency mismatch is included, then according to Eqs. (55) and (56) we estimate that |τ | ≤ 0.041 and δ ≤ 3.513 • . In view of these estimated bounds, we conservatively assume that |τ | ≤ 0.06 and δ ≤ 6 • in our security analysis.
We remark that there are several previous works on how to relax assumptions on measurement devices. In particular, the measurement-device-independent (MDI) QRNG schemes proposed in Refs. [10,11] have the potential to be resistant to all adversarial imperfections in the measurement apparatus including efficiency mismatch. However, a MDI scheme requires that the states prepared are fully characterized, which is difficult to promise in practice. There is also a dimension-witness-based QRNG scheme as developed in Ref. [12]. Such scheme does not require characterizing same basis. 3) The simultaneous photon detections at the X-and Z-bases (that is, cross clicks), which include the simultaneous detections at the early and middle time-bin modes, the simultaneous detections at the middle and late time-bin modes, and the simultaneous detections at all the three time-bin modes, are all mapped to the outcome −1 in the X-basis.
Several remarks are as follows. First, we allow Eve to manipulate the distribution of photon numbers emitted from the source, as long as Eve's manipulation satisfies the lower bound q 1,lb used in our security analysis. Second, because of the detection inefficiency or the insertion loss of the MZI in practice, the single-photon or multiphoton components of the optical pulse are not always detected. For security analysis, we make the fair-sampling assumption that the detected sub-ensemble is a fair representation of the whole ensemble. Therefore, the same as zero-photon events, the no-click events due to photon loss do not affect our security analysis but only the randomness-generation rate and latency achieved in practice. Strictly speaking, the no-click events due to photon loss does not affect the validity of our security-analysis method, but they affect the bound q 1,lb used in our model constructions which is detailed as follows. Suppose that the probability of emitting j photons from the photon source is p j , j = 0, 1, 2, . . ., and that the total photon loss in the MZI and in the detection process is ξ, where ξ ∈ [0, 1). Note that ξ is the probability that a single photon is lost in the process of measurement and detection. Then, the detector-click probability conditional on the event that j photons are emitted from the source is p(click|j) = 1 − ξ j . Therefore, the probability of the joint event that a detector clicks and j photons are emitted is p(click, j) = p j (1 − ξ j ). So, the lower bound q 1,lb for model constructions should satisfy where . When j ≥ 2, the factor f j (ξ) is a monotonically increasing function of ξ ∈ [0, 1) (by checking the positivity of the first derivative of f j (ξ) when ξ ∈ [0, 1)). Accordingly, we have where µ = j≥1 jp j is the average photon number in a pulse emitted from the source. In our experiment, µ = 0.0035 and p 1 = 3.48 × 10 −3 (see the above paragraph). So, we can set q 1,lb to be a value as large as 0.993 regardless of the specific photon loss, while in our security analysis we used the more conservative choice q 1,lb = 0.98. Third, our treatment of multiphoton events is different from the treatments in Refs. [16,[19][20][21]. These references rely on a particular post-processing map-a squashing map that assigns a double-click event uniformly randomly to be one of the two single-click events at the same measurement basis and assigns all cross-click events to no detection. After applying the squashing map, one can prove that for some particular measurements (for example, when the X-and Z-bases are mutually unbiased and when there is no detection-efficiency mismatch), the distribution of measurement outcomes according to a multiphoton state can be reproduced by the same measurements on some single-photon state. Therefore, with a squashing map the measurement outcomes according to a multiphoton state can be treated in the same way as those according to a single-photon state, and so the security analysis taking into account multiphoton events can be reduced to the security analysis in the presence of only single-photon events. In contrast, we treat the measurement outcomes according to a multiphoton state in a different way: We explicitly separate the contribution of multiphoton events from that of single-photon events when constructing the models C and Q, and for the multiphoton part we conservatively assume that the measurements are un-characterized as long as each measurement has only two outcomes (after applying our strategy for assigning simultaneous photon detections as specified in the previous paragraph). So, our treatment works for arbitrary measurements, unlike the squashing map [16,[19][20][21] which exists only for ideal measurements without any misalignment or efficiency mismatch. Another advantage is that the outcome assignments for simultaneous photon detections used by us do not consume additional random numbers as compared with the squashing map used in Refs. [16,[19][20][21].
4. Randomness Extraction: Conditional on success, run the TMPS extractor to generate k q (or k c ) random bits with error .
The first step is to set the parameters required for our analysis. To achieve the goal, we set the smoothness error to be s = 0.8 ≈ 4.34 × 10 −20 and the extractor error to be x = 0.2 ≈ 1.08 × 10 −20 . We emphasize that the positive errors s and x need to satisfy that s + x ≤ , but their choices are not unique. The splitting ratio 0.8:0.2 between s and x used by us was not optimized; instead it was chosen heuristically for the same reason as in Ref. [5]. Moreover, to achieve our goal with the TMPS extractor [22,23], the amount of quantum (or classical) s -smooth conditional min-entropy to be certified (without taking into account the probability of success) is 8516 (or 16712) bits. According to Theorem 2 (or Theorem 1) of the main text, the amount of quantum (or classical) s -smooth conditional minentropy certified in a successful QRNG instance (without taking into account the probability of success) is bounded from below by − log 2 (p). Consequently, for our goal it suffices to set the parameter p in Theorem 2 (or Theorem 1) of the main text to be 2 −8516 (or 2 −16712 ).
The second step is calibration based on the data obtained prior to the data used for randomness generation. The data sets for calibration and randomness generation are called the calibration data and the analysis data, respectively. As a result of this step, we determined a well-performing QEF F q (I, O) and its power β q (or a well-performing PEF F c (I, O) and its power β c ) to be used for each executed trial for randomness generation, and fixed the maximum number of trials n q (or n c ) that can be used for each QRNG instance, see Supplementary Note 4 A for details.
In the calibration step, we also estimated that a QRNG instance with a high probability of success requires about 0.1 s runtime (see Supplementary Note 4 A for details). Considering that we have 420 s of data for randomness generation, we thus decided ahead of time to run 4200 instances of our QRNG. We also decided to use the same QEF (or PEF) as well as the corresponding maximum number of trials n q (or n c ) for all instances. That is, when we ran 4200 instances of our QRNG, we performed the calibration step only once before processing all 420 s of data for randomness generation. In this way, the time cost for calibration would not contribute to the latency in randomness generation.
The third step consists of acquiring up to n q (or n c ) trials. After each trial m, we update the running value of the QEF product T q,m = m k=1 F q (i k , o k ) (or the running value of the PEF product T c,m = m k=1 F c (i k , o k )). Both QEFs and PEFs have the advantage that one can stop the experiment early as soon as the running value T q,m (or T c,m ) surpasses the corresponding success threshold 1/ p βq ( 2 s /2) in Theorem 2 (or 1/ p βc s in Theorem 1) of the main text. This early-stopping feature is due to the fact that the constant function F (i, o) = 1 for all inputs i and outputs o of a trial is both a PEF with any positive power for all classical models and a QEF with any positive power for all quantum models, see Refs. [1][2][3][4] for more detailed explanations. If for each possible m, T q,m (or T c,m ) fails to surpass the corresponding success threshold, the instance of our QRNG fails. When we ran the analysis, we found that each instance succeeded before the data block used for the instance was completely processed, see Supplementary Note 4 B for details. Let n q,act (or n c,act ) be the actual number of trials executed.
The fourth and last step consists of running the TMPS extractor with the trial outputs (when the instance succeeds). The extractor input is exactly n q bits (or n c bits) long and consists of the trial outputs padded with zeros to n q bits if n q,act < n q (or padded with zeros to n c bits if n c,act < n c ). The number of seed bits required for running the extractor is determined by the length of the extractor input, the number of random bits to be extracted and the extractor error as instructed in Refs. [2,4,5,24]. In each QRNG instance, to extract 8192 (or 2 × 8192) random bits with error 2 −64 in the presence of quantum (or classical) side information the number of seed bits required is 1031048 (or 1288810). In our numerical implementation of the TMPS extractor, the extraction of 8192 (or 2 × 8192) random bits took about 0.02 s (or 0.04 s) on a cluster with 30 threads. We remark that as the TMPS extractor is a strong extractor, the extractor seed could be reused for the future instance of our QRNG as long as the extractor seed is independent of the other random variables or quantum systems involved in the future instance; however, the seed bits to be reused would have a statistical distance from uniform that needs to be added in order to adjust the soundness error claimed in the future instance.

A. Calibration details
For the calibration purpose, we took 360 s of experimental data (see Supplementary Table II (1)). This calibration data is the same as that used for calibrating the imperfections in our measurement apparatus and shown in Supplementary Table I(2); however, the same data is represented in different forms in Supplementary Table I(2) and in Supplementary Table II(1). With the calibration data shown in Supplementary Table II(1), we can get an estimate of the distribution of inputs and outputs at an experimental trial. The estimated distribution is denoted by P est (I, O), which is given by the relative frequency distribution of the processed results presented in Supplementary Table II(2). Also, based on the calibration data, we observed that every 0.1 s runtime the experiment generates about 47000 trials with clicks. To determine a well-performing QEF F q (I, O) and its power β q (or a well-performing PEF F c (I, O) and its power β c ), we solved the QEF optimization of Eq. (40) (or the PEF optimization of Eq. (10)) with the parameters n = 47000, s = 0.8 ≈ 4.34 × 10 −20 , and the replacement of the distribution P exp (I, O) by the above estimated distribution P est (I, O). The solutions of these two optimization problems are presented in Supplementary Table III. Note that to obtain these solutions, both the quantum and classical models used for each experimental trial are constructed in the presence of the adversarial imperfections specified in Supplementary Note 3. classical) side information. For the estimates, we assume that the quantum device to be used for our QRNG is honest with the distribution P est (I, O) for each trial. Then, by Hoeffding's inequality [25], we found that the success probability is at least 1 − 2 −380 (or 1 − 2 −478 ) using the QEF (or PEF) in Supplementary Table III. Therefore, the goal specified for each QRNG instance can be satisfied within 0.1 s runtime with an extremely high success probability. For randomness generation, we thus decided ahead of time to use the data block obtained in 0.1 s runtime for each instance, and set the maximum number of trials allowed for the instance to an upper bound on the number of trials with clicks collected over 0.1 s runtime. Particularly, we set n q = n c = 50000.

B. Analysis results
In view of the early-stopping property of both QEFs and PEFs, when running the data analysis we found that each QRNG instance succeeded before the data block obtained in 0.1 s runtime was completely processed. Let the number of trials with clicks collected over 0.1 s runtime be n, which depends on the instance slightly but is always smaller than n q = n c = 50000. After the success event happened, we continued to process the remainder of the data observed in each 0.1 s runtime. Finally, we obtained the value for T q,n (or T c,n ) according to the description in Box. 1. By Theorem 2 (or Theorem 1) of the main text, the amount of quantum (or classical) s -smooth conditional min-entropy certified in a successful instance (without taking into account the probability of success) is at least − log 2 (p), where − log 2 (p) ≤ log 2 (T q,n )/β q +log 2 ( 2 s /2)/β q (or − log 2 (p) ≤ log 2 (T c,n )/β c +log 2 ( s )/β c ). Therefore, we set the maximum amount of certifiable quantum (or classical) s -smooth conditional min-entropy to be log 2 (T q,n )/β q + log 2 ( 2 s /2)/β q (or log 2 (T c,n )/β c + log 2 ( s )/β c ). Then we determined the (maximum) number of random bits extractable using the TMPS extractor with error x = 0.2 × 2 −64 in the presence of quantum (or classical) side information. The results for all 4200 instances are shown in Figure 4 of the main text.
In the same way as above, we can vary the soundness error and study the expected number of random bits certifiable from the data observed every 0.1 s runtime. For this, we set the smoothness error to be s = 0.8 and the extractor error to be x = 0.2 . For each value of studied, we solved the QEF optimization of Eq. (40) and PEF optimization of Eq. Then we set the expected value of the maximum amount of certifiable quantum (or classical) s -smooth conditional min-entropy to be log 2 (T q,n )/β q + log 2 ( 2 s /2)/β q (or log 2 (T c,n )/β c + log 2 ( s )/β c ). Finally, we obtained the expected number of random bits extractable using the TMPS extractor with error x = 0.2 in the presence of quantum (or classical) side information. The results are shown in Figure 3 of the main text, where the soundness error varies from 10 −5 to 10 −30 .