I know that the topic is classical and even "folklore", but many treatments make use of local coordinates and such treatments are rather messy. Could somewhere maybe provide some reference(s) to a modern treatment, using the language of modern differential geometry (with complex line bundles, connections and so on), of magnetic monopoles including a complete description, which is as explicit as possible without being too messy, of the basic Dirac monopole? If the holomorphic structure of $S^2$, which turns it into $\mathbb{C}P^1$, plays a role, I would also like that to be mentioned and made clear, if possible. Of course a combination of references would be ok too. If someone would like to describe it in an answer, that would be great of course!

I suspect I may find that somewhere in the various Cambridge lecture notes online (maybe in one of D. Tong's lecture notes), or in some articles/books by N. Manton for example. I remember it was very briefly mentioned in the Atiyah and Hitchin book on magnetic monopoles, but I kind of would like more details please. I will dig in the literature, but I suspect such references would be generally useful to others too, which is why I thought about writing this post. I will edit and report on what I find too.

Edit 1: first, there is a lot of information in Tong's lecture notes (for example, his gauge theory notes) on Dirac monopoles. What follows is a very short description which I found in the article by Gibbons and Manton called "The Moduli Space Metric for Well-separated BPS Monopoles".

A Dirac vector potential, which really means a connection on a complex line bundle, in this case defined on $\mathbb{R}^3 \setminus \{ \mathbf{0} \}$, is one which satisfies:

$$ \nabla \times \mathbf{w} = - \frac{\mathbf{r}}{r^3} $$ and $$ \mathbf{w}(-\mathbf{r}) = \mathbf{w}(\mathbf{r}). $$

I understand the first equation, which is really the Bogomolny equation. Indeed, we can think of it as $F = * d \phi$, where $F$ is the curvature of the connection $\mathbf{w}$, $*$ denotes the Hodge star with respect to the (flat) Euclidean metric on $\mathbb{R}^3 \setminus \{ \mathbf{0} \}$ and $\phi$ is a scalar potential (physically a Higgs field). In this case, $\phi = 1 / r$.

Could someone perhaps comment on the second equation please? How can you compare the values of a connection at $2$ different points? Is that using the fact that we are working inside a domain of $\mathbb{R}^3$, whose tangent bundle is naturally trivial? I guess this is what they mean please?

Edit 2: I should add that I have just stumbled on a book by Yakov Shnir entitled "Magnetic Monopoles" which seems to contain a lot of relevant material (including on multimonopoles). It definitely seems "modern", judging from the preview provided by Google.

I also found a nice detailed treatment of Dirac monopoles in the book by N. Manton and P. Sutcliffe entitled "Topological Solitons", in section 8.1. I think I will definitely benefit from reading that section. I do find the work of P. Dirac on one hand beautiful mathematically, yet difficult to interpret physically, and often his formulas have many possible interpretations. I guess I am mostly thinking about the Dirac equation, but his monopoles also have their own subtleties too!

Update: I found section 8.1 in the book by N. Manton and P. Sutcliffe called "Topological Solitons" to be exactly what I wanted. Not only is the language "modern", but they also answer many questions I had too. For example, why can't we just modify Maxwell's equations so that we have that the divergence of $B$ is (up to some factor) the magnetic charge density? Well, such a naive way of introducing magnetic monopoles has several issues which they discuss.

Dirac's work was much more subtle. He essentially notices that he could introduce a magnetic monopole mathematically using a connection on a $U(1)$ bundle over $\mathbb{R}^3 \setminus \{ \mathbf{0} \}$. One may use 2 patches to cover that region, say, using spherical coordinates, with $\theta \neq 0$ being one region and $\theta \neq \pi$ being another. In these coordinates, it looks as though the magnetic potential has a singular ray in each of the two patches (the so called Dirac strings), but these are just fake coordinate singularities. The connection only has a singularity at the origin of $\mathbb{R}^3$. The integrality of the first Chern class of the $U(1)$ bundle leads to a quantization condition on the electric/magnetic charges. Dirac (or "Monopoleon", as W. Pauli used to call him, I think) was brilliant at finding subtle loopholes in the equations which we use to describe nature.

Whether they exist or not in our universe remains to be seen (as of today's date). But their existence does not lead to contradictions with Maxwell's equations.

Finally, I wonder, purely speculatively (and I am not a physicist, so it is most probably false), whether there are black hole solutions of Einstein's equations, which are also magnetic monopoles, such that the singularity of the monopole coincides with the singularity of the black hole. Do those exist (at least mathematically)? If not, is there any fundamental reason why they cannot exist, please (mathematically and/or physically)?