Sinusoidal currents in AC circuits (I had edited this question a number of times but did not receive a satisfactory answer. So I am re-wording this question yet again. The earlier version of this question can be found here  : https://docs.google.com/document/d/1_Z6iAc_LwoTxCrevoSqpoJw5TDHvVqLuDVXktPOu2_w/ )
My whole question consists of many parts (I have added a few additional doubts) :


*

*Justifying (without skipping any logic) that current in any branch of a linear AC Circuit will be of the form $I(t)=I_0 \cos(\omega t + \phi)$. Start by writing a general set of differential equations, $n$ equations for $n$ loops and then prove that the solution will always be just a single sinusoidal term. 

*Now that the sinusoidal nature of currents is properly established and justified, we can represent the sinusoidal functions by complex numbers. But the actual current flowing in any branch will be the real part of the complex current. The imaginary component has no important role. It is just there. Now that the current is of the form $I_0 e^{i (\omega t + \phi)}$, we can use the formulas $V= IR$, $V= L \frac {dI}{dt}$  and $V= \frac Q C$ to find voltages across resistors, inductors and capacitors respectively. But why do we use $I=I_0 e^{i (\omega t + \phi)}$ for calculating voltages across inductors and capacitors? Shouldn't we only use $I=Re(I_0 e^{i (\omega t + \phi)})$ because only the real part of the complex current is the actual current flowing in the branch? If we use the real as well as the imaginary part of find the voltage across an inductor, then during differentiation the iota operator will "come down" and will be multiplied with the complex number and so the earlier-real part will become imaginary and the earlier-imaginary part will become real. Why is it okay to use complex current to find the voltage? 

*Also, while calculating power we use $S= VI^*$, where $I^*$ conjugate of complex current. Why is it not $S=VI$ or $S=Re(V)Re(I)$ ?

*Any sinusoidal quantity associated with the circuit will be of the form of $Q=Q_0 e^{i(\omega t + \phi)}$. But  when we apply Kirchhoff's Voltage Law we just use $Q=Q \angle {\phi}= Qe^{i \phi}$. Where did the $\omega t$ part in the complex exponential go ?
Thanks in advance :-).
EDIT 1:  All of my doubts will be solved if you can point me to a reference book or something which explains everything about the basics of AC circuits in proper mathematical detail with proper justification of each and every step, calculation and assumption. Does such a book even exist ?
EDIT 2: Still no answers !
 A: I'm not sure if I truly understand why you're unsure that all the branch currents should be sinusoidal in steady state. 
First, consider the simple case of just one source and linear circuit elements (and, as usual, stipulate that these linear circuit elements do not vary with time). Each branch current can be written as the output of a linear time-invariant (LTI) system with the source as the input. If you're not in agreement with this, then let me know in a comment.


*

*Now, a fundamental property of a linear time-invariant system is that there no
frequencies present in the output that are not present in the input.
Further, if there are frequencies present in the output that are not
present in the input, the system is not linear.


(If you're unsure why this is so, you might find this question and answers helpful:  Why do linear systems show sinusoidal fidelity?)
If the input to a LTI system is sinusoidal, there is just one frequency present in the input and thus, there is just one frequency present in the output of the system.
So, if the source driving a linear circuit is sinusoidal with angular frequency $\omega$, each branch current is sinusoidal with angular frequency $\omega$.
When there are two or more sources present, linearity implies that each branch current can be written as a sum of responses, one from each of the sources (Superposition theorem).
It follows that if there are two or more sinusoidal sources present, each with angular frequency $\omega$, each branch current can be expressed as a sum of sinusoids with angular frequency $\omega$.
A: This 

solutions are sinusoidal in nature because the differential equations are linear

is a fairly inaccurate capture of the true gist in which these choices are taken. A much better rendition is the following:

it is sufficient to handle only solutions which are sinusoidal in nature because, since the differential equations are linear, every other solution can be built from sinusoidal solutions.


It's a bit hard to go into more detail because the set of configurations that you're dealing with is so huge, but let's try to capture the broad-strokes essentials of the circuits you're worried about. In essence, any such system consists of:


*

*some driving term $v(t)$, which might in principle have arbitrary time dependence;

*some quantity of interest $q(t)$; and

*a linear response function that determines how the quantity of interest depends on the driving term, $q(t) = \mathcal L[v(t)]$.


The cornerstone reason why it's sufficient to deal only with sinusoidal waveforms is a theorem from Fourier analysis, which tells us that any arbitrary (sufficiently well-behaved) driving term $v(t)$ can be expressed as a linear superposition of a bunch of sinusoidal functions:
$$
v(t) = \int_{-\infty}^\infty \tilde v(\omega) e^{-i\omega t}\mathrm d\omega.
\tag 1
$$
(Here $\tilde v(\omega)$ is known as the Fourier transform of $v(t)$ and it can be calculated explicitly from $v(t)$ through a formula that's very similar to $(1)$. For now, though, this doesn't really matter ─ the only relevant thing is that such a function exists.)
The second cornerstone of the analysis is that the response function is linear, which allows us to break out the calculation of the full quantity of interest $q(t)$ into a linear combination of responses to linear driving terms:
\begin{align}
q(t)
& = \mathcal L[v(t)]
\\&  = \mathcal L\mathopen{}\left[
\int_{-\infty}^\infty \tilde v(\omega) e^{-i\omega t}\mathrm d\omega
\right] \mathclose{}
  = 
\int_{-\infty}^\infty \tilde v(\omega) 
\mathcal L\mathopen{}\left[e^{-i\omega t}\right] \mathclose{}
\mathrm d\omega
\\&  = 
\int_{-\infty}^\infty \tilde v(\omega) 
\tilde{\mathcal L}_\omega e^{-i\omega t}
\mathrm d\omega
  = 
\int_{-\infty}^\infty \tilde q(\omega) e^{-i\omega t} \mathrm d\omega
\end{align}
Moreover, it will generally be the case that if the driving term of a linear circuit is sinusoidal, the response will also be sinusoidal, so you're left directly with the Fourier transform of your quantity of interest.
A: First, a sinusoid is just an exponential in disguise:
$$\sin\theta = \frac{e^{i\theta}-e^{-i\theta}}{2i}$$
$$\cos\theta = \frac{e^{i\theta}+e^{-i\theta}}{2}$$
Second, the exponentials are eigenfunctions of the derivative and integral operators. That means if you put an exponential function in to a derivative, you get a complex exponential out.
Therefore, since a sinusoid is just a different way of writing an exponential, if you put a sinusoid into a a derivative or integral, you get another sinusoid out.
Since the functions we allow in a linear system are essentially scalar multiplications, derivatives, and integrals, this means that if we have a sinusoidal source driving a linear system, all the response waveforms of that system will also be sinusoids.
And then, as other answers have mentioned, we're able to use Fourier analysis to express any waveform we like in terms of sinusoids, so this analysis is actually useful for any source waveform, not just for pure sinusoids.
A: It is a mathematical fact that adding sinusoids of identical frequency yields another sinusoid of that frequency, regardless of the amplitudes and phases.
