\section{Introduction}
\label{intro}
Determining total internal energies and how they change with respect to some perturbation are central objectives in both quantum chemistry and condensed matter physics. Total internal energies can help identify the stability of an atomistic system, and derivatives  with respect to different variables yield a multitude of desired quantities. Recently, machine learning has been used to calculate total electronic energies for toy models \cite{snyder2012finding, ryczko2018convolutional, ryczko2019deep, ryabov2020neural}, molecular systems \cite{behler2007generalized, rupp2012fast, jager2018machine, schutt2018schnet, brockherde2017bypassing, zhou2019toward}, and solid state systems \cite{faber2015crystal, ryczko2018convolutional, mills2019extensive, ryabov2020neural, ryczko2021orbital}. Thus far, the majority of these works have focused on learning Hartree-Fock (HF) or density functional theory (DFT) computed properties. These models can make accurate and rapid predictions, but their utility is limited to the accuracy of the electronic structure method they were trained upon. When studying atomistic systems, one must properly describe electron correlation to achieve high accuracy. Such a treatment leads to computationally expensive electronic structure methods. For molecules, common techniques include M{\o}ller-Plesset perturbation theory (MP2) or coupled-cluster singles, doubles and perturbative triples (CCSD(T)). For solid-state systems, the gold standard is quantum Monte Carlo (QMC) \cite{foulkes2001quantum}.

Several studies have utilized machine learning to predict quantities computed with MP2 or CCSD(T) \cite{wilkins2019accurate, peyton2020machine, margraf2018making, schran2019automated, bogojeski2020quantum, smith2019outsmarting}. In Ref. \cite{peyton2020machine}, a novel representation, called the density tensor representation, was introduced and was used to predict accurate energies and dipoles of small molecules at the MP2 level. In Ref. \cite{bogojeski2020quantum}, the density $\Delta$-DFT approach was introduced and utilized for small molecules. This methodology used the DFT electronic density as input to a machine learning model to predict a difference in energy between CCSD(T) and DFT, allowing for a molecular dynamics simulation with quantum chemical accuracy. Training on the difference between DFT and CCSD(T), rather than the absolute value, allowed for an improved model accuracy when performing inference (predictions).

However, machine learning studies for solid-state systems with accurate electronic structure calculations are currently absent in the literature. This is due to the limitations of typical machine learning implementations, the computational cost of the electronic structure calculations, and the lack of abundant QMC data for solid-state systems. Most current machine learning approaches require large quantities of training data, making it infeasible to train on properties calculated with QMC. However, recent machine learning methodologies have eliminated this problem by focusing on scalar functions, rather than scalar quantities \cite{zhou2019toward, ryabov2020neural, ryczko2021orbital}. Namely, with voxel deep neural networks (VDNNs) it was shown that electron kinetic energies could be predicted within chemical accuracy from only 2 DFT calculations using a technique called voxel deep neural networks (VDNNs) \cite{ryczko2021orbital}.

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{methods.pdf}
    \caption{Visual representation of the methodologies used. In (a) we show how a voxel deep neural network is used to map an electron density computed via density functional theory to energy densities computed via diffusion Monte Carlo. These include the interacting kinetic energy density ($\mathcal{T}_{\text{int}}$), the electron-electron energy density ($\mathcal{V}_{\text{ee}}$), and the electron-ion energy density ($\mathcal{V}_{\text{ei}}$). The Hamiltonian $\hat{H}=T_{\text{int}} + V_{\text{ee}} + V_{\text{ei}}$,  $\Psi_T$ is the trial many body wavefunction, and $\tau$ is imaginary time. In (b) we show a 1D example of a voxel deep neural network. A model is used to map the function $g(x)$ to another function $f(x)$. For every value $x$, there is a fixed field of view that is used to construct inputs for the machine learning model. In (c) we show a visualization of how atomic contributions to the total energy are computed and a schematic of  the kernel ridge regression (KRR) training process.}
    \label{fig:fig_1}
\end{figure*}

In this report, we present and compare two methodologies to predict diffusion Monte Carlo (DMC) total energies using small datasets. The first methodology uses VDNNs to compute QMC energy densities by mapping an electron density computed via DFT to kinetic, electron-electron, and electron-ion energy densities computed with diffusion Monte Carlo (DMC). The second approach uses kernel ridge regression (KRR) with atomic environment vectors as inputs to predict atomic contributions to the DMC total energy. We first compare these methodologies for pristine graphene lattices. Afterwards, we study the generalizability of the KRR method, and predict the energy barrier associated with a Stone-Wales defect in a graphene lattice. Lastly, we move from 2D to 3D materials and use a KRR model to predict DMC total energies of liquid water along a molecular dynamics trajectory. To date, this is the first report that calculates QMC energies for a solid-state system with supervised machine learning.

\section{Methods}
\label{methods}
In Kohn-Sham DFT \cite{kohn1965self} the electron density is written as
\begin{equation}
\label{rho}
    \rho(\bm{r}) = 2\sum_{n}^{\text{occ}}\sum_{\bm{k}}w_{\bm{k}} \psi_{n, \bm{k}}^*(\bm{r})\psi_{n, \bm{k}}(\bm{r}).
\end{equation}
In Equation \ref{rho}, $n$ is the band index, $\bm{k}$ is the k-point, $w_{\bm{k}}$ is the weighting associated with the k-point, and $\psi$ is a Kohn-Sham orbital. For the graphene lattices we used a non-orthorhombic supercell with lattice vectors $a_1 = (12.310, 0, 0)$, $a_2 = (-6.155, 10.661, 0)$, and $a_3=(0, 0, 6.000)$ \AA. For liquid water, we took 32 equally spaced snapshots from a previously used dataset \cite{PhysRevLett.120.036002} that consisted of 1000 atomic configurations from 200 ps of path integral molecular dynamics with 32 water molecules. The 200 ps were run under constant atmospheric pressure and are subdivided amongst 6 replicated runs, each with a different temperature (300 - 350 K). We refer the reader to the original manuscript for more information \cite{PhysRevLett.120.036002}. The QMC electronic contribution to the energy density, as derived in Ref. \cite{krogel2013quantum} is
\begin{equation}
    \label{qmc}
    \mathcal{E}(\bm{r}) = \mathcal{T}(\bm{r}) + \mathcal{V}_{\text{ee}}(\bm{r}) + \mathcal{V}_{\text{ei}}(\bm{r}),
\end{equation}
where the many-body kinetic energy density is
\begin{equation}
\label{mbke}
    \mathcal{T}(\bm{r})=-\frac{1}{2}\sum_{i=1}^N\int d\bm{R}~\delta(\bm{r} - \bm{r}_i)\Psi^*(\bm{R})\nabla_i^2\Psi(\bm{R}),
\end{equation}
and the many-body potential energy densities are
\begin{equation}
    \label{mbpot}
    \mathcal{V}_\text{ab}(\bm{r})=\frac{1}{2}\sum_{i=1}^{N_a}\sum_{j=1}^{N_b}\int d\bm{R}~ \delta(\bm{r} - \bm{r}_i)|\Psi(\bm{R})|^2v_{\text{ab}}(\bm{r}_i, \bm{r}_j).
\end{equation}
In Equations \ref{mbke} and \ref{mbpot}, $\bm{R}\equiv\{\bm{r}_1, \bm{r}_2, \hdots, \bm{r}_N\}$ is the full many-body space, $v_{\text{ab}}$ is the interaction potential between two species, and $\Psi$ is the many-body wavefunction. For $\mathcal{V}_\text{ee}$, the interaction potential $v_\text{ee}$ is the Coulomb potential. For the potential energy density term $\mathcal{V}_\text{ei}$, the interaction potential is a norm-conserving pseudopotential of the Burkatzki-Filippi-Dolg form \cite{burkatzki2007energy}. The many-body wavefunction is constructed from a Slater determinant of Kohn-Sham orbitals and an optimized Jastrow factor. The QMC energy densities are generated in 3 stages. The first stage is a DFT calculation to obtain the Kohn-Sham orbitals needed for the Slater determinant. To perform the DFT calculation, we used Quantum Espresso \cite{giannozzi2009quantum} included with QMCPACK \cite{kim2018qmcpack} with an energy cut-off of 150 Ha. For the graphene lattices we used a $4\times4\times1$ k-point grid and for liquid water we used a $2\times2\times2$ k-point grid. The second stage is the Jastrow factor optimization. We used a one-body Jastrow factor for the nuclei and a two-body Jastrow factor for the electrons. The final stage is evaluating the energy which was done in 3 steps. This included a variational Monte Carlo (VMC) calculation followed by two DMC calculations, each with different time steps. We refer the reader to the SI for more information on the Jastrow optimization and other QMC parameters \cite{si}. All three stages were performed at each $\bm{k}$-point, and the final energy densities were obtained by averaging over all $\bm{k}$-points. All calculations were performed using QMCPACK \cite{kim2018qmcpack}, Nexus \cite{krogel2016nexus}, and a custom Python library \cite{github-qmc}. To compute the atomic contributions to the total energy, the closest energy density voxels to a particular atom (using the minimum image convention) were integrated. For graphene this yielded 50 atomic contributions per calculation and 96 atomic contributions per calculation for liquid water.

When generating a dataset for VDNNs, we used $\approx 5\times10^{5}$, $19\times19\times19$ slices of the DFT electron densities that were collected from the training configurations such that a uniform distribution of the function $\mathcal{F} = \sqrt{\mathcal{T}^2 + \mathcal{V}_{\text{ee}}^2 + \mathcal{V}_{\text{ei}}^2}$ was produced. These parameter choices were found to be optimal in a previous study on the prediction of DFT kinetic energy densities \cite{ryczko2021orbital}. Of these electron density slices, 99\% were used for training and 1\% were used for validation. When training the VDNNs, we map the electron density slices to values of $\mathcal{T},~\mathcal{V}_{\text{ee}}$, and $\mathcal{V}_{\text{ei}}$, as shown in \cref{fig:fig_1}a. We used the convolutional neural network (CNN) from Ref. \cite{ryczko2021orbital}, a batch size of 512, a learning rate of $10^{-5}$, the Adam optimizer \cite{kingma2014adam}, and layer-wise adaptive rate scaling with clipping \cite{you2017scaling}. We use ensemble models, where 5 independent networks were trained simultaneously and total energies were averaged over the models. We trained the networks for 250 epochs, and used the model that had the minimal mean squared error on the validation set. We found that the models began to overfit beyond $\approx250$ epochs. This overfitting can be attributed to the stochastic nature of the underlying training set. Beyond a certain point in training, the model begins to fit to statistical noise. We stopped training prior to this and used the model with the lowest mean squared error on the validation set. Training and inference was done across multiple nodes, each with 4 NVIDIA V100 GPUs. Our codes can be found here \cite{github-qmc}.

For learning the atomic energy contributions, we tried various atomic environment descriptors including atom centred symmetry functions (ACSF) \cite{behler2011atom}, atomic environment vectors (AEV) from the ANI models \cite{smith2018less, smith2019outsmarting, devereux2020extending},  and smooth overlap of atomic positions (SOAP) \cite{bartok2013representing}. Apart from AEV, we used the Python library dscribe \cite{dscribe} for generating the atomic environment descriptors with the default parameters. For AEV, we used the parameters which can be found in torchani \cite{gao2020torchani}. For the machine learning model, we employed Gradient Boosted Decsion Trees (GBDT), Gaussian Process Regression (GPR), artificial neural networks (ANNs), and KRR. We used the Python library scikit-learn \cite{scikit-learn} to implement all models and tried various sets of parameters for each model using grid searches. For both the graphene lattices and liquid water, we found that the KRR models had the lowest mean absolute errors (MAEs).

\begin{figure}[ht]
    \centering
    \includegraphics[width=\linewidth]{graphene_training_curve.pdf}
    \caption{Learning curves for various models along with their respective linear fits. The linear fits allows one to estimate how much data is necessary to achieve a particular accuracy. We estimate for VDNNs that we would need $\approx 400$ calculations to achieve chemical accuracy. However, for KRR with SOAP as the atomic environment vectors, we would only need $\approx 12$ calculations.}
    \label{results_1}
\end{figure}

\begin{figure}[ht]
    \centering
    \includegraphics[width=\linewidth]{sw_graphene_krr.pdf}
    \caption{The energy barrier curve associated with creating a Stone-Wales defect in a pristine graphene lattice. The error bars for QMC show the standard deviation of the calculation. The model KRR-16 has the lowest mean absolute error with 4.23 meV / atom.}
    \label{results_2}
\end{figure}

\section{Results}
\label{results}
The first task is the prediction of DMC total energy for pristine graphene lattices with random atomic perturbations. When building this dataset, we generated 8 atomic configurations where random perturbations of the atomic positions were applied to the equilibrium structure. To generate the random displacements, uniform random numbers were drawn with the range $[-0.1, 0.1]$ and were added to the pristine coordinates. 4 of these calculations were used to generate training/validation data and the remaining 4 were used for testing. For the VDNNs, we trained 4 independent ensemble models, where each model was allowed to see 1, 2, 3 or 4 training configurations. We then computed the total energies of the testing set by integrating the predicted energy densities from the ensemble models. Similar to Ref. \cite{ryczko2021orbital}, a rigid shift was applied to the predictions based on the residuals of the training set. This shift is needed to compensate for the errors from the model that accumulate when integrating on large grids. The shift can be calculated by computing the average error within an energy window. The number of energy windows depends on the structure of the residuals, and therefore differs for each model within the ensemble. We used a grid search to compute the optimal number of energy windows, and found the number of energy windows to be in the range $[8,32]$. In \cref{results_1}, we plot the learning curve for the VDNNs. The best MAE found for VDNNs was 120 meV / atom. For both KRR and VDNNs, the learning curve allows one to estimate the MAE of a model if additional calculations were available. A  linear fit was done in order to obtain these estimates. In \cref{results_1}, these estimates are labelled by the keyword ``fit.'' If one would like to have a VDNN that is capable of making predictions within chemical accuracy, one would need $\approx~400$ additional calculations. For the KRR models, we found that the SOAP representation had the lowest MAE which was 3.40 meV / atom. In comparison to the VDNNs, this is a 35x improvement. In addition, we estimate that we would need $\approx~8$ additional calculations to obtain a model capable of making DMC energy predictions within chemical accuracy. This accuracy, however, is limited to pristine graphene lattices with atomic perturbations.

\begin{figure*}[ht]
    \centering
    \includegraphics[width=\linewidth]{water_dmc.pdf}
    \caption{Top: QMC, DFT, and KRR model prediction energies for liquid water. Bottom: Residuals for DFT and the KRR model. The KRR model shows a 1.38x improvement over DFT with a MAE of 38.1 meV / molecule.}
    \label{results_3}
\end{figure*}

Given the performance of KRR in comparison to VDNNs, we concentrate on analyzing the capabilities of KRR in more complex scenarios. Namely, we are interested in studying the generalizability of the KRR model in regards to calculating the energy barrier curve associated with the introduction of a Stone-Wales defect within a graphene lattice, which was previously investigated with DFT in Ref. \cite{li2005defect}. As shown in \cref{results_2}, by rotating two atoms in the graphene lattice, one transitions from a pristine lattice ($\theta=0^{\circ}$) to a lattice with a Stone-Wales defect ($\theta=90^{\circ}$). For the training set, we started with the 3 configurations shown in \cref{results_2}. These configurations do not have any random atomic perturbations added. We then ran 16 additional DMC calculations of graphene lattice with a single Stone-Wales defect and atomic perturbations. The testing dataset consisted of 10 configurations, where two atoms were rotated such that the $nth$ configuration had an angle $\theta = n \pi/ 18,~n\in[1,9]$. For each configuration, a DFT structural relaxation was performed while holding the two rotated atoms fixed. As done in Ref. \cite{li2005defect}, the bond length between the two rotated atoms was also reduced. We first trained a KRR model ignoring the 16 additional DMC calculations and predicted the energy barrier curve. In \cref{results_2}, this is labelled by ``KRR-0.'' With only 3 total energy reference points, the KRR model does exceptionally well, with a MAE of 7.32 meV / atom. We then added in 2, 4, 8, or 16 additional DMC calculations (labelled ``KRR-2'', ``KRR-4'', ``KRR-8'', and ``KRR-16'' respectively in \cref{results_2}), trained a KRR model, and made predictions for the energy barrier. Interestingly, when comparing the models KRR-0 and KRR-[2,4,8], the MAEs increase (up to 1.8x worse with KRR-4) with the addition of data. We argue that with the addition of new data, there are two areas of chemical space that must be accurately modelled and connected in a physically meaningful way. The models KRR-[2,4,8] do not connect these two areas of chemical space well, leading to higher MAEs. However, with all 16 additional points, the MAE becomes 4.23 meV / atom (1.7x times improvement from KRR-0). The KRR-16 model learned to generalize its knowledge from the additional calculations to make accurate predictions of the energy barrier.

 We define the difference in energy between the highest peak seen in \cref{results_2} and the energy at $\theta = 0$ to be the energy barrier needed to be overcome to transition from a pristine lattice to one with a Stone-Wales defect. We find the energy barrier to be 224(46) meV / atom based on the DMC calculation. DFT underestimates the energy barrier by 17.6 meV / atom and KRR-16 underestimates the energy barrier by only 5 meV / atom.

 We now move from 2D to 3D materials and study the performance of KRR when applied to liquid water. Since there are two atom types, one KRR model was trained to predict atomic energy contributions for H, and another model was trained to predict atomic energy contributions for O. In both cases we found that the AEV representation had the lowest MAE in comparison to other representations. As described in \cref{methods}, there were a total of 32 DMC calculations that were done. Half of these were used in training, and the remaining half were used as a test set. In \cref{results_3}, these data points are highlighted under ``Training'' (in cyan) and ``Testing'' (in light red). For KRR, most of the predictions are within chemical accuracy. The MAE of the test set was 27.6 meV / molecule. In comparison to DFT, which had a MAE of 38.1 meV / molecule, the KRR model shows a 1.38x improvement.

Before concluding, we discuss the limitations and future work of both methods presented here. Firstly, both VDNNs and KRR models are limited to the data they are trained with (as is the case with all machine learning models). Secondly, for machine learning models that predict just scalar quantities (i.e. total energy), there is no error correction that must be done. However, as discussed previously, error from VDNNs accumulates when integrating, and must be accounted for. Thirdly, in VDNNs there is no dependence between neighbouring pixels that are output by the model. The dependence could be applied by outputting neighbouring pixels and averaging over the predictions or in some sort of post-processing routine. Fourthly, due to use of images, rotational invariance is not inherently built-in to the model, but must be learned during training. When constructing diversity in a dataset for VDNNs, one must consider the structure of the electron density, rather than a chemical environment. The two are inter-linked, but one must consider how the electron density images change within a dataset, and use this information when applying it to a new system. Lastly, since VDNNs must make millions of inference calls for a 3D grid, prediction time is relatively slow (on the order of minutes) compared to machine learning methods that output a scalar quantity. This is where predicting atomic contributions to the total energy excels. The KRR models take seconds to train and perform inference, and have rotational and translational symmetry built in. However, our current KRR implementation does not have automatic differentiation. Automatic differentiation with kernel methods has been done previously \cite{gardner2018gpytorch}, and future work of ours involves using these packages to obtain forces at the DMC level. In addition, this work shows that it is possible to perform supervised learning with DMC data. Given that our KRR method trains on atomic contributions to the total energy, we expect that size of datasets needed for training DNNs will be smaller by a factor of $\langle N \rangle$, which is the average number of atoms per calculation. Other future work of ours involves investigating $\Delta$ learning from DFT to DMC atomic contributions to the total energy.

\section{Conclusion}
We have used two machine learning methodologies to predict DMC total energies. The first methodology we used was VDNNs to map electron densities computed with DFT to energy densities calculated with DMC for pristine graphene lattices. After comparing these results with our second methodology, which used KRR and atomic environment vectors to predict atomic contributions to the total energy, we found that the errors from KRR were an order of magnitude lower. Given this result we then further studied the capabilities of KRR to predict atomic contributions to the energy for more complex tasks. This included predicted the energy barrier when introducing a Stone-Wales defect into a pristine graphene lattice and studying liquid water. In both cases, we found that KRR trained on DMC outperforms DFT. In addition, the KRR method requires no initial DFT calculation, but only the computation of atomic environment vectors. This allows one to rapidly infer energies at the DMC level. Since our method trains on atomic contributions to the total energy, we also conclude that dataset sizes in future studies can be reduced by a factor of $\langle N \rangle$, which is the average number of atoms per calculation. This would allow for DNNs to be used to predict DMC energies, furthering the performance showed in this report.

\section{Acknowledgements}
The Authors acknowledge Compute Canada, the Vector Institute for Artificial Intelligence, and the National Energy Research Scientific Computing Center for computational resources. KR and IT acknowledges the National Sciences and Engineering Council of Canada for funding. Work performed by JK (mentorship) was supported by the US Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.

\section{More Information on the QMC Calculations}

\subsection{Jastrow Factors}
The Jastrow factors were represented by B-splines and the coefficients were optimized using the quartic optimizer in QMCPACK \cite{kim2018qmcpack}. The meshes used in the QMC calculations were increased by 20\% and we found the average and corresponding standard deviation of the variance to energy ratio of $0.0205\pm0.0014$. The optimization was run for 100 blocks each consisting of 51200 Monte-Carlo samples.

\subsection{Calculation of the Energy Densities}
In this subsection, we provide additional details for the calculation of the energy densities. The first step is performing variational Monte Carlo (VMC) for 40 blocks with 30 warm-up steps and 10 steps used for averaging. Warm-up steps are steps taken where data is excluded from calculating quantities. The second step was DMC with a larger time step (0.02) for 20 blocks with 20 warm-up steps and 5 steps used for averaging. The final step was DMC with a smaller time step (0.01) for 200 blocks with 20 warm-up steps followed by 10 steps. The total number of Monte-Carlo samples completed for each of the final DMC calculations was 20 blocks  $\times$ 10 steps $\times$ 1024 walkers $\times$ 64 MPI processes $\approx13\times10^6$. The average and corresponding standard deviation of the total energy error across all calculations was $0.426\pm0.083$ meV.

\section{Introduction}
\label{intro}
Determining total internal energies and how they change with respect to some perturbation are central objectives in both quantum chemistry and condensed matter physics. Total internal energies can help identify the stability of an atomistic system, and derivatives  with respect to different variables yield a multitude of desired quantities. Recently, machine learning has been used to calculate total electronic energies for toy models \cite{snyder2012finding, ryczko2018convolutional, ryczko2019deep, ryabov2020neural}, molecular systems \cite{behler2007generalized, rupp2012fast, jager2018machine, schutt2018schnet, brockherde2017bypassing, zhou2019toward}, and solid state systems \cite{faber2015crystal, ryczko2018convolutional, mills2019extensive, ryabov2020neural, ryczko2021orbital}. Thus far, the majority of these works have focused on learning Hartree-Fock (HF) or density functional theory (DFT) computed properties. These models can make accurate and rapid predictions, but their utility is limited to the accuracy of the electronic structure method they were trained upon. When studying atomistic systems, one must properly describe electron correlation to achieve high accuracy. Such a treatment leads to computationally expensive electronic structure methods. For molecules, common techniques include M{\o}ller-Plesset perturbation theory (MP2) or coupled-cluster singles, doubles and perturbative triples (CCSD(T)). For solid-state systems, the gold standard is quantum Monte Carlo (QMC) \cite{foulkes2001quantum}.

Several studies have utilized machine learning to predict quantities computed with MP2 or CCSD(T) \cite{wilkins2019accurate, peyton2020machine, margraf2018making, schran2019automated, bogojeski2020quantum, smith2019outsmarting}. In Ref. \cite{peyton2020machine}, a novel representation, called the density tensor representation, was introduced and was used to predict accurate energies and dipoles of small molecules at the MP2 level. In Ref. \cite{bogojeski2020quantum}, the density $\Delta$-DFT approach was introduced and utilized for small molecules. This methodology used the DFT electronic density as input to a machine learning model to predict a difference in energy between CCSD(T) and DFT, allowing for a molecular dynamics simulation with quantum chemical accuracy. Training on the difference between DFT and CCSD(T), rather than the absolute value, allowed for an improved model accuracy when performing inference (predictions).

However, machine learning studies for solid-state systems with accurate electronic structure calculations are currently absent in the literature. This is due to the limitations of typical machine learning implementations, the computational cost of the electronic structure calculations, and the lack of abundant QMC data for solid-state systems. Most current machine learning approaches require large quantities of training data, making it infeasible to train on properties calculated with QMC. However, recent machine learning methodologies have eliminated this problem by focusing on scalar functions, rather than scalar quantities \cite{zhou2019toward, ryabov2020neural, ryczko2021orbital}. Namely, with voxel deep neural networks (VDNNs) it was shown that electron kinetic energies could be predicted within chemical accuracy from only 2 DFT calculations using a technique called voxel deep neural networks (VDNNs) \cite{ryczko2021orbital}.

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{methods.pdf}
    \caption{Visual representation of the methodologies used. In (a) we show how a voxel deep neural network is used to map an electron density computed via density functional theory to energy densities computed via diffusion Monte Carlo. These include the interacting kinetic energy density ($\mathcal{T}_{\text{int}}$), the electron-electron energy density ($\mathcal{V}_{\text{ee}}$), and the electron-ion energy density ($\mathcal{V}_{\text{ei}}$). The Hamiltonian $\hat{H}=T_{\text{int}} + V_{\text{ee}} + V_{\text{ei}}$,  $\Psi_T$ is the trial many body wavefunction, and $\tau$ is imaginary time. In (b) we show a 1D example of a voxel deep neural network. A model is used to map the function $g(x)$ to another function $f(x)$. For every value $x$, there is a fixed field of view that is used to construct inputs for the machine learning model. In (c) we show a visualization of how atomic contributions to the total energy are computed and a schematic of  the kernel ridge regression (KRR) training process.}
    \label{fig:fig_1}
\end{figure*}

In this report, we present and compare two methodologies to predict diffusion Monte Carlo (DMC) total energies using small datasets. The first methodology uses VDNNs to compute QMC energy densities by mapping an electron density computed via DFT to kinetic, electron-electron, and electron-ion energy densities computed with diffusion Monte Carlo (DMC). The second approach uses kernel ridge regression (KRR) with atomic environment vectors as inputs to predict atomic contributions to the DMC total energy. We first compare these methodologies for pristine graphene lattices. Afterwards, we study the generalizability of the KRR method, and predict the energy barrier associated with a Stone-Wales defect in a graphene lattice. Lastly, we move from 2D to 3D materials and use a KRR model to predict DMC total energies of liquid water along a molecular dynamics trajectory. To date, this is the first report that calculates QMC energies for a solid-state system with supervised machine learning.

\section{Methods}
\label{methods}
In Kohn-Sham DFT \cite{kohn1965self} the electron density is written as
\begin{equation}
\label{rho}
    \rho(\bm{r}) = 2\sum_{n}^{\text{occ}}\sum_{\bm{k}}w_{\bm{k}} \psi_{n, \bm{k}}^*(\bm{r})\psi_{n, \bm{k}}(\bm{r}).
\end{equation}
In Equation \ref{rho}, $n$ is the band index, $\bm{k}$ is the k-point, $w_{\bm{k}}$ is the weighting associated with the k-point, and $\psi$ is a Kohn-Sham orbital. For the graphene lattices we used a non-orthorhombic supercell with lattice vectors $a_1 = (12.310, 0, 0)$, $a_2 = (-6.155, 10.661, 0)$, and $a_3=(0, 0, 6.000)$ \AA. For liquid water, we took 32 equally spaced snapshots from a previously used dataset \cite{PhysRevLett.120.036002} that consisted of 1000 atomic configurations from 200 ps of path integral molecular dynamics with 32 water molecules. The 200 ps were run under constant atmospheric pressure and are subdivided amongst 6 replicated runs, each with a different temperature (300 - 350 K). We refer the reader to the original manuscript for more information \cite{PhysRevLett.120.036002}. The QMC electronic contribution to the energy density, as derived in Ref. \cite{krogel2013quantum} is
\begin{equation}
    \label{qmc}
    \mathcal{E}(\bm{r}) = \mathcal{T}(\bm{r}) + \mathcal{V}_{\text{ee}}(\bm{r}) + \mathcal{V}_{\text{ei}}(\bm{r}),
\end{equation}
where the many-body kinetic energy density is
\begin{equation}
\label{mbke}
    \mathcal{T}(\bm{r})=-\frac{1}{2}\sum_{i=1}^N\int d\bm{R}~\delta(\bm{r} - \bm{r}_i)\Psi^*(\bm{R})\nabla_i^2\Psi(\bm{R}),
\end{equation}
and the many-body potential energy densities are
\begin{equation}
    \label{mbpot}
    \mathcal{V}_\text{ab}(\bm{r})=\frac{1}{2}\sum_{i=1}^{N_a}\sum_{j=1}^{N_b}\int d\bm{R}~ \delta(\bm{r} - \bm{r}_i)|\Psi(\bm{R})|^2v_{\text{ab}}(\bm{r}_i, \bm{r}_j).
\end{equation}
In Equations \ref{mbke} and \ref{mbpot}, $\bm{R}\equiv\{\bm{r}_1, \bm{r}_2, \hdots, \bm{r}_N\}$ is the full many-body space, $v_{\text{ab}}$ is the interaction potential between two species, and $\Psi$ is the many-body wavefunction. For $\mathcal{V}_\text{ee}$, the interaction potential $v_\text{ee}$ is the Coulomb potential. For the potential energy density term $\mathcal{V}_\text{ei}$, the interaction potential is a norm-conserving pseudopotential of the Burkatzki-Filippi-Dolg form \cite{burkatzki2007energy}. The many-body wavefunction is constructed from a Slater determinant of Kohn-Sham orbitals and an optimized Jastrow factor. The QMC energy densities are generated in 3 stages. The first stage is a DFT calculation to obtain the Kohn-Sham orbitals needed for the Slater determinant. To perform the DFT calculation, we used Quantum Espresso \cite{giannozzi2009quantum} included with QMCPACK \cite{kim2018qmcpack} with an energy cut-off of 150 Ha. For the graphene lattices we used a $4\times4\times1$ k-point grid and for liquid water we used a $2\times2\times2$ k-point grid. The second stage is the Jastrow factor optimization. We used a one-body Jastrow factor for the nuclei and a two-body Jastrow factor for the electrons. The final stage is evaluating the energy which was done in 3 steps. This included a variational Monte Carlo (VMC) calculation followed by two DMC calculations, each with different time steps. We refer the reader to the SI for more information on the Jastrow optimization and other QMC parameters \cite{si}. All three stages were performed at each $\bm{k}$-point, and the final energy densities were obtained by averaging over all $\bm{k}$-points. All calculations were performed using QMCPACK \cite{kim2018qmcpack}, Nexus \cite{krogel2016nexus}, and a custom Python library \cite{github-qmc}. To compute the atomic contributions to the total energy, the closest energy density voxels to a particular atom (using the minimum image convention) were integrated. For graphene this yielded 50 atomic contributions per calculation and 96 atomic contributions per calculation for liquid water.

When generating a dataset for VDNNs, we used $\approx 5\times10^{5}$, $19\times19\times19$ slices of the DFT electron densities that were collected from the training configurations such that a uniform distribution of the function $\mathcal{F} = \sqrt{\mathcal{T}^2 + \mathcal{V}_{\text{ee}}^2 + \mathcal{V}_{\text{ei}}^2}$ was produced. These parameter choices were found to be optimal in a previous study on the prediction of DFT kinetic energy densities \cite{ryczko2021orbital}. Of these electron density slices, 99\% were used for training and 1\% were used for validation. When training the VDNNs, we map the electron density slices to values of $\mathcal{T},~\mathcal{V}_{\text{ee}}$, and $\mathcal{V}_{\text{ei}}$, as shown in \cref{fig:fig_1}a. We used the convolutional neural network (CNN) from Ref. \cite{ryczko2021orbital}, a batch size of 512, a learning rate of $10^{-5}$, the Adam optimizer \cite{kingma2014adam}, and layer-wise adaptive rate scaling with clipping \cite{you2017scaling}. We use ensemble models, where 5 independent networks were trained simultaneously and total energies were averaged over the models. We trained the networks for 250 epochs, and used the model that had the minimal mean squared error on the validation set. We found that the models began to overfit beyond $\approx250$ epochs. This overfitting can be attributed to the stochastic nature of the underlying training set. Beyond a certain point in training, the model begins to fit to statistical noise. We stopped training prior to this and used the model with the lowest mean squared error on the validation set. Training and inference was done across multiple nodes, each with 4 NVIDIA V100 GPUs. Our codes can be found here \cite{github-qmc}.

For learning the atomic energy contributions, we tried various atomic environment descriptors including atom centred symmetry functions (ACSF) \cite{behler2011atom}, atomic environment vectors (AEV) from the ANI models \cite{smith2018less, smith2019outsmarting, devereux2020extending},  and smooth overlap of atomic positions (SOAP) \cite{bartok2013representing}. Apart from AEV, we used the Python library dscribe \cite{dscribe} for generating the atomic environment descriptors with the default parameters. For AEV, we used the parameters which can be found in torchani \cite{gao2020torchani}. For the machine learning model, we employed Gradient Boosted Decsion Trees (GBDT), Gaussian Process Regression (GPR), artificial neural networks (ANNs), and KRR. We used the Python library scikit-learn \cite{scikit-learn} to implement all models and tried various sets of parameters for each model using grid searches. For both the graphene lattices and liquid water, we found that the KRR models had the lowest mean absolute errors (MAEs).

\begin{figure}[ht]
    \centering
    \includegraphics[width=\linewidth]{graphene_training_curve.pdf}
    \caption{Learning curves for various models along with their respective linear fits. The linear fits allows one to estimate how much data is necessary to achieve a particular accuracy. We estimate for VDNNs that we would need $\approx 400$ calculations to achieve chemical accuracy. However, for KRR with SOAP as the atomic environment vectors, we would only need $\approx 12$ calculations.}
    \label{results_1}
\end{figure}

\begin{figure}[ht]
    \centering
    \includegraphics[width=\linewidth]{sw_graphene_krr.pdf}
    \caption{The energy barrier curve associated with creating a Stone-Wales defect in a pristine graphene lattice. The error bars for QMC show the standard deviation of the calculation. The model KRR-16 has the lowest mean absolute error with 4.23 meV / atom.}
    \label{results_2}
\end{figure}

\section{Results}
\label{results}
The first task is the prediction of DMC total energy for pristine graphene lattices with random atomic perturbations. When building this dataset, we generated 8 atomic configurations where random perturbations of the atomic positions were applied to the equilibrium structure. To generate the random displacements, uniform random numbers were drawn with the range $[-0.1, 0.1]$ and were added to the pristine coordinates. 4 of these calculations were used to generate training/validation data and the remaining 4 were used for testing. For the VDNNs, we trained 4 independent ensemble models, where each model was allowed to see 1, 2, 3 or 4 training configurations. We then computed the total energies of the testing set by integrating the predicted energy densities from the ensemble models. Similar to Ref. \cite{ryczko2021orbital}, a rigid shift was applied to the predictions based on the residuals of the training set. This shift is needed to compensate for the errors from the model that accumulate when integrating on large grids. The shift can be calculated by computing the average error within an energy window. The number of energy windows depends on the structure of the residuals, and therefore differs for each model within the ensemble. We used a grid search to compute the optimal number of energy windows, and found the number of energy windows to be in the range $[8,32]$. In \cref{results_1}, we plot the learning curve for the VDNNs. The best MAE found for VDNNs was 120 meV / atom. For both KRR and VDNNs, the learning curve allows one to estimate the MAE of a model if additional calculations were available. A  linear fit was done in order to obtain these estimates. In \cref{results_1}, these estimates are labelled by the keyword ``fit.'' If one would like to have a VDNN that is capable of making predictions within chemical accuracy, one would need $\approx~400$ additional calculations. For the KRR models, we found that the SOAP representation had the lowest MAE which was 3.40 meV / atom. In comparison to the VDNNs, this is a 35x improvement. In addition, we estimate that we would need $\approx~8$ additional calculations to obtain a model capable of making DMC energy predictions within chemical accuracy. This accuracy, however, is limited to pristine graphene lattices with atomic perturbations.

\begin{figure*}[ht]
    \centering
    \includegraphics[width=\linewidth]{water_dmc.pdf}
    \caption{Top: QMC, DFT, and KRR model prediction energies for liquid water. Bottom: Residuals for DFT and the KRR model. The KRR model shows a 1.38x improvement over DFT with a MAE of 38.1 meV / molecule.}
    \label{results_3}
\end{figure*}

Given the performance of KRR in comparison to VDNNs, we concentrate on analyzing the capabilities of KRR in more complex scenarios. Namely, we are interested in studying the generalizability of the KRR model in regards to calculating the energy barrier curve associated with the introduction of a Stone-Wales defect within a graphene lattice, which was previously investigated with DFT in Ref. \cite{li2005defect}. As shown in \cref{results_2}, by rotating two atoms in the graphene lattice, one transitions from a pristine lattice ($\theta=0^{\circ}$) to a lattice with a Stone-Wales defect ($\theta=90^{\circ}$). For the training set, we started with the 3 configurations shown in \cref{results_2}. These configurations do not have any random atomic perturbations added. We then ran 16 additional DMC calculations of graphene lattice with a single Stone-Wales defect and atomic perturbations. The testing dataset consisted of 10 configurations, where two atoms were rotated such that the $nth$ configuration had an angle $\theta = n \pi/ 18,~n\in[1,9]$. For each configuration, a DFT structural relaxation was performed while holding the two rotated atoms fixed. As done in Ref. \cite{li2005defect}, the bond length between the two rotated atoms was also reduced. We first trained a KRR model ignoring the 16 additional DMC calculations and predicted the energy barrier curve. In \cref{results_2}, this is labelled by ``KRR-0.'' With only 3 total energy reference points, the KRR model does exceptionally well, with a MAE of 7.32 meV / atom. We then added in 2, 4, 8, or 16 additional DMC calculations (labelled ``KRR-2'', ``KRR-4'', ``KRR-8'', and ``KRR-16'' respectively in \cref{results_2}), trained a KRR model, and made predictions for the energy barrier. Interestingly, when comparing the models KRR-0 and KRR-[2,4,8], the MAEs increase (up to 1.8x worse with KRR-4) with the addition of data. We argue that with the addition of new data, there are two areas of chemical space that must be accurately modelled and connected in a physically meaningful way. The models KRR-[2,4,8] do not connect these two areas of chemical space well, leading to higher MAEs. However, with all 16 additional points, the MAE becomes 4.23 meV / atom (1.7x times improvement from KRR-0). The KRR-16 model learned to generalize its knowledge from the additional calculations to make accurate predictions of the energy barrier.

 We define the difference in energy between the highest peak seen in \cref{results_2} and the energy at $\theta = 0$ to be the energy barrier needed to be overcome to transition from a pristine lattice to one with a Stone-Wales defect. We find the energy barrier to be 224(46) meV / atom based on the DMC calculation. DFT underestimates the energy barrier by 17.6 meV / atom and KRR-16 underestimates the energy barrier by only 5 meV / atom.

 We now move from 2D to 3D materials and study the performance of KRR when applied to liquid water. Since there are two atom types, one KRR model was trained to predict atomic energy contributions for H, and another model was trained to predict atomic energy contributions for O. In both cases we found that the AEV representation had the lowest MAE in comparison to other representations. As described in \cref{methods}, there were a total of 32 DMC calculations that were done. Half of these were used in training, and the remaining half were used as a test set. In \cref{results_3}, these data points are highlighted under ``Training'' (in cyan) and ``Testing'' (in light red). For KRR, most of the predictions are within chemical accuracy. The MAE of the test set was 27.6 meV / molecule. In comparison to DFT, which had a MAE of 38.1 meV / molecule, the KRR model shows a 1.38x improvement.

Before concluding, we discuss the limitations and future work of both methods presented here. Firstly, both VDNNs and KRR models are limited to the data they are trained with (as is the case with all machine learning models). Secondly, for machine learning models that predict just scalar quantities (i.e. total energy), there is no error correction that must be done. However, as discussed previously, error from VDNNs accumulates when integrating, and must be accounted for. Thirdly, in VDNNs there is no dependence between neighbouring pixels that are output by the model. The dependence could be applied by outputting neighbouring pixels and averaging over the predictions or in some sort of post-processing routine. Fourthly, due to use of images, rotational invariance is not inherently built-in to the model, but must be learned during training. When constructing diversity in a dataset for VDNNs, one must consider the structure of the electron density, rather than a chemical environment. The two are inter-linked, but one must consider how the electron density images change within a dataset, and use this information when applying it to a new system. Lastly, since VDNNs must make millions of inference calls for a 3D grid, prediction time is relatively slow (on the order of minutes) compared to machine learning methods that output a scalar quantity. This is where predicting atomic contributions to the total energy excels. The KRR models take seconds to train and perform inference, and have rotational and translational symmetry built in. However, our current KRR implementation does not have automatic differentiation. Automatic differentiation with kernel methods has been done previously \cite{gardner2018gpytorch}, and future work of ours involves using these packages to obtain forces at the DMC level. In addition, this work shows that it is possible to perform supervised learning with DMC data. Given that our KRR method trains on atomic contributions to the total energy, we expect that size of datasets needed for training DNNs will be smaller by a factor of $\langle N \rangle$, which is the average number of atoms per calculation. Other future work of ours involves investigating $\Delta$ learning from DFT to DMC atomic contributions to the total energy.

\section{Conclusion}
We have used two machine learning methodologies to predict DMC total energies. The first methodology we used was VDNNs to map electron densities computed with DFT to energy densities calculated with DMC for pristine graphene lattices. After comparing these results with our second methodology, which used KRR and atomic environment vectors to predict atomic contributions to the total energy, we found that the errors from KRR were an order of magnitude lower. Given this result we then further studied the capabilities of KRR to predict atomic contributions to the energy for more complex tasks. This included predicted the energy barrier when introducing a Stone-Wales defect into a pristine graphene lattice and studying liquid water. In both cases, we found that KRR trained on DMC outperforms DFT. In addition, the KRR method requires no initial DFT calculation, but only the computation of atomic environment vectors. This allows one to rapidly infer energies at the DMC level. Since our method trains on atomic contributions to the total energy, we also conclude that dataset sizes in future studies can be reduced by a factor of $\langle N \rangle$, which is the average number of atoms per calculation. This would allow for DNNs to be used to predict DMC energies, furthering the performance showed in this report.

\section{Acknowledgements}
The Authors acknowledge Compute Canada, the Vector Institute for Artificial Intelligence, and the National Energy Research Scientific Computing Center for computational resources. KR and IT acknowledges the National Sciences and Engineering Council of Canada for funding. Work performed by JK (mentorship) was supported by the US Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.

\section{More Information on the QMC Calculations}

\subsection{Jastrow Factors}
The Jastrow factors were represented by B-splines and the coefficients were optimized using the quartic optimizer in QMCPACK \cite{kim2018qmcpack}. The meshes used in the QMC calculations were increased by 20\% and we found the average and corresponding standard deviation of the variance to energy ratio of $0.0205\pm0.0014$. The optimization was run for 100 blocks each consisting of 51200 Monte-Carlo samples.

\subsection{Calculation of the Energy Densities}
In this subsection, we provide additional details for the calculation of the energy densities. The first step is performing variational Monte Carlo (VMC) for 40 blocks with 30 warm-up steps and 10 steps used for averaging. Warm-up steps are steps taken where data is excluded from calculating quantities. The second step was DMC with a larger time step (0.02) for 20 blocks with 20 warm-up steps and 5 steps used for averaging. The final step was DMC with a smaller time step (0.01) for 200 blocks with 20 warm-up steps followed by 10 steps. The total number of Monte-Carlo samples completed for each of the final DMC calculations was 20 blocks  $\times$ 10 steps $\times$ 1024 walkers $\times$ 64 MPI processes $\approx13\times10^6$. The average and corresponding standard deviation of the total energy error across all calculations was $0.426\pm0.083$ meV.


\section{Introduction}
Quantum information science is an interdisciplinary field where computer science meets physics, and where long lasting principles of information transfer, cryptography and secret sharing are challenged \cite{NielsenBook}.  Quantum random number generators (QRNG) are mature devices that directly benefit from the properties of  quantum systems in order to produce truly random numbers \cite{Herrero-Collantes2017}. By harnessing the fundamental randomness from the outputs of quantum mechanical measurements it is possible to generate sequences of fundamentally unpredictable random numbers. A popular implementation for QRNGs is based on performing measurements on single-photons. Typically, in these photonic QRNGs, one degree of freedom of a photon is used to prepare a superposition state of the form $\ket{\psi} = \alpha\ket{0} + \beta\ket{1}$, where $\left|\alpha\right|^2 + \left|\beta\right|^2 = 1$, which is then projected onto the basis formed by $\ket{0}$ and $\ket{1}$. By letting one of the basis vectors represent a binary ``0'' and the other a binary ``1'', we can utilize the fact that measurements of an equally weighted superposition state yield uniformly distributed random measurement outcomes to generate a random sequence by performing repeated measurements.

One type of photonic QRNG is known as a branching-path generator \cite{Herrero-Collantes2017, Jennewein2020, Stefanov2000, Wang2006, Naruse2015}, and is typically implemented either using polarized light and polarizing beamsplitters, or with unpolarized light and conventional beamsplitters. Regardless of the scheme, it is of utmost importance that the beamsplitters' splitting ratio is $50:50$ for transmitted and reflected light respectively, in order to guarantee an unbiased sequence. In practice, due to manufacturing limitations, the splitting ratio will not be perfectly balanced, and this leads to non-uniformity of the bit distribution of the generated sequence. Furthermore, the splitting ratio may also change over time, due to temperature fluctuations and other environmental factors, as well as component degradation. Branching-path generators may also suffer from detector imbalances, where two identical detectors used at each of the beamsplitter's outputs have slightly different detection efficiencies that also lead to a bias in the sequence \cite{Herrero-Collantes2017}. Such a bias reduces the unpredictability of the sequence, and also so does its usability, not only for cryptography, but also for uses like Monte-Carlo simulations and gambling. For cryptography specifically, the security is directly proportional to the randomness of the secret key, and therefore a low-entropy sequence increases the probability that an attacker is able to successfully guess the key, which thus reduces the effective strength of the cipher \cite{Goldwasser1984}.

In this article, we demonstrate a tunable-ratio branching-path photonic QRNG by implementing a dynamically tunable beamsplitter using a fiber-optical Sagnac interferometer. The tunable beamsplitting function is implemented with an active phase modulator placed in a Sagnac loop working as a single-photon router \cite{Alarcon2020}. By modulating the phase of a light packet propagating in the clockwise direction vs the counter-clockwise one, we are able to tune the output probabilities from the Sagnac interferometer through a change in the amplitude of the driving signal fed to the phase modulator. This allows us to change the entropy of the generated sequence as we wish. Both outputs are then detected through a time-division multiplexing scheme to allow detection with only one single-photon detector. {\color{black}Interferometers in a Sagnac-type configuration have previously been used to modulate polarization \cite{Agnesi2019,Avesani2021} and intensity \cite{Roberts2018} in order to prepare states for quantum key distribution experiments. While building on the same principle, our system is the first demonstration of interferometric routing of single photons for quantum random number generation.} In this work we are able to generate a continuous stream of random numbers over 36 hours of continuous operation at an average rate of 131 kbit/s. The generation rate and entropy are highly stable over time due to the intrinsic phase stability of the Sagnac interferometer configuration. We finally demonstrate the randomness of the sequence by passing it successfully through the stringent statistical test suites NIST 800-22 \cite{Rukhin2010} and Dieharder \cite{R.Brown2011} aimed at assessing random number generators. A critical issue in QRNGs, especially the ones aimed at cryptographic applications, is whether the user can trust the device. A solution to this problem is performing independent real-time testing of the internal devices of the QRNG, usually through the user changing either the prepared quantum states, or a measurement parameter \cite{Ma2016}. {\color{black} In order to test the device and achieve high confidence that the generated sequence is private, a high interferometric visibility is necessary \cite{Carine2020}}. The dynamic tunability provided by our setup is therefore useful to tackle this issue, and is particularly appealing to recently demonstrated measurement-device-independent QRNGs \cite{Nie2016, Carine2020}. Our results show that Sagnac interferometers can be successfully employed as tunable beamsplitters for the construction of quantum random number generators, thus opening up for implementation of state-of-the-art measurement-device independent quantum random number generators.

\section{The experiment}

{\color{black}The principle for a QRNG scheme using a TBS is based on a source of strongly attenuated pulses, a beamsplitter with a variable splitting ratio and single-photon detection of the two outputs of the beamsplitter (\cref{fig:simpleexpsetup}).} A detection at one output is assigned the bit ``0'', while ``1'' is assigned to the other one. The entire system works synchronously, driven by Field Programmable Gate Array (FPGA) electronics, which also performs the data acquisition. {\color{black} If both outputs are simultaneously detected (either from two photons generated together, or from noise sources), the FPGA discards that result in order to remove the effects of after-pulsing and multiphoton events}. Similarly, if no photons are detected for a generated pulse (due to absorption) no result is recorded. The generated bits are then streamed from the FPGA to a personal computer for randomness extraction and storage. The bitstream then passes through statistical tests to certify its randomness.

\begin{figure}[t!]
    \centering
    \includegraphics[clip, trim=0cm 3.7cm 0cm 0cm, width=0.4\textwidth, angle=0]{figures/expsetupsimplified.pdf}
    \caption{\label{fig:simpleexpsetup} Principle of a tunable-ratio beamsplitter for quantum random number generation. Light from a photon source is split using a tunable beamsplitter (TBS). Detectors (D1, D2) placed at the outputs of the TBS register photon events and record a binary ``0'' or a binary ``1'' depending on whether the photon impinges on D1 or D2. The choice of assignment of ``0''/``1'' to D1/D2 is arbitrary.}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[clip, trim=0cm 3cm 0cm 0cm, width=\textwidth, angle=0]{figures/expsetup_labeled_sagnac.pdf}
    \caption{\label{fig:expsetup}Experimental setup. The CW-laser is first fed to an intensity modulator (IM) that chops the laser into 20 ns pulses, which are then attenuated using a variable attenuator (ATT) to weak coherent pulses (WCP) with an average photon probability of approximately $10$ photons per pulse. The WCPs are then fed to a Sagnac interferometer that acts as a tunable beamsplitter where the superposition $\ket{\psi}= 1/\sqrt{2}\left(\ket{\textrm{CW}} + ie^{i\phi}\ket{\textrm{CCW}} \right)$ is prepared, and CW and CCW denote the photon state of propagation inside the Sagnac loop. At the output of the Sagnac loop, we use a {\color{black} 150 meter} fiber-optic delay line ($L_{TB}$) to add {\color{black}$\approx 750$} ns of delay to one term, effectively time-bin multiplexing the two outputs from the Sagnac interferometer into an \emph{early} and a \emph{late} time-bin for detection using one single-photon detector.}
\end{figure}

The key aspect in our experiment is the use of a tunable beamsplitter (\cref{fig:simpleexpsetup}). Although in principle it is possible to create one by adding extra attenuators in each output of a conventional beamsplitter, or by using a half-wave plate in tandem with a polarizing beamsplitter, these options do not provide possibilities of fast tuning. Furthermore, in the case of the first option, extra losses are added to the system. Here we pursue an interferometric approach to implement tunability. For instance, a Mach-Zehnder interferometer can realize a tunable beamsplitter, with the tunability provided by setting the relative phase $\phi$ between the interferometer arms \cite{Ma2011}. The output probability of the single-photon is then given by a sinusoidal function of $\phi$. However, the drawback of using conventional Mach-Zehnder interferometers is that they suffer from intrinsic phase instability, and therefore require either mechanical isolation (i.e. dampened tables, environmental protection, etc) \cite{Micuda2014} or active stabilization \cite{Xavier2011,Reynaud1992}. In our work we follow the approach of \cite{Alarcon2020} where we employ a Sagnac interferometer to provide the same functionality, with the added benefit of intrinsic phase stability. In our Sagnac interferometer a high-speed telecom fiber-pigtailed phase modulator is placed to give a relative phase shift $\phi$ between the two propagation directions (\cref{fig:expsetup}). This is done by synchronizing the phase shift with the wave packet, such the phase shift is only applied in one direction by the modulator.

We employ a continuous wave (CW) distributed feedback telecom laser diode, centered at $\lambda$ = 1546.92 nm, connected to variable optical attenuators (ATT) to produce weak coherent states. {\color{black}The FPGA board generates electrical pulses at a repetition rate of $250$ kHz in order to drive an intensity modulator which chops the CW input light into optical pulses.} The electronic pulses are shaped to $20$ ns width, as limited by the employed electronics driver before being fed to the intensity modulator. The attenuated optical pulses pass through an optical circulator before entering one of the input ports of a bidirectional $50:50$ fiber coupler. The output ports form the Sagnac loop, together with a phase modulator and a 50 m long optical fiber spool ($L_s$) that works as a delay line (\cref{fig:expsetup}). The delay line is needed to ensure there is sufficient time separation for the phase modulation signal to act on the wave packet propagating in only in one direction and not in the opposing one, thus creating the necessary relative phase shift to generate the different output probabilities in the interferometer outputs \cite{Alarcon2020}. {\color{black}Not shown in \cref{fig:expsetup} are two manual polarization controllers added in the Sagnac loop, where one is used to align the pulses' polarization with the input polarization of the phase modulator, while the second one is used to ensure that the two counter-propagating terms have the same polarization at the recombining beamsplitter.}

Inside the Sagnac loop the attenuated pulses can be represented as a coherent superposition of the two counter-propagating paths as $\ket{\psi} = 1/\sqrt{2}(\ket{\textrm{CW}} + ie^{i\phi}\ket{\textrm{CCW}})$ where $\ket{\textrm{CW}}$ and $\ket{\textrm{CCW}}$ represent the clockwise and counter-clockwise propagation paths respectively. The detection probabilities at the outputs A and B of the Sagnac loop are proportional to $\cos^2(\frac{\phi}{2})$ and $\sin^2(\frac{\phi}{2})$. For voltages corresponding to $\phi =\pi/2, 3\pi/2$ the Sagnac works as a beamsplitter with a perfectly balanced splitting ratio. The phase shift $\phi$ is proportional to the modulator's driving voltage, thus allowing us to precisely control the splitting ratio between 0 and 1 in the ideal case. We fine-tune the desired ratio by applying a small offset $\delta V_\phi$ to the theoretical voltage value $V_\phi$. Therefore, the voltage we apply to achieve a balanced splitting ratio is $\overline{V}_{\pi/2} = V_{\pi/2} + \delta V_{\pi/2}$. {\color{black} Due to insertion losses and propagation losses in the time-bin multiplexer, a balanced splitting ratio between the early and late time bin does not imply a 50:50 probability for a photon to emerge at the outputs A:B of the Sagnac interferometer.} The two outputs of the Sagnac (denoted A and B in \cref{fig:expsetup}) are then routed to the detection system (with output A passing through the circulator to decouple it from the attenuated pulses coming from the source, propagating in to the Sagnac loop). We opted to use one SPD and mapped each output path to an \textit{early} or a \textit{late} time-bin, with both time slots detected sequentially by the SPD. The time-bin separation is created by a 750 ns time delay (a $\approx$150 m optical fiber spool) {\color{black}placed at the A output path} from the Sagnac loop. The two paths are then recombined with a 50:50 bidirectional fiber coupler, and one of its outputs is connected to the SPD (idQuantique id210), operating in Geiger mode with 10\% detection efficiency, 2.5 ns gate width and approximately $10^{-5}$
dark count probability per gate. The SPD is externally gated by the FPGA by sending pairs of pulses, each separated by the time-bin delay, that are synchronized with the signal sent to the intensity modulator that generates the photon pulses. {\color{black}Each of these pulses correspond to a gate in the \emph{early} and \emph{late} bin respectively.} The output signals from the SPD are also read by the FPGA for data acquisition. We set the deadtime of the SPD to 500 ns, so that it is shorter than the time separation between the \emph{early} and \emph{late} time-bins. We then let a detection of a photon in the \emph{early} time-bin represent a ``0'' and a detection in the \emph{late} time-bin represent a ``1''.

\section{Results}

First we demonstrate the variable splitting ratio of the Sagnac interferometer acting as a tunable beamsplitter. We set the attenuators such that we maximize the detection probability at the detector for both time-bins, while still operating in a linear regime (i.e. the detector is not saturating).
This is observed from \cref{fig:countsphasevoltage}, where the number of detection events per integration time fits well with the theoretical prediction that output A [B] from the Sagnac loop should follow a $\cos^2(\frac{\phi}{2})$ [$\sin^2(\frac{\phi}{2})$] distribution. We also observe from these results that the Sagnac interferometer is operating correctly, and complementary output behaviour of the two outputs is observed, as expected.

\begin{figure}[t!]
    \centering

    \subfloat[]{\label{fig:countsphasevoltage}\includegraphics[width=.49\textwidth]{figures/voltage_vs_counts-eps-converted-to.pdf}}
    \hspace{1em}
    \subfloat[]{\label{fig:entropysweep}\includegraphics[width=.47\textwidth]{figures/entropy_vs_voltage_with_inset-eps-converted-to.pdf}}

    \caption{\label{fig:extinctionentropy} a) Characterization of detected photon events measured in the \emph{early} and \emph{late} time-bin, as a function of applied voltage to the phase modulator. The error bars are calculated assuming Poissonian photon statistics ($\sigma = \sqrt{N}$ for a photon count of $N$). Also plotted is a fit of the theoretical detection probabilities that are proportional to $\cos^2$ [$\sin^2$] for the \emph{early} [\emph{late}] time-bin, respectively. b) Raw entropy (before extraction) as a function of applied voltage to the phase modulator Also shown in the inset is a plot of the entropy versus applied voltage to the phase modulator with a smaller step size. }

\end{figure}

We then wish to optimize the operating point of the Sagnac loop for optimal random number generation. We scan the driving voltage of the phase modulator and measure the Shannon entropy of the generated bits when considered as $8$-bit symbols. The Shannon entropy $\mathcal{H}$ is defined as $\mathcal{H} = -\sum_{i=0}^{255}p_i\log_2{(p_i)}$, where $p_i$ is the probability of finding the $i_{th}$ 8-bit sequence. We increase the voltage $V_\phi$ applied to the phase modulator in the range $0\leq V_\phi \leq 4.2$ in steps of $200$ mV and measure the entropy, as can be seen in \cref{fig:entropysweep}. An inset is shown in \cref{fig:entropysweep} where we also sweep the voltage between $2-2.3$ volts in steps of $20$ mV, and likewise measure $\mathcal{H}$ to identify the optimal operating point where the entropy is maximized. It can be clearly seen that the maximum obtained entropy (7.98 bits/byte) occurs for an applied driving voltage of around $2.15$ V.

\begin{figure}[t!]
    \centering
    \subfloat[]{\label{fig:entropytime}\includegraphics[clip, trim=0cm 5.1cm 0cm 0cm, width=.48\textwidth]{figures/bitrate_entropy_over_time-eps-converted-to.pdf}}
    \hfill
    \subfloat[]{\label{fig:bitratetime}\includegraphics[clip, trim=0cm 0cm 0cm 5.28cm, width=.48\textwidth]{figures/bitrate_entropy_over_time-eps-converted-to.pdf}}

    \caption{\label{fig:entropybitrate} Characterization of the stablity of the raw bitrate  a) and entropy stability b) over a period of 36 hours. Each of the data points represents the average over 30 minutes.}
\end{figure}

\begin{figure}[t!]
    \centering
    \subfloat[]{\label{fig:nisttable}\includegraphics[clip, trim=0.4cm 0cm 1.2cm 0.5cm, width=.5\textwidth]{figures/nist_with_double_siglev-eps-converted-to.pdf}}
    \subfloat[]{\label{fig:diehardertable}\includegraphics[clip, trim=0cm 0cm 1.3cm 0.5cm, width=.5\textwidth]{figures/dieharder-eps-converted-to.pdf}}

    \caption{\label{fig:statisticalresults} Results from statistical testing of randomness from NIST STS 800-22\cite{Rukhin2010}, a) and Dieharder\cite{R.Brown2011} b). In \cref{fig:nisttable} we plot the proportion of tests that have passed, along with the confidence interval that the proportion must fall within. For \cref{fig:diehardertable}, each colored vertical line corresponds to a p-value from the tests. Lines in green are p-values that fall outside of the confidence bound $10^{-2} < p < 10^0$. The red lines are the resulting p-value of a Kolmogorov-Smirnov test done on test results that yield multiple p-values. }
\end{figure}

With the entropy maximized, we then proceed to acquire a long continuous sequence of bits to feed to the statistical randomness tests, and also to demonstrate the robustness of our implementation. The raw bit sequence is buffered in the FPGA in blocks of 32 kilobytes and then transmitted over an Ethernet connection to a personal computer for storage. Before assessing the sequence, it is passed through a Toeplitz randomness extraction procedure in order to remove any remaining non-uniformities from the raw data. To extract randomness, we split the generated sequence into $N$ sequences of length $n$, and multiply each subsequence with an $n\times m $ binary Toeplitz matrix which yields an output of $m$ hashed bits \cite{Qi2017}. We choose to use $n=400$, and as defined in \cite{Qi2017}, we obtain $m = \mathcal{H}_\mathrm{min} -2\log_2\varepsilon$ where $\mathcal{H}_\mathrm{min}$ is the min-entropy of 8-bit strings which is defined as {\color{black}$\mathcal{H}_\mathrm{min} = -\log_2 \left[ \max_{x\in\left\{0, 1\right\}^8} \mathrm{Pr}\left\{X=x\right\}  \right] $. We measure a min-entropy of $7.7451$ from the acquired sequence before extraction.} We choose a security parameter $\varepsilon=2^{-100}$ which is derived from the leftover hash lemma \cite{Ma2013}. We seed the Toeplitz matrix $T$ at the beginning of the extraction procedure with $n+m-1$ bits taken from the raw sequence, as it is sufficient to only specify the first row and first column. The matrix $T$ can be reused for each of the subsequences, and we obtain the extracted sequence of bits by concatenating the results of the multiplications $Tn_1, Tn_2, Tn_3,...,T_N$, where the subscript represents the subsequence index \cite{Qi2017, Ma2013}. The extraction procedure is lossy, and for our choice of parameters, we are able to achieve an efficiency of $\approx 50\%$. {\color{black} After extraction we calculate the min-entropy of the sequence to be $7.9982$ bits.}

While recording the raw data, we also measure the entropy of each streamed 32 kilobyte block, and as can be seen in \cref{fig:entropytime}, we achieve an average entropy of $7.9925 \pm 0.0014$ bits/byte over $36$ hours, which demonstrates exceptional phase stability of the QRNG. Likewise, we also measure the raw bitrate when streaming the data, and we achieve an average bitrate of $131.8 \pm 2.3$ kbps (\cref{fig:bitratetime}). The bitrate is mainly limited by the delays in the electronics that amplify the signals from the FPGA in order to drive the phase modulator, and in the triggering electronics which is necessary to synchronize the generation of pulses with the detection. By improving the driver electronics, the bitrate is only limited by the separation of the time bins.

After generating $2.2$ gigabytes of data, and extracting it,
we pass the sequence through the NIST 800-22 Statistical Test Suite\cite{Rukhin2010} (NIST test) and the Dieharder\cite{R.Brown2011} test suite. Both test suites subject the random sequence to a number of statistical tests aimed at assessing the randomness of any given sequence. We run the NIST tests with $1000$ bitstreams of length $1$ Mb. To deem any given test as passed, we form a confidence bound as defined in \cite{Rukhin2010} for the proportion of bitstreams that must yield a p-value  $p \geq 10^{-2}$ using the formula $\hat{p} \pm 3\sqrt{\frac{\alpha(1-\alpha)}{n}}$ where $\alpha=10^{-2}$, $n=1000$ is the number of bitstreams and $\hat{p}=1-\alpha$. This gives us a confidence interval of $0.980 \leq r < 0.999$, and if the ratio of tests that have passed falls outside this range, there is evidence that the sequence is not random \cite{Rukhin2010}. \Cref{fig:nisttable} shows the proportions of tests that pass, for any given test category, along with the significance level.

The Dieharder test suite yields a verdict ("Fail", "Weak", "Pass") for each test, and a test is deemed to pass if the p-value associated with it is in the interval $10^{-2} < p < 1-10^{-2}$. As can be seen in \cref{fig:diehardertable}, there are parameterized tests where each test yields multiple p-values (one for each value of the parameter). We plot all outputted p-values, and compute a resulting p-value from a Kolmogorov-Smirnov test of all the p-values. As is well-known, the p-values should be uniformly distributed under the null hypothesis that the sequence is indeed random, and the Kolmogorov-Smirnov test is used to eliminate false rejection of the null hypothesis (Type-I error) \cite{Hung1997}.
For three of the tests, we observe that the resulting p-values fall outside the confidence bound, and for those values of the parameter, the test verdict is "weak". Since the Dieharder suite is a statistical test, it is probable that there will be occasional "weak" results. If the generator truly is weak, then it will certainly fail other tests \cite{R.Brown2011}.

\begin{figure}[t!]
    \centering

    \includegraphics[clip, trim=0cm 0cm 0cm 0cm, width=.6\textwidth]{figures/partial_mdi-eps-converted-to.pdf}

    \caption{\label{fig:partial_mdi} { \color{black} Achieved visibilities for the states $\ket{\Psi}$ and $\ket{\Phi}$ that are expected to yield deterministic measurement results, together with the visibilities for the balanced state $\ket{\Omega}$ used for randomness generation. } }
\end{figure}

{\color{black}To further demonstrate that the tunability can be used as a building block in measurement-device independent protocols we run a prepare-and-measure type protocol where we use the TBS to prepare either of the orthogonal states $\ket{\Psi}$ and $\ket{\Phi}$ which are expected to yield deterministic measurement results, and a state $\ket{\Omega}$ which yields random (but balanced) measurement outcomes. We prepare the states by applying different voltages to the PM inside the Sagnac for each state. From \cref{fig:countsphasevoltage} we determine the voltages to be  $V_\Psi = 0 \textrm{ V}$,  $V_\Phi = 4.2 \textrm{ V}$ and  $V_\Omega = 2.15 \textrm{ V}$.  For each acquired 32 kilobyte block of data we then randomly choose between preparing either of the states with a probability $\mathrm{Pr}\left\{ \ket{\Psi}\right\}  = \mathrm{Pr}\left\{\ket{\Phi}\right\} = 0.005 $ and $\mathrm{Pr}\left\{\ket{\Omega}\right\} = 0.99$, where $\mathrm{Pr}\left\{ \ket{\bullet}\right\} $ is the probability to prepare a certain quantum state, in order to demonstrate the switching of states necessary to implement measurement-device independent QRNG protocols. In \cref{fig:partial_mdi} we show the recorded visibilities for each of the states as a function of time, where the visibility $\nu_x \in \left[0,1\right]$ for the state $x\in\left\{\ket{\Psi},\ket{\Phi}, \ket{\Omega}\right\}$  is defined $\nu_{x} = \left| \frac{n_\mathrm{early} - n_\mathrm{late}}{n_\mathrm{early} + n_\mathrm{late}} \right| $ where $n_\mathrm{early}$ and $n_\mathrm{late}$ are the detector counts in the early and late bin respectively. We are able to achieve a mean visibility of $\nu_\Psi = 0.997 \pm 3.3\cdot 10^{-4}$,  $\nu_\Phi = 0.990 \pm 3.1\cdot 10^{-3}$ and $\nu_\Omega = 0.004 \cdot 10^{-2} \pm 7.7 \cdot 10^{-3}$ for each of the three states without subtracting dark counts. }

\section{Conclusions}
We have proposed and experimentally demonstrated a photonic quantum random number generator based on a dynamically tunable beamsplitter. The random numbers are generated depending on which output from the beamsplitter the single-photons are detected. By tuning the beamsplitting ratio we are able to maximize the generated entropy of the random sequence. This also helps maximize the entropy in spite of differences in losses following the two outputs. We implemented the tunable beamsplitter through a Sagnac interferometer, in which the tunability is given by a relative phase between the two counter-propagating paths, applied with an electro-optical telecom phase modulator. This technique allows ultra-fast tunability, due to the short response time of the phase modulator.

We are able to generate a continuous stream of random numbers over 36 hours of continuous operation, showing an average raw generation rate of 131 kbit/s and an average generated entropy of 7.9925 bits per byte in the raw sequence, thus demonstrating a stable source of entropy. {\color{black} As the bitrate is mainly limited by internal delays in the employed driving electronics, by improving the driver circuitry, a considerable bitrate gain can be achieved. By reducing these delays, we would also reduce the lower limits on the time-bin spacing and thus the repetition rate can be increased.} Our setup is robust and built using only off-the-shelf telecom fiber components, thus also showing the practicability of our setup.  Our results can aid in future designs of dynamically tunable quantum random number generators, and serve as a foundation for measurement-device quantum random number generation due to ease of changing the beamsplitting ratio, through the appropriate selection of the applied phase in the Sagnac loop.

{\color{black}Our system also demonstrates that the tunable beamsplitter enables realization of measurement-device independent quantum random number generators thanks to the high visibility when measuring orthogonal states, as well the fact that the visibilities show minimal degradation over time.}

\section*{Acknowledgments}
We acknowledge financial support from CENIIT Link\"{o}ping University, the Swedish Research Council (VR 2017-04470), QuantERA grant SECRET (VR 2019-00392), the Knut and Alice Wallenberg Foundation through the Wallenberg Center for Quantum Technology (WACQT).

\section*{Author Contributions}
J.A. and A.A. conceived the work together with G.B.X. J.A. developed the electronic control system, and A.A. developed the optical system. Data analysis was done by J.A. with assistance from G.B.X. The work was coordinated by G.B.X. All authors contributed to the writing of the manuscript.

\section*{Competing Interests}
The authors declare no competing interests.

\section*{References}
\printbibliography[heading=none]

\end{document}