  We derive and analytically solve renormalization group (RG) equations of
gauge invariant non-local Wilson line operators which resum logarithms for
event shape observables $\tau$ at subleading power in the $\tau\ll 1$
expansion. These equations involve a class of universal jet and soft functions
arising through operator mixing, which we call $\theta$-jet and $\theta$-soft
functions. An illustrative example involving these operators is introduced
which captures the generic features of subleading power resummation, allowing
us to derive the structure of the RG to all orders in $\alpha_s$, and provide
field theory definitions of all ingredients. As a simple application, we use
this to obtain an analytic leading logarithmic result for the subleading power
resummed thrust spectrum for $H\to gg$ in pure glue QCD. This resummation
determines the nature of the double logarithmic series at subleading power,
which we find is still governed by the cusp anomalous dimension. We check our
result by performing an analytic calculation up to ${\cal O}(\alpha_s^3)$.
Consistency of the subleading power RG relates subleading power anomalous
dimensions, constrains the form of the $\theta$-soft and $\theta$-jet
functions, and implies an exponentiation of higher order loop corrections in
the subleading power collinear limit. Our results provide a path for carrying
out systematic resummation at subleading power for collider observables.
