  We use the gradient flow for the renormalization of the Polyakov loop in
various representations. Using 2+1 flavor QCD with highly improved staggered
quarks and lattices with temporal extents of $N_\tau=6$, $8$, $10$ and $12$ we
calculate the renormalized Polyakov loop in many representations including
fundamental, sextet, adjoint, decuplet, 15-plet, 24-plet and 27-plet. This
approach allows for the calculations of the renormalized Polyakov loops over a
large temperature range from $T=116$ MeV up to $T=815$ MeV, with small errors
not only for the Polyakov loop in fundamental representation, but also for the
Polyakov loops in higher representations. We compare our results with standard
renormalization schemes and discuss the Casimir scaling of the Polyakov loops.
