Unrotate a square polygon with Python: all 3 axes
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)