I'm not certain that this is the intended reason from the progenitor of the deviance measure, but it seems likely that the reason for the factor of two is to ensure that the deviance corresponds to an appropriate scaled version of the sum-of-squares in the normal case. That is, in the case where you have a GLM using a normal distribution with variance $\sigma^2$ (i.e., standard regression) you get deviance values:
$$\begin{equation} \begin{aligned}
D_{TOT} = 2(\hat{\ell}_{s} - \hat{\ell}_0) &= \frac{1}{\sigma^2} \sum_{i=1}^n (y_i - \bar{y})^2 = \frac{1}{\sigma^2} \cdot SS_{TOT}, \\[6pt]
D_{REG} = 2(\hat{\ell}_{p} - \hat{\ell}_0) &= \frac{1}{\sigma^2} \sum_{i=1}^n (\hat{y}_i - \bar{y})^2 = \frac{1}{\sigma^2} \cdot SS_{REG}, \\[6pt]
D_{RES} = 2(\hat{\ell}_{s} - \hat{\ell}_{p}) &= \frac{1}{\sigma^2} \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \frac{1}{\sigma^2} \cdot SS_{RES}, \\[6pt]
\end{aligned} \end{equation}$$
where the log-likelihoods under the null model, actual model, and saturated model, are given respectively by:
$$\begin{equation} \begin{aligned}
\hat{\ell}_0 = \ln p(y|\hat{\theta}_0) \quad \quad \quad
\hat{\ell}_p = \ln p(y|\hat{\theta}_p) \quad \quad \quad
\hat{\ell}_s = \ln p(y|\hat{\theta}_s).
\end{aligned} \end{equation}$$
If you go through the algebraic derivation of these results you will see that the factor of two in the deviance cancels out with the $\frac{1}{2}$ in the exponential for the normal density. Presumably the reason to have this factor is to define the deviance measures in a way that naturally extends the sums-of-squares in standard Gaussian regression.
