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öberl, IMA Journal of numerical Analysis, 2019] that approximates the velocity in an H(div)-conforming space and introduces a new stress-like variable for the approximation of the gradient of the velocity within the function space H(curl 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 ψ and its approximation ψ h are elements of H(curl) and the H(curl)-conforming Nédélec finite element space, respectively. We conclude with an error analysis revealing optimal convergence rates for the error of the discrete velocity u h = curl(ψ 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.Stokes equations, Hellan-Herrmann-Johnson, Stream function formulation, incompressible flows.represent the vector valued and matrix-valued versions of). This notation is extended in an obvious fashion to other function spaces as needed.Depending on the type of the function, the gradient ∇ is to be understood from the context as an operator that results in either a vector whose components are [∇φ] i = ∂ i φ (where ∂ i is the partial derivative ∂/∂x i ) for φ ∈ C ∞ c ′ (Ω, R) or a matrix whose entries are [∇φ] ij = ∂ j φ i for φ ∈ C ∞ c ′ (Ω, R d ). Similarly, the "curl" is given as any of the following