Dirichlet/Chebotarev's Theorem for natural/analytic density
First, here is some perspective on why most accounts only treat Dirichlet density: the difference between proving a theorem about primes using Dirichlet density or natural density is essentially about the difference between proving $L$-functions are nonzero at $s = 1$ or on the whole line ${\rm Re}(s) = 1$. Dirichlet density is concerned only with limiting behavior as $s \rightarrow 1^+$, while natural density is (not obviously!) related to behavior on the line ${\rm Re}(s) = 1$ because of results like the Wiener-Ikehara Tauberian theorem in its original form or in the simplified form later found by D. H. Newman. When you know Dirichlet $L$-functions $L(s,\chi)$ are nonvanishing on ${\rm Re}(s) = 1$, you can conclude that the logarithmic derivatives $L'(s,\chi)/L(s,\chi)$ are holomorphic on this line (except for a simple pole at $s = 1$ when $\chi$ is trivial) and feeding a suitable linear combination of such logarithmic derivatives into the Wiener-Ikehara theorem or Newman's version of it lets natural density theorems on primes or prime ideals drop out. In the simplest case, a proof of the prime number theorem $\pi(x) \sim x/\log x$, methods using complex analysis involve showing $\zeta(s) \not= 0$ on the whole line ${\rm Re}(s) = 1$ excluding the pole at $s = 1$.
To address your reference questions, the natural density version of Dirichlet's theorem can be found as Proposition 6 and Theorem 1 on pp. 148-150 of Geometric and Analytic Number Theory by Hlawka, Schoißengeier, and Taschner. Lemma 2 on pp. 145-146 shows Dirichlet $L$-functions are nonvanishing on the line ${\rm Re}(s) = 1$.
As for the Chebotarev density theorem, let $L/K$ be a Galois extension of number fields, $G = {\rm Gal}(L/K)$, and $C$ be a conjugacy class in $G$. Picking $g \in C$, the Chebotarev density theorem for $C$ in $G$ follows from the Chebotarev density theorem for $g$ in the cyclic Galois group ${\rm Gal}(L/L^g)$, either using Dirichlet density for both cases or natural density for both cases. (Note $g$ generates the group ${\rm Gal}(L/L^g)$.) This reduction step, which requires a change in ground field from $K$ to $L^g$, is due to Deuring and is used in all current proofs of the theorem because of the way it simplifies what has to be done: you are reduced to the case of a cyclic extension. (Chebotarev did not change the ground field in his own proof, and in fact he only proved his theorem for Galois extensions of $\mathbf Q$. This is evident from the title of his thesis: "Determination of the density of the set of prime numbers belonging to a given class of substitutions." If you can read Russian or German, the thesis is available at mathnet.ru here or the Göttingen digital archive here) At this point there is not much difference between working with the cyclic extension $L/L^g$ and a general finite abelian extension $L/M$, where Chebotarev's theroem says: for each $\sigma \in {\rm Gal}(L/M)$, the natural density of the set of prime ideals in $M$ unramified in $L$ with Frobenius element $\sigma$ is $1/[L:M]$.
Just as Dirichlet's theorem (with Dirichlet density or natural density) is proved using $L$-functions of characters of $(\mathbf Z/m\mathbf Z)^\times \cong {\rm Gal}(\mathbf Q(\zeta_m)/\mathbf Q)$, proving Chebotarev's density theorem (with Dirichlet density or natural density) for the abelian extension $L/M$ uses $L$-functions of characters of $H = {\rm Gal}(L/M)$. These $L$-functions are at first defined and analytic in the open half-plane ${\rm Re}(s) > 1$, and you need to prove properties analogous to what is done with Dirichlet $L$-functions for Dirichlet's theorem.
First suppose $L/M$ is a cyclotomic extension, say $L = M(\zeta_m)$. For each character $\chi$ of ${\rm Gal}(M(\zeta_m)/M)$, we want to show $L(s,\chi)$ extends analytically to ${\rm Re}(s) \geq 1$ except for a simple pole at $s = 1$ for the trivial character. For nontrivial Dirichlet characters, a simple technique shows their $L$-functions converge for ${\rm Re}(s) > 0$. In the case of a nontrivial character $\chi$ of ${\rm Gal}(M(\zeta_m)/M)$, a harder argument shows that when the series $L(s,\chi)$ as a sum over ideals of $\mathcal O_M$ is rewritten as a series over positive integers $\sum_{n \geq 1} a_n/n^s$ , it converges for ${\rm Re}(s) > 1 - 1/[M(\zeta_m):M]$ and thus is analytic there, so in particular on the line ${\rm Re}(s) = 1$. (This step involves viewing ${\rm Gal}(M(\zeta_m)/M)$ as a quotient group of a generalized ideal class group, and in particular uses surjectivity of the Artin map so that a nontrivial character of that cyclotomic Galois group can be viewed as a nontrivial character of a generalized ideal class group of $M$.) The same type of argument made with the trivial character shows the zeta-function of every number field $F$, initially defined for ${\rm Re}(s) > 1$, extends analytically to ${\rm Re}(s) > 1 - 1/[F:\mathbf Q]$ except for a simple pole at $s = 1$.
To prove a nonvanishing theorem for the $L$-function of $\chi$, use the argument found in Serre's Course in Arithmetic in its account of Dirichlet's theorem: show the product of $L(s,\chi)$ over all characters $\chi$ of ${\rm Gal}(M(\zeta_m)/M)$ is $\zeta_{M(\zeta_m)}(s)$ up to finitely many Euler factors $1/(1 - 1/{\rm N}(\mathfrak p)^s)$, which by their very definition are all nonvanishing on ${\rm Re}(s) > 0$, so the product of $L(s,\chi)$ over all nontrivial characters of ${\rm Gal}(M(\zeta_m)/M)$ is, up to finitely many Euler factors, $\zeta_{M(\zeta_m)}(s)/\zeta_M(s)$. Since zeta-functions of number fields all have a simple pole at $s = 1$, that ratio is holomorphic and nonzero at $s = 1$, so $L(s,\chi)$ for nontrivial $\chi$ are all nonzero at $s = 1$ because we know they are all holomorphic at $s = 1$. This is enough to get Chebotarev's theorem for $M(\zeta_m)/M$ using Dirichlet density by the same argument Dirichlet used for his own theorem (the case $M = \mathbf Q$). It's also true that all zeta-functions of number fields are nonvanishing on ${\rm Re}(s) = 1$, so $L(s,\chi)$ for nontrivial $\chi$ are all nonzero on this whole line because we know they are holomorphic on this line. This is enough to get Chebotarev's theorem for $M(\zeta_m)/M$ using natural density by the same argument used to prove Dirichlet's theorem with natural density.
Suppose $L/M$ is a general abelian extension. A deduction of this case from the cyclotomic case using Dirichlet density for both or natural density for both, is in the appendix (p. 18) of Lenstra and Stevenhagen's article "Chebotarev and his density theorem," which is on Lenstra's webpage here. Lenstra and Stevenhagen work with Dirichlet density (see p. 9), but the reasoning in their appendix works in exactly the same way using natural density if you have proved the cyclotomic case with natural density. Just follow their instructions to use $\varliminf$ and $\varlimsup$ in the definition of natural density instead of Dirichlet density.
This only gives the Chebotarev density theorem as an asymptotic statement, not with an error term. For most applications, just as with Dirichlet's theorem, the infinitude of the set of prime ideals is all that's needed, so the asymptotic statement without an error term is sufficient.
You can find a proof in Chapter VIII of Cassels-Fröhlich (eds.): Algebraic number theory (Academic Press, 1967). The chapter was written by Heilbronn, and his proof uses the Wiener-Ikehara tauberian theorem. See Theorem 4 there.
See also Lagarias-Odlyzko: Effective versions of the Chebotarev density theorem, In: Algebraic number fields: L-functions and Galois properties, pp. 409–464, Academic Press, London, 1977.
There is a nice treatment of this and discussions of the error term in Serre's book "Lectures on $N_X(p)$". He also gives various applications but I don't think he proves it.
I should note that the natural density version holds over number fields but not over global function fields. Serre discusses this in his book.