As Richard Hardy comments, exponential smoothing methods like Holt-Winters typically do not allow for external regressors.
On the one hand, extending the state space framework in which
forecast::ets() is fitted to incorporate external regressors should be straightforward. You may want to look at Forecasting with Exponential Smoothing - The State Space Approach by Hyndman, Koehler, Ord & Snyder. However, I do not know of any implementation, so you would need to code this up by yourself.
On the other hand, you could take inspiration from the way
forecast::auto.arima() actually models external regressors. Contrary to first impressions, this is not an ARIMAX model, but a regression with ARIMA errors.
This suggests a way forward: you could simply regress your time series against your regressors (using, say,
lm()), and then run your favorite exponential smoothing model fitting algorithm (
forecast::ets()) on the residuals from this original regression. This would do something different from a true state space model, but it might be worth looking at.
(And if your Holt-Winters model gives you better results than
auto.arima() with regressors, you may want to investigate your model and/or data a bit more. Perhaps your regressors are in reality not all that relevant. Or you have something else in your data that it hard to handle for
auto.arima(), e.g., multiple-seasonalities, missing data, long seasonal cycles or something similar.)