There are two modes of operation for masers. For unsaturated masers, the gain is exponential with the path length, and for saturated ones, the gain only grows linearly with path length. The masers we find in the cosmos are typically saturated. The following is working toward understanding why there are two modes in masers. We begin with the expression for specific intensity from emissivity: \begin{aligned}{dI_\nu\over dz}&=j_\nu-\alpha_\nu I_\nu\\ &={h\nu\over4\pi}\phi(\nu)n_2A_{21}-{h\nu\over4\pi} \phi(\nu)[n_1\bot-n_2\bto]I_\nu\\ \end{aligned} In principle, the line-profile function governing spontaneous emission and the line-profile governing stimulated emission might not have to be the same, but evidence seems to suggest they are. Let's assume they are, and rewrite this: \begin{aligned}{dI_\nu\over dz}&={h\nu\over4\pi}\phi(\nu){n_2\ato\over g_2}\cdot g_2 -{h\nu\over4\pi}\phi(\nu)g_2\bto\left[{n_1\over g_1}-{n_2\over g_2}\right] \end{aligned} Defining $N_1\equiv{n_1\over g_1}$, $N_2\equiv{n_2\over g_2}$, $A\equiv \ato g_2$, and $B=\bto g_2$, our equation looks like: $${dI_\nu\over dz}={h\nu\over4\pi}\phi(\nu)[N_2A+B(N_2-N_1)I_\nu]$$ Now we'll integrate over frequency. $\phi(\nu)$ is a sharply peaked function, so we'll treat it as a $\delta$-function. Then substituting $I\equiv\int{\phi(\nu)I_\nu d\nu}$, we have: $${\int{I_\nu d\nu}\over\Delta\nu}\equiv\int{I_\nu\phi(\nu)d\nu}$$ which is true by definition of $\Delta\nu$: $\int{I_\nu d\nu}=I\Delta\nu$. So now we have: $${dI\over dz}={h\nu\over4\pi\Delta\nu}[(N_2-N_1)BI+N_2A]$$ This is an equation with three unknowns. To close this system, we'll use the equations of statistical equilibrium, which say: $${dN_2\over dt}=-N_2BJ+N_1BJ-N_2A+R_2(N-N_1-N_2)-\Gamma_2N_2$$ where $J\equiv\int{I_\nu\phi(\nu)d\nu}$, and $J_\nu=\inv{4\pi}\int{I_\nu d\Omega}$. $R_2$ is the pumping rate. It describes the rate at which molecules not in state 1 or 2 (counted by $N-N_1-N_2$) are radiatively or collisionally knocked into state 2. $\Gamma_2$ is the loss rate of molecules in state 2 into any state other than state 1. In the above equation, we've neglected two terms: collisional excitation $1\to2$ and $2\to1$. These are crucial terms, being tightly related to local thermal equilibrium. However, we will neglect them to simplify analysis. We also have a similar equation for state 1: $${dN_1\over dt}=N_2BJ-N_1BJ+N_2A+R_1(N-N_1-N_2)-\Gamma_1N_1$$ We'll further simply matters by setting the loss-rates for the two populations equal to each other ($\Gamma_1=\Gamma_2=\Gamma$). Now let's solve for the steady-state solution (${dN\over dt}=0$). $${d(N_2-N_1)\over dt}=-(N_2-N_1)2BJ-2N_2A+(R_2-R_1)(N-N_{12})-\Gamma(N_2-N_1)$$ where $N_{12}=N_1+N_2$. Two terms in this equation are acting to reduce the population inversion of state 2 with respect to state 1: $(N_1-N_1)2BJ$ and $\Gamma(N_2-N_1)$. In the {\it unsaturated regime} where $BJ\ll\Gamma$, stimulated emission is a minor perturbation to the inverse. In this case we can solve the system of equations: \begin{aligned}{d\Delta N\over dt}={d(N_2-N_1)\over dt}&=-2N_2A+(R_2-R_1)(N-N_{12}) -\Gamma\Delta N=0\\ {dN_{12}\over dt}&=(R_2+R_1)(N-N_{12})-\Gamma(N_1+N_2)=0\\ \end{aligned} Then solving for $N-N_{12}$: $$N-N_{12}={\Gamma(N_1+N_2)\over R_1+R_2}={\Gamma N_{12}\over R_1+R_2}$$ Substituting this into the difference equation, which reads: $$-2N_2A={R_2-R_1\over R_1+R_2}\Gamma N_{12}-\Gamma\Delta N=0$$ the using $N_2={N_{12}+\Delta N\over2}$ (this is just an identity), we get: $$\underbrace{\left(1+{\Gamma\over A}\right)}_{\equiv2\beta}\Delta N =N_{12}\left[\underbrace{{R_2-R_1\over R_1+R_2}{\Gamma\over A}}_{\equiv \inv{\alpha}}-1\right]$$ Thus: $$\Delta N={N_{12}\over2\beta}\left[\inv{\alpha}-1\right]={N_{12}(1-\alpha)\over 1\alpha\beta}$$ Now we'll define $S\equiv{\alpha\beta\over1-\alpha}$. $S$ is meant to connote the source function here. The reason we might expect $S$ to be related to the source function is that $S_\nu\propto{N_2\over N_2-N_1}\propto{N_2\over \Delta N}$, so $\Delta N\propto {N_{12}\over2S}$. This is why we use an $S$ here. Getting back to our original equation, we have: $${4\pi\Delta\nu\over h\nu A}{dI\over dz}=\Delta N\underbrace{BI\over A}_{ \equiv\mathfrak{I}}+N_2$$ Using our newly defined $S$, this becomes: $${4\pi\Delta\nu\over h\nu B}{d\mathfrak{I}\over dz}={N_{12}\over2S}\mathfrak{I}+ \hf(N_{12}+\Delta N)$$ Finally, defining an unsaturated gain length $L\equiv{4\pi\Delta\nu\over Bh\nu} {2S\over N_{12}}$, we have: $$L=\mathfrak{I}+S+\hf$$ Then using $ds={dz\over L}$, we have: $${d\mathfrak{I}\over ds}=\mathfrak{I}+\hf+S$$ To order of magnitude: \begin{aligned}\mathfrak{I}&\equiv{B\over A}I={B\over A}B_\nu(T_{bright})\\ &={c^2\over2h\nu^3}{2h\nu^3\over c^2}\inv{e^{h\nu\over kT_{bright}}-1}\\ &=\inv{e^{h\nu\over kT_{bright}}-1}\\ \end{aligned} Now since $1\ll{kT_{bright}\over h\nu}$ (typical maser wavelengths are in mm and $T_{bright}$ is typically several K), we can throw away our $\hf$ in our equation for ${d\mathfrak{I}\over ds}$: $${d\mathfrak{I}\over ds}=\mathfrak{I}+S\Rightarrow \int{d\mathfrak{I}\over\mathfrak{I}+S}=\int{ds}\Rightarrow \mathfrak{I}+S=De^s$$ Choosing $\mathfrak{I}=0, s=0$ (we're assuming there is no background source), then $D=S$, so: $$\boxed{\mathfrak{I}=S(e^s-1)}$$ If $s\ll1$, we have $\mathfrak{I}=Ss$ from spontaneous emission, and this is the saturated case. If $s\gg1$, the $\mathfrak{I}=Se^s$ from stimulated emission, and this is the unsaturated case. Earlier we threw away the $BJ$ term, but we could have done all of this including that term, and the algebra would have been the same. If we do this, we ge the answer: $$\boxed{{d\mathfrak{I}\over ds}={\beta(\mathfrak{I}+\hf)\over\beta+\mathfrak{J}}+S}$$ where $\mathfrak{J}$ is the non-dimensionalized, integrated flux, $\mathfrak{J}={BJ\over A}$.