With such a large data set, it's possible that the "statistically significant" deviation from proportional hazards (PH) has no practical importance. This is similar to the problems with testing for normality in linear regression in a large data set. With enough data points, even a trivial deviation from the assumption might be "statistically significant."
The first thing to do is to examine the plot of smoothed, scaled Schoenfeld residuals over time. I'd recommend using the original plot functions from the standard R survival package, as at least prior versions of the survminer package had a coding error that made its plots pretty much useless. Therneau and Grambsch explain what to look for in Section 6.5, as explained on this page.
Fundamentally, you want to see how large the variation in the estimated Cox regression coefficient over time, $\hat\beta(t)$ from the smoothed plot, is with respect to the absolute value of the single value, $\hat\beta$, returned by the Cox model under the PH assumption. If the variation is relatively small (admittedly a judgment call), then you can assume that PH holds well enough, explain your reasoning in your report, and show the plot in supplemental information so that reviewers and readers can judge for themselves.
If the variation in $\hat\beta(t)$ is fairly large, particularly if the variation has a systematic trend with time, then either the model is mis-specified in terms of associations of continuous predictors with outcome or the PH assumption doesn't hold very well. If the PH assumption doesn't hold, then there is no single hazard ratio (HR) that can describe the predictor's association with outcome over all time periods. Having a change in HR over time might be clinically meaningful; apply your understanding of the subject matter in that case.
Even if PH doesn't hold, the single coefficient estimate $\hat\beta$ represents a sort of event-weighted average value that might be informative. AdamO shows the rationale on this page, with reference to the literature. Again, whether to accept that value is a judgment call that must be based on your understanding of the subject matter.
Alternatively, survival data that violate PH can sometimes be fit well with a parametric accelerated failure time (AFT) model. You should examine AFT models other than a Weibull model, as a Weibull model also effectively assumes PH. Log-normal or log-logistic models are examples worth examining.
AFT models can be a good choice, as some find AFT models to have a simpler intuitive interpretation than PH models. In AFT models, predictors effectively speed up or slow down the time scale. That's pretty easy to think about. The concept of hazard that underlies PH models can be trickier to grasp and often confuses those who are new to survival analysis. See this recent question, for example.
In response to new graph showing smoothed residuals
There is an extremely large drop-off in the time-varying coefficient up to about 10 years. It's far greater than the associated confidence intervals. The initial coefficient value of ~1.4 is equivalent to a hazard ratio of 4; the value of ~0.4 is equivalent to a hazard ratio of 1.5. It's reasonably constant thereafter.
The violation of PH for this variable is thus fairly substantial, not just representing your very large number of observations. There will be no single hazard ratio that applies over all times. The single coefficient reported by coxph() will provide an event-averaged estimate, if you wish some single value, but (in any event) follow the recommendation for robust error estimates made by @AdamO in a comment on your question.
I would recommend reporting your Cox model and then showing this graph to show how the estimated Cox regression coefficient (and thus the hazard ratio) changes over time. That's a direct way to show what you have found.
