One of the hallmarks of short and intense laser-pulse-induced reactions at surfaces is that their yield depends superlinearly on fluence. In Figure we show the desorption yield of H 2 from a fully covered (4×4) cell, as a function of fluence for fluences in the range of 10-60 J/m 2 . Note that the training data is for fluence of 60 J/m 2 , so using larger fluences would result in an extrapolation regime for our neural network potential. We observe a superlinear behavior, with an exponent of about 1.8. One can also observe the onset of a saturation regime for fluences larger than 40 J/m 2 . Before saturation, the exponent is around 1.9. Fitting experimental yields before saturation results in a slightly larger exponent of 2.8, and saturation is reached only at around 100 J/m 2 . Further, in ref 35, no desorption for coverages lower than 0.5 ML was found. In Figure we show, for a (4×4) cell, desorption yields for H 2 and D 2 as a function of coverage with either H or D atoms. Coverages between 1/8 (2 adsorbate atoms per supercell) and 1 (16 atoms per supercell) are considered. We again obtain qualitative agreement with experiments, but desorption is possible even for coverages of about 0.375 ML (or 6 adsorbate atoms per (4×4) cell). Note that previous theoretical studies did not simulate coverage dependence as previous analytical potentials by construction could not simulate it, and in DFT-based dynamics, these simulations are expensive, especially when low yields are expected. In ref 35 the authors proposed that the results can be explained by the coverage dependence of the binding energy and/or the electronic friction coefficient. In our methodology, the electronic friction coefficient does not depend on the coverage, but as can be seen in We have also performed molecular dynamics simulations for different cell sizes, from (2×1) up to (8×8), always at full coverage. Results are shown in Figure . Clearly, simulations using smaller unit cells severely underestimate desorption yields. One can rationalize this result with a reduced configurational space, leading to desorption, in smaller unit cells as shown in Figure . This shows that better agreement of theoretical studies with experimental isotope ratios and desorption yields prior to ref 17, e.g. in ref 7, can be due to the fact that in ref 7, a single moving H 2 /D 2 in a (2×2) cell was considered and the somewhat too low desorption energy (1.042 eV compared to the effective barrier of 1.35 eV in ref 35) was compensated with a reduced configurational space. Note that even the AIMDEF methodology using already, for DFT, very large (4×4) unit cells underestimates yields as much as 50%. Only with very large unit cells such as (6×6) and (8×8) we obtain convergence in yields. Therefore, true first principles agreement with experiment should be achievable only with the methodology presented in this manuscript, but trained on DFT data calculated with more accurate exchange-correlation functional. Note also that in this particular case DFT data corresponds to only 3-layer slab that does not provide fully converged results, and all Ru atoms were fixed neglecting energy exchange with the surface.