  In this paper, we study the error in first order Sobolev norm in the
approximation of solutions to linear parabolic PDEs. We use a Monte Carlo Euler
scheme obtained from combining the Feynman--Kac representation with a Euler
discretization of the underlying stochastic process. We derive approximation
rates depending on the time-discretization, the number of Monte Carlo
simulations, and the dimension. In particular, we show that the Monte Carlo
Euler scheme breaks the curse of dimensionality with respect to the first order
Sobolev norm. Our argument is based on new estimates on the weak error of the
Euler approximation of a diffusion process together with its derivative with
respect to the initial condition. As a consequence, we obtain that neural
networks are able to approximate solutions of linear parabolic PDEs in first
order Sobolev norm without the curse of dimensionality if the coefficients of
the PDEs admit an efficient approximation with neural networks.
