The hypothetical presence of a nonzero electric dipole moment of the neutron would imply a violation of charge-parity symmetry (CP) which is absent from the Standard Model. Toward the characterization of potential BSM sources of CP violation, we study effective interactions in which the high-energy BSM fields are integrated out, leaving a basis of higher-dimensional, local QCD operators. Since these are composite operators, any attempt at renormalization involves some sort of operator mixing, which at hadronic scales must be studied nonperturbatively on the lattice. Naturally, the operator mixing is parametrized solely by the lattice spacing, producing power and logarithmic divergences as the regulator is removed. We use the gradient flow to reparametrize the operator mixing in terms of the flow time, so that at fixed flow times, each continuum limit may be taken independently, taking the physical $t\rightarrow0^+$ limit only at the end. Additionally, since operators are multiplicatively renormalized at positive flow time, it is trivial to identify the physical renormalized operators once the mixing matrix is computed. To that end, the divergent operator mixing may be calculated in perturbation theory to arbitrary order, constraining the high-energy matching to nonperturbative results. Moreover, since the continuum limit is possible at each positive flow time (identified with the inverse-squared renormalization scales), we may renormalize our operator bases over a wide range of hadronic scales as well. We study the renormalization of the quark chromoelectric dipole moment (qCEDM), presenting both nonperturbative and next-to-leading-order perturbative results on the power-divergent mixing with the pseudoscalar density with some perturbative results on the leading logarithms.