I will only address the question to know why $\mathbb{R}^{4}$ admits some non-standard, "exotic", differentiable structure. The question to know why there are infinitely many (and even uncountably many) examples simply requires an extension of the same kind of techniques.

There are several ways to construct examples of exotic $\mathbb{R}^{4}$. All use some deep results of topological nature due to Freedman and some deep results of differentiable nature due to Donaldson. I don't know if the results of Freedman have any physical interpretation. The Yang-Mills theory appears in Donaldson's results, on the differentiable side (one needs a differentiable structure to write partial differential equations).

Here is a sketch of one of the standard construction. The existence of exotic $\mathbb{R}^{4}$ will be a consequence of results on compact four manifolds. Let $M$ be an oriented compact differentiable manifold. The most important invariant of $M$ is its intersection form on $H^{2}(M, \mathbb{Z})$, it is an integral quadratic form unimodular by Poincaré duality (i.e. $H^{2}(M, \mathbb{Z})$ is a self-dual lattice).

Examples: $M = \mathbb{CP}^{2}$, $H^{2}$ is of dimension 1, intersection form $x^{2}$, of signature (+). $M = S^{2} \times S^{2} = \mathbb{CP}^{1} \times \mathbb{CP}^{1}$, $H^{2}$ is of dimension 2, intersection form $xy$ of signature (+-), we have a standard hyperbolic lattice $H$.

Real quadratic forms are classified by their signature but integral quadratic forms are more complicated: they can have the same signature but still not to be equivalent over the integers. Example: $x^{2}-y^{2}$ and $xy$ are two quadratic forms of signature (+-) but there are not equivalent over the integers (one needs to divide by 2 to find a change of coordinates passing from one to the other). Other example: the quadratic form $E_{8}$ is unimodular of signature (++++++++) but is not isomorphic over the integers to the standard (=sum of squares) quadratic form of signature (++++++++).

A basic question is; which integral unimodular quadratic form can appear as intersection form of a oriented compact 4-manifold? It happens that the answer is very different if one considers topological manifolds or differentiable manifolds. Roughly, Freedman says that in the topological case everything can be realized and Donaldson says that in the differentiable case there are some obstructions.

A precise result of Freedman relevant for the exotic $\mathbb{R}^{4}$ is the following, if the intersection form of $M$ is a direct sum of quadratic forms, then this decomposition can be realized at the level of topology, $M$ can be written has a connected sum of topological manifolds corresponding to the decomposition of the quadratic form.

Example: let $M$ be a K3 surface, it is a differentiable 4-manifold. Its intersection form is of signature (3+,19-) and is in fact isomorphic to the direct sum of minus two copies of $E_{8}$ and of 3 hyperbolic planes $H$. By a refinement of Freedman's result, one can write $M$ as a topological connected sum of a topological manifold $N$ of intersection form minus two copies of $E_{8}$ and of three copies of the differential manifold $S^{2} \times S^{2}$.

A precise result of Donaldson relevant for the exotic $\mathbb{R}^{4}$ is : let $M$ be a differentiable compact oriented manifold of intersection form positive definite. Then this intersection form is the standard one. In particular, $E_{8}$ or the sum of two copies of $E_{8}$ are forbidden in the differentiable world. Ampplying this to the K3 example shows that $M$ is not a differentiable manifold. As the total K3 is a differentiable manifold, the 3 copies of $S^{2} \times S^{2}$ are differentiable manifolds, the "neck" defining the connected some is an exotic differentiable structure on $S^{3} \times \mathbb{R}$. With some work, we can pass from $S^{3} \times \mathbb{R}$ to $\mathbb{R}^{4}$.

The statement of the Donaldson's result seems to have nothing to do with Yang-Mills theory. Yet, the proof uses it. The key point is that because the manifold is differentiable we can do some analysis on it, in particular partial differential equations as the Yang-Mills equations (equations of motion of the classical Yang-Mills theory). In four dimensions, it is possible to construct solutions of Yang-Mills equations, which are of the second order, starting from a first order equation called the ASD (=anti self-dual) equation, which is of first order (the point is that in four dimensions, the 2-form curvature can be decompose in a self-dual part and a anti-self-dual part. The ASD equation is simply self-dual part =0). A solution of this equation is called an instanton (such particular solutions of Yang-Mills equations extremize a natural BPS like bound).

Now the idea of Donaldson is to use the moduli spaces of instantons to study the 4-manifold $M$. More precisely, he takes for gauge group $SU(2)$ and considers instantons of instanton number 1 (the instanton number is a topological invariant of the instanton, roughly minus the second Chern class of the corresponding $SU(2)$-principal bundle). Such instanton is a configuration of the gauge field which can be more or less localized on $M$. A relatively localized instanton looks like an extended particles, with 4 parameters to describe its position and 1 parameter to describe its size (an instanton of instanton number k would generically look like k extended particles). So the corresponding moduli space is of dimension 5. This space is not compact because an instanton can shrink to a zero size point. As this point can be any point of $M$, $M$ is in fact a natural boundary of (a compactification of) the moduli space of instantons. So if one studies the moduli space, one should be able to obtain informations on $M$.

The moduli space is smooth everywhere except at some points where the ASD connection is reducible, i.e. when the $SU(2)$ principle bundle $E$ decomposes as $L \oplus L^{-1}$ with $L$ a complex line bundle. As $c_{1}(L)^{2}=- c_{2}(E)=1$, this happens precisely in $m$ points where $m$ is the number of pure squares in the intersection form of M. As by hypothesis, the intersection form of M is positive definite, it is standard if and only $m$ is equal to the rank $n$ of the $H^{2}$. A local study shows that near a singularity the moduli space is diffeomorphic to $\mathbb{CP}^{3}/S^{1}$ (there is a $S^{1}$ group of automorphisms) which is a cone over $\mathbb{CP}^{2}$. Cutting the ends of these cones containing the singular points of the instanton moduli space, we obtain a smooth manifold of dimension 5 with boundary $M$ and the union of $m$ copies of $\mathbb{CP}^{2}$, i.e. we have constructed a cobordism between $M$ and $m$ copies of $\mathbb{CP}^{2}$. But it is known that the signature intersection form of four manifolds is invariant under cobordism. As the signature of $M$ is $n$ and the signature of the union of $m$ copies of $\mathbb{CP}^{2}$ is $m$, we conclude that $n=m$.