What is the most efficient way to finely "slice" layers though objects with python
A more straight-forward approach.
Not sure this is &quot;most efficient&quot; but it's relatively quick and is a straight forward load mesh into bmesh &amp; single bisect cut per rep
Now it's clearer the result desired, re just the edges of the bisect cuts have re-jigged answers from
Python: Bisect mesh into n parts, without separate &#39;LOOSE&#39;
How to extract side-view outline (e.g. top view) of a 3D object to 2D surface?
to create objects using the result of the bisect operator.
Notes.

Suprisingly it was quicker to constantly reload the bmesh rather than use recycling the same

For example sake have simply run from (0, 0, -1) to (0, 0, 1) local coordinates to test on sphere / icosphere / cube etc.

See Replace matrix @ vector list comprehensions with something more efficient  re a fast way to convert between spaces and calculate the required limits.

To convert the bmesh to global (or another space) use bm.transform(matrix)

To use a mesh with modifiers and shapes applied use bm.from_object(ob, depsgraph)

Other commented out code re attempts to squeeze out a bit more speed, have AFAIK left those that ran the quickest.


Test script.
import bpy
import bmesh
from mathutils import Vector
from bpy import context 


def slicer(ob, start, end, cuts):
    #slices = [] could instead return unlinked objects
    axis = end - start
    dv = axis / cuts
    #axis.normalize()
    mw = ob.matrix_world
    bm0 = bm = bmesh.new()
    
    # transform to world coords
    # bm.transform(mw)
    # mesh = bpy.data.meshes.new(&quot;Slice&quot;)

    for i in range(cuts + 1):
        bm.from_mesh(ob.data)
        #bm.from_object(ob, dg) # use modified mesh
        #bm = bm0.copy()
        #bm.transform(mw) # make global
        plane_co = start + i * dv
       
        cut = bmesh.ops.bisect_plane(
                bm,
                plane_co=plane_co,
                plane_no=axis,
                geom=bm.verts[:] + bm.faces[:] + bm.edges[:],
                clear_inner=True,
                clear_outer=True,
                )[&quot;geom_cut&quot;]
        if not cut:
            bm.clear()
            continue
        me = bpy.data.meshes.new(f&quot;Slice{i}&quot;)
        # me = mesh.copy()
        bm.to_mesh(me)
        '''
        # only slightly slower
        me.from_pydata(
            [v.co for v in cut if hasattr(v, &quot;co&quot;)],
            [(e.verts[0].index, e.verts[1].index) for e in cut if hasattr(e, &quot;verts&quot;)],
            []
            )
        '''

        slice = bpy.data.objects.new(f&quot;Slice{i}&quot;, me)
        slice.matrix_world = mw
        context.collection.objects.link(slice)
        bm.clear()        
        #bm.free() # if using bm = bm0.copy()
    #bm.free()
    
slicer(
        context.object, 
        Vector((0, 0, -1)), 
        Vector((0, 0, 1)), 
        200)

Timing it.
Using the decorator method
REPS = 10
import time
def time_it(func):
    from builtins import print
    def wrapper(*arg, **kw):
        t1 = time.time()
        for i in range(REPS):
            ret = func(*arg, **kw)
        t2 = time.time()
        print(func.__name__, (t2 - t1))
        return ret
    return wrapper
@time_it
def slicer(....

10 Test runs of 10 reps making 200 cuts in default UV sphere.  Using the command for your script of
slicer(0.01, (0,0,1)) # 

which should be 200 cuts on default sphere at origin.
To time the two scripts, after removing prints from yours  to make things relatively equal (as mentioned didn't calculate the bounds of mesh, however as shown in link provided using numpy or other methods, it's hardly rate determining) consistently came up with a speed ratio, this answer to yours of around 8 to 11.
However an issue was yours is making 2500 objects, whereas mine makes 2000 over 10 runs, which will also account for some time differential when fixed.
Not seeing the glitch.

Unfortunately there seems to be a glitch where sometimes layers are
lost or only partially formed - I think it's a bug in the bisect
function, not my code - eg bottom photo: the icosphere is missing the
middle layer if the slice would be at exactly the midline of the
icosphere .

No issue with code above run on icosphere,

Left using 2 cuts (prob should be segments as it's 1 middle cut) and on right with 200 and middle edge ring selected