The three-loop anomalous dimension and the four-loop $\beta$-function for ${\cal N}=1$ SQED regularized by higher derivatives

For ${\cal N}=1$ SQED with $N_f$ flavors regularized by higher derivatives in the general $\xi$-gauge we calculate the three-loop anomalous dimension of the matter superfields defined in terms of the bare coupling constant and demonstrate its gauge independence. After this the four-loop $\beta$-function defined in terms of the bare coupling constant is obtained with the help of the NSVZ equation, which is valid for these renormalization group functions in all loops. Next, we calculate the three-loop anomalous dimension and the four-loop $\beta$-function defined in terms of the renormalized coupling constant for an arbitrary subtraction scheme supplementing the higher derivative regularization. Then we construct a renormalization prescription for which the results coincide with the ones in the $\overline{\mbox{DR}}$-scheme and describe all NSVZ schemes in the considered approximation. Also we demonstrate the existence of a subtraction scheme in which the anomalous dimension does not depend on $N_f$, while the $\beta$-function contains only terms of the first order in $N_f$. This scheme is obtained with the help of a finite renormalization compatible with a structure of quantum corrections and is NSVZ. The existence of such an NSVZ scheme is also proved in all loops.