Multiple rotations do not add up
My first thought was that the vectors from the origin of a regular unit tetrahedron are normal to the plane tangent to the sphere where the vector intercepts, and we must be able to recover the rotation from that.
We see that Wikipedia gives the rather ugly values of $( 0, 0, 1 )$, $(\sqrt{\frac{8}{9}},         0, -1/3 )$, $( -\sqrt{\frac{2}{9}}, -\sqrt{\frac{2}{3}}, -\frac{1}{3} )$, and $( -\sqrt{\frac{2}{9}}, \sqrt{\frac{2}{3}}, -\frac{1}{3} )$ for these vertices, so I wrote some Python to draw it.
D = bpy.data
from math import sqrt, radians
from mathutils import Vector

radius = 1

bpy.ops.mesh.primitive_solid_add(source='4', size=radius)

tetrahedron = C.active_object


empties = []

theta = radians(120)
vectors = [ Vector( (   0, 0, 0, 1 ) ),
            Vector( ( theta, sqrt(8/9),         0, -1/3 ) ),
            Vector( ( theta,-sqrt(2/9), -sqrt(2/3), -1/3 ) ),
            Vector( ( theta,-sqrt(2/9),  sqrt(2/3), -1/3 ) ),]
hack=[ 0, 2, 3, 1]


for n in range(4):
    bpy.ops.object.empty_add(type='SINGLE_ARROW', align='WORLD', location=(0, 0, 0))
    empties.append(C.active_object)

    empties[n].rotation_mode= 'AXIS_ANGLE'
    empties[n].rotation_axis_angle = vectors[n]
    empties[n].location = radius * Vector(vectors[hack[n]][1:4])


for n in range(4):
    bpy.ops.mesh.primitive_monkey_add(size=0.5, enter_editmode=False, align='WORLD', 
        location=radius * Vector(vectors[hack[n]][1:4]) )
    monkey = C.active_object
    monkey.rotation_mode= 'AXIS_ANGLE'
    monkey.rotation_axis_angle = vectors[n]