There is quite some literature on survival analysis under population heterogeneity, but like you notice yourself, I rarely see such models being used - or even considered - in applied research. I'll give some brief intuition, and hopefully others can add mathematically-heavier explanations if that's what you're looking for.
So let's say our observations look like this - quarter of the sample is "cured" and does not experience the event during observation period, while the rest are "susceptible" and slowly dying/failing/divorcing:

The entire blue fraction will be censored at $t=50$, and you can have additional censoring in both groups as usual (not seen in the plot).
Now, Cox PH models assume that survival time is independent of censoring. However, if you fit such a model to the pictured sample, the censored group will have much higher time-to-event a different hazard function because it contains the entire "cured" fraction. Simulations (e.g. 2) show that the HR estimate can be strongly biased under these conditions.
The problem could be solved by excluding the "cured" observations, but those usually are not distinguishable from censored, but "susceptible" individuals. Another analysis option is to only consider the binary outcome event vs. no event, but this again means excluding all censored individuals.
A conceptually simple solution is to treat the sample as a mixture of "cured" and "susceptible" distributions with weights $p$ and $1-p$, and fit a mixture model with survival $S(t) = 1$ for the "cured" fraction and some decreasing $S(t)$ for the "susceptible" fraction. These are known as mixture cure models; there are also non-mixture solutions, but I'm unsure if those are popular in practice.
Some references, in the order from most-layman to most mathematical:
Cure Models as a Useful Statistical Tool for Analyzing Survival
Mixture and non-mixture cure fraction models based on the generalized modified Weibull distribution with an application to gastric cancer data
What Cure Models Can Teach us About Genome-Wide Survival Analysis
A General Approach for Cure Models in Survival Analysis
EDIT
Given the good comments by @JarleTufto below, I should clarify this answer. I do not propose that merely censoring at the end of observation period is bad, or that the hazard after final $t$ is somehow relevant, but that problems are caused by underlying heterogeneity in population - i.e. fractions with different hazard functions $h(t)$ are censored differently. The corresponding Cox model assumption is stated as:
$$h(t|covariates) = h(t|C_i > t, covariates)$$
(e.g. Relaxing the independent censoring assumption in the Cox proportional hazards model using multiple imputation)
Let's take the simulation code from the answer below, and add a "cured" fraction with $h(t)=0$ (also added seed and increased sample size for easier reproducibility):
set.seed(1234)
# simulated survival times from the model
n <- 5000
x <- rnorm(n)
beta <- 0.5
# variation of inversion method
u <- runif(n)
eventtime <- 1/(1 + exp(-beta*x)*log(u)) - 1
eventtime[eventtime < 0] <- Inf
# simulate independent right censoring points
censoringtime <- runif(n,0,20)
# compute the observed data, that is, the censoring indicator and
# the time of whichever event comes first
delta <- eventtime < censoringtime
time <- pmin(eventtime, censoringtime)
# add "cured" fraction, censored at max observation time:
n2 <- 1000
time2 <- c(time, rep(20, n2))
delta2 <- c(delta, rep(FALSE, n2))
x2 <- c(x, rnorm(n2))
So the observed (event or censoring) timepoints look like:
qplot(time2)

Using only the "susceptible" fraction, we get the expected estimate of $\beta=0.5$:
# Fit the cox proportional hazards model
library(survival)
model <- coxph(Surv(time,delta) ~ x)
model
coef exp(coef) se(coef) z p
x 0.5067 1.6598 0.0198 25.6 <2e-16
Likelihood ratio test=661 on 1 df, p=0
n= 5000, number of events= 2923
However, using the full sample, the HR is underestimated:
model2 <- coxph(Surv(time2,delta2) ~ x2)
model2
coef exp(coef) se(coef) z p
x2 0.4192 1.5207 0.0192 21.8 <2e-16
Likelihood ratio test=478 on 1 df, p=0
n= 6000, number of events= 2923