In this article, following the paradigm of bias–variance trade-off philosophy, we derive parametrix expansions of order two, based on the Euler–Maruyama scheme with random partitions, for the purpose of constructing an unbiased simulation method for multidimensional stochastic differential equations. These formulas lead to Monte Carlo simulation methods which can be easily parallelized. The second order method proposed here requires further regularity of coefficients in comparison with the first order method but achieves finite moments even when Poisson sampling is used for the partitions, in contrast to Andersson and Kohatsu-Higa (2017). Moreover, using an exponential scaling technique one achieves an unbiased simulation method which resembles a space importance sampling technique which significantly improves the efficiency of the proposed method. A hint of how to derive higher order expansions is also presented.