Alright, I found the solution with some help of the answer of batFINGER here.
There's surprisingly little information on the internet about manipulating face normals at the vertex positions, confusingly called 'split normals' or 'loop normals'.
To work with split normals with python you first need to use the method calc_normals_split() on the mesh or else you won't get accurate vector values for the split normals or loop normals.
It took me hours finding that out.
In stead of setting x scale to -1 and applying transforms at the object level, you should edit the vertices, face normals and split normals at the mesh level.
code:
import bpy
get active object and mesh
ob = bpy.context.object
me = ob.data
Very important to include line below for working with custom normals!
me.calc_normals_split()
me.use_auto_smooth = True
go through each polygon, reverse second to last loop,
multiply x vector of loops with -1, store mirrored vector in
mirror_split_normals for later use.
mirror_split_normals = []
for iter_poly in me.polygons:
reverse_normals_indices = [iter_poly.loop_indices[0]] + [i for i in reversed(iter_poly.loop_indices[1:])]
for i in reverse_normals_indices:
mirror_split_normals.append((me.loops[i].normal[0] * -1, me.loops[i].normal[1], me.loops[i].normal[2]))
mirror vertices in x by multiplying x by -1
for v in me.vertices:
v.co = v.co[0] * -1, v.co[1], v.co[2]
flip face normals and apply mirrored split normals
me.flip_normals()
#me.update()
me.normals_split_custom_set(mirror_split_normals)
Original mesh:
Mirrored mesh:

abnormal-addon. You may find some help in its code base or its faeatures – Reigen Feb 22 '21 at 06:33