How to deal with non-proportional risk? Let's start with a toy example. We have an RCT with 2 treatment arms (Treatment A vs. Treatment B), and one outcome.
When we analyse survival curves according to the treatment arm, we found the following scenario:

It seems that the proportionality-hazard assumption is violated, and we therefore cannot use a Cox Regression Analysis to analyse the association between treatment arm and risk of outcome in this scenario.
How can we deal with this situation?
Is there any different model than Cox Regression analysis which are better suited to deal with such a situation?
 A: Your example shows that Treatment A has better survival until about Time = 3500, when Treatment B starts to become superior. In this case there is no single answer to which Treatment is superior: it depends on the time at which you evaluate outcomes.
A simple accelerated failure time (AFT) model, suggested in another answer, wouldn't help here. A standard AFT model assumes that a change in a covariate (here, treatment) leads to a compression or dilation of the time scale. Under such an AFT model these survival curves thus also wouldn't cross.
You need to decide what aspect of the survival curve is most important to you.  Example 7.2 of Klein and Moeschberger, 2nd edition covers an example like this, and shows how modifications of the log-rank test can be used to differentially weight early versus late events.
A standard Cox model would return a type of event-weighted average log-hazard, probably weighted toward early events when there is censoring. The methods in the R coxphw package might provide a less-biased estimate of the overall hazard ratio, if that is of interest.
A parametric model that goes beyond a simple AFT setup could be used. A standard form for an AFT model is:
$$\log T = \beta' X + \sigma W $$
where $T$ is event time, $X$ is the set of covariates with associated coefficients $\beta$, $W$ is a defined probability distribution (e.g., standard normal for a log-normal model), and $\sigma$ is a fixed scale factor. You can allow $\sigma$ also to be a function of the covariates, for example via the tools in the R flexsurv package, which could model the crossing curves. But then it's no longer an AFT model, just a parametric description of the results, and that doesn't address the fundamental question of which survival times are most important.
A: I'm wondering if an accelerated-life model would be more appropriate here.
We would take the log life length to be a linear function of the explanatory variables,
$$ \mathrm{log}(Y_i) = X_i \beta + \varepsilon_i$$
or, equivalently,
$$ Y_i = \exp(X_i \beta) U_i\,
$$
where $U_i=\exp(\varepsilon_i)$. The conditional survival function for $Y_i$ given explanatory variables $X_i=x_i$ is then
$$
S_{Y_i|X_i}(y|x_i)=S_U(y/\exp(x_i \beta))\,.
$$
In other words, the survival function has a fixed form, but its argument is scaled by the values of the explanatory variables. For example, if $x_i \beta=0$ and $x_j \beta=\log2$, the effect is that individual $j$ ages twice as fast as individual $i$.
The downside is that this is usually a fully parametric model, so we need to choose a distribution for $U$ (typically log-logistic or Weibull).
