Not very efficient, but you may find it useful for some experiments. I perused the code from the link you provided (kuba's), although there are better alternatives in the answers. ClearAll[spiral, genTri, mp]; spiral[n_?OddQ] := Nest[ With[{d = Length@#, l = #[[-1, -1]]}, Composition[ Insert[#, l + 3 d + 2 + Range[d + 2], -1] &, Insert[#\[Transpose], l + 2 d + 1 + Range[d + 1], 1]\[Transpose] &, Insert[#, l + d + Range[d + 1, 1, -1], 1] &, Insert[#\[Transpose], l + Range[d, 1, -1], -1]\[Transpose] & ][#]] &, {{1}}, (n - 1)/2]; genTri[n_] := genTri[n] = IntegerQ[(-1 + Sqrt[1 + 8 #])/2] & /@ Range@n mp[n_] := mp[n] = Image[Unitize[spiral[n] /. Thread[Flatten@Position[genTri[n^2], True] -> 0]]] Erosion[mp[601], 1] ![Mathematica graphics](https://i.sstatic.net/xPwHy.png) There they are, your 17 arms Erosion[mp[1001], 1] ![Mathematica graphics](https://i.sstatic.net/Q9z9z.png) For the "stability" of the number of arms a reasonable condition (rule of dumb thumb) is that the density of triangular numbers remain almost constant in each "layer" of the spiral, and that is what effectively seems to happen: r[n_] := Range[(2 n + 1)^2 + 1, (2 n + 3)^2] f[n_] := Count[IntegerQ[(-1 + Sqrt[1 + 8 #])/2] & /@ r[n], True] ListLinePlot[f /@ Range[200]] ![Mathematica graphics](https://i.sstatic.net/nCSGn.png) The maths are left as an exercise :)