Density matrix formalism The density matrix $\hat{\rho}$ is often introduced in textbooks as a mathematical convenience that allows us to describe quantum systems in which there is some level of missing information.
$\hat{\rho} = \sum_{i=1}^N c_i \rvert\psi_i\rangle \langle \psi_i \rvert$
I have two questions regarding density matrices.
First, it is clear that a generic pure state $\rvert \psi \rangle$ belongs to some Hilbert space $\mathcal{H}$. But what mathematical space do density matrices belong to? It is clear that the expression $\hat{\rho} \in  \mathcal{H}$ is incorrect, as a mixed state cannot possibly be described in terms of a state in a Hilbert space. Can we think of the space of all possible density matrices (of a given dimension) as a metric space? Does this space have the topological properties of a manifold?
Secondly, can the density matrix be considered 'physical'? For example, if we take a single photon described in the Fock basis (neglecting polarization), can the fundamental description of that photon ever be the completely mixed state $\frac{1}{2}(\rvert 0 \rangle \langle 0 \rvert+\rvert 1 \rangle \langle 1 \rvert)$? Or is this only a reflection of some ignorance on the part of the experimentalist, and in reality the photon must be described by a pure state?
 A: I'd just like to add a point about the last question:

Secondly, can the density matrix be considered 'physical'? For example, if we take a single photon described in the Fock basis (neglecting polarization), can the fundamental description of that photon ever be the completely mixed state $\frac{1}{2}(\lvert 0\rangle \langle 0 \lvert + \lvert1\rangle \langle 1 \lvert)$? Or is this only a reflection of some ignorance on the part of the experimentalist, and in reality the photon must be described by a pure state?

If one assumes that the photon is not correlated with any other system, then a mixed state can only arise as a consequence of experimental ignorance. However, there are many situations in which a mixed state unavoidably arises due to the fundamental quantum uncertainty. In particular, if the photon (or any quantum system) is entangled with another system, then the state of the photon alone is necessarily mixed. For example, consider a pair of light modes in the state
$$ \lvert \psi\rangle = \frac{1}{\sqrt{2}}\left ( \lvert 0 \rangle \lvert 0 \rangle + \lvert 1 \rangle \lvert 1 \rangle \right ), $$
If you can only make measurements on one of the modes, then your outcomes will be equivalent to those obtained from the mixed state
$$\rho = \frac{1}{2}(\lvert 0\rangle \langle 0 \lvert + \lvert1\rangle \langle 1 \lvert)$$
which gives totally random outcomes. The fluctuations in the measurement results come from the fact that the global system has been prepared in a pure quantum state $\psi$ which does not have a definite number of photons in either mode individually. This uncertainty is an unavoidable quantum effect and occurs even when the observer has the maximum possible information about the individual mode in question. In this sense the density matrix is a more complete description of physical systems than pure states alone.
When the observer has access to both modes then it is possible to devise measurements which have definite outcomes, since the global system is in a pure state. For example, a measurement which checks whether both modes always have the same number of photons in each individual experimental run will always yield the answer "yes" with 100% certainty. 
A: A mixed state is mathematically represented by a bounded, positive trace-class operator with unit trace $\rho : \cal H \to \cal H$. 
Here $\cal H$ denotes the complex Hilbert space of the system (it may be nonseparable). The set of mixed states $S(\cal H)$ is a convex body in the complex linear space of trace class operators $B_1(\cal H)$ which is a two-side $*$-ideal of the $C^*$-algebra of bounded operators $B(\cal H)$. 
Convex means that if $\rho_1,\rho_2 \in S(\cal H)$ then a convex combination of them, i.e. $p\rho_1 + q\rho_2$ if $p,q\in [0,1]$ with $p+q=1$, satisfies $p\rho_1 + q\rho_2 \in S(\cal H)$. 
two-side $*$-ideal means that linear combinations of elements of $B_1(\cal H)$
belong to that space (the set is a subspace), the adjoint of an element of $B_1(\cal H)$ stays in that space as well and $AB, BA \in B_1(\cal H)$ if $A\in B_1(\cal H)$ and $B \in B(\cal H)$.
I stress that, instead, the subset of states $S(\cal H)\subset B_1(\cal H)$ is not a vector space since only convex combinations are allowed therein.
The extremal elements of $S(\cal H)$, namely the elements which cannot be decomposed as a nontrivial convex combinations of other elements, are all of the pure states. They are of the form $|\psi \rangle \langle \psi|$ for some unit vector of $\cal H$. (Notice that, since phases are physically irrelevant the operators $|\psi \rangle \langle \psi|$ biunivocally determine the pure states, i.e. $|\psi\rangle$ up to a phase.)
The space $B_1(\cal H)$ and thus the set $S(\cal H)$ admits at least three relevant normed topologies induced by corresponding norms. One is the standard operator norm $||T||= \sup_{||x||=1}||Tx||$ and the remaining ones are:
$$||T||_1 = || \sqrt{T^*T} ||\qquad \mbox{the trace norm}$$
$$||T||_2 = \sqrt{||T^*T||} \qquad \mbox{the Hilbert-Schmidt norm}\:.$$
It is possible to prove that:
$$||T|| \leq ||T||_2 \leq ||T||_1 \quad \mbox{if $T\in B_1(\cal H)$.}$$
Moreover, it turns out that  $B_1(\cal H)$ is a Banach space with respect to $||\cdot||_1$ (it is not closed with respect the other two topologies, in particular, the closure with respect to $||\cdot||$ coincides to the ideal of compact operators $B_\infty(\cal H)$). 
$S(\cal H)$ is closed with respect to $||\cdot ||_1$ and, more strongly, it is a complete metric space with respect to the distance $d_1(\rho,\rho'):= ||\rho-\rho'||_1$.
When $dim(\cal H)$ is finite the three topologies coincide (though the norms do not), as a general result on finite dimensional Banach spaces.
Concerning your last question, there are many viewpoints. My opinion is that a density matrix is physical exactly as pure states are. It is disputable whether or not a mixed state encompasses a sort of physical ignorance, since there is no way to distinguish between "classical probability"  and "quantum probability" in a quantum mixture as soon as the mixture is created. See my question Classical and quantum probabilities in density matrices and, in particular Luboš Motl's answer.
See also my answer to Why is the application of probability in QM fundamentally different from application of probability in other areas?
ADDENDUM. In finite dimension, barring the trivial case $dim({\cal H})=2$ where the structure of the space of the states is pictured by the Poincaré-Bloch ball as a manifold with boundary, $S(\cal H)$ has a structure which generalizes that of a manifold with boundary. A stratified space. Roughly speaking, it is not a manifold but is the union of (Riemannian) manifolds with different dimension (depending on the range of the operators) and the  intersections are not smooth. When the dimension of $\cal H$ is infinite, one should deal with the notion of infinite dimensional manifold and things become much more complicated.
A: Nobody seems to have answered directly one of the question so: 
$\rho$ lives in
${\cal H}\otimes {\cal H}^*$, with ${\cal H}^*$ the dual of ${\cal H}$.  
Alternatively, if kets live in ${\cal H}$ then bras live in ${\cal H}^*$ so that $\rho$ is an object that contains linear combinations of the form
$\vert \psi_j\rangle\langle \phi_k\vert$; typically $\vert \psi_j\rangle \in {\cal H}$ and $\langle \phi_k\vert\in {\cal H}^*$.  
Not every linear combination of the type $\vert \psi_j\rangle\langle \phi_k\vert$ can enter in $\rho$ so $\rho$ describes a physical state.  In a fixed orthonormal basis, the diagonal entries of $\rho$ are classical probabilities and so must be real non-negative and sum to $1$, while $\rho$ itself must be hermitian.
A: The density matrix is an operator that picks up ''classical indetermination'', I mean, indetermination coming from the experiment. For example, you may have a population of spin $\frac 1 2$ particles, and you know that a half of this population is in state $S_Z = \frac 1 2$ and the other half is in $S_z = -\frac{1}{2}$. The corresponding $\hat{\rho}$ is
$$\hat{\rho} =  \frac 1 2 (|+\rangle\langle+| +|-\rangle\langle-|) $$
This not means that the state of each particle is
$$|\psi\rangle = \frac{1}{\sqrt{2}}(|+\rangle + |-\rangle)$$
This indetermination is not a quantum indetermination, because you know that a particle is in $\frac 1 2$ state or $\frac 1 2$ . That's why density matrix is useful: it includes effects of 'classical' indetermination.
The density matrix is not a 'genuine' operator since its temporal evolution doesn't follow Heisenberg equation. It satisfies Von Neumann's equation.
In your example, the density operator represents that with 50% of probability, there are 0 photons (that is, the system is in the ground state $|0\rangle$) or 1 photon (state $|1\rangle$). This means that one of the states $|0\rangle$ or $|1\rangle$ has been prepared before the experiment.
A: 
But what mathematical space do density matrices belong to? It is clear that the expression $\hat{ρ}\in\mathcal{H}$ is incorrect, as a mixed state cannot possibly be described in terms of a state in a Hilbert space.

A bounded linear operator on $\mathcal{H}$ is of trace class iff it has finite trace independent of choice of basis (cf. also nuclear operator). A density matrix is a trace class positive linear operator $\rho:\mathcal{H}\to\mathcal{H}$ with trace $1$. I will assume that $\mathcal{H}$ is separable.
The eigenfunction decomposition theorem (Hilbert-Schmidt theorem) implies that it has the form $\rho = \sum_k p_k|\psi_k\rangle\langle\psi_k|$, where $p_k>0$ are the nonzero eigenvalues and $\sum_k p_k = 1$, so in that sense can be described in terms of the pure states in Hilbert space. Actually, we could even define the terminology like so: a density matrix is a "state", and in the particular case that it is a projection operator to a linear subspace of $\mathcal{H}$, it is a "pure state". So we can equivalently consider pure states as special cases of density matrices, rather than the usual way of considering density matrices a generalization of pure states.

Can we think of the space of all possible density matrices (of a given dimension) as a metric space?

Yes, in fact they form a normed vector space with the trace norm / nuclear norm. Moreover, there's a closely related notion of Hilbert-Schmidt operators that form a Hilbert space via the trace inner product.

Does this space have the topological properties of a manifold?

I don't know of interesting results along those lines for density matrices, but I think I've seen a subclass of Hilbert-Schmidt operators being given a Riemannian manifold structure based on the trace inner product, but I'm not aware of the details.

Secondly, can the density matrix be considered 'physical'? For example, if we take a single photon described in the Fock basis (neglecting polarization), can the fundamental description of that photon ever be the completely mixed state... Or is this only a reflection of some ignorance on the part of the experimentalist, and in reality the photon must be described by a pure state?

A density matrix is a reflection of ignorance. I don't think this makes it 'unphysical'. Whether they're 'non-fundamental' depends on what you mean by 'fundamental', esp. in light of the fact that we can consider pure states as just special cases of density matrices.
