Run a polynomial regression without combinations of the features
Tbh, this seems like very weird idea to me. A lot of strength of poly features come out of interaction features.
I'm pretty sure there's nothing built-in, but

you can take a look at PolynomialFeatures pre-processor, and make your own with modification. You'll need to "revert" logic of interaction_only parameter
just build your pre-processor, your case is quite simple, .fit will do smth like concat(a, a**2)
I am running a polynomial regression for order p. To make it simple, we use order p = 2 in this question.
Suppose we have X with two feature x1, x2 and y. And I am trying to run a polynomial regression of
y = ε + α + β1•x1 + β2•x2 + β3•x1^2 + β4•x2^2
I find that the sklearn have a sklearn.preprocessing.PolynomialFeatures. However, if I use order p = 2 and it automatically gives the combination of features. It will result in a regression of:
y = ε + α + β1•x1 + β2•x2 + β3•x1^2 + β4•x2^2 + β5•x1x2
However, I do not want the combination of the features, i.e. things like x1x2. Is there any package that can do the polynomial regress as I wished?
Thanks!