# SCRIPT: BH-DRIFT-SYM-02
# Claim: for L s = -Δs + 2<X,∇s> + Q s with an EXACT drift X = ∇ψ, the weighted similarity
# e^{-ψ} L e^{ψ} = -Δ + (Q + |∇ψ|² - Δψ) is self-adjoint (no first-order term).
# Kill condition: any first-order term surviving, or a different zeroth-order term.
# Second claim: for a NON-exact drift (X = (-y, x), a rotation) no function ψ removes it:
# the conjugated first-order coefficient is 2(X - ∇ψ), zero only if X is a gradient (curl X = 0).
import sympy as sp
x, y = sp.symbols('x y')
psi, f, Q = sp.Function('psi')(x, y), sp.Function('f')(x, y), sp.Function('Q')(x, y)
lap = lambda u: sp.diff(u, x, 2) + sp.diff(u, y, 2)
grad = lambda u: (sp.diff(u, x), sp.diff(u, y))
def L(s, X):
    gs = grad(s); return -lap(s) + 2*(X[0]*gs[0] + X[1]*gs[1]) + Q*s
Xg = grad(psi)
conj = sp.expand(sp.simplify(sp.exp(-psi)*L(sp.exp(psi)*f, Xg)))
target = sp.expand(-lap(f) + (Q + Xg[0]**2 + Xg[1]**2 - lap(psi))*f)
print('exact drift: e^{-psi} L e^{psi} - [-Δ + (Q+|∇psi|^2-Δpsi)] =', sp.simplify(conj - target))
# rotational drift: conjugate with arbitrary psi, read first-order coefficients
Xr = (-y, x)
conj_r = sp.expand(sp.simplify(sp.exp(-psi)*L(sp.exp(psi)*f, Xr)))
c_fx = sp.simplify(conj_r.coeff(sp.Derivative(f, x)))
c_fy = sp.simplify(conj_r.coeff(sp.Derivative(f, y)))
print('rotational drift, first-order coefficients after conjugation:', c_fx, ',', c_fy)
print('curl of (-y,x) =', sp.diff(x, x) - sp.diff(-y, y), ' (nonzero => no psi can remove both)')