wilmott.com

oops -- my bad

It could happen to a bishop.

What does Fichera say?

This brings us back to the favourite topic of you two, characteristics!

What does Fichera say?

Fichera is only for degenerate elliptic/parabolic pde.
Goursat pde seems to be neither. But how do I classify it as ....?

This brings us back to the favourite topic of you two, characteristics!

It's a term from an elliptic pde, I see as adding some anisotropic terms to dffusion. Maybe they are hiding there.
But characteristic direction is important.

Looks as if this problem has everyone stumped..

Maybe the answer is that your problem is always ill-posed. 

Suppose there was a time-independent solution [$]u(x,y)[$], which is then a solution to a 2D hyperbolic problem. After 45 degree coordinate rotation to [$](\xi,\eta)[$], the characteristics leaving the origin are [$]\xi = \pm \eta[$]. Then, rotating back, aren't the characteristics in the original coordinates just the coordinate axes?  

Now, in general, to make any hyperbolic problem well-posed, you have to specify data along curves which are non-characteristic. 

But, if the horizontal and vertical lines are exactly the characteristics here, then you are asking for illegal boundary conditions. 

This argument suggests the time-independent problem is ill-posed. It makes one suspicious about the time-dependent problem.

Not to mention that it has aspects of solving heat eqn in the wrong direction!

Alan,
I disagree with most of what you have written in this case. I would first forget the 45 degrees rotation. It does not give insights, at least not to me.

The (time-independent) Goursat problem is well-posed (check Riemann method + function)

https://www.encyclopediaofmath.org/index.php/Goursat_problem

https://mathworld.wolfram.com/GoursatProblem.html

Here's Exhibit C (existence and uniqueness)

http://www.kjm-math.org/article_65161_8fed2a0e71a8af7d76b734e82401f3b2.pdf

And numerically, it works well as I have already mentioned. If it were ill-posed you would see it in no time in the numerical solution. 
So, I need to make the pure maths look respectable.

I was brainstorming on

1. Separation of variables (Sturm-Liouville)
2. Wirtinger dervivatives (complex derivatives).

Not to mention that it has aspects of solving heat eqn in the wrong direction!

This is a completely different kettle of fish. Very few analogies here. Heat is parabolic/elliptic.

If ill-posed, a proof is needed.

Eg rho positive, initial data f(x-y). Soln of form g(x-y, t).

Alan,
I disagree with most of what you have written in this case. I would first forget the 45 degrees rotation. It does not give insights, at least not to me.

The (time-independent) Goursat problem is well-posed (check Riemann method + function)

https://www.encyclopediaofmath.org/index.php/Goursat_problem

https://mathworld.wolfram.com/GoursatProblem.html

Here's Exhibit C (existence and uniqueness)

http://www.kjm-math.org/article_65161_8fed2a0e71a8af7d76b734e82401f3b2.pdf

And numerically, it works well as I have already mentioned. If it were ill-posed you would see it in no time in the numerical solution. 
So, I need to make the pure maths look respectable.

Fair enough. I looked up some of that Goursat problem stuff and see that of two things I said, one was right and one was wrong. The right thing was the rotation argument and so, indeed, the characteristics of [$]u_{xy} = 0[$] are lines parallel to the coordinate axes. The wrong thing was saying that you would get an ill-posed problem if you specified data along those characteristics. Now I see that you can. 

For example, just specify [$]u(0,y) = g(y)[$] for [$]0 \le y \le 1[$] and  [$]u(x,0) = f(x)[$] for [$]0 \le x \le 1[$], with compatibility at the origin. Then, the pde solution anywhere in [$]0 \le (x,y) \le 1[$] (your square) is given by [$]u(x,y) = f(x) + g(y) - u(0,0)[$], where [$]u(0,0) = f(0) = g(0)[$]. 

Geometrically, given [$](x,y)[$] you just draw a horizontal line and vertical line to hit the axes, sum the two values you hit, and subtract the (common) corner value. By symmetry, I will guess you can pick any two adjacent edges of the square and specify [$]u[$] there. Then, to get the solution at [$](x,y)[$] you again just draw a horizontal and vertical to hit the two edges you have chosen for the prescribed data, again sum the two values and subtract the corner value. Note the values on the non-chosen edges are determined by the values on the chosen edges.

What this means for the full time-dependent problem, I don't know. A wild guess for that one might be that it makes sense if you specify [$]u(t,0,y) = g(t,y)[$] and [$]u(t,x,0) = f(t,x)[$] or, similarly on two adjacent edges. But what about an initial value? Arbitrary as long as compatible with the chosen edge values at [$]t=0[$]?

The last paragraph assumes rho positive, and the usual 'initial' -> 'terminal' (t -> T-t)  when rho is negative. 
  

Nice answer!
I've be looking in Polyanin an Zaitsev where they use  self-similarity ansatz [$]\xi = xy[$] for the sine-Gordon equation. In our case the time-dependent pde

[$]\frac{\partial u}{\partial t} = \rho  \frac{\partial^2 u}{\partial x \partial y}[$], [$]0 \lt x,y \lt 1[$]

 now becomes

[$]\frac{\partial u}{\partial t} = \rho (\xi \frac{\partial^2 u}{\partial \xi^2} + \frac{\partial u}{\partial \xi})[$], [$]0 \lt \xi \lt 1[$].

Looks like a degenerate parabolic pde for [$]\rho \gt 0[$]. But for [$]\rho \lt 0[$]?

When [$]\xi = 0[$] the pde becomes

[$]\frac{\partial u}{\partial t} = \rho \frac{\partial u}{\partial \xi}[$].

BTW if it makes life easier you can use

[$]\frac{\partial u}{\partial t} = \rho  \frac{\partial^2 u}{\partial x \partial y} + b_1\frac{\partial u}{\partial x} + b_2\frac{\partial u}{\partial y} + cu + f(x,t)[$].

Sigh...