Singular fit with simplest random structure in lmer (lme4), is a Bayesian approach the only option?

Question

I'm running a mixed model with the lmer function from the lme4 package in R and ran into some issues with singular fits. I get the warning message 'singular fit', which is specified as 'convergence code 0; 1 optimizer warnings; 0 lme4 warnings' in summary(model).

Just as background information, I recently re-installed my computer (MacOS) and I can't remember that I had the same issue before. Thus, something might have changed, which could explain that the issue emerged (see also this thread here, where somebody only got the message after some unknown changes, possibly updates).

Based on this thread here, I first thought that the random effect structure might be too complex for my dataset and simplified the model. However, I still got the same warning after building a model only including the random intercept for subject and nothing else: lmer(y ~ 1 + (1|Subject))

Therefore, I tried to reconstruct my dataset by only using the means and standard deviations from all subjects, and then simulate new data based on these values. My idea was to test whether the issue arises in similar datasets, and I also wanted to have a dataset that is easier to share here.

Running the same model as above on some of these new (partly random) datasets, I still have the same issue with the singular fit. However, I don't get it all the time (see below).

Therefore, I wanted to ask if anybody could give me a good explanation for this phenomenon.

Reproducible Example:

First, I defined a function to get n normally distributed, random values per subject, based on the mean and sd provided for each subject.

get_subject_df <- function(means, sds, samples_per_ind){
  ind_vars = vector(length = length(means) * samples_per_ind)  
  for(i in 1:length(means)){
    vals <- rnorm(n = samples_per_ind, mean = means[[i]], sd = sds[[i]])
    irows <-  ((i - 1) * samples_per_ind + 1):(i*samples_per_ind)
    ind_vars[irows] <- vals
    names(ind_vars)[irows] <- paste0("ID_", i)
  }
  df = data.frame(ID = names(ind_vars), vals = ind_vars)
  return(df)
}

Here, I create lists with means and sds for all subjects based on real data

subject_means <- c(9.70,  9.89,  9.72, 10.66,  9.83, 10.32, 10.35, 10.93, 14.36, 12.38, 
                   10.91, 10.38,  8.71, 10.21, 12.28, 11.42, 11.90, 11.41,  9.69,  
                   9.10, 10.20, 10.49, 10.11,  9.16, 10.74)
subject_sds   <- c(3.28, 3.12, 4.91, 4.49, 3.08, 3.67, 4.08, 3.99, 6.10, 5.32, 3.04, 
                   4.69, 3.54, 4.63, 5.47, 3.46, 6.05, 5.94, 3.27, 4.18, 3.41, 4.21, 
                   4.59, 3.28, 4.15)

I ran the model y ~ 1 + (1|Subject) several times using different values for set.seed() and singular fit issues arises with some but not all models. (E.g. when the value for set.seed() is 2, 4, 8, there is an issue but not with 1, 3, 5, ....). Here's a model with the issue (note that the sd for the random intercept for ID is 0.)

set.seed(2)
df_1 <- get_subject_df(means=subject_means, sds=subject_sds, samples_per_ind=5)
lme4::lmer(vals ~ 1 + (1|ID), data = df_1)

singular fit 
Linear mixed model fit by REML ['lmerMod'] 
Formula: vals ~ 1 + (1 | ID)   
Data: df_1 

REML criterion at convergence: 751.7272 

Random effects:
  Groups   Name        Std.Dev.
  ID       (Intercept) 0.000
  Residual             4.917   
 Number of obs: 125, groups:  ID, 25 

Fixed Effects: 
 (Intercept)  
      10.64 

convergence code 0; 1 optimizer warnings; 0 lme4 warnings

Models where values for all subjects have the same standard deviation (here 2) don't seem to have that problem (I tried all values between 1 and 1000 in set.seed(), and never had an issue with singular fit).

set.seed(1)
df_2 <- get_subject_df(means=subject_means, sds=rep(2, length(subject_means)),
                         samples_per_ind=5)
lme4::lmer(vals ~ 1 + (1|ID), data=df_2)

Linear mixed model fit by REML ['lmerMod']
Formula: vals ~ 1 + (1 | ID)
   Data: df_2
REML criterion at convergence: 524.5715
Random effects:
 Groups   Name        Std.Dev.
 ID       (Intercept) 1.134   
 Residual             1.766   
Number of obs: 125, groups:  ID, 25
Fixed Effects:
(Intercept)  
      10.81  

---

Edit:
I went through the links provided by Ben Bolker in the comments and also through the other sources that were linked within these threads. If I understand it correctly, there are the following (potential) options to deal with such singular models (as summarized here):

1. Simplify the random structure. Not possible in my case as I already used the simplest random effect structure that is possible and an intercept-only model.
2. Fit the random effect as fixed effect. However, it seems that this suggestion applies to scenarios where the number of random-effect levels is very low (<5), which is not the case here (25 levels).
3. Use a Bayesian approach.

Thus, my open questions is whether my only option is to use a Bayesian approach?