set_index(..., append=True) act as with append=False with 'Dimensions without coordinates'

What happened:

I get into this strange behaviour when trying to recreate a stacked (MultiIndex) coordinate using set_index(...,append=True).

Since it is not possible to save Dataset to netCDF or Zarr containing stacked / MultiIndex coordinates, before writing to disk I used reset_index(<stacked_coordinate>). When reading such data, I need to use set_index(.., append=True) to recreate such stacked coordinate.

What you expected to happen:

I would expect that set_index(..., append=True) would recreate the MultiIndex stacked coordinate. However, this does not occur if the dimension coordinate specified within set_index() is a 'dimension without coordinate'.
In such situation, set_index(..., append=True) behaves as set_index(, append=False).

Minimal Complete Verifiable Example:

import xarray as xr 
import numpy as np

### Create Datasets 
arr1 = np.random.rand(4, 5).reshape(4,5)
arr2 = np.random.rand(4, 5).reshape(4,5)
da1 = xr.DataArray(arr1, 
                   dims=['nodes','time'],
                   coords={"time": [1,2,3,4,5],
                           "nodes": [1,2,3,4]},
                   name='var1')
da2 = xr.DataArray(arr2, 
                   dims=['nodes','time'],
                   coords={"time": [1,2,3,4,5],
                           "nodes": [1,2,3,4]},
                   name='var2')
ds_unstacked = xr.Dataset({'var1':da1,'var2':da2})
print(ds_unstacked) 

# - Stack variables across a new dimension 
da_stacked = ds_unstacked.to_stacked_array(new_dim="variables", variable_dim='variable', 
                                           sample_dims=['nodes','time'], name="Stacked_Variables") 
ds_stacked = da_stacked.to_dataset()

# - Look at the stacked MultiIndex coordinate 'variables'
print(ds_stacked) 
print(da_stacked.variables.indexes)

### Remove MultiIndex (to save Dataset to netCDF/Zarr, ...)
ds_stacked_disk = ds_stacked.reset_index('variables')
print(ds_stacked_disk)
 
### Try to recreate MultiIndex  
print(ds_stacked_disk.set_index(variables=['variable'], append=False)) # GOOD ! Replace 'variable'  coordinate with 'variables'
print(ds_stacked_disk.set_index(variables=['variable'], append=True)) # BUG ! Do not create the expected MultiIndex !  

### Current workaround to obtain a MultiIndex stacked coordinate
tmp_ds = ds_stacked_disk.assign_coords(variables=(np.arange(0,2)))
ds_stacked1 = tmp_ds.set_index(variables=['variable'], append=True)   
print(ds_stacked1)  # But with level 0 - 'variables_level_0'

### Unstack back  
# - If the BUG is solved, no need to specify the level argument  
ds_stacked1['Stacked_Variables'].to_unstacked_dataset(dim='variables', level='variable')

Environment:

Start by running the minimal example and inspect Dataset.set_index with append=True after reset_index('variables'). Trace the set_index path for a dimension without a coordinate. Done means append=True recreates the expected MultiIndex, so the stacked data can be unstacked without the temporary coordinate workaround.