  We develop two unfitted cut finite element methods for the Stokes equations
based on $\mathbf{H}^{\text{div}}$-conforming finite elements which exhibit
optimal convergence order for the velocity, pointwise divergence-free velocity
fields, and well-posed linear systems, independently of the position of the
boundary relative to the computational mesh. The first method is based on the
Brezzi-Douglas-Marini (BDM) elements and involves interior penalty terms to
enforce tangential continuity of the velocity at interior edges in the mesh.
The second method is a 3-field formulation involving the vorticity, velocity,
and pressure and uses the Raviart-Thomas (RT) space for the velocity. We
present mixed ghost penalty stabilization terms for both methods so that the
resulting discrete problems are stable and the divergence-free property of the
$\mathbf{H}^{\text{div}}$-conforming elements is preserved also on unfitted
meshes. In both methods boundary conditions are imposed weakly. We show that
imposing Dirichlet boundary conditions weakly introduces additional challenges;
1) The divergence-free property of the RT and the BDM finite elemens may be
lost depending on how the normal component of the velocity field at the
boundary is imposed. 2) Pressure robustness is affected by how well the
boundary conditions are satisfied and may not hold even if the
incompressibility condition holds pointwise. We study two approaches of weakly
imposing the normal component of the velocity at the boundary; we either use a
penalty parameter and Nitsche's method or a Lagrange multiplier method. We show
that appropriate conditions on the velocity space has to be imposed when
Nitsche's method is used. Pressure robustness can hold with both approaches by
reducing the error at the boundary but this impacts the condition numbers of
linear systems, independent of if the mesh is fitted or unfitted to the
boundary.
