here is my elegant implementation

l[c_]:=TakeList[Range@Sum[k!,{k,c}],Range@c!][[c-1]];

T[x_]:=Graph[(F=Flatten)@Table[MapThread[#->#2&,{Sort@F@Table[l@i,i],l[i+1]}],{i,2,x+1}]];



T@3   

which returns

but if your Mathematica version doesn't support TakeList here is another way

s[x_] := Sum[k!,{k,x}];

z[y_] := Partition[Range[s@y+1,s[y+1]],1+y];

v[n_] := Table[{Flatten[z[n-1]][[i]]->z[n][[i,j]]},{i,n!},{j,n+1}];

tree[t_] := Graph[Flatten[Array[v@#&,t],3]];



tree@3    

tree@6    

Update 2: a more streamlined version for 2D graphs:

ClearAll[g]

g = GraphComputation`ExpressionGraph[ConstantArray[x, #], VertexLabels -> None] &;

Examples:

g[Range[2, 4]]

SetProperty[g[{3, 1, 3, 1, 2, 1, 4}],

  {GraphLayout -> "RadialEmbedding", EdgeShapeFunction -> "Line"}]

Original answer:

ClearAll[f]

f[g_: Graph][n_List, o : OptionsPattern] := g[UndirectedEdge @@@ EdgeList@

    GraphComputation`ExpressionGraph[ConstantArray[x, n]], 

     o, GraphLayout -> {"BalloonEmbedding"}, ImageSize -> Large]

f[g_: Graph][n_Integer, o : OptionsPattern] := f[g][Range[n], o]

Examples:

f[6]

f[6, GraphLayout -> {"RadialEmbedding"}]

g1 = f[Graph3D][6]

g2 = f[Graph3D][6, GraphLayout -> {"RadialEmbedding"}]

Use a list for number of vertices on each layer as the argument:

f[{3, 5, 2, 4}, GraphLayout -> {"RadialEmbedding", "LayerSizeFunction" -> (# &)}]

Update: ... to take the Graph3D and somehow improve the layout on the high-n layers: There seems to be lots of wasted space.

One way to change the box ratios without distorting the vertex shapes is to modify the VertexCoordinates using ScalingTransform:

SetProperty[g1, VertexCoordinates -> ScalingTransform[{1, 1, 3}][GraphEmbedding@g1]]

SetProperty[g2, VertexCoordinates -> ScalingTransform[{1, 1, 3}][GraphEmbedding@g2]]

Or add the suboption "LayerSizeFunction" in "RadialEmbedding":

g3 = f[Graph3D][6, GraphLayout -> {"RadialEmbedding", "LayerSizeFunction" -> (# &)}];

SetProperty[g3, VertexCoordinates -> ScalingTransform[{1, 1, 3}][GraphEmbedding@g3]]

SetProperty[#, VertexCoordinates -> ScalingTransform[{1, 1, 3}][GraphEmbedding@#]] &@

 f[Graph3D][{3, 5, 2, 4}, GraphLayout -> {"RadialEmbedding", "LayerSizeFunction"->(#&)}]

IGraph/M already has this built-in as IGSymmetricTree. You can specify the number of branches at each level.

IGSymmetricTree[

 Range[2, 4],

 DirectedEdges -> True,

 GraphLayout -> "LayeredEmbedding"

]

The implementation is mostly in C (not Mathematica) and not from igraph this time. This is simply easier to implement procedurally, for which C is a good fit. This is why I did not do it in pure WL.

Here's another structure, with a different number of branches at each level.

IGSymmetricTree[{5, 4, 3, 2}]

I don't know if you find this elegant. But I give it a try.

maxdepth = 5;

Graph[

 Transpose[{

   Join @@ MapIndexed[

     {x, i} [Function] Join @@ Transpose[ConstantArray[Range @@ x, (i + 1)]], 

     Join[{{1}}, Partition[Accumulate[Range[maxdepth - 1]!], 2, 1] + ConstantArray[{1, 0}, maxdepth - 2]]

     ],

   Range[2, Total[Range[maxdepth]!]]

   }],

 DirectedEdges -> True,

 GraphLayout -> "BalloonEmbedding"

 ]

Edit

Out of curiosity, I adapted the algorithm above to produce also other symmetric trees.

SymmetricTree[branchlist_?(VectorQ[#, IntegerQ] &)] := 

 Module[{levelnodecounts}, 

  levelnodecounts = FoldList[#1 #2 &, 1, branchlist];

  Graph[Transpose[{

     Join @@ MapIndexed[

       {x, i} [Function] Join @@ Transpose[ConstantArray[Range @@ x, branchlist[[i[[1]]]]]],

       Join[

        {{1}},

        Partition[Accumulate[Most[levelnodecounts]], 2, 1] + ConstantArray[{1, 0}, Length[branchlist] - 1]

        ]

       ],

     Range[2, 1 + Total[Rest[levelnodecounts]]]}], 

   DirectedEdges -> True

   ]

  ]

Regarding speed, it seems to be on par with IGSymmetricTree. Of course, I cannot provide such a detailed user interface as Szabolcs so I would suggest to use IGraphM whenever possible.

Edit 2

Adapting my (slow) code for fractal trees, here is another way to embedd the tree:

BoccoliEmbedding[branchlist_] := 

 Module[{data0, data, θ, stem, thickness, s1, s2, f, F},

  θ = Pi/4.;

  s1 = 1/GoldenRatio // N;

  s2 = 1/GoldenRatio // N;

  stem = {0., 0., 1.};

  thickness = 0.15;

  data0 = {Join[

     {{0., 0., 0.}},

     {stem},

     {{thickness, 1., 0.}},

     Table[

      RotationMatrix[2. k Pi/branchlist[[1]], {0, 0, 1}].{Cos[θ], 0.,Sin[θ]}, 

      {k, 0, branchlist[[1]] - 1}]

     ]

    };

  f = {U, n} [Function] Table[

     Join[

      {U[[1]] + U[[2]]},

      {U[[i]]},

      {s2 U[[3]]},

      Dot[

       s1 Table[RotationMatrix[2. Pi j/n, U[[2]]].U[[i]], {j, 0, n - 1}], 

       RotationMatrix[{U[[i]], U[[2]]}]

       ]

      ],

     {i, 4, Length[U]}];

  F = {data, n} [Function] Join @@ (f[#, n] & /@ data);

  data = Join @@ FoldList[F, data0, Join[Rest[branchlist], {1}]];

  data[[All, 1]] + data[[All, 2]]

  ];

And this is how we apply it:

b = Range[2, 7];

plot = Graph[

  EdgeList[SymmetricTree[b]],

  VertexCoordinates -> BoccoliEmbedding[b]

  ]