How to get to original T-Pose coordinates after bpy.ops.pose.armature_apply?

Question

very related to this question:

PoseBone local rotation values to global with axis changed

For the addon we're building we do the following steps:

1. Create Mesh and armature from third party system
2. Change the pose of the mesh using the standard Blender rigging
3. Send bone rotations in global space to a third party system
4. Third party system returns new geometry
5. We substitute the mesh with the new geometry
6. Set new armature rest pose as the current state
7. Iterate between 2 to 6 until the pose looks like we want
8. Export to file

The problem we find is that while iterating between 2 to 6 there are moments where the geometry coming back does not correspond to what we expect. As the third party system needs the rotations coming from the rest pose we need to accumulate the rotations to output what the system expects. Currently we do the following:

        for pose_bone in context.object.pose.bones:
        # We'll translate from bone coords to avatar coords
        # using the matrix_local matrix in the actual bone
        # (not the pose bone) which contains the rest pose matrix

        # From matrix_local we'll extract the rotation
        rot = pose_bone.bone.matrix_local.decompose()[1]
        rot = rot.to_matrix().to_3x3()

        # Apply the rotation to the direction axis angle
        axis_angle = pose_bone.rotation_quaternion.to_axis_angle()
        new_axis = rot @ axis_angle[0]
        new_axis.normalize()

        # Get the axis angle to matrix
        axis_to_matrix: Matrix = Matrix.Rotation(
            axis_angle[1], 3, new_axis
        )

        # Rotate that axis with the accumulated bone matrix
        bone_accumulated_rotations[pose_bone] = (
            axis_to_matrix @
            bone_accumulated_rotations[pose_bone]
        )

        accumulated_axis_angle = (
            bone_accumulated_rotations[pose_bone]
            .to_quaternion()
            .to_axis_angle()
        )

        # Get the actual axis angle multiplying the magnitude
        # with the unitary vector calculated earlier
        result = (
            accumulated_axis_angle[0]
            .normalized() *
            accumulated_axis_angle[1]
        )

        # Assign the rotations to the locations avatar-core expects
        # on the pose numpy array
        bone_idx: int = get_bone_index(skeleton_info, pose_bone)

        pose_np[bone_idx*3] = result[0]
        pose_np[bone_idx*3+1] = result[1]
        pose_np[bone_idx*3+2] = result[2]

    mesh: Mesh = avatar_obj.data
    mesh.vertices.foreach_set("co", verts.ravel())
    mesh.update()

    # Set rest pose to armature current pose
    bpy.ops.pose.armature_apply(selected=False)

Is there anything else we should take into account when accumulating those rotations?

Thanks!