The polarization of the Λ particle offers the unique opportunity to study the hydrodynamic gradients in the Quark-Gluon Plasma formed in heavyion collisions. However, the theoretical formula commonly used to calculate polarization is only a linear order expansion in thermal vorticity and neglects higher-order corrections. Here, I present an exact calculation to all orders in (constant) thermal vorticity at global equilibrium, obtaining the analytic form of the spin density matrix and the polarization vector for massive particles of any spin. Finally, I extend these results to local equilibrium and assess their phenomenological impact by numerically calculating the polarization vector in a 3+1 hydrodynamic simulation.