We present results from an ongoing project concerned with the computation of $\mathcal{O}(1/m_Q)$ and $\mathcal{O}(1/m_Q^2)$ relativistic corrections to the static potential. These corrections are extracted from Wilson loops with two chromo-field insertions. We use gradient flow, which allows to renormalize the inserted fields and leads to a significantly improved signal-to-noise ratio, providing access to loops with large spatial and temporal extents.