Use stat_function with facets

Prep Steps

library(tidyverse)
penguins <- palmerpenguins::penguins |> 
  filter(!is.na(sex))

penguin_colors <- c("#E69F00", "#009E73", "#0072B2")
names(penguin_colors) <- unique(penguins$species) 

params_penguins <- penguins |> 
  summarize(
    mean_body_mass_g = mean(body_mass_g),
    sd_body_mass_g   = sd(body_mass_g),
    .by              = species
  )

facet_plot <- penguins |> 
  ggplot() +
  geom_density(
    aes(x = body_mass_g, fill = species),
    col = NA
  ) +
  facet_wrap(vars(species), ncol = 1) +
  labs(
    x = 'Body weight (in g)',
    y = element_blank(),
    title = 'Penguin weights are somewhat normally distributed'
  ) +
  scale_fill_manual(values = penguin_colors) +
  scale_y_continuous(
    breaks = c(0, 6e-4),
    labels = scales::label_comma(),
    expand = expansion(mult = c(0, 0.005))
  ) +
  theme(legend.position = 'none')

Add lines to each facet

range_body_mass <- range(penguins$body_mass_g)
params_penguins
## # A tibble: 3 × 3
##   species   mean_body_mass_g sd_body_mass_g
##   <fct>                <dbl>          <dbl>
## 1 Adelie               3706.           459.
## 2 Gentoo               5092.           501.
## 3 Chinstrap            3733.           384.

computed_vals <- params_penguins |> 
  mutate(
    # wrap vector in list() to save whole vector into the cell
    x = list( 
      seq(range_body_mass[1], range_body_mass[2], 1)
    ),
    y = pmap(
      list(x = x, mean = mean_body_mass_g, sd = sd_body_mass_g),
      ~{
        l <- list(...)
        dnorm(l$x, mean = l$mean, sd = l$sd)
      }
    )
  ) |> 
  select(species, x, y)
computed_vals
## # A tibble: 3 × 3
##   species   x             y            
##   <fct>     <list>        <list>       
## 1 Adelie    <dbl [3,601]> <dbl [3,601]>
## 2 Gentoo    <dbl [3,601]> <dbl [3,601]>
## 3 Chinstrap <dbl [3,601]> <dbl [3,601]>

coords_penguins <- computed_vals |> 
  unnest(cols = c(x, y))
coords_penguins |> print(n = 3)
## # A tibble: 10,803 × 3
##   species     x         y
##   <fct>   <dbl>     <dbl>
## 1 Adelie   2700 0.0000784
## 2 Adelie   2701 0.0000788
## 3 Adelie   2702 0.0000791
## # ℹ 10,800 more rows

facet_plot +
   geom_line(
    data = coords_penguins,
    aes(x = x, y = y),
    linetype = 2,
    linewidth = 1.25
  )