How code a three-levels cox model for hierarchical data in R?

Question

I have data I would like to analyze using a Cox proportional hazards model. The dataset contains an outcome variable (e.g. alive/dead), a follow-up time variable, quantitative and qualitative covariables, and nested cluster variables : patients (with multiple patients clustered within physicians), physicians (with multiple physicians clustered within regions) and regions. Data look like this :

set.seed(666)
n <- 1000 
outcome <- sample(c(0,1), n, replace=TRUE)
follow_time <- round(rnorm(n, mean=2, sd=1), 2)*100 
covariate1 <- sample(c("A","B","C"), n, replace=TRUE)
covariate2 <- rnorm(n)
covariate3 <- sample(c("X","Y","Z"), n, replace=TRUE)
covariate4 <- rnorm(n)
patient_id <- c(1:1000)
physician_id <- rep(1:200, each = 5)
region_id <- rep(1:10, each = 100)

df <- data.frame(outcome, follow_time, covariate1, covariate2, covariate3, covariate4,
  patient_id, physician_id, region_id)

outcome follow_time covariate1 covariate2 covariate3  covariate4 patient_id
       1         299          B  1.1179509          Z  0.30611075          1
       1         218          B  0.2552784          Z  0.57186034          2
       1         259          C  0.2763426          Z -1.55426489          3
       0          23          B  0.2477036          Y -0.68620755          4
       0         205          B  0.7430005          X -2.00498162          5
       1         399          B -1.1510293          X  0.06120172          6
       1         204          A -1.0451297          Y  0.10244376          7
       0         100          C -0.6770400          Z -1.38118875          8
       0         263          C  0.2681764          Z  0.04140385          9
       1         130          A -0.7490197          Z  0.99470583         10
 physician_id region_id
            1         1
            1         1
            1         1
            1         1
            1         1
            2         1
            2         1
            2         1
            2         1
            2         1

The data has a multilevel structure due to clustering of patients within physicians and physicians within regions. I would like to model the effect of theses three levels. I have no problem for a two-levels model, but i can't find any explicite ressource explaining how code this.

Articles I've read briefly mention that function coxme and packages like and frailtypack allow fitting three-multilevel Cox models in R. However, all examples and documentation seem focused on two-level structures only.

I'm uncertain how to properly specify the random effects to account for the three hierarchical levels in my data. Could someone provide a simple code example or steps to fit a three-level Cox model in R? I'm particularly unsure about how to define the frailties across levels.