I'm trying to linearize the following PDE in an attempt to find an analytical solution:$$\frac{\sigma_I}{r} \frac{\partial}{\partial r} \left(r \frac{\partial z(r,t)}{\partial r} \right)=\rho_D gz(r,t)+\rho_T\left(\frac{\partial z(r,t)}{\partial t}\right)^2-2 \mu_T \frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial z(r,t)}{\partial r} \right) \frac{\partial z(r,t)}{\partial t}\tag{1}$$ which is applicable on $r_m\leq r\lt\infty$. The initial condition is $z(r,0)=0$, the outer boundary condition is $z(\infty,t)=0$, and the inner boundary condition varies with time and essentially works to drive the change. Under steady state, equation 1 reduces to:$$\frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial z}{\partial r} \right)-\frac{\rho_D g}{\sigma_I}z=\frac{\partial^2 z}{\partial r^2}+\frac{1}{r} \frac{\partial z}{\partial r}-\frac{z}{\lambda^2}=0\tag{2}$$Which has a solution in terms of the modified Bessel functions of order zero:$$z=AI_0\left(\frac{r}{\lambda}\right)+BK_0\left(\frac{r}{\lambda}\right)\tag{3}$$I know that as $r$ goes to infinty, $z$ goes to zero so $A=0$ since $I_0$ is unbounded at large argument values. Technically, the inner boundary condition varies with time which means $B$ varies with time but I will treat it as a constant for now since I'm considering a steady state solution. I then substitute $z(r,t)=BK_0\left(\frac{r}{\lambda}\right)+\varepsilon w(r,t)$. I'm going to drop the argument of the Bessel fucntion so it's a bit neater:$$\frac{\partial z}{\partial r}=-BK_1+\varepsilon\frac{\partial w}{\partial r}\tag{4}$$ $$\frac{\partial^2 z}{\partial r^2}=BK_0+B\frac{1}{r}K_1+\varepsilon\frac{\partial^2 w}{\partial r^2}\tag{5}$$ $$\frac{\partial z}{\partial t}=\varepsilon\frac{\partial w}{\partial t}\tag{6}$$ Expanding equation 1 and substituting in equations 4-6 then gives:$$\sigma_I\left(\frac{\partial^2 z}{\partial r^2}+\frac{1}{r} \frac{\partial z}{\partial r}\right)=\rho_D gz+\rho_T\left(\frac{\partial z}{\partial t}\right)^2-2 \mu_T \left(\frac{\partial^2 z}{\partial r^2}+\frac{1}{r} \frac{\partial z}{\partial r}\right) \frac{\partial z}{\partial t}\tag{7}$$$$\sigma_I\left(BK_0+B\frac{1}{r}K_1+\varepsilon\frac{\partial^2 w}{\partial r^2}+\frac{1}{r} \left(-BK_1+\varepsilon\frac{\partial w}{\partial r}\right)\right)=\rho_D g(BK_0+\varepsilon w)+\rho_T\left(\varepsilon\frac{\partial w}{\partial t}\right)^2-2 \mu_T \left(BK_0+B\frac{1}{r}K_1+\varepsilon\frac{\partial^2 w}{\partial r^2}+\frac{1}{r} \left(-BK_1+\varepsilon\frac{\partial w}{\partial r}\right)\right) \varepsilon\frac{\partial w}{\partial t}\tag{8}$$$$\sigma_IBK_0+ \varepsilon\sigma_I\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial w}{\partial r}\right)=\rho_D g(BK_0+\varepsilon w)+\varepsilon^2\rho_T\left(\frac{\partial w}{\partial t}\right)^2 -\varepsilon2 \mu_TBK_0\frac{\partial w}{\partial t}-\varepsilon^22 \mu_T\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial w}{\partial r}\right)\frac{\partial w}{\partial t} \tag{9}$$ From here, I've seen two methods: (1) take the derivative with respect to $\varepsilon$ and then let $\varepsilon$ go to zero, or (2) just drop the $\varepsilon^2$ terms and keep the $\varepsilon$ terms. Applying each to equation 9 gives me equations 10 and 11 respectively after rearranging some terms:$$-B\frac{2\mu_T}{\sigma_I}K_0\frac{\partial w}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial w}{\partial r}\right)-\frac{w}{\lambda^2}\tag{10}$$ $$-B\frac{2\mu_T}{\sigma_I}K_0\frac{\partial w}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial w}{\partial r}\right)-\frac{w}{\lambda^2}+BK_0\left(1-\frac{1}{\lambda^2}\right)\tag{11}$$ where $\lambda^2=\frac{\sigma_I}{\rho_Dg}$ as before. And thus I have two possible solutions for the linearization.

I have two main concerns: (1) Did I even do this correctly?, and (2) which is the proper linearization if either?

I posted this in MSE here, but didn't get any response. Any suggestions or help would be greatly appreciated, thanks!