We present preliminary results for the Renormalization Group (RG) running of the complete basis of $\Delta F=2$ four-fermion
operators in QCD with $N_f=3$ dynamical massless flavours. We use O(a)-improved Wilson fermions in a mixed action setup, with chirally rotated Schrödinger functional ($\chi$SF) boundary conditions for the valence quarks and Schrödinger functional (SF) boundary conditions for the sea quarks.
The RG evolution operators are evaluated non-perturbatively via the matrix step-scaling functions (matrix SSF) using a SF coupling from the perturbative region down to $\sim4$GeV and a Gradient Flow (GF) coupling from $\sim4$GeV down to $\sim 250$MeV. The perturbative running is computed through a novel approach that extends the usual computations in the literature relying on consequences of the Poincaré-Dulac theorem.
