Return to Answer

Remark

Adding PrecisionGoal -> 3, AccuracyGoal -> 3 to NDSolveValue will reduce the timing to about 277.04 seconds! But the beginning part of $S_1$ (for about $t<0.2$) changes. (Not widely different, but quite observable. ) Anyway, if the accuracy of transient behavior isn't important for you, consider adding these options.

V[ϕ_, θ_] := 
  0.045360921162651446/Sqrt[1.00001 + Cos[ϕ] Sin[θ]] + 
   0.045360921162651446/Sqrt[1.00001 - Sin[θ] Sin[π/6 - ϕ]] + 
   0.045360921162651446/Sqrt[1.00001 - Sin[θ] Sin[π/6 + ϕ]];

Clear[points, grid, domain]
points@ϕ = points@θ = 100;
eps = 0(*Pi/10^2*);
domain@ϕ = {0, 2 Pi}; domain@θ = {eps, Pi - eps};
(grid@# = Array[# &, points@#, domain@#]) & /@ {ϕ, θ};

difforder = 2;
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = 
  pdetoode[{Subscript[S, 0], Subscript[S, 1]}[ϕ, θ, t], t, 
   grid /@ {ϕ, θ}, difforder, {True, False}];
μ0 = 10^-2; μ1 = 5 10^-2;
Unevaluated[
   eq = {μ0 (D[Subscript[S, 0], {θ, 2}] + 
          D[Subscript[S, 0], {ϕ, 2}]) - 
       D[Subscript[S, 0], 
        t] == -V[ϕ, θ] + (D[Subscript[S, 0], θ]^2 + 
         Csc[θ]^2 D[Subscript[S, 0], ϕ]^2), μ1 (D[Subscript[S, 
           1], {θ, 2}] + D[Subscript[S, 1], {ϕ, 2}]) - 
       D[Subscript[S, 1], 
        t] == -(-Cot[θ] D[Subscript[S, 0], θ] + 
          D[Subscript[S, 0], {θ, 2}] + 
          Csc[θ]^2 D[Subscript[S, 0], {ϕ, 2}]) + 
       2 (D[Subscript[S, 0], θ] D[Subscript[S, 1], θ] + 
          Csc[θ]^2 D[Subscript[S, 0], ϕ] D[Subscript[S, 1], ϕ])};
   ic = {Subscript[S, 0] == Sin[θ]^2, Subscript[S, 1] == Sin[θ]^2} /. 
     t -> 0;
   bc = {D[Subscript[S, 0], θ] == 0, 
      D[Subscript[S, 1], θ] == 
       0} /. {{θ -> domain[θ][[1]]}, {θ -> 
        domain[θ][[-1]]}}] /. {Subscript[S, 0] -> 
    Subscript[S, 0][ϕ, θ, t], 
   Subscript[S, 1] -> Subscript[S, 1][ϕ, θ, t]};
del = #[[2 ;; -2]] &;
ode = Map[del, ptoofunc@eq, {2}] // Quiet;

odebc = ptoofunc@diffbc[t][bc];
odeic = ptoofunc@ic;

varlst = Outer[#@##2 &, {Subscript[S, 0], Subscript[S, 1]}, grid@ϕ, 
   grid@θ];

tend = 10; 

showStatus[status_]:=LinkWrite[$ParentLink,
  SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],ToString[status]]];
clearStatus[]:=showStatus[""];
clearStatus[]
jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]]
 sollst = NDSolveValue[{ode, odeic, odebc}, varlst, {t, 0, tend}, jianshi[t], 
   Method -> {"EquationSimplification" -> "MassMatrix"}(*,SolveDelayed->
   True*)]; // AbsoluteTiming
(* {251.958, Null} *)

{sols0, sols1} = rebuild[#, grid /@ {ϕ, θ}, -1] & /@ sollst;

piclst1 = Table[
   SphericalPlot3D[sols1[f, th, t], {th, #3, #4}, {f, #1, #2}, PlotPoints -> 50] & @@
     Flatten[domain /@ {ϕ, θ}], {t, 0, tend, tend/25}];

ListAnimate[piclst1]

V[ϕ_, θ_] := 
  0.045360921162651446/Sqrt[1.00001 + Cos[ϕ] Sin[θ]] + 
   0.045360921162651446/Sqrt[1.00001 - Sin[θ] Sin[π/6 - ϕ]] + 
   0.045360921162651446/Sqrt[1.00001 - Sin[θ] Sin[π/6 + ϕ]];

Clear[points, grid, domain]
points@ϕ = points@θ = 200;
eps = 0(*Pi/10^2*);
domain@ϕ = {0, 2 Pi}; domain@θ = {eps, Pi - eps};
(grid@# = Array[# &, points@#, domain@#]) & /@ {ϕ, θ};

difforder = 2;
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = 
  pdetoode[{Subscript[S, 0], Subscript[S, 1]}[ϕ, θ, t], t, 
   grid /@ {ϕ, θ}, difforder, {True, False}];
μ0 = 10^-2; μ1 = 10^-2;
Unevaluated[
   eq = {μ0 (D[Subscript[S, 0], {θ, 2}] + 
          D[Subscript[S, 0], {ϕ, 2}]) - 
       D[Subscript[S, 0], 
        t] == -V[ϕ, θ] + (D[Subscript[S, 0], θ]^2 + 
         Csc[θ]^2 D[Subscript[S, 0], ϕ]^2), μ1 (D[Subscript[S, 
           1], {θ, 2}] + D[Subscript[S, 1], {ϕ, 2}]) - 
       D[Subscript[S, 1], 
        t] == -(-Cot[θ] D[Subscript[S, 0], θ] + 
          D[Subscript[S, 0], {θ, 2}] + 
          Csc[θ]^2 D[Subscript[S, 0], {ϕ, 2}]) + 
       2 (D[Subscript[S, 0], θ] D[Subscript[S, 1], θ] + 
          Csc[θ]^2 D[Subscript[S, 0], ϕ] D[Subscript[S, 1], ϕ])};
   ic = {Subscript[S, 0] == Sin[θ]^2, Subscript[S, 1] == Sin[θ]^2} /. 
     t -> 0;
   bc = {D[Subscript[S, 0], θ] == 0, 
      D[Subscript[S, 1], θ] == 
       0} /. {{θ -> domain[θ][[1]]}, {θ -> 
        domain[θ][[-1]]}}] /. {Subscript[S, 0] -> 
    Subscript[S, 0][ϕ, θ, t], 
   Subscript[S, 1] -> Subscript[S, 1][ϕ, θ, t]};
del = #[[2 ;; -2]] &;
ode = Map[del, ptoofunc@eq, {2}] // Quiet;

odebc = ptoofunc@diffbc[t][bc];
odeic = ptoofunc@ic;

varlst = Outer[#@##2 &, {Subscript[S, 0], Subscript[S, 1]}, grid@ϕ, 
   grid@θ];

tend = 10; 

showStatus[status_]:=LinkWrite[$ParentLink,
  SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],ToString[status]]];
clearStatus[]:=showStatus[""];
clearStatus[]
jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]]
 sollst = NDSolveValue[{ode, odeic, odebc}, varlst, {t, 0, tend}, jianshi[t], 
   Method -> {"EquationSimplification" -> "MassMatrix"}(*,SolveDelayed->
   True*)]; // AbsoluteTiming
(* {1984.24, Null} *)

{sols0, sols1} = rebuild[#, grid /@ {ϕ, θ}, -1] & /@ sollst;

piclst1 = Table[
   SphericalPlot3D[sols1[f, th, t], 
     {th, #3, #4}, {f, #1, #2}, PlotPoints -> 50] & @@
       Flatten[domain /@ {ϕ, θ}], {t, 0, tend, tend/25}];

ListAnimate[
 Show[#, PlotRange -> {{-1.2, 1.2}, {-1.2, 1.2}, {-3.2, 3.2}}] & /@ piclst1]

I've solved $S_0$ and $S_1$ all in one, but since firsteq only involves $S0$ (and further trial suggests it's easier to solve than secondeq), you may consider solve firsteq first.

I've solved $S_0$ and $S_1$ all in one, but since firsteq only involves $S_0$ (and further trial suggests it's easier to solve than secondeq), you may consider solve firsteq first.