It seems reasonable that the mesh needs to be very fine only near the boundary, where the initial datum is discontinuous. You may try the option MeshRefinementFunction.

The purpose to make the input exact instead of floating-point is good. In your case it could be done more easily by giving exact parameters at the start. Instead of floating point values, as in

\[Eta] = 5.; 
\[Kappa] = .75; 
Tmax = 5.;

I would suggest

\[Eta] = 5; 
\[Kappa] = 3/4; 
Tmax = 5;

This way there is no need for Rationalize. Even more: with your original floating-point \[Eta] = 5. the triangle definition

\[CapitalOmega] = 
  Polygon[Rationalize[{{0, 0}, {\[Eta], 
      0}, {\[Eta]/2, (\[Eta]   Sqrt[3])/2}}, 
    0]]

actually decreases the precision, because \[Eta] Sqrt[3])/2 becomes first floating-point, and then Rationalize freezes the approximation instead of the exact value.

Hi Cameron, I agree with you that your formulation of the problem and your PDEs look correct. If you plot your solution over time at one of the boundary points such as P(0,0,t), it looks like Mathematica is setting the boundary conditions to be P(0,0,t)*e^(-t) with the coefficient (-1) being exact. This is quite annoying, but you can remedy it by explicitly setting the boundary condition to P(0,0,t)*e^(-10000*t).

bC = Rationalize[
   DirichletCondition[
    P[x1, x2, t] == Exp[-10000  t] 1/((Sqrt[3]   \[Eta]^2)/4), True], 
   0]; 

In the study of linear PDEs this discontinuity (setting P(x,y,0)=const. on the domain and P(x,y,t)=0 on the boundary) isn't much of an issue and it's standard to do, so it's frustrating that Mathematica makes such a big deal about it! At least with this fix applied you can see that the derivative of the total probability mass approaches -infinity near t=0, so that hints at some potential numerical difficulties.

Update 2: My formulation of the differential equation matches that produced automatically by Mathematica

eqn = ItoProcess[{\[DifferentialD]x1[
      t] == \[Kappa] (x1opt - 
        x1[t])  \[DifferentialD]t + \[DifferentialD]w1[
       t], \[DifferentialD]x2[
      t] == \[Kappa] (x2opt - 
        x2[t])   \[DifferentialD]t + \[DifferentialD]w2[t]}, {x1[t], 
   x2[t]}, {{x1, x2}, {x10, x20}}, 
  t, {w1 \[Distributed] WienerProcess[], 
   w2 \[Distributed] WienerProcess[]}]

eqn["KolmogorovForwardEquation"] // TraditionalForm 

I am totally ignorant on numerical methods for PDEs, but I wonder what accuracy you can expect with discontinuous initial data.

That's a great point regarding rationalising the parameters only, I'll be instituting that.

I'll have to think a bit more about how to square your point about spatial derivatives with the fact that mass evidently does decrease within the domain. I'll get back when I have a decisive answer, but one thing to note is that the flux of probability is in the opposite direction to the gradient of $$p$$

Thanks David, this is really great input! I think this along with the comments above practically solves the problem!