It was (apparently) indeed the dispersion that was wrong. The function

norm_loglik2 <- function(mod, newdata) {
  sum(dnorm(
    newdata$y,
    predict(mod, newdata, type = "response"),
    sqrt(mod$deviance / nobs(mod)),
    log = TRUE
  ))
}

yields the same likelihood as logLik .