As I indicated above, $\log{N(d)} \asymp d^2$ reduces to a problem about irreducible polynomials, the very likely affirmative answer to which would prove the lower bound $n(h,d) \geq (h-o(1))d^2$, as $d \to_h \infty$, for the logarithm of the number of points in $\mathbb{G}_m(\bar{\mathbb{Q}})$ of degree $d$ and height $< h$. Northcott's upper bound is $n(h,d) \leq (h+\log{2})d^2 + O(d)$.
I nonetheless leave below the remarks I had made.
I do find Will's model to be a convincing evidence that the logarithmic asymptotic count should likewise be $\log{N(d)} = (1+o(1))d^2$. Below I suggest an explanation for the apparent discrepancy in the two models as far as the distribution of points of very small height is concerned. Taking into account this and the (Schanuel-Schmidt)-Masser-Vaaler theorem, I would think that the statistical parallels between number and function fields are much too precise to think of any other answer.
Let $n(h,d)$ be the logarithm of the number of points of degree $d$ and logarithmic height $< h$. Northcott's estimate, $n(h,d) \leq (h+\log{2})(d^2+d) + O(d)$, is valid also in the function field case (where the term $\log{2}$, which arises from the triangle inequality at the archimedean places, is not needed); but there, by an extension of Will's argument, we have a precise asymptotic: Assuming as we may that $h \in \frac{\log{q}}{d} \mathbb{Z}^{> 0}$, then in fact $n(h,d) = h(d^2+d) - 1 + o(1)$ as $\max(h,d) \to \infty$. In other words: in our function field model, Northcott's bound is always sharp.
For the number field case, Masser and Vaaler, extending previous work of Schanuel and Schmidt, consider the count for $d$ fixed and $h \to \infty$; they prove in particular that $n(h,d) = h(d^2+d) + k(d) + o(1)$, with an explicit term $k(d) \asymp d\log{d}$. (Their result is more precise than this, and it applies more generally to relative extensions over a fixed number field). Hence, in this regime, Northcott's logarithmic bound is again sharp.
In $\overline{\mathbb{Q}}^{\times}$, however, in contrast to the situation for our function field model, Northcott's prediction breaks down completely when we take $h = c/d$, where $0 < c < \infty$ is a constant: $\exp(hd^2)$ is then exponential in $d$, whereas we expect $\# \{ \alpha \mid h(\alpha) < c/d, \, [\mathbb{Q}(\alpha):\mathbb{Q}] \leq d \}$ to be polynomially bounded, perhaps even by $\kappa(c)d^{1+\exp(c)}$. (This looks like an interesting problem; does such a bound appear anywhere in the literature?)
Here, I believe, is the explanation for this apparent discrepancy (exponential versus polynomial) in the distribution of points of very small height - the torsion points in particular. There are two things to keep in mind. The first is the distinction between general Weil heights and canonical (dynamical, normalized, Neron-Tate) heights: they differ by a bounded amount and therefore they produce comparable asymptotics - certainly the same rate of growth - for large heights; whereas the distribution of small points is a very fine intrinsic property of the latter heights. The second point to keep in mind is that our particular function field model here is isotrivial (constant, in fact); this accounts for the profusion of torsion (as well small non-torsion) points. An analogy here is to think of a constant elliptic curve over a complex function field: its points of zero canonical height are all the constant sections (an uncountable set), whereas for non-isotrivial elliptic curves they only comprise of the countable set of torsion points.
In either setting, think of our height $h(\cdot)$ as the canonical dynamical height - the global Tate limit attached by Call and Silverman to a dynamical system on $\mathbb{P}^1$ - of the iteration $z \mapsto z^2$. The difference is that in the function field case, this iteration is isotrivial. Isotriviality has no counterpart in the number field setting, and a more accurate function field model in our problem would be to consider the dynamical height $\hat{h}_f$ attached to a non-isotrivial map of $f : \mathbb{P}^1 \to \mathbb{P}^1$ over $\mathbb{F}_q(t)$; or, if you prefer, the canonical height on a non-isotrivial elliptic curve over $\mathbb{F}_q(t)$. In those cases, as in the arithmetic setting, the number of points of degree $\leq d$ and height $< c/d$ should be similarly bounded polynomially in $d$, and the discrepancy in the count will not occur.
In any case, $\hat{h}_f - h = O_f(1)$, and this implies that for $h \gg_ f 0$ sufficiently large, the logarithmic count $n_f$ for $\hat{h}_f$ still satisfies $n_f(h,d) \asymp h(d^2+d)$. (A moment of reflection will show that the proportionality constant should still be $1$, although this could be difficult to prove; all we know at this point is that $n_f(h,d)$ is locked between two positive multiples of $h(d^2+d)$.) But for any fixed positive $h > 0$, we may construct enough points of height $< h$ by taking inverse images under $f^{\circ N}$ of points of height $< h \cdot (\deg{f})^N$, which will be large enough as soon as $N \gg 0$. Hence, for any fixed $h > 0$, we still know that $n_f(h,d) \asymp h(d^2+d)$.
This now should hold in the number field setting too.
Let me finish by adding one more remark about the small point analogy. It is often said that the crux of Lehmer's problem, stating $h(\alpha) > c/d$ unless $\alpha$ is torsion, is in the presence of archimedean places, and that the question is pointlessly trivial in the function field setting. This of course is true for the constant function field model considered in Will's answer, and it is true more generally for polarized dynamical systems with everywhere potential good reduction over a global function field (e.g., abelian varieties with everywhere potential good reduction). However, the real analogy arises when we have degenerations; for instance, for the case of the canonical height on a non-isotrivial elliptic curve over $\mathbb{F}_q(t)$, where I believe that Lehmer's problem would be every bit as difficult as the original one for $\mathbb{G}_m$ over $\bar{\mathbb{Q}}$.
If this is correct then, from a point of view of algebraic dynamics, the essential difficulty in Lehmer's problem is not so much in the failure of the ultrametric triangle inequality in number fields, as it is in the presence of places of degenerating - or chaotic - dynamics (though of course the two points are related). After all, it is a common wisdom in arithmetic geometry that archimedean fibres (or complex spaces) should be regarded as totally degenerate: think of the Tate curve, Mumford's theory of $p$-adic Shottky groups, or the non-triviality of Julia sets in Berkovich analytic spaces (they are always trivial for iterations with good reduction over a non-archimedean local field).
In view of this, all the subtleties in the distribution of small points in $\mathbb{G}_m(\bar{\mathbb{Q}})$ appear to be present also in a non-isotrivial dynamical system on $\mathbb{P}^1$ over $\mathbb{F}_q(t)$, or in a non-isotrivial elliptic curve over a global function field. But we saw that in such a dynamical system, $n_f(h,d) \asymp h(d^2+d)$ for any fixed $h > 0$ (or for any fixed $d$). So the same should persist in the number field setting too - presumably, with proportionality constant equal to $1$.
To sum up, as soon as the height is bounded away from zero, Northcott's bound is pretty sharp. This is something I did not expect at the time of asking this question.
