There's probably a cleaner solution, but for the time being, you can use:

    import bpy
    from math import radians
    from mathutils import Vector, Matrix
    
    obj = bpy.context.active_object
    mesh = obj.data
    
    
    def get_layers():
        """
        return: tuple of 20 booleans. True if layer is active layer,False if not
        Allows the mesh to be created on the active layer. There is only 1 active layer
        """
        active_layer = bpy.context.scene.active_layer
        layers = [False] * 20
        layers[active_layer] = True
        return tuple(layers)
    
    
    def create_light_at_face():
        """
        layers: tuple. 20 booleans where exaclty one boolen is True
        return: None
        creates an area light the has the same orientation as the normal at the center of each face
        """
        mat_world = obj.matrix_world
        up = Vector((0,0,1))
        
        for poly in mesh.polygons:
            co = mat_world * Vector(poly.center)
            
            forward = poly.normal.copy()
            forward.rotate(mat_world)
            right = forward.cross(up).normalized() # Vector.length closer to 1.0
            up = right.cross(forward)
        
            rot = Matrix((right, up, -forward)).transposed().normalized().to_4x4()
            mat = Matrix.Translation(co) * rot
        
            bpy.ops.object.lamp_add(type='AREA', view_align=False, location=co, layers=get_layers())
            # scale lamp up a little fo we can see the area
            bpy.context.object.data.size = 3
            bpy.context.object.matrix_world = mat
    
    create_light_at_face()

Also tried eulers, but did not work reliably (some lamps pointing inwards).

As ideasman_42, suggested, `to_track_quat()` can be used to do the job, too:

            quat = poly.normal.to_track_quat('-Z', 'Y')
            loc = Matrix.Translation(poly.center)
            mat = mat_world * loc * quat.to_matrix().to_4x4()
            bpy.context.object.matrix_world = mat

It's not much code less, but more understandable I hope.

Another variant:

            loc = mat_world * Vector(poly.center)
            # use loc as co for lamp_add operator here
            quat = poly.normal.to_track_quat('-Z', 'Y')
            quat.rotate(mat_world)
            bpy.context.object.matrix_world *= quat.to_matrix().to_4x4()