In which mathematical space do magnetic fields live in? A magnetic field may be described as a mapping $\mathbb{R}^3 \to \mathbb{R}^3$. Therefore, it is a function, and apparently can be described as an element in a function space. Is there some particular (perhaps Hilbert) space that magnetic fields belong to?
 A: $\def\vA{{\vec A}}\def\vB{\vec{B}}\def\vH{{\vec H}}\def\vS{{\vec S}}\def\vpsi{{\vec\psi}}\def\div{{\rm div}}\def\grad{{\rm grad}}\def\rot{{\rm rot}}\def\rmr{{\rm r}}\def\pd{\partial}\def\nR{{\mathbb{R}}}\def\ltag#1{\tag{#1}\label{#1}}$
It strongly depends on what you want to model. If you want to include dipoles in your model you must switch to distribution spaces. But, let us keep in a subspace of $L^2(\Omega)$.
The space $L^2(\Omega)$ is too big to give equations like
\begin{align}
\div\vB &= 0.
\end{align}
sense. On the other hand one often has discontinuities of $\mu_\rmr$ at some boundary layers. There the tangent component of $\vB$ can be discontinuous and $\div$ in the classical context does not make sense.
Therefore one must switch to generalized derivatives. The integration by parts rule
\begin{align}
\int_{\Omega} \psi\cdot \div(\vB) d V &= \int_{\pd\Omega} \psi\; \vB\cdot d\vA
- \int_{\Omega} \grad(\psi) \cdot \vB dV\ltag{weakDiff}
\end{align}
is used to transfer the differentation from $\vB$ to test functions $\psi$ and one uses the right hand side to define the generalized $\div$ operator.
If there is a $L^2$ function "$b$" such that
\begin{align}
\int_{\Omega} \psi b\; d V &= \int_{\pd\Omega} \psi\; \vB\cdot d\vA
- \int_{\Omega} \grad(\psi) \cdot \vB dV
\end{align}
for all smooth $\psi$ with bounded support in $\bar\Omega$ then we define the generalized divergence operator for $\vB$ as $\div\vB:=b$ and we say that for $\vB$ there exists the weak divergence on $\Omega$. The space $H(\div,\Omega)$ is the set of all $L^2$-functions on $\Omega$ for which there exists the weak divergence.
Note, that the surface integral in \ref{weakDiff} must be treated in special manner, since the surface is actually a set of zero measure in $\nR^3$. If we would only use $L^2(\bar \Omega)$ then we could modify $\vB$ on $\partial\Omega$ arbitrarily. Therefore, the $\vB$-field on the boundary is actually to be understood as the trace of the field on the boundary. For the trace to exist there are additional restrictions on the boundary of $\Omega$ (piecewise $C^1$, edges not too "sharp"). Thus, not all domains $\Omega$ are really appropriate for modeling of magnetic fields.

Note, that the situation is similar for Ampere's law
\begin{align}
\rot\vH = \vS
\end{align}
where we need $\rot\vH$. The normal component of the field strength can be discontinuous at boundary layers. In this case one uses the equation
\begin{align}
\int_{\Omega}\vpsi\cdot \rot(\vH)dV = \int_{\pd\Omega}\vpsi\cdot(d\vA\times\vH) + \int_{\Omega} \vH\cdot \rot(\vpsi) dV
\end{align}
to define a generalized curl operator $\rot\vH$. The set of fields for which $\rot\vH\in L^2(\Omega)$ is named $H(\rot,\Omega)$.
Interestingly, $\vB$ and $\vH$ live in different spaces and $\mu$ must mediate between these spaces $\mu:\vH\in H(\rot,\Omega)\mapsto \vB=(\mu\vH)\in H(\div,\Omega)$.
