How do I include a quadratic time trend to a plm fixed effects model in R?

Question

I am working with an unbalanced panel of 16 entities and 38 years. So far, I've simply used the plm() twoways specification for my panel and I clustered heteroscedastic standard errors at the state level. Now, I'd like to include a state-specific quadratic time trend.

So far, I've tried to do it with a variable stating the years of the observations for each state from 1980 to 2019 in a monochronic order and just included this variable as explaining variable:

fixed_trend <- plm(X ~ Y + Z + time + time^2, data = df, model = "within", effect = "twoways")
summary(fixed_trend, vcovHC(fixed_trend, type="HC3", cluster="group"))

The output of this regression, however, doesn't show the trends. It only gives me the coefficients of my explaining variables (in this example, hence, Y and Z). Can you tell me what I did wrong?