The answer guessed in the question is correct, but it would be helpful to know what the constant $a$ is. Deriving the solution is a bit tedious, but it is the only way to answer the question, so here we go. I won't show every gory detail, but enough to show how to do the problem.
Make the change of variables $s = r \mu$ and $p = r \sqrt{1 - \mu^2}$ where $s$ is in the direction of the observer and $p$ is the so-called impact parameter. Using the chain rule, we find
$$\frac{\partial I(r,\mu)}{\partial r} = \frac{\partial I(s, p)}{\partial s} \frac{\partial s}{\partial r} + \frac{\partial I(s, p)}{\partial p} \frac{\partial p}{\partial r} = \mu \frac{\partial I(s, p)}{\partial s} + \sqrt{1 - \mu^2} \frac{\partial I(s, p)}{\partial p},$$
$$\frac{\partial I(r,\mu)}{\partial \mu} = \frac{\partial I(s, p)}{\partial s} \frac{\partial s}{\partial \mu} + \frac{\partial I(s, p)}{\partial p} \frac{\partial p}{\partial \mu} = r \frac{\partial I(s, p)}{\partial s} - \frac{r \mu}{\sqrt{1 - \mu^2}} \frac{\partial I(s, p)}{\partial p}.$$
If we substitute this into the radiative transfer equation (the first one in the question), we get
$$\frac{\partial I(s, p)}{\partial s} + \rho \kappa I(s, p) = 0,$$
which has solution
$$I(s, p) = f(p) \exp(-\rho \kappa s). \qquad \mathrm{\left[ \frac{erg}{s \, cm^2 \, sr} \right]}$$
We can convert this back to $r$ and $\mu$. Doing so, we find
$$I(r, \mu) = f\left( r \sqrt{1 - \mu^2} \right) \exp(-\rho \kappa r \mu). \qquad \mathrm{\left[ \frac{erg}{s \, cm^2 \, sr} \right]}$$
Surprisingly, $f(\cdot)$ can be any function of $r \sqrt{1 - \mu^2}$ and still satisfy the radiative transfer equation.
Although we now have the solution to the radiative transfer equation, the form of the solution is not particularly helpful. We have found an infinite number of solutions due to the generic function $f(\cdot)$. To try to alleviate this, we expand $f(\cdot)$ in a Laurent series to get
$$I(r, \mu) = \sum_{l = -\infty}^\infty a_l \left( r \sqrt{1 - \mu^2} \right)^l \exp(-\rho \kappa r \mu). \qquad \mathrm{\left[ \frac{erg}{s \, cm^2 \, sr} \right]}$$
We chose a Laurent series because we know the problem has a singularity at $r = 0$, a feature automatically handled by a Laurent series. Since the domain of the problem is $r > R$, this singularity never enters the problem. However, the Laurent series is helpful nonetheless since $R$ can be arbitrarily small.
If we apply the boundary condition, we find
$$I(R, \mu) = \frac{\delta(\mu - 1)}{2\pi} I_R = \sum_{m = -\infty}^\infty a_m \left( R \sqrt{1 - \mu^2} \right)^m \exp(-\rho \kappa R \mu),$$
which implies that
$$\exp(-\rho \kappa R \mu) = \frac{\delta(\mu - 1)}{2\pi} \frac{I_R}{\sum_{m = -\infty}^\infty a_m \left( R \sqrt{1 - \mu^2} \right)^m}.$$
This can be substituted into our solution as follows
\begin{align}
I(r, \mu) &= \sum_{l = -\infty}^\infty a_l \left( r \sqrt{1 - \mu^2} \right)^l \exp(-\rho \kappa r \mu) \\
&= \sum_{l = -\infty}^\infty a_l \left( r \sqrt{1 - \mu^2} \right)^l \exp(-\rho \kappa [r - R] \mu) \exp(-\rho \kappa R \mu) \\
&= \frac{\delta(\mu - 1)}{2\pi} \frac{\sum_{l = -\infty}^\infty a_l \left( r \sqrt{1 - \mu^2} \right)^l}{\sum_{m = -\infty}^\infty a_m \left( R \sqrt{1 - \mu^2} \right)^m} I_R \exp(-\rho \kappa [r - R] \mu).
\end{align}
Because of the Dirac delta, this expression diverges unless all terms in the Laurent series vanish except for one. That is, we let $a_m = 0$ for all $m$ except for $m = -n$ (the minus sign is for convenience). Applying this condition, we find
$$I(r, \mu) = \frac{\delta(\mu - 1)}{2\pi} \left( \frac{R}{r} \right)^n I_R \exp(-\rho \kappa [r - R] \mu). \qquad \mathrm{\left[ \frac{erg}{s \, cm^2 \, sr} \right]}$$
We are getting closer to a solution, but we still need to determine $n$. It turns out only one value of $n$ conserves energy. To start, we compute the energy density
$$u(r) = \frac{2\pi}{c} \int_0^1 I(r, \mu) \, d\mu = \frac{I_R}{c} \left( \frac{R}{r} \right)^n \exp(-\rho \kappa [r - R]). \qquad \mathrm{\left[ \frac{erg}{cm^3} \right]}$$
We can convert this to specific energy with
$$e(r) = \frac{u(r)}{\rho} = \frac{I_R}{c \rho} \left( \frac{R}{r} \right)^n \exp(-\rho \kappa [r - R]) \qquad \mathrm{\left[ \frac{erg}{g} \right]}$$
and energy fluence with
$$\Psi(r) = \frac{e(r)}{\kappa} = \frac{I_R}{c \rho \kappa} \left( \frac{R}{r} \right)^n \exp(-\rho \kappa [r - R]). \qquad \mathrm{\left[ \frac{erg}{cm^2} \right]}$$
We can think of energy fluence as the total energy $E$ [erg] hitting a surface area $A(r) = 4 \pi r^2$ [cm$^2$] attenuated by the exponential term. Writing fluence this way, we find
$$\Psi(r) = \frac{E}{4 \pi r^2} \exp(-\rho \kappa [r - R]). \qquad \mathrm{\left[ \frac{erg}{cm^2} \right]}$$
Comparing the previous two equations, we require $n = 2$ and find that the total energy is
$$E = \frac{4 \pi R^2 I_R}{c \rho \kappa}. \qquad \mathrm{\left[ erg \right]}$$
To confirm that this conserves energy, we can integrate the energy density $u(r)$ over the domain. Doing so also gives the total energy
$$E = 4 \pi \int_R^\infty u(r) r^2 dr = \frac{4 \pi R^2 I_R}{c} \exp(\rho \kappa R) \int_R^\infty \exp(-\rho \kappa r) \, dr = \frac{4 \pi R^2 I_R}{c \rho \kappa}. \qquad \mathrm{\left[ erg \right]}$$
Hence, $n = 2$ is the only solution that conserves energy. Finally, to answer the question, the energy density is
$$u(r) = \frac{I_R}{c} \left( \frac{R}{r} \right)^2 \exp(-\rho \kappa [r - R]). \qquad \mathrm{\left[ \frac{erg}{cm^3} \right]}$$
