In limma paper in 2015 (https://academic.oup.com/nar/article/43/7/e47/2414268), it stated that it accommodates unequal variances of different treatment arms through two ways under the section titled "Variance models allow for unequal variability". My question is that should we assume Limma deals with unequal variability automatically or should we set some parameter(s) based on estimation using our data? Is there any example if it's the latter case?

I've edited the title of your question because you already know that limma does not require equal variances. The question is what limma does by default and how different group variances should be requested is desired.

The use of weights and the ability to model global parameters allow limma to incorporate unequal variances in a number of ways. One way is through estimating a mean-variance trend, which can either be incorporated into the empirical Bayes procedure as mentioned above or used to generate observation weights (10). A recent development is the ability to estimate precision weights associated with treatment groups or more generally with any given set of covariates. More generally again, the mean-variance trend can be estimated in a treatment-specific way, combining the two types of heteroscedasticity mentioned above.

The mean-variance trend is available as part of eBayes as one might imagine from reading that it can 'be incorporated into the empirical Bayes procedure'. From ?eBayes:

trend: logical, should an intensity-trend be allowed for the prior
variance? Default is that the prior variance is constant.

I've edited the title of your question because you already know that limma does not

requireequal variances. The question is what limma does by default and how different group variances should be requested is desired.