General Result
If the linear PDE is well-posed, it is enough to determine the stability of your numerical scheme by only studying the numerical discretization of the homogenous PDE.
The Details...
As long as the inhomogenous function $f$ isn't too nasty (i.e. the problem is well-posed), it suffices determine the numerical stability of your scheme by just working with the homogenous pde. To see why, let's consider the following linear inhomogenous differential equation defined on
\begin{align}
u_t = Du + f(x,t)
\end{align}
where $D$ is the linear differential operator and $f(x,t)$ is the inhomogenous forcing term. Let $u^n_j = u(j\Delta x ,n\Delta t)$ denote the discretization of $u$. Using the forward Euler method to discretize this equation in time, the discretized version of the PDE becomes
$$\frac{u^{n+1}_j - u^n_j}{\Delta t} = Au_j^n + f^n_j$$
where $A$ is the matrix which approximates the linear differential operator $D$ to some desired level of accuracy. Furthermore, we can represent the grid functions above using Fourier analysis. We have
$$u_j^n = \int_{-\frac{\pi}{\Delta x}}^{\frac{\pi}{\Delta x}} \hat{u}(\xi)e^{i\xi\Delta x j}d\xi, \quad f^n_j = \int_{-\frac{\pi}{\Delta x}}^{\frac{\pi}{\Delta x}}\hat{f}(\xi)e^{i\xi \Delta xj}d\xi.$$
Plugging in this expression into the discrete PDE above, we obtain
$$\int_{-\frac{\pi}{\Delta x}}^{\frac{\pi}{\Delta x}} \left(\frac{\hat{u}^{n+1} - \hat{u}}{\Delta t} - g(\xi,\Delta x)\hat{u}^n - \hat{f}^n\right)e^{i\xi\Delta x j}d\xi = 0$$
where $g(\xi,\Delta x)$ is the amplification factor resulting from applying the matrix $A$ to the eigenvector $e^{i\xi\Delta x j}$.
Since the entire integral is equal to zero, we must have that the integrand itself is zero and so
$$\hat{u}^{n+1} = \hat{u}^n\left(1 + \Delta t g(\xi,\Delta x)\right) + \Delta t \hat{f}^n.$$
This implies that
$$\hat{u}^n = \left(1 + \Delta t g(\xi,\Delta x)\right)^n\hat{u}^0 + \Delta t \sum_{k=0}^{n-1}\hat{f}^k(1+g(\xi,\Delta x))^{n-1-k}$$
Now, if we assume that $|1+\Delta tg(\xi,\Delta x)| \le 1$ where $|\cdot|$ is the complex modulus, we can multiply the above equation by its complex conjugate $\bar{\hat{u}}^n$, perform some somewhat tedious algebraic manipulation, and apply the triangle inequality to obtain that
$$|\hat{u}|^2\le |\hat{u}^0|^2 + 16\Delta t|\hat{u}^0|\sum_{k=0}^{n-1}|\hat{f}^k| + \Delta t^2\sum_{k=0}^{n-1}|\hat{f}^k|^2$$
Integrating both sides of this equation from $\frac{-\pi}{\Delta x}$ to $\frac{\pi}{\Delta x}$ and applying the Cauchy-Schwarz inequality then gives
$$||\hat{u}^n||_{2}^2 \le ||\hat{u}^0||_2^2 + 16\Delta t||\hat{u}^0||_2\sum_{k=0}^{n-1}||\hat{f}^k||_2 + \Delta t^2\sum_{k=0}^{n-1}||\hat{f}^k||_2^2$$
Therefore, thanks to Parseval's theorem, we may conclude that the solution is stable so long as $n\Delta t \sim \mathcal{O}(1)$, $f^n_j \in l_2$ ($f$ isn't too nasty) for each $n$ and our initial data $u^0_j\in l_2$. Since we only required that the modulus $|1+\Delta tg(\xi,\Delta x)|\le 1$, which is the same requirement needed for the homogenous problem, it is enough to deduce stability just by studying the numerical discretization of the homogenous problem. We applied the forward Euler time discretization to arrive at this conclusion, but you can easily apply this result to other time discretizations by applying slight modifications.
Best Answer
The concept of stability is, in and of itself, meaningless. I mean this in the sense that stability is trivial to achieve. We always have the possibility to set all the finite difference coefficients to zero, giving a numerical solution that is zero. This is stable no matter the definition of stability we use. However, as a concept it is clearly void of any information.
Stability is only meaningful if we also account for some notion of accuracy. What this notion should entail is strongly problem dependent. For example, consider the linear PDE $$ u_t + a u_x = 0, $$ posed on a periodic domain with a periodic initial condition $u(0,x) = u_0$. Here $a$ is some constant. The solution satisfies the estimate $\| u \| = \| u_0 \|$ for all times. In this example I am considering the $L^2$ norm. A relevant stability concept for this problem would be something like $\| \mathbf{u}^{n+1} \|_h \leq \| \mathbf{u}^n \|_h$. The discrete norm is chosen in a way that approximates the norm used in the continuous setting. We would aim for $\| \mathbf{u}^{n+1} \|_h = \| \mathbf{u}^n \|_h + \mathcal{O}(\Delta t^p, \Delta x^q)$ where the ordo-term is non-positive and $p$ and $q$ are constants that describe the desired accuracy of the numerical scheme. If we can't obtain such an estimate, then we should not expect this accuracy.
On the other hand, suppose that we are instead solving the PDE $$ u_t + a u_x = c u, $$ where $c$ is another constant. The solution now satisfies $\| u \| = e^{ct} \| u_0 \|$ for all times and a relevant concept may be something like $\| \mathbf{u}^{n+1} \|_h \leq K \| \mathbf{u}^n \|_h$, where $K \approx e^{c \Delta t}$. The symbol $\rho$ you provide in the question is such an approximation. The sign of $c$ obviously has a strong influence on the nature of this stability concept.
Thus, to answer your questions, we should bother about different stability definitions because they only make sense in the context of the problem we are solving. The Neumann or the 'regular' stability concepts should not be thought of as more or less flexible - after all, what is the point in allowing for exponential growth if the exact solution is decaying? Rather, what the Neumann analysis allows us to do is to obtain some specific information about the numerical solution that a 'regular' (say, based on energy estimates for example) stability analysis may not provide, and vice versa.