Adding mortality rate (/1000 patients-years) + 95CI in a tbl_regression (gtsummary) from Cox regression mode?

Question

I used the following code:

library(survival)
library(gtsummary)
library(dplyr)

df <- data.frame(
  time = sample(1:10, 20, replace = TRUE),
  status = sample(0:1, 20, replace = TRUE),
  age = sample(40:80, 20, replace = TRUE),
  sex = sample(c("male", "female"), 20, replace = TRUE))

coxph(Surv(time, status) ~ age + sex, data = df) %>%
  tbl_regression(exponentiate = TRUE) %>%
  add_n(location = "level") %>% 
  add_nevent(location = "level") %>%
  modify_header(exposure ~ "**Total follow-up (Person years)**") 

To get this great {gtsummary} table:

How can I add another column with the mortality rate per 100 patients-years and the corresponding 95% CI ? Ie, with the formula: n_event/total followup*1000, as we can see in this example:

Thanks for your help!

Edit: my attemp below

coxph(Surv(time, status) ~ age + sex, data = df) %>%
      tbl_regression(exponentiate = TRUE) %>%
      add_n(location = "level") %>%
      add_nevent(location = "level") %>%
      bold_p() %>%
      modify_header(exposure ~ "**Total follow-up (years)**") %>%
      modify_table_body(~ .x %>%
                          mutate(
                            #add mortality rate
                            mortality_rate = case_when(
                              !is.na(stat_nevent) & !is.na(exposure) & exposure > 0 ~ 
                                round((stat_nevent / exposure) * 1000, 2),
                              TRUE ~ NA_real_),
                            #add 95%CI
                            ci_mortality_rate = case_when(
                              !is.na(stat_nevent) & !is.na(exposure) & exposure > 0 & stat_nevent > 0 ~ {
                                rate = (stat_nevent / exposure) * 1000
                                se = sqrt(stat_nevent) / exposure * 1000
                                lower = round(rate - 1.96 * se, 2)
                                upper = round(rate + 1.96 * se, 2)
                                sprintf("%.2f - %.2f", lower, upper)
                              },
                              TRUE ~ NA_character_))) %>%
      modify_header(mortality_rate ~ "**Mortality Rate (/1000 person-years)**",
                    ci_mortality_rate ~ "**95% CI of Mortality Rate**")

It give the mortality rate and the 95% CI computed with Poisson.

Anybody can confirm that this is correct?