What is the distribution of primes modulo $n$?

Let $n\geq 2$ and let $k$ be "considerably larger" than $n$ (like some large multiple of $n$). Then for each $i$ such that $0<i<n$ and $\gcd(i,n)=1$ let's define $$c_i=\left|\{p_j\;|\; p_j\equiv i \mod n,\;\mbox{where $p_j$ is the $j$-th prime, $1\leq j\leq k$}\}\right|$$ so $c_i$ represents how many of the first $k$ primes are congruent to $i$ modulo $n$.

What can we say about $c_i$s, that is about the distribution of the first $k$ primes modulo $n$?

I thought the distribution will be seemingly random, and that is mostly true - for example $c_i$s are always very close together. But there are observable non-random patterns. For example for $n=3$, for various $k$s I've tried (up to $10^6$) I always got $c_2>c_1$. If this were just some kind of random discrepancy due to distribution of small primes, it would eventually vanish for large $k$s, which I don't observe.

Dirichlet's theorem on primes in arithmetic progression tells us that the proportion of primes will be the same, for values of $i$ that are coprime to $n$ (and 0 otherwise).

However, Chebyshev's bias tells us that numerically there are more primes with a quadratic non-residue, than a residue, when you're counting primes up to $N$.

I once tried to answer the particular case of primes $4n\pm 1$ to myself. (I think it will strongly help to read the answers by Raymond, Raymond and Greg. At the end of the last answer there's also a link to the chat, where we continued the discussion.)

Here is how far I got with an explicit formula for the number of primes of the form $4n+3$ below $x$, $\pi^*(x;4,3)$, expressed in terms of (sums of) sums of Riemann's $R$ functions over roots of Riemann's $\zeta$ resp. Dirichlet $\beta$ function:

\begin{align*} \Pi^*(x;4,3) &= \pi^*(x;4,3) + \tfrac12 \sum_{\substack{b\pmod 4 \\ b^2\equiv 3\pmod 4}} \pi^*(x^{1/2};4,b) + \tfrac13 \sum_{\substack{c\pmod q \\ c^3\equiv 3\pmod 4}} \pi^*(x^{1/3};4,c) + \cdots \\ \end{align*} Then I try to complete things by adding several up

\begin{align*} \Pi^*(x;4,3) &= \tfrac11\pi^*(x;4,3) + \tfrac13 \pi^*(x^{1/3};4,3) + \cdots \\ \tfrac12\Pi^*(x^{1/2};4,3) &= \tfrac12\pi^*(x^{1/2};4,3) + \tfrac16 \pi^*(x^{1/6};4,3) + \cdots \\ \tfrac14\Pi^*(x^{1/4};4,3) &= \tfrac14\pi^*(x^{1/4};4,3) + \tfrac1{12} \pi^*(x^{1/12};4,3) + \cdots \\ &\vdots&\\ \hline\\ \tag{1}\sum_{k=0}^\infty 2^{-k}\Pi^*(x^{2^{-k}};4,3)&=\sum_{m=1}^\infty \tfrac1m \pi^*(x^{1/m};4,3) \end{align*} Using Möbuis inversion I'll get

\begin{align*} \pi^*(x;4,3)&=\sum_{m=1}^\infty \tfrac{\mu(m)}m\sum_{k=0}^\infty 2^{-k}\Pi^*(x^{2^{-k}/m};4,3)\\ \tag{2}&=\sum_{k=0}^\infty 2^{-k}\sum_{m=0}^\infty \tfrac{\mu(m)}m\Pi^*(x^{2^{-k}/m};4,3) \end{align*} Now I use

\begin{align*} \Pi^*(x^{2^{-k}};4,3)&=\frac1{\phi(4)} \sum_{\chi\pmod 4} \overline{\chi(3)}\Pi^*(x^{2^{-k}},\chi)\\ \tag{3}&=\frac12 \left( \Pi^*(x^{2^{-k}},\chi_1)- \Pi^*(x^{2^{-k}},\chi_2) \right) \end{align*} and then

\begin{align*} \tag{$4_1$}\Pi^*(x^{2^{-k}},\chi_k)&=\operatorname{li}(x^{1/2^{k}})-\sum_{\rho_\zeta} \operatorname{li}(x^{\rho_\zeta/2^k})\text{ if $k=1$}\\ \tag{$4_2$}&=\phantom{\operatorname{li}(x^{1/2^{k}})}-\sum_{\rho_\beta} \operatorname{li}(x^{\rho_\beta/2^k})\text{ if $k=2$}\\ \end{align*} which gives

\begin{align*} \tag{3'}\Pi^*(x^{2^{-k}};4,3)&=\frac12 \left( \operatorname{li}(x^{1/2^{k}})-\sum_{\rho_\zeta} \operatorname{li}(x^{\rho_\zeta/2^k}) +\sum_{\rho_\beta} \operatorname{li}(x^{\rho_\beta/2^k}) \right) \end{align*} so finally

\begin{align*} \pi^*(x;4,3)&=\sum_{k=0}^\infty 2^{-k}\sum_{m=0}^\infty \tfrac{\mu(m)}m\frac12 \left( \operatorname{li}(x^{1/2^{k}})-\sum_{\rho_\zeta} \operatorname{li}(x^{\rho_\zeta/2^k}) +\sum_{\rho_\beta} \operatorname{li}(x^{\rho_\beta/2^k}) \right)\\ \tag{5}&=\sum_{k=0}^\infty 2^{-k-1}\left( \operatorname{R}(x^{1/2^{k}})-\sum_{\rho_\zeta} \operatorname{R}(x^{\rho_\zeta/2^k}) +\sum_{\rho_\beta} \operatorname{R}(x^{\rho_\beta/2^k}) \right) \end{align*}

I would be very, very glad to read your opinion...