We also point out that the double-peak structure, which is not visible in the experiments, persists and remains visible even for the simulation with 3,000 snapshots. It becomes more pronounced when thinner gaussians are used. This feature is reminiscent of the double-peak structure observed on the x-ray emission spectra in liquid water 56 . Fig. shows the spectra constructed from either one IE (of the central water) or 5 IEs (of the cluster of 5 waters). As one can see, for the small QM (5 waters), the spectrum changes qualitatively when the IEs of all 5 waters are used, leading to a much smaller red shift relative to the gas-phase peak than the calculation using the IEs of 1 central water (the difference between the two calculations is 0.25 eV). This illustrates the limitations of the electrostatic embedding, which apparently is not able to correctly describe the waters on the QM/MM boundary. Using the QM of 20 waters, results in a larger red shift (by about 0.5 eV), consistent with the benchmark calculations for a single snapshot (Fig. ). Moreover, the calculations with the large QM yield the same spectra (in terms of the peak maximum and its width) whether using the IEs of one or five waters, suggesting that the boundary is sufficiently far in this case. The spectrum computed with the IEs of 5 waters is more smooth, as in such calculation the sampling is effectively improved by a factor of 5. Hence, we construct our best spectrum by using the IEs of the five water molecules embedded in the QM cluster of 20 waters (with the rest of waters described by point charges). The effect of triples is evaluated using the MLCC3 method with 20 CCSD waters and the triples operator restricted to orbitals within the first solvation shell of the central water molecule. Fig. shows our best spectrum and compares it with the best available experimental spectra . We note that the reported spectrum of Olivieri et al. 12 is for zero bias and, therefore, does not include the correction due to the interface. The corrected values from this experiment are given in Table , which also lists the computed peak positions and widths for different protocols and compares them with the experimental values. To extract the value of the peak maxima and the width from the computed spectra, we convoluted the computed spectra with broad gaussians (FWHM=0.5 eV) to remove the noise due to finite sampling. The effect of this broadening is analyzed in the SI: Figure in the SI shows the sample spectra obtained by convoluting the raw data with gaussians of different widths. As one can see, using broad gaussians does not introduce noticeable change in the band width.