  We perform a global fit of the available polarized Semi-Inclusive Deep
Inelastic Scattering (SIDIS), polarized pion-induced Drell-Yan (DY) and
$W^\pm/Z$ boson production data at N$^3$LO and NNLO accuracy of the Transverse
Momentum Dependent (TMD) evolution, and extract the Sivers function for $u$,
$d$, $s$ and for sea quarks. The Qiu-Sterman function is determined in a model
independent way via the operator product expansion from the extracted Sivers
function. The analysis is supplemented by additional studies, such as the
estimation of applicability region, the impact of the unpolarized
distributions' uncertainties, the universality of the Sivers functions,
positivity constraints, the significance of the sign-change relation, and the
comparison with the existing extractions
