How can an empty be added at intersection of edges using Python?
Getting the location is indeed trivial, even more with mathutils.geometry.intersect_line_line().

Getting rotation here is not so hard either. Let's call your two edges e1 and e2 (note that this is theoretical code, assuming you have your vertices array and empty object already at hand):

# We get 'average' intersection point - in case the two edges do not actually exactly intersect.
inter = mathutils.geometry.intersect_line_line(vertices[e1.v1].co, vertices[e1.v2].co, vertices[e2.v1].co, vertices[e2.v2].co)
vi = (inter[0] + inter[1]) / 2

# Then, we get our first axis.
ax = (vertices[e1.v1].co - vi).normalized()

# We get our third axis by cross product of first axis and second edge.
az = (ax.cross(vertices[e2.v1].co - vi)).normalized()

# Then second axis is crossproduct of the first and third ones.
# Note that if you know for sure your two edges are orthogonal,
# You can skip that and rather create ay exactly as we did with ax.
ay = (az.cross(ax)).normalized()

# Those three axes allow us to create a rotation matrix directly.
# However, we need to put them as columns, not rows. Since matrix
# creation in Blender uses rows, we'll have to transpose it.
mrot = mathutils.Matrix((ax, ay, az))
mrot.transpose()

# And now, we just have to make this a 4d matrix and add translation component.
mat = mrot.to_4x4()
mat.translation = vi

# And set this a base matrix of our empty:
empty.matrix_basis = mat