Gradient flow step-scaling function for SU(3) with twelve flavors

We calculate the step scaling function, the lattice analog of the renormalization group $\beta$-function, for an SU(3) gauge theory with twelve flavors. The gauge coupling of this system runs very slowly, which is reflected in a small step scaling function, making numerical simulations particularly challenging. We present a detailed analysis including the study of systematic effects of our extensive data set generated with twelve dynamical flavors using the Symanzik gauge action and three times stout smeared M\"obius domain wall fermions. Using up to $32^4$ volumes, we calculate renormalized couplings for different gradient flow schemes and determine the step-scaling $\beta$ function for a scale change $s=2$ on up to five different lattice volume pairs. Our preferred analysis is fully $O(a^2)$ Symanzik improved and uses Zeuthen flow combined with the Symanzik operator. We find an infrared fixed point within the range $5.2 \le g_c^2 \le 6.4$ in the $c=0.250$ finite volume gradient flow scheme. We account for systematic effects by calculating the step-scaling function based on alternative flows (Wilson or Symanzik) as well as operators (Wilson plaquette, clover) and also explore the effects of the perturbative tree-level improvement.