Does bmesh guarantee order of vertex in a edge according to the time they were added I would like to represent directions in a single edge in a mesh similar to a curve. 

However I cannot use curves themselves due to performance issues when dealing with huge amount of curves and not able to have more than 2 endpoints. Linked:
Creating/Handling pathfinding large world waypoints in Blender

So I am thinking to represent the order of vertices in an edge as direction. Then I would simply either expose the direction via a gl draw or in the UI panel and also allow to invert the direction by deleting the edge and adding it again in opposite order.

For instance, I add a place using the following code

    vertex1 = bm.verts.new((-0.5, -0.5, 0.0)) # 0
    vertex2 = bm.verts.new( (0.5, -0.5, 0.0)) # 1
    vertex3 = bm.verts.new( (0.5, 0.5, 0.0) ) # 2
    vertex4 = bm.verts.new( (-0.5, 0.5, 0.0)) # 3

    bm.verts.index_update()

    # How to create edges    
    bm.edges.new( (vertex1, vertex2) ) # 0 
    bm.edges.new( (vertex2, vertex3) ) # 1
    bm.edges.new( (vertex3, vertex4) ) # 2
    bm.edges.new( (vertex4, vertex1) ) # 3


Later when I access the bmesh again, using 

    obj = bpy.context.edit_object
    me = obj.data 
    bm = bmesh.from_edit_mesh(me)

    # are the vertices stored in edge guarentee'd to be
    # in order at the time they were added
    bm.edges[0].verts[0] == bm.verts[0]
    bm.edges[0].verts[1] == bm.verts[1]


I have tested it and it does seem correct for few of my test cases. However, I would like to be certain.
The order of both vertex indices stored in an edge (me.edges[#].verts), or the vertices references respectively (bm.edges[#].verts) should not change unless you make it so (by flipping an edge I guess).

That is because the order is important. Think of a line of 3 connected edges: if you flipped the vertices referenced by the middle edge, it would still look like the same edge. But both vertices would need to be switched (as in coordinates) - but that doesn't work because they are used by other edges.

(Unscientific) stress test:

import bpy

for area in bpy.context.screen.areas:
    if area.type == 'VIEW_3D':
        break
else:
    raise Exception()

ctx = bpy.context.copy()
ctx["area"] = area
ctx["region"] = area.regions[-1]
bpy.ops.ed.undo_push(ctx)

ob = bpy.context.object
me = ob.data
bpy.ops.object.mode_set(mode='EDIT')
print()
print()
t = {}
for e in me.edges:
    v1, v2 = e.vertices
    t[(v1, v2)] = "{: 8.4f} {: 8.4f} {: 8.4f} | {: 8.4f} {: 8.4f} {: 8.4f}".format(*(me.vertices[v1].co[:] + me.vertices[v2].co[:]))

bpy.ops.mesh.sort_elements(type='RANDOMIZE')

bpy.ops.ed.undo_push(ctx)
me.update()
bpy.ops.transform.translate(value=(0,0,2))
bpy.ops.object.editmode_toggle()
bpy.ops.object.editmode_toggle()
bpy.ops.ed.undo(ctx)
bpy.ops.mesh.sort_elements(type='RANDOMIZE')
bpy.ops.mesh.flip_normals()

for e in me.edges:
    v1, v2 = e.vertices
    if t[(v1, v2)] != "{: 8.4f} {: 8.4f} {: 8.4f} | {: 8.4f} {: 8.4f} {: 8.4f}".format(*(me.vertices[v1].co[:] + me.vertices[v2].co[:])):
        raise Exception()


The order of geometry elements is not guaranteed on the other hand (me..index / bm..index), although it is practically stable, as long as you add and remove in stack principle (last in, first out).

If a vertex / edge / face other than the one with the highest index is deleted, all elements with a higher index will shift to fill the gap.

There are also mesh operations to sort mesh elements based on properties like distance to camera, or randomize the indices even.

You could use a custom data layer to store an explicit order (e.g. bm.verts.layers.int), but beware of the value propagation: if you extrude a vertex for instance, that custom data will be carried over to the new vertex.

This patch might be interesting for you: https://developer.blender.org/T42865