The distance between strings is in fact an infrared cut-off for the divergent string tension. To see this, note that in cylindrical coordinates $(r, \theta, z)$,
$$\Phi \simeq
\begin{cases}
0 &\text{ for } r=0\,,\\
f_a\ e^{i\theta} &\text{ for } r \gg\delta \,,
\end{cases}$$
is a static, straight axionic string stretched along the $z$-axis. With this, one can compute the string tension to be
$$ \mu = \int d^2x \left[ |\nabla\Phi|^2 + \frac{\lambda}8 \left(|\Phi|^2-f_a^2\right)^2\right] \simeq \int_{r>\delta} d^2x\ \ f_a^2|\nabla e^{i\theta}|^2 = 2\pi f_a^2\ln \left(\frac L\delta\right)\,,$$
where $L$ is an infrared cut-off, which is then interpreted as the distance between strings.
By definition, the string tension $\mu$ measures the energy per unit length of the string. If $\xi$ is a parameter measuring the average length of a string present in a volume $L^3$, then the energy density of strings is given by
$$ \rho = \xi \frac{\mu(L)}{L^2} \,.$$
Now, in the background of an expanding radiation-dominated universe with Hubble parameter $H = 1/2t$, since long strings extend across the horizon, the infrared cutoff can be taken to be a Hubble horizon, $L \sim 1/H \sim t$. The parameter $\xi$ is of order $1$. Thus,
$$ \rho \sim \frac{\mu(L)}{L^2}\Bigg|_{L = t} = \frac{2\pi f_a^2}{t^2} \ln \left(\frac t\delta\right)\,.$$