We can compute a, e first to show, that there are no such events we are looking for. This code is only working minimal example I been able to derive from the code above: (Constants)

The last picture is the expression Papsis[t] - Roche[2000]] used for event detection. There is no time where Papsis[t] - Roche[2000]]<=0, and therefore there are no such events.

We can compute a, e first to show, that there are no such events we are looking for. This code is only working minimal example I been able to derive from the code above:

The last picture is the expression Papsis[t] - Roche[2000]] used for event detection. There is no time where Papsis[t] - Roche[2000]]<=0, and therefore there are no such events.

Update 1. In a case of initial data a[0] == 1*au, e[0] == 0.3 there are several events and we can detect and plot solution with events as follows (here we use a[t] in astronomical unit, and t in Myr)

au = UnitConvert[Quantity[1, "AstronomicalUnit"], "Meters"] // 
    QuantityMagnitude // Rationalize;
c = UnitConvert[Quantity[1, "SpeedOfLight"], "MetersPerSecond"] // 
    QuantityMagnitude // Rationalize;
Qpr = 1;
Lsun = UnitConvert[Quantity[1, "SolarLuminosity"], "Watts"] // 
    QuantityMagnitude // Rationalize;
Rsun = UnitConvert[Quantity[1, "SolarRadius"], "Meters"] // 
    QuantityMagnitude // Rationalize;
Msun = UnitConvert[Quantity[1, "SolarMass"], "Kilograms"] // 
    QuantityMagnitude // Rationalize;
G = UnitConvert[
     Quantity[1, "GravitationalConstant"], ("Meters"^2*"Newtons")/
      "Kilograms"^2] // QuantityMagnitude // Rationalize;
year = UnitConvert[Quantity[1, "Years"], "Seconds"] // 
    QuantityMagnitude // Rationalize;
Myr = year*10^6;
Gyr = year*10^9;
Mwd = 6/10*Msun;
Cst = 127/100;
U = 10^17;


(*Functions*)

L[t_] := (326/100*Lsun*(Mwd/(6/10*Msun)))/(1/10 + t)^(118/100);
Roche[dens_] := (65/100*Cst*
     Rsun*(Mwd/(6/10*Msun))^(1/3))/(dens/3000)^(-1/3);
Papsis[t_] := au a[t]*(1 - e[t]);

(*Radiative Drag*)

RDdadtR = -((1/
       c^2)*((3*L[t]*Qpr*(2 + 3*e[t]^2))/(16*Pi*2000*Rast*
         au a[t]*(1. - e[t]^2)^(3/2))));

RDdedtR = -((1/
       c^2)*((15*L[t]*e[t])/(32*Pi*Rast*2000*au^2 a[t]^2*
         Sqrt[1. - e[t]^2])));


Do[tmax[r] = 9 Gyr/Myr; 
 sol[r] = NDSolve[{Derivative[1][a][t] == (Myr/au) RDdadtR /. {Rast ->
          r}, Derivative[1][e][t] == Myr RDdedtR /. {Rast -> r}, 
      a[0] == 1, e[0] == 3/10, 
      WhenEvent[{Evaluate[Papsis[t] <= Roche[2000]]}, {tmax[r] = t, 
        "StopIntegration"}]}, {a[t], e[t]}, {t, 0, 9 Gyr/Myr}, 
     AccuracyGoal -> 20, PrecisionGoal -> 20, 
     WorkingPrecision -> 35][[1]] // Quiet; 
 tm[r] = If[tmax[r] < 9 Gyr/Myr, tmax[r], 9 Gyr/Myr]; 
 p[r] = {Plot[a[t] /. sol[r], {t, 0, tm[r]}, 
    FrameLabel -> {"Time, Myr", "Semi-major Axis, AU"}, Frame -> True,
     PlotLabel -> Row[{"Rast = ", r 1.}], PlotRange -> All], 
   Plot[e[t] /. sol[r], {t, 0, tm[r]}, 
    FrameLabel -> {"Time, Myr", "Eccentricity"}, Frame -> True, 
    PlotLabel -> Row[{"tmax = ", N[tm[r] 1., 3]}], 
    PlotRange -> All]};, {r, 10^-3, 1/100, 10^-3}]
Table[p[r], {r, 10^-3, 1/100, 10^-3}]