Appendix \[sec:alg\] we design an integration algorithm, which is consistent with the time reversibility of the Langevin Dynamics (LD) in the macroscopic limit, [Eq. (\[eq:mac1\])]{}. Below we report the results of simulations, obtained with this numerical scheme.

As described in [Sec. \[sec:mic\]]{}, we determined the parameters of [Eq. (\[eq:lgv2\])]{} for the shear current $P_{yx}$ and the two models of athermal noise, $\epsilon_\mathrm{ex}(t)$ and $\epsilon_\mathrm{sh}(t)$, from our MD calculations, obtained at $\gamma = 1$ r.u. [Table \[tab:stat\]]{} summarizes results of our simulations. Statistics of the original MD data are reproduced very accurately and precisely by our numerical LD model up to the third-order statistics.

A notable discrepancy between the simulations is observed in [Table \[tab:stat\]]{} only for the excess kurtosis. This should be expected, because the parameters of [Eq. (\[eq:lgv2\])]{} were obtained by fitting solely the first three cumulants of the MD data, and the two LD models, considered here, differ by their fourth cumulants in principle. Such level of precision, however, is beyond the statistical and numerical resolution of our simulations. Indeed, the magnitude of the excess kurtosis for the MD data is comparable to its standard deviation. Its theoretical value in our LD models is $17\cdot10^{-5}$ for the case 