Fundamental misunderstanding with crs and map projections in R

The CRS used in your picture looks like EPSG:3184, which is the appropriate longitudinal slice of the Universal Transverse Mercator system. As long as you tell coord_sf this is the co-ordinate reference system that you want to replicate, it will automatically convert your map to that system for plotting. You can do this as follows:

library(rnaturalearth)
library(rnaturalearthdata)
library(ggplot2)
library(sf)

north <- 72.5; south <- 59; east <- -52; west <- -45

extent <- st_sf(a = 1:2, crs = 4326,
            geom = st_sfc(st_point(c(east, south)), 
                          st_point(c(west, north)))) |>
            st_transform(crs = 'EPSG:3184') |>
            st_bbox()

lat <- lapply(c(60, 62, 64, 66, 68), 
            function(x) cbind(seq(-70, -40, 0.1), x)) |>
  st_multilinestring() |> st_sfc(crs = 4326) |> st_sf() |>
  st_transform(crs = 'EPSG:3184') 

lon <- lapply(c(-60, -55, -50, -45, -40), 
              function(x) cbind(x, seq(56, 74, 0.2))) |>
  st_multilinestring() |> st_sfc(crs = 4326) |> st_sf() |>
  st_transform(crs = 'EPSG:3184') 
  
ne_countries(scale = "large", country = "Greenland",returnclass = "sf") |>
  ggplot() + 
  geom_sf(fill = 'white') +
  geom_sf(data = lat, linetype = 2, color = 'gray50') +
  geom_sf(data = lon, color = 'gray', linewidth = 0.2) +
  coord_sf(crs = 'EPSG:3184', xlim = extent[c(1, 3)], ylim = extent[c(2, 4)]) +
  scale_x_continuous(breaks = -5 * 1:10) +
  theme_bw() +
  theme(panel.background = element_rect(fill = '#f4f6f5'),
        panel.grid = element_blank())