  We derive the two-loop QCD helicity amplitudes for the processes $l q \to l
qg$ ($l \bar q \to l \bar q g$) and $l g \to l q\bar q$, which are the partonic
reactions yielding $(2+1)$-jet final states in deep inelastic lepton nucleon
scattering. The amplitudes are obtained by analytic continuation of the known
helicity amplitudes for $e^+e^- \to q\bar q g$. We separate the infrared
divergent and finite parts of the amplitudes using Catani's infrared
factorization formula. The analytic results for the finite parts of the
amplitudes are expressed in terms of one- and two-dimensional harmonic
polylogarithms. To evaluate these functions numerically, we list in detail the
non-trivial (and kinematic region dependent) variable transformations one needs
to perform.
