Consider a GAM model, expressed in mgcv just to fix ideas:
my_model <- gam(y ~ ti(x1)+ti(x2) + ti(x1, x2), method= "REML")
The model is linear in the parameters, right? Each smooth is a linear combination of basis functions, which are independent of the data set (unless I use bs = "ad"). Thus the model is linear in the parameters, which are the coefficients of the basis functions. Right? And this should be true, whether or not there are interaction terms - it doesn't really matter. The only exception are adaptive smooths, because in that case the coefficients of the basis functions are themselves functions of the covariates (x1 and x2, in my example). Correct?