Combining FEs in feols.fit

Question

When using fixest's feols it is simple to add combinations of fixed-effects with ^ as in z ~ x1 + x2 | fe_1^fe_2.

However, how could I implement this type of interaction in feols.fit, which requires matrices/vectors instead of a formula?

Assuming that z, fe_1 and fe_2 are now vectors and that a matrix X0 contains the variables x1 and x2, I tried

feols.fit(y = z,
          X = X0,
          fixef_df = fe_1^fe_2)

but it didn't work.

Thanks