Liquid density as a function of pressure and temperature, how to model experimental data of $\rho(p,T)$?

I would do this very differently. Your Eqns. 2 and 3 suggest that a more appropriate form of functionality to start with would be $$\rho=\rho(T_0,P_0)\exp{[-\alpha(T-T_0)+\beta(P-P_0)]}\tag{1}$$or, equivalently, $$\ln{\rho}=\ln{\rho(T_0,P_0)}-\alpha(T-T_0)+\beta(P-P_0)\tag{2}$$ That is, I would work with the natural log of $\rho$ rather than $\rho$ itself. This should also go a long way to more nearly linearizing the behavior with temperature and pressure, and should also help to minimize the percentage errors in density (rather than the absolute errors).

As it stands, my Eqn. 2 represents a 5 parameter fit to the density data ($T_0$, $P_0$, $\rho(T_0,P_0)$, $\alpha$, and $\beta$). I would start out by doing a least squares minimization fit to all the data using this functionality. This would involve reasonable initial guesses for the 5 parameters appropriate, say, to the middle of the range of data.

If the residual errors were found to be too large using the 5 parameters, I would then add a cross-coupling parameter and quadratic terms involving $(T-T_0)^2$ and $(P-P_0)^2$ to the functionality:$$\ln{\rho}=\ln{\rho(T_0,P_0)}-\alpha(T-T_0)+\beta(P-P_0)+\gamma (T-T_0)(P-P_0)+\delta (T-T_0)^2+\epsilon (P-P_0)^2\tag{3}$$

Then, I would re-implement the least squares fit to find a new optimized set.