  We develop a cut finite element method for the Bernoulli free boundary
problem. The free boundary, represented by an approximate signed distance
function on a fixed background mesh, is allowed to intersect elements in an
arbitrary fashion. This leads to so called cut elements in the vicinity of the
boundary. To obtain a stable method, stabilization terms is added in the
vicinity of the cut elements penalizing the gradient jumps across element
sides. The stabilization also ensures good conditioning of the resulting
discrete system. We develop a method for shape optimization based on moving the
distance function along a velocity field which is computed as the $H^1$ Riesz
representation of the shape derivative. We show that the velocity field is the
solution to an interface problem and we prove an a priori error estimate of
optimal order, given the limited regularity of the velocity field across the
interface, for the the velocity field in the $H^1$ norm. Finally, we present
illustrating numerical results.
