comparison of two different mediators for the models with the same DV and IVs

Question

Based on my theoretical arguments, I have two competing mediators (M1 and M2), each of which possibly mediates the relationship between my independent (X) and dependent (Y) variable. I implemented two separate mediation analyses using the lavaan package in R with bootstrapping technique for the same dataset: one with the mediator M1 (X>M1>Y), and one with the mediator M2 (X>M2>Y). The coefficient estimate for the mediating effect (i.e. the indirect effect a1*b1) of M2 (β=0.199) is higher than that of M1 (β=0.170). How can I find out if the mediating effects of these two mediators differ in a statistically significant way? Through which method or procedure can I assess this?

I have tried to find some information on comparing the coefficients of different mediators for the same statistical model, but what I could get was related only to comparing the coefficients of different independent variables for the same or different models but not mediators.Although a very similar question was asked some time ago in this platform, surprisingly it has not been answered.

As can be seen from the R codes, my model is actually a moderated mediation model, but at the moment I am only interested in the formal comparison of the mediation effects, not the effects of moderated mediation.

I would very much appreciate if someone could help me resolve this issue. Thanks a lot in advance.

Here's my code:

# For Mediator1 (M1)

ModMediation1 <- ' M1 ~ a1*RD + a2*PD + a3*RDXPD + a4*FA +
              a5*FS + a6*FG + a7*T
              
              Perf ~ b1*M1 + b2*FA + b3*FS + b4*FG + b5*T +
              c2*PD + c3*RDXPD
              
              # indirect effect
              IndEff := a1*b1
              
              # index of moderated mediation 
              IndModMed := a3*b1
              
            '

ModMediation1_fit <- lavaan::sem(ModMediation1, data = p2_df, se = "bootstrap", bootstrap = 10000)

summary(ModMediation1_fit, fit.measures = TRUE, rsq = TRUE, standardized = TRUE, ci = TRUE)

parameterestimates(ModMediation1_fit, boot.ci.type = "bca.simple", standardized = TRUE, level = 0.95)

# For Mediator2 (M2)

ModMediation2 <- ' M2 ~ a1*RD + a2*PD + a3*RDXPD + a4*FA +
              a5*FS + a6*FG + a7*T
              
              Perf ~ b1*M2 + b2*FA + b3*FS + b4*FG + b5*T +
              c2*PD + c3*RDXPD
              
              # indirect effect
              IndEff := a1*b1
              
              # index of moderated mediation 
              IndModMed := a3*b1
              
            '


ModMediation2_fit <- lavaan::sem(ModMediation2, data = p2_df, se = "bootstrap", bootstrap = 10000)

summary(ModMediation2_fit, fit.measures = TRUE, rsq = TRUE, standardized = TRUE, ci = TRUE)

parameterestimates(ModMediation2_fit, boot.ci.type = "bca.simple", standardized = TRUE, level = 0.95)

PS: I posted this question some time ago in StackOverflow, which was apparently not a right platform.