Perturbative gradient flow coupling of the twisted Eguchi-Kawai model with the numerical stochastic perturbation theory

The gradient flow scheme has emerged as a prominent nonperturbative renormalization scheme on the lattice, where flow time is introduced to define the renormalization scale. In this study we perturbatively compute the gradient flow coupling for the SU($N$) Yang-Mills theory in the large-$N$ limit in terms of the lattice bare coupling up to three-loop order. This is achieved by combining the twisted Eguchi-Kawai model with the numerical stochastic perturbation theory. We analyze the flow time dependence of the perturbative coefficients to determine the perturbative beta function coefficients, successfully computing the one-loop coefficient in the large-$N$ limit using three matrix sizes $N=289,441,529$. However, the higher-order coefficients are affected by large statistical errors. We also explore the potential for reducing these statistical errors through variance reduction combined with the large-$N$ factorization property of the SU($N$) Yang-Mills theory, and estimate the required number of samples for the precise determination of the higher-order coefficients.