Here's the algorithm I'm using:
Pastebin -- Tetrahexacontree Voxel Volume Raymarching
P.S. I'm using slab method ray box intersection. It's more or less the same code (GLSL + 3D) as this Scratchpixel - A Minimal Ray-Tracer. I'm also having to use dynamic epsilon values to slightly enlarge bounding boxes for intersection to make it work correctly for nearest side and farthest side ray bounding box intersections to get correct "t" values inside or outside of bounding box. Additionally, sometimes the bounding box intersections fail, so the code has to handle that gracefully or GLSL shader threads will get stuck in infinite loops.
Tree nodes (including the root node) are stored in an array of nodes. "tree root node index" is an integer variable containing the index of the root node in said array. Many array index based tree implementations have the root at the first element in the array (index zero), but in my implementation the root node can be any node in the array of nodes.The reason I've implemented things this way is so that if the tree's bounds need to be expanded (voxel volume(s) need to be added to the tree outside of tree's current bounds [in C++ code]), the tree doesn't need to be rebuilt, but instead tree's bounds are expanded accordingly and a free/unused node in the array of nodes (which is a dynamic size data structure like a C++ std::vector of tree nodes) becomes the new root node and the old root node becomes one of new root node's children.
P.S. Disregard the part I set as strike though text in the previous reply. That was left over code from the octree version and I had to re-implement tree bounds size increase by rebuilding the tetrahexacontree due to integer coordinate powers of two alignment restrictions if the old root was made a child of a new root. So "tree root node index" is zero (always first element of array of nodes) and "tree root node index" is no longer passed to the tetrahexacontree voxel volumes raymarching shader via a uniform buffer object field (just hard coded now).