Reporting multiple p-values with stargazer

Question

I would like to export the p-values of coefficients from lm model. I also perform a randomisation inference test as I am interested in a treatment effect. I also have p-values on the same coefficients from the test. I would like to report the treatment coefficients with both p-values underneath it to enable comparison which to show validity of the randomisation procedure based on which was collected my real data.

The problem is that the p-values can be different and point out to different significance level or no significance in either case. Ideally I would like to have a star character for each type of p-value. Stargazer can report starts based on p-values of the coefficients, but I am not able to have stars on both types of p-values, for now the stars are the same for both p-values even when their values are different.

This is my code to illustrate the problem I am having:

library(stargazer)
library(dplyr)
library(lmtest)
library(sandwich)
library(fastDummies)


set.seed(123)

ngroup = 4

nrep = 15

treat <- rep(c("control", "treat1", "treat2", "treat3"), each = nrep)

Y <- rnorm(n = 60, mean = 50, sd = 10)
                 
data <- cbind(Y, treat)

data <- as.data.frame(data)

data <- data %>%
mutate(Y = as.numeric (Y))

data <- fastDummies::dummy_cols(data, select_columns = "treat")

### lm model
m1 <- lm(Y ~ treat_treat1 + treat_treat2 + treat_treat3, data = data)

### extract the p-values from m1
p_m1 <- coeftest(m1)[,4]

## p-vlues generated through ritest (just as a example)

p_values <- c(0.1, 0.001, 0.05, 0.1)

stargazer(m1, se = list(p_m1), p = list(p_values), report=('vcs*p*'), omit.stat = c("adj.rsq", "f"), out = "text")

Any ideas how I can report stars for both types of p-values in this case? Maybe there could be an alternative package to stargazer that can also do this and I am just not aware.