We have done a theoretical study of the γd → π 0 ηd reaction starting with a realistic model for the γN → π 0 ηN reaction that reproduces cross sections and polarization observables at low energies and involves the γN → ∆(1700) → η∆(1232) → ηπ 0 N process. For the coherent reaction in the deuteron we considered the impulse approximation together with the rescattering of the π 0 and the η on a different nucleon than the one where they are produced. We found this second mechanism very important since it helps share between two nucleons the otherwise large momentum transfer of the reaction. We find a natural explanation, tied to the dynamics of our model, for the shift of the η − d mass distribution to lower invariant masses, and of the π 0 − d mass distribution to larger invariant masses, compared to a phase space calculation. We also study theoretical uncertainties related to the large momenta of the deuteron wave function involved in the process. Striking differences are found with the experimental angular distribution and a discussion is done around this issue.