We present a robust approach to incorporating models for the time-variable broadening of the pulse profile due to scattering in the ionized interstellar medium into profile-domain pulsar timing analysis. We use this approach to simultaneously estimate temporal variations in both the dispersion measure (DM) and scattering, together with a model for the pulse profile that includes smooth evolution as a function of frequency, and the pulsars timing model. We show that fixing the scattering timescales when forming time-of-arrival estimates, as has been suggested in the context of traditional pulsar timing analysis, can significantly underestimate the uncertainties in both DM, and the arrival time of the pulse, leading to bias in the timing parameters. We apply our method using a new, publicly available, GPU accelerated code, both to simulations, and observations of the millisecond pulsar PSR J1643$-$1224. This pulsar is known to exhibit significant scattering variability compared to typical millisecond pulsars, and we find including low-frequency ($< 1$ GHz) data without a model for these scattering variations leads to significant periodic structure in the DM, and also biases the astrometric parameters at the $4sigma$ level, for example, changing proper motion in right ascension by $0.50 pm 0.12$. If low frequency observations are to be included when significant scattering variations are present, we conclude it is necessary to not just model those variations, but also to sample the parameters that describe the variations simultaneously with all other parameters in the model, a task for which profile domain pulsar timing is ideally suited.
Download