Why do populations only change in second order of the driving field? In the field of quantum optics when solving master equations it is well know that the populations1 are constants to linear order in the driving field. I.e. weak driving fields will only affect the coherences off the system.
I am interested in how this statement can be made more precise and what the most general system is that it applies too. So far I could only find derivations for example systems (e.g. two level system) and imprecise statements of this notion.
To be a bit more formal: Let there be a general multi-level system with unspecified couplings between the levels and a loss channel from each level. The system is driven by a harmonic field of given frequency and intensity and couples to (some of the) levels via a dipole approximation Hamiltonian. For such a system, under what conditions do the populations1 only change in second order of the driving field and how can this be proven in general?

As an example a typical Master equation equation for a two level system is
$$\dot{\rho} = \frac{1}{i \hbar} \left[H,\rho\right] - \frac{\gamma}{2} \left(\sigma^+ \sigma^- \rho + \rho \sigma^+ \sigma^- - 2 \sigma^- \rho \sigma^+ \right)$$
where $\sigma^\pm$ are the raising and lowering operators for the two level system, $\rho$ is the system density matrix in the interaction picture and $H=\hbar \Omega (\sigma^+ + \sigma^-)$ is the laser driving Hamiltonian. $\Omega$ is proportional to the driving field and the dipole matrix element.
If one solves this system (e.g. most simply finds the steady state) one will see that the populations1 are not affected by the driving in linear order of $\Omega$. Note that the question is not about an example, but about formulating this notion in its most general form and how to prove it.

1 The diagonal elements of the density matrix.
 A: I will consider a general multi-level system.
Like in your example, let us write $\rho$ for the density matrix of the system in the interaction picture.
That means that the populations of the system are* the diagonal entries of $\rho$,
$$ p_n(t) = \langle n | \rho(t) | n \rangle . $$
Generically, the system Hamiltonian (in the interaction picture) can be split into two parts:
$$ H = H^0 + K . $$
$H^0$ is the diagonal part.
For example, for the two-level system it is the Lamb/Stark shift proportional to $\sigma^z$, which was neglected in your example.
On the other hand, $K$ is the off-diagonal part, i.e. the driving (in your example, the $\sigma^+$ and $\sigma^-$ terms).
Remark: $H$ is often time-dependent in practice. Time-dependence would not change the following much, and I want to keep the notation simple.
The time evolution of the density matrix is then
$$ \dot \rho(t) = -\frac{\mathrm i}{\hbar} [H^0 + K, \rho(t)] + \hat D \rho(t) , $$
where $\hat D$ is the dissipative part which is not terribly important for this question.
Using this, we can calculate how the populations evolve in time:
$$ \dot p_n(t) = \langle n | [K, \rho(t)] | n \rangle + \langle n | \hat D \rho(t) | n \rangle . \tag{1} $$
Note that the first term is zero, because $\langle n | [H^0, \rho(t)] | n \rangle = E^0_n\, p_n(t) - p_n(t)\, E^0_n = 0$, where $E^0_n$ are the eigenvalues of $H^0$.

Assume now that the density matrix is diagonal initially:
$$ \rho(0) = \sum_n p_n(0)\, |n \rangle\!\langle n| .$$
In that case,
$$ \dot p_n(0) = \langle n | \hat D \rho(0) | n \rangle $$
does not depend on $K$.
The first term is zero, because $\rho(0) | n \rangle = p_n(0) | n \rangle$ and we can use the same trick that we used earlier for $\langle n | [H^0, \rho(t)] | n \rangle$.
Dissipation does its work immediately, but the driving will only change the populations after a short while, when $\rho(\Delta t)$ is not diagonal any more:
$$ \rho(\Delta t) = \rho(0) + \left( -\frac{\mathrm i}{\hbar} [K, \rho(0)] + \hat D\rho(0) \right) \Delta t $$
so that
$$ \dot p_n(\Delta t) = -\frac{\mathrm i}{\hbar} \langle n | [K, [K, \rho(0)]] | n \rangle \Delta t + \cdots . $$
You see that the effect is quite general: The driving field only changes the populations at second order.
And its easy to understand, why: (1) tells us that the change of the populations due to $K$ is proportional to the coherences.
If the system is initially diagonal, the driving has to create these coherences first.

If the system is not initially diagonal, the whole statement is not true at all.
Consider your example without the dissipative part, that can be solved easily.
For $\rho(0) = \scriptstyle\begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}$, you get
$$ \rho(t) = \frac 1 2 \begin{pmatrix} 1 - \cos(2\Omega t) & -\mathrm i \sin(2\Omega t) \\ \mathrm i \sin(2\Omega t) & 1 + \cos(2\Omega t) \end{pmatrix} , $$
the populations change in second order.
After a $\pi/2$-pulse, the system is in the state $\scriptstyle \frac 1 2 \begin{pmatrix} 1 & -\mathrm i \\ \mathrm i & 1 \end{pmatrix}$.
Let us for example see what happens if we start in that state, $\rho(0) = \scriptstyle \frac 1 2 \begin{pmatrix} 1 & -\mathrm i \\ \mathrm i & 1 \end{pmatrix}$:
$$ \rho(t) = \frac 1 2 \begin{pmatrix} 1+\sin(2\Omega t) & -\mathrm i \cos(2\Omega t) \\ \mathrm i \cos(2\Omega t) & 1-\sin(2\Omega t) \end{pmatrix} . $$
The populations change in the first order.

*Side note: This is only true in the interaction picture. Without the interaction picture, the populations are $p_n(t) = \langle E_n(t) | \rho(t) | E_n(t) \rangle$, where $|E_n(t)\rangle$ are the - in general time-dependent - eigenvectors of the system Hamiltonian.)
A: The question is actually not affected by the fact that the system is dissipative (except the disspation is much stronger than the driving, but then the question makes no sense whatsoever), so you need not the full master equation approach but could deal with Hamiltonian dynamics. 
That said, first you have to recall how your master equation was derived. The 'Bible' in quantum optics is the Cohen-Tannoudjji - Atom-photon interactions, where the derivation of the quantum optical master equation is described in Chapter IV-B. It starts by introducing the full Hamiltonian
$$ H = H_a + H_r + V ,$$
where the parts are the atomic Hamiltonian $H_a$, the Hamiltonian of the radiation field $H_r$ and the interaction Hamiltonian $V$ between them.
The Liouville-von Neuman- equation of the full system, characterized by a density operator $\rho_{a+r}$ ,is given in the interaction picture by
$$
\frac{\partial}{\partial t}\tilde\rho_{a+r}(t) = \frac{1}{\mathrm i \hbar} [\tilde V(t) ,\tilde\rho_{a+r}(t) ] ,
$$
where we have already the transformed quantities denoted by $\tilde{}$. From this equation the master equation is then derived by taking several assumptions, and by tracing out the degrees of freedom you're not interested in, so that you obtain the density operator of a subsystem only. (That could be either $\tilde\rho_{a}$ when you treat the radiation field as a reservoir as in spontaneous emission, e.g.; or $\tilde\rho_{r}$ in the case of the micromaser experiments of the Haroche group.)
In any case, you always start with the full exact model. Formally, you solve it by integrating the last equation, where you get
$$
\rho_{a+r}(t) = \rho_{a+r}(0) + \frac{1}{\mathrm i \hbar} \int_0^t \mathrm d t' [\tilde V(t') , \tilde \rho_{a+r}(t') ].
$$
Substituting this back into the Liouville-von Neumann equation gives
$$
\frac{\partial}{\partial t}\tilde\rho_{a+r}(t) = \frac{1}{\mathrm i \hbar} [\tilde V(t) ,\tilde\rho_{a+r}(0) ] - \frac{1}{\hbar^2} \int_0^t \mathrm d t' [\tilde V(t) ,[\tilde V(t') ,\tilde \rho_{a+r}(t') ]].
$$ 
this equation is still exact. Integrating again, we obtain a second order equation,
$$
\rho_{a+r}(t) = \rho_{a+r}(0) + \frac{1}{\mathrm i \hbar} \int_0^t \mathrm d t' [\tilde V(t') , \tilde \rho_{a+r}(t') ] - \frac{1}{\hbar^2} \int_0^t \mathrm d t' \int_0^{t'} \mathrm d t'' [\tilde V(t') ,[\tilde V(t'') ,\tilde \rho_{a+r}(t'') ]].
$$
The last steps I've made here could be iterated further. The next step would give a commutator of the third order, etc. If you look closely here, it reminds us of the derivation of the Dyson series, or the derivation of the transition amplitudes in time-dependent perturbation theory. It does not remind us, it is the same terms that are appearing in all three approaches. 
In the derivation of the master equation, the next step is the Born approximation, where you neglect the higher-than-second order couplings (that are implicitly contained in $\tilde \rho_{a+r}(t')$ and $\tilde \rho_{a+r}(t'')$) by replacing $\rho_{a+r}(t')$ and $\tilde \rho_{a+r}(t'')$ by $ \tilde \rho_{a+r}(0)$, so that we obtain
$$
\rho_{a+r}(t) = \rho_{a+r}(0) + \frac{1}{\mathrm i \hbar} \int_0^t \mathrm d t' [\tilde V(t') , \tilde \rho_{a+r}(0) ] - \frac{1}{\hbar^2} \int_0^t \mathrm d t' \int_0^{t'} \mathrm d t'' [\tilde V(t') ,[\tilde V(t'') ,\tilde \rho_{a+r}(0) ]].
$$
Here we threw away higher-order interaction terms that were before implicitly contained.
Next steps would be now the trace over the bath variables, the Markov approximation, etc. But for the point I'm trying to make, it is just important to realize that taking the Born approximation in this way is equivalent to second-order perturbation theory, see also e.g. The Markov master equations and the Fermi golden rule, Alicki, R. Int J Theor Phys (1977) 16: 351. 
Therefore I continue my reasoning with the terms in perturbation theory.
In perturbation theory, the proportionality of the driving field enters via the order in the terms. Usually, the first order is considered where you have terms of the form (in the interaction picture, of course)
$
\langle \psi_f | V | \psi_i \rangle
$, where $| \psi_i \rangle $ and $| \psi_f \rangle $ are the initial and final state, and $V$ is the interaction potential, that presumably is proportional to you driving field. 
In second order, you have terms of the form 
$
\sum_k \langle \psi_f | V | \psi_k \rangle \langle \psi_k | V | \psi_i \rangle 
$. where the sum is over all intermediate states $| \psi_k \rangle $.
So these are the interaction terms where the transition would be second order in the driving field.

Coming to your main question: If you want the populations to change only in second order, you need to make the first order vanish. That is, you choose an interaction and the states such, that the first order matrix elements are zero. An example for this is a dipole radiation field interacting on an electric quadrupole transition. 
An example is the S-D transition in Hydrogen-like atoms.

As you want to know how to find the conditions and how to prove this in general, you need to follow the steps in deriving your Master equation from perturbation theory, where you then identify the symmetries of the atomic operators in relation to your interaction Hamiltonian. But this really depends on the specifics of the system you consider.
An example of this approach a quick search gave is arXiv:1608.04163
