The root cause, as you no doubt have guessed, is that a part of the Alaska archipelago crosses the antimeridian.
sf::st_crop() will be your friend here, and you can apply it either in original CRS (your limits look metric and the plot has the look of the Web Mercator, which is the CRS of the popular Kontur dataset, but I speculate here) or in WGS 84 (= decimal degrees).
But you have to set the limits in units your object understands = to apply a degree based crop you must have your object in degree based CRS. Or use a meters crop for a metric CRS. A slight complication is that metric coordinates are so unfamiliar, and decimal degrees so familiar, that {ggplot2} likes to annotate even metric plots by decimal degrees. A confusion we must live with.
To get your object from one CRS to another sf::st_transform() is your friend.
Consider using a more friendly projection. It's getting late for me, so I haven't done the boundary multipolygon overlay. Just code and image of plot, too long running for reprex. You might also consider binning population to create a discrete scale.
require(ggplot2)
require(sf)
# https://data.humdata.org/dataset/kontur-population-united-states-of-america?
boundaries <- st_read("/Users/ro/projects/demo/kontur_boundaries_US_20220407.gpkg")
population <- st_read("/Users/ro/projects/demo/kontur_population_US_20220630.gpkg")
# identify source data projection
st_crs(boundaries)[1]
st_crs(population)[1]
# Mercator type projections are distorted
# consider changing to NAD83, Alaska Albers
# https://epsg.io/3338
boundaries <- st_transform(boundaries, 3338)
population <- st_transform(population, 3338)
# subset out alaska boundary
alaska <- boundaries[which(boundaries$name_en == "Alaska"),]
# subset out alaska population
ak_pop <- population[st_intersects(alaska,population)[1][[1]],]
# plot
ggplot(ak_pop) +
geom_sf(aes(fill = population)) +
scale_fill_continuous(type = "viridis") +
theme_minimal()