Rotate by a given relative quarternion
I figured it out. The code is simply:
rot = normal.rotation_difference(target_vector)
obj.rotation_quaternion = rot

BUT I had transformed my normal vector the the global coordinate system like this:
world_matrix = obj.matrix_world 
normal = obj_face.normal           
normal = world_matrix.to_3x3() @ normal  

It worked if I skipped the last line with the matrix multiplication. I'm haven't wrapped my head around why this did the trick, as I just tried it out of intuition (and desperation). So I guess I have something to think about over the weekend :-)