This is a very crude attempt to explain the growth rate of $13$,  can be almost certainly be improved upon.

In the situation where $20$ is the diameter, there are $18$ moves (or monoid generators), with three moves for each face, i.e. two quarter-turns and on half-turn.

For each position except for the starting position, moving the same face as on the previous move will not yield a new position, so that reduces the number of moves to new positions to $15$.

I think the main reason why we get $13$ rather than $15$ is that moves of opposite faces commute. In approximately $1/5$ of the positions, the penultimate move commutes with the last move. In those cases, moving the opposite face to the final move will not give a new position (but unfortunately I think the $1/5$ is an overestimate).

In the remaining $4/5$ of the positions, moving the opposite face will give a position that could also have been reached by moving the opposite face first. So the number of new positions reached in this way should be divided by two.

That yields an estimate of $12/4 + 4(12 + 3/2)/5 = 66/5$, which is just slightly more than $13$.

But on reflection, I think the proportion of positions where the final two moves commute is closer to $1/10$ than $1/5$, giving the estimate $12/10 + 9(12+3/2)/5 = 267/20 = 13.35$.

Of course these are all overestimates, because there are other group relations that reduce the number of positions, but I believe that I have identified the principal relations involved that reduce the growth rate from that of the free monoid.

*Added later* After writing the above, I realized that my estimate is calculating the growth rate of the group with the presentation
$$\langle a_i, b_i\ (1 \le i \le 6) \mid a_i^4= 1,a_i^2=b_i\ (1 \le i\le 6), a_1a_4=a_4a_1,a_2a_5=a_5a_2,a_3a_6=a_6a_3 \rangle,$$

and that I could do that exactly on the computer. So I did, and the principal term in the growth function came out as (a multiple of) $(13.3485)^n$, which is surprsingly (to me) close to my second crude estimate of $13.35$.

I did this computation in Magma - here is the code, followed by the Taylor series expansion of the growth function $f$.

    > G := Group<a,a2,b,b2,c,c2,d,d2,e,e2,f,f2 | a^4,b^4,c^4,d^4,e^4,f^4,
    >                          a^2=a2,b^2=b2,c^2=c2,d^2=d2, e^2=e2,f^2=f2,
    >                          a*d=d*a, b*e=e*b, c*f=f*c >;
    > A := AutomaticGroup(G);
    > f := GrowthFunction(A);
    > PZ<x> := FunctionField(IntegerRing());
    > PZ!f;
      (-9*x^2 - 6*x - 1)/(18*x^2 + 12*x - 1)
    > PR := PolynomialRing( RealField(6) );
    > [ r[1]^-1 : r in Roots(PR!Denominator(f)) ];
      [ -1.34847, 13.3485 ]

    > LR<x> := LaurentSeriesRing(Z);
    > LR!f;
    1 + 18*x + 243*x^2 + 3240*x^3 + 43254*x^4 + 577368*x^5 + 7706988*x^6 + 
      102876480*x^7 + 1373243544*x^8 + 18330699168*x^9 + 244686773808*x^10 + 
      3266193870720*x^11 + 43598688377184*x^12 + 581975750199168*x^13 + 
      7768485393179328*x^14 + 103697388221736960*x^15 + 1384201395738071424*x^16 +
      18476969736848122368*x^17 + 246639261965462754048*x^18 + 
      3292256598848819251200*x^19 + O(x^20)