Overlaying shapefile over multiple raster plots in R

I have multiple NetCDF data and I need to plot only one component of each data. I have achieved this in plotting all the data. But now I am unable to overlay my shapefile over the plots created. Is there any way to overlay the shapefile over the maps I have plotted using levelplot? I tried it with spplot but still, I am not getting the shapefile over my maps. I also tried using sp.polygons but can still not make it. Here is a snippet of my code and the generated plot.

library(ncdf4)
library(raster)
library(rgdal)
library(rasterVis)
library(RColorBrewer)
library(colorRamps)
path<- setwd('D:/DATA/2015/MIXING_RATIO/')

# Set the path to the shapefile
shapefile_path <- "C:/Users/IITM/Downloads/India_State/India_State/India_State.shp"

# Load the shapefile
shapefile <- readOGR(shapefile_path)
# Get a list of all netCDF files in the directory
nc_files <- list.files(path, pattern = ".nc", full.names = TRUE)
raster_list <- list()
# Loop through each netCDF file
for (nc_file in nc_files) {
  # Read the netCDF data
  data <- nc_open(nc_file)
  
  # Extract the necessary variables
  lon <- ncvar_get(data, "lon")
  lat <- ncvar_get(data, "lat")
  dust_MR <- ncvar_get(data, "DU001")
  
  # Plot the first level slice of the data
  dust1.slice1 <- dust_MR[, , 60]
  
  # Create a raster object
  r <- raster(t(dust1.slice1), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat), crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
  
  # Flip the raster vertically
  r <- flip(r, direction = 'y')
  
  # Crop the raster based on the shapefile extent
  r_crop <- crop(r, extent(shapefile))
  
  # Mask the raster based on the shapefile geometry
  r_mask <- mask(x = r_crop, mask = shapefile)
  
  # Add the raster to the list
  raster_list[[nc_file]] <- r_mask
  
  # Close the netCDF file
  nc_close(data)
}


# Stack the rasters
raster_stack <- stack(raster_list)
crs(raster_stack)
DUST<- stack(mget(rep("raster_stack")))
names(DUST)
raster_names<- gsub('D..DATA.2015.MIXING_RATIO.MERRA2_400.inst3_3d_aer_Nv.|.SUB.nc'," ",names(DUST))
raster_names
raster_dates <- as.Date(trimws(raster_names), format = "%Y%m%d")
raster_dates
txtdates <- format(raster_dates, "%d%b%Y")

plot(raster_stack, nc = 4)
# Create a single plot using levelplot
p<-levelplot(raster_stack,names.attr=txtdates,
          margin = FALSE,colorkey = list(space = "right"), main = "Dust Mixing Ratio (bin 1)(Kg/Kg)")
p + layer(sp.polygons(polygon, lwd=0.8, col='darkgray'))

Thank you for any suggestions.

This is another version of the answer by Robert Hijmans using rasterVis. The code includes a function, llines.SpatVector, that will be included in rasterVis after next update of the lattice package (more information here).

library(terra)
library(latticeExtra)
library(rasterVis)

llines.SpatVector <- function(x, ...) {
  xy <- crds(x, list=TRUE)
  names(xy) <- c("x", "y")
  lattice::llines(xy, ...)
}

r <- rast(system.file("ex/elev.tif", package="terra")) * 1:2
v <- vect(system.file("ex/lux.shp", package="terra"))
lns <- as.lines(v)

levelplot(r) + layer(llines(lns))