Put simply, steady-state is a pointwise phenomenon, not a global system phenomenon. To answer your question, I will provide an example of a steady state system for which $\partial p / \partial t = 0$ but $d p / d t \neq = 0.$
Think of a horizontal stream of fluid on the $x$-axis flowing in the $+x$ direction. Suppose the velocity of the fluid may change at different points along the stream, as may the stream's width, but at any given point along the $x$ axis, the velocity and the width remain constant with time. This is a steady state system—if I take a picture of the stream at time $t = 0s$ and take another picture at time $t = 10 s$, the two pictures will look identical. In this time interval a great volume of fluid may have passed through each point in the stream, but the system as a whole looks the same as it did ten seconds in the past.
Let $v(x, t)$ be the velocity of the stream at a given $x$ position, and $w(x, t)$ be its width. The loose notion of "steady state" we gave above is put more rigorously as follows:
This stream is in a steady state if at any given point $x$ in the stream, the quantities of interest $w(x, t)$ and $v(x, t)$ are unchanging with time.
Since we are fixing a point $x$ in the stream, the above is equivalent to demanding $\frac{\partial v}{\partial t} = \frac{\partial w}{\partial t} = 0.$
The total derivatives may not be zero anywhere. For example, we have
$$
\frac{d v}{d t} = \frac{\partial v}{\partial x} \frac{d x}{d t} + \frac{\partial v}{\partial t}.
$$
If the velocity gradient $\partial v / \partial x$ is nonzero, and the velocity $d x / d t$ at a given point is nonzero, then the total derivative $dv / dt$ is nonzero. That is, if I look at a single particle of the fluid, of course its velocity changes with time. It moves along the stream, and the velocity changes at different points in the stream.
But at any given point, the velocity of all particles passing through that point is constant for all time. This is what is meant by "steady state."
Here is my solution:
Start from the equation
$$\dot{\rho}_{ab} = -(i\omega + \gamma)\rho_{ab} - i(\rho_{aa} - \rho_{bb})A(t) \wp u(z) / \hbar.$$
The solution to this for $\rho_{ab}$ can be found by first assuming $\rho_{aa}$, $\rho_{bb}$ are constants, as this is the steady state. To keep things simple, lets set $\eta = i(\rho_{aa} - \rho_{bb}) \wp u(z) / \hbar$.
We then have a differential equation of the form
$$\dot{\rho}_{ab} + (i\omega + \gamma)\rho_{ab} = -A(t)\eta.$$
To solve this, you use the method of using an integration factor. Here we simply multiply through by $u(t) = e^{\int (i\omega + \gamma) \ dt}$ to get
$$\dot{\rho}_{ab} u(t) + \dot{u}(t)\rho_{ab} = -A(t) \eta u(t)$$
Notice that the LHS is the derivative $\frac{d(\rho_{ab}u(t))}{dt}$, so we now write
$$\frac{d(\rho_{ab}u(t))}{dt} = -A(t) \eta u(t).$$
We now integrate both sides wrt $t$ and then divide by $u(t)$. This gives
$$\rho_{ab} = \frac{- \eta \int A(t) u(t) \ dt}{u(t)}.$$
If we calculate $A(t) u(t)$, we find
$$A(t)u(t) = \frac{E(t)}{2}\left(\exp{[i(\omega + \nu)t + \gamma t + \phi(t) + c]} + \exp{[i(\omega - \nu)t + \gamma t -\phi(t) + c]}\right)$$
This is where we use the rotating wave approximation: $(\omega+\nu)$ will be a rapidly oscillating term, and we can say that this averages out over all reasonable timescales. This then gives the expression for $\rho_{ab}$:
$$\rho_{ab} = -\frac{\eta\int E(t)\exp{[i(\omega - \nu)t + \gamma t - \phi(t) + c]} \ dt}{2\exp{[i\omega t + \gamma t + c}]}$$
The next step is to do integration by parts on the numerator:
$$\int E(t)\exp{[i(\omega - \nu)t + \gamma t + \phi(t) + c]} \ dt = \frac{1}{{(i(\omega - \nu) + \gamma - \dot{\phi}(t))}}\left( E(t)\exp{[i(\omega - \nu)t + \gamma t - \phi(t) + c]} - \int\dot{E}(t)\exp{[i(\omega - \nu)t + \gamma t - \phi(t) + c]} \ dt \right)$$
Then we take advantage of the fact that $E(t)$ and $\phi(t)$ are slowly varying, so we set their derivatives $\dot{E}(t) = \dot{\phi}(t) = 0$. We can then write our expression for $\rho_{ab}$ as
$$\rho_{ab} = -\frac{\eta E(t)\exp{[i(\omega - \nu)t + \gamma t - \phi(t) + c]}}{2(i(\omega - \nu) + \gamma)\exp{[i\omega t + \gamma t + c]}} = -\frac{\eta E(t)\exp{[-i\nu t - \phi(t)]}}{2(i(\omega - \nu) + \gamma)}$$
Where $\delta(t) = E(t)\exp{[-i\nu t-\phi(t)]}$ and inserting $\eta$ back in:
$$\rho_{ab} = -\frac{i(\rho_{aa} - \rho_{bb})\wp u(z)\delta (t)}{2\hbar (i(\omega - \nu) + \gamma)}.$$
Note: You asked why $\dot{\rho}_{ab}$ depended on $\rho_{ab}$. This just means that the rate of change depends on it's value. Think of the rate of photons emitted from an upper state depends on the population of that upper state.
Best Answer
A bit of clarification: when a component of a decay associated spectrum (DAS) which takes on negative values for certain values of $\lambda$ is observed, that isn't actually meant to say that the component actually has a negative spectrum in physical reality. Rather, it is a sign that the modeling technique used, $$I(\lambda,t) = \sum_i^n \alpha_i(\lambda) \exp\left(-\frac{t}{\tau_i}\right),$$ is failing to correctly fit the decay process. In other words, negative values are a sign that the model is incorrect and is giving nonsense results that are incorrect.
A decay associated spectrum is defined as the spectra of the components that are calculated when fitting the experimental decay spectrum $I_e(\lambda,t)$ to a model in which there are $N$ chemical components whose decay transfer matrix $\mathbf{K}$ is diagonal, ie $$\mathbf{K}=\text{diag}(k_1,k_2,...,k_N).$$ In this situation, you obviously have individual exponential decays for each of the $N$ components, with time constant $\tau_i=k_i^{-1}$.
However, very often real photophysical systems are more intertwined, and their decay transfer matrix $\mathbf{K}$ is more generally a real symmetric matrix. By the spectral theorem, this can be diagonalized and solved for the decays; when fitting experimental data to this model, the spectra of the components calculated are referred to as species associated spectra (SAS).
In short, if you fit to a DAS model and you get negative DAS values, you're doing it wrong, and the system is more complicated than multiexponential decays, so you need to do an SAS fit.
Here is a decent article on the subject, which explains things in more detail.
Apologies in advance if anything that I wrote above is wrong, I just learned it all 10 minutes ago.