Eliminating NDSolve Instability

pde = {D[Q[t, ϕ, θ], t] == -Cos[θ]*D[Q[t, ϕ, θ], ϕ] + Sin[θ]*D[Q[t, ϕ, θ], θ, ϕ]/2,
Q[t, 0, θ] == Q[t, 2*Pi, θ], Q[t, ϕ, 0] == 0,
Q[0, ϕ, θ] == ((1 - Cos[ϕ]*Sin[θ])^2 - 1)/4};
sol = NDSolveValue[pde, Q, {t, 0, 5}, {ϕ, 0, 2*Pi}, {θ, 0, Pi}, MaxStepSize -> {.01, .25, .25}]

generates the message

Warning: scaled local spatial error estimate of 7560.422997509231 at t = 5. in the direction of independent variable \[Phi] is much greater than the prescribed error tolerance. Grid spacing with 27 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options. >>

and exhibits rapid numerically unstable growth. By t = 5, sol is as large as 4.23527*10^8, even though the analytical solution is known to be less than one in absolute value.

The warning message not withstanding, increasing resolution in ϕ or θ increases the instability growth rate. I have tried a number of Methods, including "Adams", "BDF", "ExplicitRungeKutta", "ImplicitRungeKutta", and "Extrapolaton" without substantive improvement.

I would appreciate advice on how to obtain a stable solution.

Note: Although this particular problem is solvable analytically, it is meant as a prototype for more general problems involving D[Q[t, ϕ, θ], θ, ϕ] that are not solvable analytically.