Why is the geometric solution of minimizing error via orthogonality called a "least squares" solution?

You write “the best approximate solution $\hat x$”, but the best by what criterion? What you describe is the solution that minimises $\left\|A\hat x-b\right\|$. For the euclidean norm, but that is a choice. Now minimising $\left\|A\hat x-b\right\|$ is the same as minimising $\left\|A\hat x-b\right\|^2=\sum_i (A\hat x-b)_i^2$ and those are the squares to which “least squares” refers. So you are indeed solving a least squares problem. Note that if you see $Ax=b$ as a system of equations then $(A\hat x-b)_i$ is the error for the $i$-th equation.

I'll focus on this piece:

Or is there a more fundamental reason that the geometric condition of orthogonality is itself an expression of "least squares," independent of the calculus perspective?

Note that for any $\mathbf x$, the error vector $\mathbf e(\mathbf x) = \mathbf b - A \mathbf x$ can always be decomposed into two components: one parallel to (i.e., an element of) the column space of $A$ and the other perpendicular. That is, $$ \mathbf e(\mathbf x) = \mathbf e_{\|}(\mathbf x) + \mathbf e_{\perp}. $$ Note that the parallel component depends on our choice of $\mathbf x$, but the perpendicular component does not. Indeed: for any attainable error vectors $\mathbf e(\mathbf x_1), \mathbf e(\mathbf x_2),$ the difference $\mathbf e(\mathbf x_2) - \mathbf e(\mathbf x_1)$ is always parallel to the column space.

With that, we can see by the Pythagorean theorem (generalized to $n$-dimensional space) that the sum of squares error corresponding to a given $\mathbf x$ is given by $$ \|\mathbf e(\mathbf x)\|^2 = \|\mathbf e_{\|}(\mathbf x)\|^2 + \|\mathbf e_{\perp}\|^2 $$ So, finding $\mathbf x$ such that $\mathbf e(\mathbf x)$ is orthogonal to the column space of $A$ (equivalently, such that $\mathbf e_\| = 0$) is equivalent to minimizing the sum-of-squares error.

Yes, it is called "least squares" because it solves the problem of minimizing the sum of the squares of the errors. Geometrically, as you say we are finding a "closest possible vector", and "closest" is defined in terms of the Euclidean norm, which is the square root of the sum of the squares of the coordinates.

According to https://mathshistory.st-andrews.ac.uk/Miller/mathword/m/, the term "Method of least squares" (in French: "la Méthode des moindres quarrés") goes back to Legendre in 1805. This was before the development of the modern machinery of linear algebra. So Legendre and Gauss would not be using vectors and matrices as such.

We want to make the error $\mathbf e=\mathbf b-A\mathbf x$ as close as possible to zero, so we can minimize the squared error $e^2=\mathbf e\cdot\mathbf e$. We can actually solve for this using calculus by setting $\frac{\partial}{\partial \mathbf x}e^2=0$. Working with vector derivatives makes me confused, so I will be using Einstein summation convention. Everytime you see two repeated indices, there is an implicit sum. Everytime you see an unpaired index you are viewing it as vector component. For example, $x_ix_i\equiv \sum_{i=1}^n x_ix_i$ and $y_i=2x_i$ means that the i-th component of $y$ is twice the i-th component of $x$. Using this notation we get \begin{align} \frac{\partial}{\partial x_i}e^2&=\frac{\partial}{\partial x_i}(b_j-A_{jk}x_k)(b_j-A_{jl}x_l)\\ &=\frac{\partial}{\partial x_i}\left[b_jb_j-2b_jA_{jk}x_k+x_kA_{jk}A_{jl}x_l\right] \end{align} The $b_jb_j$ term drops out, so let's look at the second and third term individually. Here $\delta_{ij}$ is the Kronecker delta. Second term: \begin{align} \frac{\partial}{\partial x_i}\left[-2b_jA_{jk}x_k\right]&=-2b_jA_{jk}\delta_{ik}\\ &=-2b_jA_{ji}\\ &\equiv-2A^T\mathbf b. \end{align} Third term:

\begin{align} \frac{\partial}{\partial x_i}\left[x_kA_{jk}A_{jl}x_l\right]&=\delta_{ik}A_{jk}A_{jl}x_l+x_kA_{jk}A_{jl}\delta_{il}\\ &=A_{ji}A_{jl}x_l+x_kA_{jk}A_{ji}\\ &=2A_{ji}A_{jl}x_l\\ &\equiv2A^TA\mathbf x. \end{align} Combined, they give \begin{align} \frac{\partial}{\partial \mathbf x}e^2&=-2A^T\mathbf b+2A^T\mathbf x\\ &=2A^T(A\mathbf x -\mathbf b)=0 \end{align}

Einstein notation may seem confusing at first, but all the individual steps make total sense. Converting back to vector notation is actually the hard part.