regress_out polynomial effect

I was wondering if someone who is familiar with sc.pp.regress_out could confirm the following:
I would like to regress out nonlinear effect, e.g. ~1 + a + a^2 + a^3, where a is non-categorical variable.

I have looked at the code of regress_out: https://github.com/theislab/scanpy/blob/8fe1cf9cb6309fa0e91aa5cfd9ed7580e9d5b2ad/scanpy/preprocessing/_simple.py#L677
It seems that the code performs the fitting for all specified variables at once, but I am not sure:
https://github.com/theislab/scanpy/blob/8fe1cf9cb6309fa0e91aa5cfd9ed7580e9d5b2ad/scanpy/preprocessing/_simple.py#L701
If the design passed to GLM is combined of all keys passed to the function then I could just create the necessary columns a, a^2, a^3 and pass this as keys.

Can someone confirm if I understand this correctly and passing the polynomial columns will do the fitting of a polynom?

Read the linked scanpy/preprocessing/_simple.py regress_out implementation, especially where the GLM design is assembled from the supplied keys. Confirm whether polynomial columns can be passed together and whether the behavior needs a documentation clarification; done means an authoritative answer or documented usage guidance.