Survival analysis: using time varying covariates for resolving violated PH p1 is a factor
p2 is a continuous variable
Model
m = Survival_function ~ p1 + p2

Both covariates - p1 and p2 - violate the proportional hazards assumption. I am interested in showing p1's effect on survival, p2's PH violation should just be "somehow managed" to say that the showable result was adjusted to p2.
rstanarm survival package for R offers modelling time-varying effects of covariates. Example of time-varying effect is given below.

Question
This package does not allow testing for PH assumption (there is no specific test for this purpose but plotting the effects of simple time-varying models can be a possible alternative?). If I specify the covariates in the model as "time-varying" ("tve()"), does it automatically "resolve" the PH violation (the results of the model are OK for publishing)?
m = Survival_function ~ tve(p1) + tve(p2)

My own understanding is that if I show a plot like above for p1, then its related PH violation would be resolved as I show a time-varying effect. Is the PH of p2 similarly "resolved"?
Or am I wrong here? The time-varying effect is similar to using restricted cubic splines in classical statistics Cox models that sometimes resolves PH and sometimes not.
 A: What you're proposing is to use time-varying effects or cofficients of covariates. To avoid confusion, restrict the phrase "time-varying covariates" to situations in which the values of covariates change over time. See the R time-varying vignette for the distinction.
Before you proceed to time-varying coefficients, try to improve the model to reduce the violation of proportional hazards (PH). Your model presently is restricted to a linear association of p2 with log hazard, and only includes two predictors. A more flexible fit of p2 and incorporation of other outcome-associated predictors might get around the PH violation. Alternatively, your data might better be handled with an accelerated-failure-time model that doesn't require PH (but comes with its own assumptions).
If you do end up needing to use time-varying coefficients in your model, there doesn't seem to be any easy computational test for PH thereafter. Although the R survival package used to report results for a cox.zph() test of PH on a model including the tt() term used to handle time-varying coefficients, that calculation was incorrect and is no longer performed starting with version 3.1-2 of the survival package. In the time-varying vignette linked above, the test for modeling PH properly thus is a visual comparison of shapes between the modeled time-dependence function and the smoothed estimates of variance-scaled Schoenfeld residuals over time from the model without the time-varying coefficient.
I'm not familiar with the rstanarm* package, but in a quick look at its manual I didn't see that the package provides those residuals for evaluation. The plot you show from the rstanarm documentation is only a comparison between the time-dependence of simulated data and what was returned by the function, so you can't really "show a plot like above for p1" for a model of empirical data. You might find a way to do a visual comparison like that in the survival package documentation if you can find a way to extract the correct residuals from the rstanarm model.
Nothing "automatically" assures that results will be "OK for publishing"; you need to be ready to document to the reader that your model is good enough for its purposes. If you choose to use a Bayesian survival model, make sure that you understand the principles very well and are prepared to defend the modeling to a skeptical reviewer. Advantages of Bayesian survival modeling have long been appreciated but still aren't what one typically sees in the literature, so your reviewers and readers might need more than the usual guidance to understand what you're doing.

*As I understand it, the Bayesian modeling used by rstanarm isn't strictly a Cox model. Unlike a standard Cox model, it seems to fit the baseline hazard directly as part of the model.
