There seems to be a problem with the starting equation and the equation you ended up with. Your intuition is on the right track. The equation you wrote down is actually the advection equation not the diffusion equation. The radiative transfer equation is the just the continuity equation for light transport.
To recap, the continuity equation keeps track of all the gains and losses of a conserved quantity $u$. In 1D
$$
\frac{\partial u}{\partial t} = \frac{\partial \mathcal{F}}{\partial x} + s(t, x, u)
$$
Where $\mathcal{F}$ is the flux of the conserved quantity and $s$ is the source function which maybe depend on the other variables.
To recover the radiative transfer equation in the limit of zero scattering, let,
$$
u = n_{\gamma} \\
\mathcal{F} = \hat{n} \cdot I_{\gamma} \\
s = - k_{\gamma, a} I_{\gamma} + j_{\gamma}
$$
- $n_{\gamma}$ $[\text{eV}^{-1} \text{cm}^{-3}]$ is the number density of photons in spectral interval $(E, E + dE)$. This is related to the photon flux by the following, where $n$ is the refractive index and $c$ is the speed of light,
$$n_{\gamma} dE = \frac{n}{c} dE \int I_{\gamma} d\Omega$$
$\hat{n}\cdot I_{\gamma}$ $[\text{eV}^{-1} \text{cm}^{-2} \text{Sr}^{-1} \text{s}^{-1}]$ is the photon flux in spectral interval $(E, E+dE)$, transported along direction $\hat{n}$, across an element of area, confined to element of solid angle per unit time. This is the specific intensity that Chandrasekhar refers to, but with different units (photons rather than energy).
$k_{\gamma, a}$ $[cm^{-1}]$ is the absorption coefficient.
$j_{\gamma}$ $[\text{eV}^{-1} \text{cm}^{-3} \text{Sr}^{-1} \text{s}^{-1}]$ is the emission coefficient.
It is not necessary to define the units of the quantities exactly as I have done here. You could consider using integrated values over photon energy and solid angle to remove the dependence on those dimensions.
If you apply the same logic using the nabla operator, you can recover the original equation (in the limit of zero scattering).
$$
\frac{1}{c}\frac{\partial I_{\gamma}}{\partial t} = -\hat{n}\cdot\nabla I_{\gamma} -k_{\gamma, a} I_{\gamma} + j_{\gamma}
$$
For your application you can probably assume local thermodynamic equilibrium and steady-state, which will give the classic radiative transfer equation from Chandrasekhar’s book,
$$
\hat{n}\cdot\frac{1}{k_{\gamma, a}} \nabla I_{\gamma} = -I_{\gamma}+ B_{\gamma}(T)
$$
The source term is derived as follows,
$$
s = -B_{\text{ein}} n_{\gamma} + j_{\gamma}
$$
where $B_{\text{ein}}$ is the Einstein B-coefficient. This is related to absorption coefficient by,
$$
B_{\text{ein}} = \frac{c}{n}k_{a, \gamma}
$$
and by using the definition from above,
$$
n_{\gamma} = \frac{n}{c} I_{\gamma}
$$
we arrive at
$$
s = - k_{\gamma, a} I_{\gamma} + j_{\gamma}
$$