Custom Split Normals
Elaborating on this answer and here
Test script. Run in object mode. Switch to edit mode to see results.
import bpy
context = bpy.context
ob = context.object
me = ob.data
visualize in edit mode
me.show_normal_vertex = me.show_normal_loop = True
me.use_auto_smooth = True
normal custom verts on each axis
me.normals_split_custom_set([(0, 0, 0) for l in me.loops])
end up with custom split normals in each of the corner directions. Blue are the vert normals
# all custom split normals pointing up.
normals = [(0, 0, 1) for v in me.vertices]
# make csn's all face up.
me.normals_split_custom_set_from_vertices(normals)
all custom split normals facing up (0, 0, 1)
I am not aware if custom split normals can be accessed from bmesh.