As mentioned by Joel in the comments, I just posted this week on arxiv a paper that explains the relationship between some of these points using the geometric Satake correspondence (arxiv.org/abs/2310.00855).
The idea is to construct a ring isomorphism between $H^\bullet(G(n,m))$ and a suitable quotient of the ring of symmetric polynomials $\mathbb{Z}[x_1,\ldots, x_n]^{S_n}$ which sends Schubert cycles to Schur polynomials by first constructing it as a representation isomorphism and then upgrading to an isomorphism of rings with bases. Moreover, this can be generalized in a straightforward way to the torus-equivariant setting, with double Schur polynomials (you can find the details of this in the paper). Let me explain roughly how this is done in the non-equivariant setting.
First, the boson-fermion correspondence usually concerns an isomorphism between the semi-infinite wedge space (or fermionic Fock space), which is a representation of the Lie algebra $\mathfrak{gl}_\infty$ of $\mathbb{Z}\times \mathbb{Z}$ matrices, with the algebra $\Lambda$ of symmetric functions in infinitely many variables. I will however discuss a variant of it involving an isomorphism between the $n^{th}$ wedge representation of the Lie algebra $\mathfrak{gl}_\infty^+$ of $\mathbb{N}\times \mathbb{N}$ matrices and the algebra of symmetric polynomials in $n$ variables.
Define $\mathfrak{gl}_\infty^+$ to be the Lie algebra of $\mathbb{N}\times \mathbb{N}$ matrices with finitely many nonzero diagonals. It acts naturally on any vector space with a basis indexed by $\mathbb{N}$, in particular on $\mathbb{C}[x]$ (using the basis $1,x,x^2,\ldots$). This induces an action on $\bigwedge^n \mathbb{C}[x]$, which we can identify with the vector space of skew-symmetric polynomials in $n$ variables, which I'll denote $\Lambda_n^\text{sgn}$. Moreover, multiplication by the Vandermonde determininant $\delta=\Pi_{i<j} (x_i-x_j)$ defines an isomorphism between the vector space of symmetric polynomials in $n$ variables (which I'll denote $\Lambda_n$) and the skew-symmetric polynomials $\Lambda_n^\text{sgn}$.
The "finite $n$" boson-fermion correspondence is then simply the isomorphism
$$\bigwedge^n \mathbb{C}[x] = \Lambda_n^\text{sgn} \underset{\cdot \delta}{\overset{\sim}{\longleftarrow}}\Lambda_n$$
(the left hand side would be called by physicists the fermionic Fock space and the right hand side the bosonic Fock space).
Moreover, an important feature of this boson-fermion correspondence is that it is equivariant for the action of the Heisenberg subalgebra. I'll discuss only the action of the negative part of the Heisenberg subalgebra, which is the abelian subalgebra $\mathfrak{h}^-$ consisting of lower triangular matrices whose entries are constant on diagonals. Identifying $\mathfrak{gl}_\infty^+$ with its image inside $\text{End}_\mathbb{C}(\mathbb{C}[x])$, these are exactly the "multiplication by $p$" operators for $p\in \mathbb{C}[x]$. Alternatively, $\mathfrak{h}^-$ is the centralizer of
$$f=\begin{pmatrix}
0 & 0 & 0 & 0 & \cdots\\
1 & 0 & 0 & 0 & \cdots\\
0 & 1 & 0 & 0 & \cdots\\
0 & 0 & 1 & 0 & \cdots\\
\vdots & \vdots & \vdots & \vdots & \ddots
\end{pmatrix}$$
(the "multiplication by $x$" operator).
If we make $\mathfrak{h}^-$ acts on $\Lambda$ by identifying $\mathfrak{h}^-=\mathbb{C}[x]$ as above and letting $p\in \mathbb{C}[x]$ act via multiplication by $p(x_1)+\ldots + p(x_n)$, then it is clear from its construction that the boson-fermion correspondence is $\mathfrak{h}^-$-equivariant.
It is also clear from the definition that the image of the weight basis elements $x^{\lambda_1+n-1}\wedge x^{\lambda_2+n-2}\wedge \cdots \wedge x^{\lambda_n}$ for $\lambda_1 \ge \cdots \ge \lambda_n \ge 0$ under the boson-fermion correspondence are the Schur polynomials, defined as
$$s_\lambda = \frac{\sum_{\sigma \in S_n} \text{sgn}(\sigma) x_{\sigma(1)}^{\lambda_1+n-1} x_{\sigma(2)}^{\lambda_2+n-2}\cdots x_{\sigma(n)}^{\lambda_n}}{\sum_{\sigma \in S_n} \text{sgn}(\sigma) x_{\sigma(1)}^{n-1} x_{\sigma(2)}^{n-2}\cdots x_{\sigma(n)}^{0}}=\frac{\sum_{\sigma \in S_n} \text{sgn}(\sigma) x_{\sigma(1)}^{\lambda_1+n-1} x_{\sigma(2)}^{\lambda_2+n-2}\cdots x_{\sigma(n)}^{\lambda_n}}{\delta}.$$
To get a representation of $\mathfrak{gl}_m$ rather than $\mathfrak{gl}_\infty^+$, we have to take a quotient. Precisely, let $J_m$ be the ideal of $\mathbb{C}[x_1, \ldots, x_n]$ generated by $x_i^m$ for $i=1, \ldots, n$, $I_m^\text{sgn}=J_m \cap \Lambda_n^\text{sgn}$ and $I_m=\{p \in\Lambda_n\mid \delta p \in I_m^\text{sgn}\}$. Then $I_m$ is spanned by Schur polynomials indexed by partitions with some part larger than $m-n$ and, if we let $\Lambda_{n,m}=\Lambda_n/I_m$, $\Lambda_{n,m}^\text{sgn}=\Lambda_n^\text{sgn}/I^\text{sgn}_m$, then the boson-fermion correspondence induces an isomorphism
$$\bigwedge^n \mathbb{C}[x]/(x^m) = \Lambda_{n,m}^\text{sgn} \underset{\cdot \delta}{\overset{\sim}{\longleftarrow}}\Lambda_{n,m}.$$
The representation of $\mathfrak{gl}_m$ on the left hand side (which comes from identifying $\mathbb{C}[x]/(x^m)$ with $\mathbb{C}^m$ via the basis $1,x,\ldots, x^{m-1}$) can therefore be transferred to a representation of $\mathfrak{gl}_m$ on the right hand side.
Moreover, we can define $\mathfrak{h}^-_m$ to be the abelian subalgebra of lower triangular matrices constant on the diagonals. Just as for $\mathfrak{h}^-$, we can identify $\mathfrak{h}^-_m$ with $\mathbb{C}[x]/(x^m)$ and then $p\in \mathbb{C}[x]/(x^m)$ act on $\Lambda_{n,m}$ via multiplication by $p(x_1)+\ldots+p(x_n)$.
We now turn to $H^\bullet(G(n,m))$. The affine Grassmannian $Gr=GL_m(\mathbb{C}(\!(t)\!))/GL_m(\mathbb{C}[\![t]\!])$ can be identified with the set of $\mathbb{C}[\![t]\!]$-lattices inside $\mathbb{C}(\!(t)\!)^n$. The $GL_m(\mathbb{C}[\![t]\!])$-orbits are indexed by coweights of $GL_m$ (which are the same as weights because $GL_m$ is self dual) an the geometric Satake correspondence tell us that the intersection cohomology of the closure of the orbit indexed by $\mu$ is isomorphic to the representation of $GL_m$ of highest weight $\mu$. When $\mu=(\underbrace{1, \ldots, 1}_{n}, 0, \ldots, 0)$, the orbit consists of the lattices sandwiched between $\mathbb{C}[\![t]\!]^m$ and $t\mathbb{C}[\![t]\!]^m$ and of index $n$ inside $\mathbb{C}[\![t]\!]^m$, which is isomorphic to $G(n,m)$. We get therefore that $H^\bullet(G(n,m))$ is isomorphic to the $n^{th}$ wedge representation of $\mathfrak{gl}_m$. The general theory of geometric Satake also tell us that the Schubert cycles form a weight basis for the action of the maximal torus.
Putting all this together, we have an isomorphism of representations of $GL_m$
$$\Lambda_{n,m} \cong H^\bullet(G(n,m))$$
which sends Schur polynomials to Schubert cycles up to scalar since both are weight bases and the weights are multiplicity-free. To see that this is in fact an isomorphism of rings, we can use the action of $\mathfrak{h}_m^-$. We already know that any element of $\mathfrak{h}_m^-$ act on $\Lambda_{n,m}$ by multiplication by some element of $\Lambda_{n,m}$. The corresponding statement is also true for $H^\bullet(G(n,m))$ (this follows from identifying $\mathfrak{h}^-_m$ with the primitive elements in $H^\bullet(Gr)$, as explained in Ginzburg's paper Perverse sheaves on a loop
group and Langlands’ duality). Therefore, we have a commutative diagram
where the two vertical arrows are given by the action on $1$ and are therefore ring homomorphism. Since the left vertical arrow is surjective (as power-sum symmetric functions generate $\Lambda_n$), it follows that $\Lambda_{n,m}\to H^\bullet(G(n,m))$ is itself a ring homomorphisms.
Showing that the scalar by which Schur polynomials and Schubert cycles differ under this isomorphism is actually $1$ is a bit subtle. By working over $\mathbb{Z}$ we can see that it has to be $\pm 1$, and then one can show that it is $1$ by using the structure constants for both bases are nonnegative.
This essentially explains the connection between your points 3 and 4, and we can also easily relate this story to 5 by noting that the wedge representation admits a unique (up to scalar) nondegenerate bilinear form satisfying $(Mv,w)+(v,-M^Tw)=0$ for $M\in \mathfrak{gl}_m$, for which the weight basis is orthonormal (and is the unique orthonormal basis up to signs if we work over $\mathbb{Z}$). The relationship with 2 is essentially Weyl's character formula.
