We propose a new method to determine quark masses using ratios of the vacuum-expectation values (vevs) of flowed quark bilinear operators. They can be expressed as functions of the flow time t and the \( \overline{\textrm{MS}} \) quark mass \( \overline{m} \) , which can then be determined by matching with the corresponding lattice results. Motivated by this, we evaluate these vevs perturbatively through next-to-leading order in the strong coupling. We provide the results as expansions in the limits of small and large \( {\overline{m}}^2t \) . To this end, we develop a new expansion technique based on the Laplace transform. Additionally, we present numerical results with the exact mass dependence over a wide range of \( {\overline{m}}^2t \) . We discuss the expected perturbative precision for the mass determination based on our next-to-leading order perturbative calculations, and possible non-perturbative corrections.