  We proposed a recipe to systematically calculate Feynman integrals containing
linear propagators using the auxiliary mass flow method. The key of the recipe
is to introduce a quadratic term for each linear propagator and then using
differential equations to get rid of their effects. As an application, we
calculated all master integrals of vacuum integrals containing a gauge link up
to four loops, and we checked the results by nontrivial dimensional recurrence
relations.
