It follows that $$r_2=0.$$ The general ansatz linear in $\hat{M}$, of order $a_0^2$, has 10 coefficients. Equating it with 0 gives a system of 16 homogeneous equations. The equations have a unique solution given by all the constants equal to 0.

In the third order, $r_3=\frac{1}{6}\rho^{3} -\frac{1}{6}r_1^{3}$ + terms of order $a_0^{3}$ coming from the use of relations $\widetilde{\mathcal{R}}$ ($\lambda=1/2$) in $(\rho+ \frac{1}{2}\rho^{2} -r_1 -\frac{1}{2}r_1^{2})$ + terms of order $a_0^{3}$ coming from the expansion of $\hat{M}_{k0}$ in $r_1$. Summing up all the contributions, we find $$\begin{aligned}
r_3 &=& -\frac{ia_0^3}{24}\cdot 1\otimes x_i \Bigl( 3p_0^3\otimes p_i -3p_i\otimes p_0^3 +3p_0\otimes p_ip_0^2 -3p_ip_0^2\otimes p_0 +2p_0^2\otimes p_ip_0 -2p_ip_0\otimes p_0^2 \nonumber \\
&& -2p_ip_jp_0\otimes p_j +2p_j\otimes p_ip_jp_0 -2p_ip_j\otimes p_jp_0 +2p_jp_0\otimes p_ip_j -2p_ip_0\otimes p_jp_j +2p_jp_j\otimes p_ip_0 \nonumber \\
&& -2p_0\otimes p_ip_jp_j +2p_ip_jp_j\otimes p_0 \Bigr).\end{aligned}$$ As for the first and second order, the obtained expression is compared with the most general ansatz linear in $\hat{M}$ and of order $a_0^{3}$, which in this case contains 28 constants. We get a system of 35 equations for the 28 coefficients. Solving the equations leaves 3 coefficients undetermined, giving us a three-parameter family of solutions for $r_3$ $$\begin{aligned}
r_3 &=& -\frac{ia_0^3}{24} \Bigl( -\alpha_1 M_{ij}p_jp_0\otimes p_i +\beta_1p_i\otimes M_{ij}p_jp_0 -\alpha_2M_{ij}p_j\otimes p_ip_0 +\alpha_2 p_ip_0\otimes M_{ij}p_j 