Unrotate a square polygon with Python: all 3 axes I have a rotated square polygon. And I want to 'unrotate' it so that its normal axes would match world axes.

I'm able to unrotate it on the Z axis by using the face normal and rotation difference (I don't care about translation at this point):
quaternion = Vector((0,0,1)).rotation_difference(face_normal)
rotation = ( +quaternion ).to_matrix().inverted()
obj.matrix_world = rotation.to_4x4()


But for the love of god I can't figure out how to also add face  X and Y rotation to the mix.
Desired result:
Use the edges.
The bmesh module has some methods for faces to calculate the tangent vector of the face based on its edges. For example sake will use the tangent vector created from two longest non connected edges.
Using the face normal, the face tangent and the cross product of the two can create a rotation matrix from the three orthogonal vectors.
The axes are the columns of a matrix, so a rotation matrix with unit vector axes x, y, and z could be
Matrix((x, y, z)).transposed()

if you recall from linear algebra, the transpose of an orthogonal matrix is also its inverse.
Test script, aligns face zero based on the matrix created as above.
import bpy
import bmesh
from mathutils import Matrix, Quaternion

ob = bpy.context.object
me = ob.data
bm = bmesh.from_edit_mesh(me)
face = bm.faces[0]

 
n = face.normal
t = face.calc_tangent_edge_pair().normalized()
c = face.calc_center_median()

M = (
    Matrix.Translation(c) @ 
    Matrix((t.cross(n).normalized(), t, n)).to_4x4() @
    Matrix.Translation(-c)
    )

bmesh.ops.transform(
    bm,
    matrix=M,
    verts=bm.verts,
    )
bmesh.update_edit_mesh(me)