Slightly different than what you asked for, but have you tried to combine all dfs and doing a faceted plot?

library(ggplot2)
library(patchwork)

# Create different datasets for each plot
df1 <- expand.grid(x = seq(300, 800, length.out = 50), y = seq(300, 600, length.out = 50))
df1$z <- with(df1, dnorm(x, mean = 500, sd = 50) * dnorm(y, mean = 400, sd = 50))

df2 <- expand.grid(x = seq(300, 800, length.out = 50), y = seq(300, 600, length.out = 50))
df2$z <- with(df2, dnorm(x, mean = 600, sd = 50) * dnorm(y, mean = 450, sd = 50))

df3 <- expand.grid(x = seq(300, 800, length.out = 50), y = seq(300, 600, length.out = 50))
df3$z <- with(df3, dnorm(x, mean = 550, sd = 50) * dnorm(y, mean = 500, sd = 50))

df4 <- expand.grid(x = seq(300, 800, length.out = 50), y = seq(300, 600, length.out = 50))
df4$z <- with(df4, dnorm(x, mean = 650, sd = 50) * dnorm(y, mean = 350, sd = 50))

# Compute global min and max for z-values across all datasets
min_z <- min(c(df1$z, df2$z, df3$z, df4$z), na.rm = TRUE)
max_z <- max(c(df1$z, df2$z, df3$z, df4$z), na.rm = TRUE)

df.grouped <- dplyr::bind_rows(list(df1=df1, df2=df2, df3=df3, df4=df4), .id = 'source')

head(df.grouped)

ggplot(df.grouped, aes(x, y, fill = z)) +
  geom_raster() +
  scale_fill_viridis_c(limits = c(min_z, max_z)) + 
  labs(y = "Excitation Wavelength / nm",
       x = "Emission Wavelength / nm") +
  facet_wrap(~source, scales = "free")+
  theme_classic()+
  theme(strip.text = element_blank())