We specified Bayesian hierarchical Poisson log-linear models to investigate the association of COVID-19 deaths and NO2 and PM2.5 independently. The LSOA specific standardised mortality ratio is known to be an unstable estimator with high variance when the number of expected deaths is small. To overcome this problem, we used a well-established hierarchical framework, specifying spatially structured and unstructured random effects, so that the model borrows strength from the other areas across the entire study region, as well as from the neighbouring ones (Best et al., 2005, Wakefield et al., 2000, Wakefield, 2006). We model these random effects using a re-parametrisation of the Besag-York-Molliè conditional autoregressive prior distribution (Besag et al., 1991, Simpson et al., 2017). We fitted four models including: 1) each pollutant (model 1), 2) each pollutant and the spatial autocorrelation term (model 2), 3) each pollutant and all confounders (model 3) and 4) each pollutant, the spatial autocorrelation term and all confounders (model 4). All models were adjusted for age, sex and ethnicity using indirect standardisation; we used the English population as the standard population to calculate the rates. We do not report results from the joint analysis including both pollutants since they are highly correlated (Supplemental Material Figure S5).