Well, you sucked me in -- clever boy (or girl).

I think the [$]\beta[$] notation is a distraction that overly complicates things.

If we simply write

[$] dS_t = \sqrt{V_t} S_t dB_t[$], where [$]V_t[$] is some diffusion (time-homogeneous or not, driven by a Brownian uncorrelated to [$]B_t[$]), then it seems to me that

(M2) [$]C_t(t,T,S_t,V_t) = \left< \frac{S_t \, e^{-d_1^2/2}}{2 \sqrt{2 \pi \, U(t,T)}} \, U_t \right> = - S_t V_t \left< \frac{e^{-d_1^2/2}}{2 \sqrt{2 \pi \, U(t,T)}} \right>[$],

since [$]U(t,T) = \int_t^T V_s \, ds[$].

Proof: just differentiate [$]C_{BS}[$] w.r.t. t inside the brackets, so no different than a proof of my earlier (*) under deterministic volatility.

To implement (M2), you just simulate many V-path's, all starting at [$]t[$] with the same [$]V_t[$], evaluate [$]U(t,T)[$] for each such path, and take the MC average of the final expression in brackets. It's not much different than using mixing to evaluate [$]C[$] itself, where you do the same thing but in that case the "expression in brackets" is the BS formula. So, each MC run can yield both [$]C[$] and [$]C_t[$] with little additional work, as you have to evaluate [$]U(t,T)[$] and [$]d_1[$] for the BS formula anyway.

Again, if doubtful, I'd say try it -- but I have noticed you are reluctant for some reason to post numerical examples.