IMHO solving for differential Faraday rotation should be included in one way or another in the solutions for the phase self-cal.
I think (hope?) that with enough bandwidth something like "CommonFaradayRotation" (doesn't exist yet in BBS) should be an useful parameter.
The questions is how to solve for it, very likely the BBS solver is not stable solving for CommonFaradayRotation even when providing significant bandwidth (due to local minima in Chi^2 space). I've done some very easily testing a few years ago on 3C295 using the full bandwidth and that did not look very stable.
dFR is only large on the longest baselines (if the direction independent selfcal is run at 20-30 arcsec resolution the effects is likely too small and can be neglected, this becomes more important for longer baselines, those are giving the larger dFR signal). So any dFR solve will need the long baselines. Not clear if this is achievable in pre-factor.
An option would be to solve it on the brightest source(s) in the DDE part and then assume it is constant across the FoV (this might be quite a good approximation, unlike dTEC). But again the question is -how- to reliable solve for it. I can imagine that you can first solve RR, LL phases only then do a manual dFR fit in python on the solutions (like in Clock-TEC fitting, dFR this is already implemented) and then refit with that as initial guesses in BBS.
Lots of testing will need to be done before I would trust any results.
IMHO solving for differential Faraday rotation should be included in one way or another in the solutions for the phase self-cal. I think (hope?) that with enough bandwidth something like "CommonFaradayRotation" (doesn't exist yet in BBS) should be an useful parameter.
Andreas