  We introduce a new discretization for the stream function formulation of the
incompressible Stokes equations in two and three space dimensions. The method
is strongly related to the Hellan-Herrmann-Johnson method and is based on the
recently discovered mass conserving mixed stress formulation [J.
Gopalakrishnan, P.L. Lederer, J. Sch\"oberl, IMA Journal of numerical Analysis,
2019] that approximates the velocity in an $H(\operatorname{div})$-conforming
space and introduces a new stress-like variable for the approximation of the
gradient of the velocity within the function space
$H(\operatorname{curl}\operatorname{div})$. The properties of the (discrete) de
Rham complex allows to extend this method to a stream function formulation in
two and three space dimensions.
  We present a detailed stability analysis in the continuous and the discrete
setting where the stream function $\psi$ and its approximation $\psi_h$ are
elements of $H(\operatorname{curl})$ and the
$H(\operatorname{curl})$-conforming N\'ed\'elec finite element space,
respectively. We conclude with an error analysis revealing optimal convergence
rates for the error of the discrete velocity $u_h =
\operatorname{curl}(\psi_h)$ measured in a discrete $H^1$-norm. We present
numerical examples to validate our findings and discuss structure-preserving
properties such as pressure-robustness.
