In $k$-means, how is it NP-hard if the dimensionality of the data is at least $2$ ($d\geq 2$)? In $k$-means, how is it NP-hard if the dimensionality of the data is at least $2$ ($d\geq 2$)? Can someone justify or give reasons to this statement?
Any guidance would be appreciated.

First, note that it's not terribly surprising that the problem is NP-hard: there are a lot of possible clusterings; the usual algorithm doesn't claim to find the optimum solution; the usual algorithm is purely local and so might well miss a global optimum in some weird data setup.
To be NP-hard, the problem only has to be hard in some non-negligible fraction of datasets; in particular, it might still be easy for any data set that genuinely has $k$ clusters but hard for some data sets that don't really have $k$ clusters. Also, $k$ is part of the input, so we're allowing $k$ to grow with $n$; the fixed-$k$ problem isn't NP-hard.
The paper linked in @jbowman's comment shows how to construct planar point layouts where lots of the interpoint differences are equal and the number of clusters is roughly $n/2$. The question of which points optimally go in the same cluster then becomes a global rather than a local problem and local algorithms aren't going to work well.  For the precise details of how this gets reduced to the known-hard problem 'planar 3SAT' you'll need to read the paper. From a data analysis point of view it's not really a problem, since it's not a setting where you'd be using $k$-means, but it does illustrate why it might be hard to prove useful results about $k$-means even in settings where it's useful.
One of the links from the paper is to this, which looks at how to modify the usual algorithm so it has probabilistic guarantees of getting close to the optimal clustering in polynomial time even for large $k$, assuming there really is a good $k$-cluster solution.
