Let's assume the following symbolic equations
$$E=\sqrt{p^2+m^2}+\frac{e d_1(p)}{r}+\frac{e^2d_2(p)}{r^2}$$
being $d_{i}$ some generic function of $p$. I would like to find solutions to this equation in terms of $p$ and for small values of the constant $e$. How can I do this using Mathematica?
I'll make an example of the functions $d_i(p)$:
d1[p_] = Sin[p]; d2[p_] = p^2 - 7; Series-expansion of the solution for small $e$: I'll use $E=$F because the symbol E is already in use. The order of the series-expansion is set to 2 here but can be increased at will,
Assuming[F > 0, AsymptoticSolve[F == Sqrt[p^2 + m^2] + (e d1[p])/r + (e^2 d2[p])/r^2, {p, Sqrt[F^2 - m^2]}, {e, 0, 2}] // FullSimplify] {{p -> (1/( 4 ((F - m) (F + m))^(3/2) r^2))(-e^2 (m^2 + 4 F (F - m) (F + m) (-7 + F^2 - m^2)) + 4 (F^2 - m^2)^2 r^2 + e^2 m^2 Cos[2 Sqrt[F^2 - m^2]] + 4 e F (-F^2 r + m^2 r + e F Sqrt[F^2 - m^2] Cos[Sqrt[F^2 - m^2]]) Sin[Sqrt[ F^2 - m^2]])}}
This is the solution that, for $e=0$, passes through $p_0=\sqrt{F^2-m^2}$. There is also the other branch that passes through $p_0=-\sqrt{F^2-m^2}$; choose the branch in the AsymptoticSolve input.
For the more complex problem in the comments, it's sufficient to extend the assumptions a bit:
Assuming[ma > 0 && mb > 0 && F > ma + mb, AsymptoticSolve[ F == Sqrt[p^2 + ma^2] + Sqrt[p^2 + mb^2] + (e d1[p])/r + (e^2 d2[p])/r^2, {p, Sqrt[(F - ma - mb) (F + ma - mb) (F - ma + mb) (F + ma + mb)]/(2 F)}, {e, 0, 2}] // Refine] (lengthy output)
AsymptoticSolve is an experimental function and may have unexpected trouble. I'm away from the computer and cannot try, the weather is splendid in Europe! $\endgroup$ Just an idea how to solve the problem without specifying d1[p],d2[p] (Perturbation theory, MMA version <12 )
First assume the solution p can be approximated by p -> p0 + p1 edr + p2 edr^2 + p3 edr^3 (* edr=e/r *)
Substituting into the equation and series expansion
Clear[d1,d2] ser = Series[-F + Sqrt[p^2 + m^2] + (edr d1[p]) + (edr^2 d2[p]) /. p -> p0 + p1 edr + p2 edr^2 + p3 edr^3, {edr, 0, 2}] // Normal (*-F + Sqrt[m^2 + p0^2] + edr ((p0 p1)/Sqrt[m^2 + p0^2] + d1[p0]) +edr^2 ((m^2 p1^2)/(2 (m^2 + p0^2)^(3/2)) + (m^2 p0 p2)/(m^2 + p0^2)^(3/2) + (p0^3 p2)/(m^2 + p0^2)^(3/2) + d1[p0] + p1 Derivative[1][d1][p0])*) The series coefficients can be solved for p0,p1,p2
CoefficientList[ser, edr] // FullSimplify[#, F > m > 0] & Simplify[Solve[% == 0, {p0, p1, p2}], F >m> 0] (*{{ p0 -> -Sqrt[F^2 - m^2], p1 -> (F d1[-Sqrt[F^2 - m^2]])/Sqrt[F^2 - m^2], p2 -> (1/(2 (F^2 - m^2)^2))(m^2 Sqrt[F^2 - m^2] d1[-Sqrt[F^2 - m^2]]^2 +2 F (F^2 - m^2)^(3/2) d2[-Sqrt[F^2 - m^2]] +2 F^2 (F^2 - m^2) d1[-Sqrt[F^2 - m^2]]Derivative[1][d1][-Sqrt[F^2 - m^2]]) }, { p0 -> Sqrt[F^2 - m^2], p1 -> -((F d1[Sqrt[F^2 - m^2]])/Sqrt[F^2 - m^2]), p2 -> -(1/(2 (F^2 - m^2)^2))(m^2 Sqrt[F^2 - m^2] d1[Sqrt[F^2 - m^2]]^2 +2 F (F^2 - m^2)^(3/2) d2[Sqrt[F^2 - m^2]] +2 F^2 (-F^2 + m^2) d1[Sqrt[F^2 - m^2]] Derivative[1][d1][Sqrt[F^2 - m^2]]) }}*) 
AsymptoticSolve does in a more automated way. $\endgroup$ d1[p],d2[p] and for users with MMA versions<12 $\endgroup$