I'm just making this up on the fly. There's got to be at least one paper out there on this, or sections in Kalman filtering books.
Supposing that Δtk is distributed normally, does this just come out naturally as disturbance? i.e. could I account for this with my Kalman gain
Nope.
Take the continuous-time model of the system:
$$\begin{aligned}
\dot x(t) &= \mathbf A_c x(t) + \mathbf B_c u(t) \\
y(t) &= \mathbf C x(t)
\end{aligned} \tag 1$$
For a sample at $t = T_s k + \Delta t$ and for a $\dot x(t)$ that doesn't change significantly over a span of $\Delta t$, $y_k$ is close to $y_k \simeq \mathbf C \left(x(T_s k) + \Delta t\ \dot x(T_s k)\right)$ or
$$y_k \simeq \mathbf C \left(x(T_s k) + \Delta t\ \mathbf A_c\ x(T_s k)\right) \tag 2$$.
That condition on $\dot x(t)$ can be satisfied if the eigenvalues of $\mathbf A_c$ are significantly smaller than $\frac{1}{\Delta t}$.
So the effect of $\Delta t$ is multiplicative.
A good start on this, if your $\mathbf A_c$ meets the eigenvalue criterion above and if the pdf of $T_s$ is Gaussian with deviation $\sigma_{T_s}$, is to make an unscented Kalman, where everything is "normal" Kalman filtering except that at each step you compute
$$\mathbf R_k = \mathbf R_m + \sigma_{T_s}^2 \mathbf C \mathbf A_c x_{k-1} x_{k-1}^T \mathbf A_c^T \mathbf C^T, \tag 3$$
where $\mathbf R_m$ is the "actual" measurement noise and $\sigma_{T_s}^2 \mathbf C \mathbf A_c x_{k-1} x_{k-1}^T \mathbf A_c^T \mathbf C^T$ is an estimate of the measurement noise added by the variation in the sampling time.
Note that you could, at the cost of longer settling times, just take $\mathbf R$ as constant, and equal to some worst-case version of (3) assuming an as-high-as-possible $\dot x(t)$. This would get you back to a plain-old Kalman though.
Suppose that Δtk is not distributed normally, and instead has an asymmetric distribution, possibly even pareto. Is there any way to account for this?
Yes, but. To the extent that the probability of $\dot x(t)$ changes much over the worst-case $\Delta t$, the approxmation in (2) holds -- it's just that the resulting measurement error has a distribution that's close to the distribution of $\Delta t$. However, you now have noise that's not just multiplicative, but which is non-Gaussian.
If your filtering is still close enough using a Gaussian assumption for your non-Gaussian noise, then (3) should still work. If not, then you'd have to go to something fancier, like Baysian filters or particle filters.
