Is there an option to maintain factor levels when save a `stars` object to a tif file?

I have a raster in R as a stars object with an attribute that is a factor. When I save the raster to a tif file using write_stars, the factor levels are lost and only the integers are stored. I read somewhere that what I need is to create a Raster Attribute Table. Can this be done using stars?

Here is a reprex:

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

m <- matrix(letters[1:9], nrow = 3, ncol = 3)
dim(m) <- c(x = 3, y = 3) # named dim
s <- st_as_stars(m)

s[["A1"]] <- factor(s[["A1"]])

plot(s)

write_stars(s, "reprex.tif")

^{Created on 2024-01-31 with reprex v2.1.0}

And this is how the legend of the tif looks in QGIS:

At the end it was not an R issue but a QGIS issue. The information is read by QGIS as it is stored in the xml:

<PAMDataset>
  <PAMRasterBand band="1">
    <CategoryNames>
      <Category></Category>
      <Category>a</Category>
      <Category>b</Category>
      <Category>c</Category>
      <Category>d</Category>
      <Category>e</Category>
      <Category>f</Category>
      <Category>g</Category>
      <Category>h</Category>
      <Category>i</Category>
    </CategoryNames>
  </PAMRasterBand>
</PAMDataset>

The only problem is that QGIS does not have support to visualize this in the legend, see this issue. Here a screenshot of this information in the properties of the raster in QGIS:

Currently, it seems that the only solution is using some plugins or some workarounds with python...