Alexis Devilez

Nicolas Bonod
email: nicolas.bonod@fresnel.fr

Brian Stout

Davy Gérard

Jérôme Wenger

Hervé Rigneault

Evgeny Popov

Three-dimensional subwavelength confinement of light with dielectric microspheres

Keywords:

Dielectric microspheres are shown to be capable of confining light in a three-dimensional region of subwavelength dimensions when they are illuminated by tightly focused Gaussian beams. We show that a simple configuration, not involving resonances, permits one to reach an effective volume as small as 0.6 (λ/n) 3 . It is shown that this threedimensional confinement arises from interferences between the field scattered by the sphere and the high angular components of the incident Gaussian beam passing aside the sphere.

Introduction

Strong concentration of light in a single spot of subwavelength dimensions is widely performed by metallic structures such as gratings, pinholes, tips or nanoparticles [START_REF] Sentenac | Subdiffraction Light Focusing on a Grating Substrate[END_REF][START_REF] Rigneault | Enhancement of single-molecule fluorescence detection in subwavelength apertures[END_REF][START_REF] Anger | Enhancement and Quenching of Single-Molecule Fluorescence[END_REF][START_REF] Nie | Probing single molecules and single nanoparticles by surface-enhanced Raman scattering[END_REF], which take advantage of electromagnetic resonances on metals such as surface plasmon modes. However, losses in metals and technically challenging nanofabrication processes may limit their interest for applications in ultramicroscopy, spectroscopy and optical data storage. Resonant dielectric structures such as gratings or photonic crystal cavities offer an alternative with almost no intrinsic absorption [START_REF] Noda | Seeking the ultimate nanolaser[END_REF][START_REF] Laroche | Optical resonances in one-dimensional dielectric nanorod arrays: field-induced fluorescence enhancement[END_REF][START_REF] Marinica | Bound states in the continuum in photonics[END_REF]. The quality factors and field enhancement of these resonance structures vary from tens to thousands or more.

In the search for simpler dielectric structures to concentrate light at the nanoscale, several groups have considered the use of dielectric spheres of micrometer dimensions [START_REF] Chen | Photonic nanojet enhancement of backscattering of light by nanoparticles: A potential novel visible-light ultramicroscopy technique[END_REF][START_REF] Li | Optical analysis of nanoparticles via enhanced backscattering facilitated by 3-D photonic nanojets[END_REF][START_REF] Lecler | Properties of a 3D photonic jet[END_REF][START_REF] Itagi | Optics of photonic nanojets[END_REF][START_REF] Kofler | Axially symmetric focusing as a cuspoid diffraction catastrophe: scalar and vector cases and comparison with the theory of Mie[END_REF][START_REF] Ferrand | Direct imaging of photonic nanojets[END_REF][START_REF] Devilez | Spectral analysis of three-dimensional photonic jets[END_REF]. When the microsphere is illuminated by a plane wave, the so-called "photonic nanojet" beam that emerges from the sphere has subwavelength transverse dimensions and low divergence, which makes it fruitful for applications in dry laser cleaning [START_REF] Mosbacher | Optical field enhancement effects in laser-assisted particle removal[END_REF][START_REF] Luk'yanchuk | Three-dimensional effects in dry laser cleaning[END_REF], nanopatterning [START_REF] Piglmayer | Laser-induced surface patterning by means of microspheres[END_REF][START_REF] Mcleod | Subwavelength direct-write nanopatterning using optically trapped microspheres[END_REF][START_REF] Pereira | Laser-fabricated porous alumina membranes for the preparation of metal nanodot arrays[END_REF], Raman spectroscopy [START_REF] Yi | Enhanced Raman scattering by self-assembled silica spherical microparticles[END_REF][START_REF] Kasim | Near-field Raman imaging using optically trapped dielectric microsphere[END_REF] and optical data storage [START_REF] Kong | Photonic nanojet-enabled optical data storage[END_REF]. However, due to its large dimension along the optical axis (typically 2-3 µm), the photonic nanojet does not provide three-dimensional subwavelength light confinement. For applications requiring high transverse and longitudinal resolutions, the classical photonic nanojet is unsuitable, as performs no better than the focusing obtained from a classical microscope objective with a high numerical aperture.

Quite surprisingly, it has recently been shown that a single microsphere illuminated by a tightly focused Gaussian beam can outperform classical microscope systems and significantly enhance the fluorescence emission from a single molecule [START_REF] Gérard | Strong electromagnetic confinement near dielectric microspheres to enhance single-molecule fluorescence[END_REF]. In that case, strong confinement of light, on the order of (λ/n) 3 , with a non-resonant dielectric structure was clearly demonstrated, but the physical origin of this effect and its implications remained untreated.

In this paper, we investigate the light confinement produced by a dielectric microsphere illuminated by a tightly focused Gaussian beam. Our simulations employ a rigorous Lorentz-Mie theory [START_REF] Stout | Light diffraction by three-dimensional object: differential theory[END_REF], and concern the experimentally relevant configuration of a 2 µm diameter latex sphere of refractive index n s = 1.6 surrounded by water (n = 1.33) as used in Ref. [START_REF] Gérard | Strong electromagnetic confinement near dielectric microspheres to enhance single-molecule fluorescence[END_REF]. The incident beam parameters are chosen to approximate the properties of a perfectly corrected objective with high numerical aperture NA ≈ 1.1: vacuum wavelength λ = 633 nm, 300 nm transverse waist (half width at 1/e²), and 1.25 µm longitudinal half width at 1/e². A circular polarization is chosen to simplify the numerical simulations. This beam is simulated using the first-order Davis coefficients [START_REF] Davis | Theory of electromagnetic beams[END_REF], the translation-addition theorem acting on the incident beam coefficients is used to tune the position d of the incident beam focus with respect to the sphere center [START_REF] Stout | A transfert matrix approach to local field calculations in multiplescattering problems[END_REF].

Three dimensional subwavelength confinement of light

Figure 1(a) presents the electric field intensity map when the incident beam focus is at the position resulting in the optimal concentration of light and strongest field intensity, denoted I max , behind the sphere. This corresponds to setting the d parameter to 1.62 µm. In this case, a central spot with subwavelength dimensions along the three spatial directions is present close to the sphere surface. Its intensity is one order of magnitude higher than the intensity in the other lobes (note the logarithmic scale in Fig. 1). The thin black line represents the I max /e² intensity contour, which is a useful guide when comparing the intensity map with the microsphere (Fig. 1(a)) to the incident Gaussian beam without the sphere (Fig. 1(b)).

Table 1 compares the transverse and longitudinal waist defined at I max /e² and denoted respectively w xy and w z for the incident beam and the beam focused by the microsphere for d = 1.62 µm. Let us remark that with the add of the microsphere, the longitudinal waist is no more symmetric with respect to the I max position and that the longitudinal waist is calculated outside the sphere, after the I max position. The incident field is further confined by the sphere both longitudinally and transversally. It is apparent that the longitudinal modification of the beam is more spectacular since the maximum intensity has both been enhanced and moved toward the sphere surface. The effective volume defined by π 3/2 w xy 2 w z /2 is reduced by one order of magnitude, and is approximately 0.6 (λ/n) 3 . It must be stressed that this light confinement permits to decrease by one order of magnitude the number of probed molecules in fluorescence spectroscopy [START_REF] Gérard | Strong electromagnetic confinement near dielectric microspheres to enhance single-molecule fluorescence[END_REF] and that this strong confinement of light is obtained with a very simple structure not involving resonances.

Incident beam

Beam + microsphere Ratio to incident beam Transverse waist w xy 300 nm 230 nm 77 %

Longitudinal waist w z 1250 nm 420 nm 33 %

Effective volume V 0.6 µm 3 ≈ 6 (λ/n) 3 0.06 µm 3 ≈ 0.6 (λ/n) 3 10 % Table 1. Summary of the characteristics widths at I max /e² corresponding to the intensity maps displayed in Fig. 1. The volume is derived for the incident beam as V = π 3/2 w xy 2 w z , and for the focus by the microsphere V = π 3/2 w xy 2 w z /2 (let us recall that only the beam exiting the microsphere is considered here).

In order to understand the physical mechanism of this confinement, we have employed a simulation method that allows the calculation of the incident and scattered field separately [START_REF] Stout | Light diffraction by three-dimensional object: differential theory[END_REF]. Figure 1(c) presents the map of the scattered field intensity. It shows an elongated region with high intensity in the shadow side of the sphere. This scattered field is very close to a classical "photonic nanojet" beam: it has both narrow lateral extent and large longitudinal extent. Only the coherent sum of the incident beam plus the scattered field is able to reach the strong confinement in the three directions (Fig. 1(a)), which highlights the relation between interference and the subwavelength confinement.

Discussion

To further understand the phenomenon, two numerical experiments are undertaken. First, the incident focused Gaussian beam is replaced by a plane wave: the incident angular spectrum being reduced to zero. In that case, in contrast to Gaussian beams, the maps of the scattered and total field intensities are very similar and present a large longitudinal extent (Fig. 2), which leads to an effective focal volume of xxx(λ/n) 3 . We conclude that non-null angular contributions present in the incident beam are needed to reach an all 3axis subwavelength confinement. The effective volume is displayed in Fig. 3 as a function of the numerical aperture. Let us note that, for each numerical aperture value, the distance d between the center of the sphere and the focus of the incident beam has to be tuned to optimize the maximum of intensity I max behind the sphere. The effective volume decreases as a function of the numerical aperture and a numerical aperture above ~0.8 is needed to obtain an effective volume below (λ/n) 3 . It results that the 3-dimensional confinement directly depends on the high angular contributions present in the incident highly focused Gaussian beam, and that a confinement below (λ/n) 3 can be achieved with a numerical aperture higher than 0.8. A second numerical experiment is carried out by tuning the distance d between the center of the sphere and the waist of the incident beam for an numerical aperture equal to 1. The movie of Fig. 4 illustrates the behavior of the high intensity region behind the sphere as d is decreased from 7 to 0 µm. In the first frame, d is set at 7 µm and the focalization of the incident beam can be clearly observed at the center of the screen. Then, as the sphere approaches the focalization area, the field is fully affected by the presence of the sphere: the intensity and the volume of the incident spot decrease to the benefit of the field focalized by the microsphere. When d equals ~2.2 µm, one can observe a dark spot due to a destructive interference between incident and scattered fields, which creates the longitudinal confinement. The near field behind the sphere is then highly enhanced. When the sphere is further approached to the focalization area, interference ceases to play an important role since there are no longer incident field components passing aside the sphere towards the region just behind it. Only one spot, arising from the collimation of the incident field by the sphere is then observed, but the maximum of intensity is located inside the sphere and presents no interest for microscopy applications.

Conclusion

This study demonstrates that the three-dimensional volume reduction observed when illuminating a microsphere with a tightly focused Gaussian beam is due the interplay of two different contributions: a well-known collimation of the incident field by the microsphere, combined with interferences between the field focused by the dielectric microsphere and the incident field propagating aside the sphere. It has been shown that a confinement behind the sphere below (λ/n) 3 requires a tightly focused beam with a numerical aperture higher than 0.8 together with a proper focusing of the incident beam respectively to the sphere. A further increase of the numerical aperture up to 1 permits a strong confinement of light behind the sphere down to 0.6 (λ/n) 3 . It is important to remark that this property was obtained using only dielectric materials and without invoking resonance phenomenon. Its utility has already been demonstrated in enhancing the fluorescence signal of single-molecules [START_REF] Gérard | Strong electromagnetic confinement near dielectric microspheres to enhance single-molecule fluorescence[END_REF] and should prove of considerable interest in Raman spectroscopy, laser nano-patterning and microscopy.

Acknowlegments

This work has been funded by the grant PEPS "NANODRILL" of the Centre National de la Recherche Scientifique

Fig. 1 .

 1 Fig. 1. (a) Total electric field intensity map for d = 1.62 µm. The white circle represents the microsphere section, the thin black line represents the I max /e² intensity contour. (b) and (c) display the intensity maps of the incident and scattered fields. The coherent summation of these two fields leads to the total intensity presented in (a).

Fig. 2 .Fig. 3 .

 23 Fig.2. Electric field intensity maps of (a) the scattered field and (b) the total field in the same conditions than previously but the sphere is illuminated by a plane wave. The incident intensity is normalized to be unitary per surface unit.

Fig. 4 .

 4 Fig. 4. (2000KB) Movie of the electric field intensity map for d decreasing from d = 7 µm to d = 0 µm and N.A. ≈ 1. The white circle represents the sphere cross-section.


Keywords: Integration, intermittency, unit commitment, optimization techniques, priority list, genetic algorithm. I

Integrating renewable energy generators with existing electric power systems is essential in developing renewable energy generation in electric power systems in islands and remote areas. The high penetration of renewable energy sources (RES) and the possibility of working in a connected mode with the existing electricity grid are significant problems that must be overcome. With the many RES implementations, a proper integration method is needed between RES plants, which have intermittent characteristics, and existing power generation systems so that the electrical system remains stable, continuous, and reliable. One of the power plant integration schemes is unit commitment (UC), a daily scheduling method to achieve optimization goals while still meeting various operational, environmental, and technical constraints. UC method development is needed to maintain system stability, resilience, and minimize operating costs, especially for small systems for electrical development. Therefore, it is necessary to develop a new optimization method in the UC scheme for microgrid power systems.

This paper discusses the development of an optimization method to integrate power generation from intermittent renewable energy using a UC scheme in a microgrid system. The algorithm used is a combination of improved priority list techniques and genetic algorithms, where the improved priority list is run first for scheduling optimization and then a combination of the algorithm priority list and genetic algorithm for economic dispatch solutions to obtain the lowest operational costs. The generating unit's capacity to operate per capability segment is also taken into account in the operation scheduling and affects the number of iterations in the genetic algorithm's operation. This method has been simulated in the electrical system in Eastern Indonesia, which is a developing power system and has an intermittent characteristic of RES power plants.

Tenggara has increased significantly. For this reason, the government of the Republic of Indonesia has implemented an electricity system development program by utilizing all available resources. One of the agenda is to develop and utilize local energy potential using renewable energy sources (RES) [START_REF]Indonesia Electricity Supply Business Plan 2016-2025[END_REF]. That is needed in many regions and islands in Indonesia that require electricity at an affordable price and low carbon intensity [START_REF] Banjar-Nahor | Study on Renewable Penetration Limits in a Typical Indonesian Islanded Microgrid Considering the Impact of Variable Renewables Integration and the Empowering Flexibility on Grid Stability[END_REF][START_REF] Zafeiratou | Sustainable island power system -Scenario analysis for Crete under the energy trilemma index[END_REF]. Efforts to provide electricity infrastructure in developing countries like Indonesia, especially for islands and remote areas, are generally carried out by building a microgrid system with diesel generators [START_REF] Rendroyoko | Optimizing Unit Commitment Schemes for Variable RES Power Plant Integration in Microgrid Systems[END_REF][START_REF] Irena | Transforming small-island power systems: Technical planning studies for the integration[END_REF]. A microgrid system consists of several generator units, electrical loads, and other supporting equipment to manage operations independently by an intelligent control system [START_REF] Hand | Renewable Electricity Futures Study (Entire Report)[END_REF][START_REF] Fossati | Unit Commitment and Economic Dispatch in Micro grid[END_REF]. Generally, most island microgrid power systems rely on diesel generators for their electricity supply. These generators are usually very flexible, with high lean capability, short start-up/shutdown times, and low minimum technicality, but their operating costs tend to be high [START_REF] Sigrist | Island Power Systems[END_REF][START_REF] Crossa | Benchmarking island power systems: Results, challenges, and solutions for Long Term Sustainability[END_REF]. However, nowadays, along with the development of RES generation technology and the incessant reduction in fuel oil consumption programs, the use of RES has been intensified. In operation, renewable energy power plants are usually interconnected with other generating units in the existing microgrid power system. Moreover, in some cases, the RES generator can cause instability in the power system due to its intermittency characteristics [START_REF] Rendroyoko | Optimizing Unit Commitment Schemes for Variable RES Power Plant Integration in Microgrid Systems[END_REF]. As the RES penetration of many microgrid systems increases in the future, the need to study its impact on power system operation will be increasingly needed. Therefore, a strategy is required so that the power system to be robust and flexible enough to deal with RES generators' uncertainty and fluctuating power supply [START_REF] Jo | Improved Genetic Algorithm-Based Unit Commitment Considering Uncertainty Integration Method[END_REF].

One method that can be used to create operational flexibility for the power system to vary the supply of RES generators is the unit commitment (UC) method [START_REF] Kazarlis | A genetic algorithm solution to unit commitment problem[END_REF][START_REF] Salam | Unit Commitment Solution Methods[END_REF][START_REF] Anjos | Unit Commitment in Electric Energy System[END_REF]. The UC method regulates the electric power system's operation by determining the operation and production schedule of electricity for each Vincent Debusschere is an Assoc. Professor at the Université Grenoble Alpes, CNRS, Grenoble INP, G2Elab, France (e-mail: vincent.debusschere@grenoble-inp.fr).

Deddy P Koesrindartoto is an Assistant Professor at the School of Business and Management, Institute Technology of Bandung, Bandung, Indonesia (email: deddypri@sbm-itb.ac.id)

A

generating unit in the electric power system, including the RES, to meet electricity load needs with minimal costs. In its operation, UC must meet several system demands, reserve requirements, and other limitations optimally and keep operating costs low. On the other side, changes in the power system load also cause variability and uncertainty, impacting UC settings. Considering the conditions of variability and uncertainty that may occur on the generation side and the load side of the power system, UC must be required to regulate and ensure the availability of continuous, and reliable generators to meet the load requirements during the operating period.

The problem of UC in microgrids systems is quite different from large electric power systems [START_REF] Fossati | Unit Commitment and Economic Dispatch in Micro grid[END_REF] [START_REF] Liang | Unit commitment in microgrids by improved genetic algorithm[END_REF]. Microgrids usually consist of thermal generators, diesel engines, micro-hydro, and some variable renewable sources such as wind turbines (WT) and solar photovoltaics (PV). The uncertainty of the variable RES generator's primary energy availability can cause instability in the microgrid power system. The variable RES generator characteristics are very different from conventional generators, which will complicate UC problems' resolution.

The UC scheme requires a technique to optimize the problems that arise between objective functions and constraints. The optimization technique is the process of minimizing or maximizing an objective function while still fulfilling the constraints. The main objective is to minimize total operating costs. Several optimization techniques have been developed to find a suitable solution for implementing the UC: Priority list, Dynamic programming, Lagrangian relaxation, Genetic algorithm (GA), etc [START_REF] Kazarlis | A genetic algorithm solution to unit commitment problem[END_REF][START_REF] Swarup | Unit Commitment Solution Methodology Using Genetic Algorithm[END_REF][START_REF] Abujarad | Recent approaches of unit commitment in the presence of intermittent renewable energy resources: A review[END_REF]. From several previous studies on the implementation of UC, many researchers use genetic algorithms combined with other techniques to solve UC problems [START_REF] Sarjiya | Unit Commitment Solution Using Genetic Algorithm Based on Priority List Approach[END_REF][START_REF] Madraswala | Genetic Algorithm Solution to Unit Commitment Problem[END_REF]. The GA technique has been cited as the most successful evolutionary programming method, and it is suitable for optimizations with wide variations and uncertainties, such as the intermittency characteristics of RES generators. The GA technique is used because it is flexible, has the right convergence level, and is characterized by multi-search-nosingle-point.

This paper is focused on discussing the use of the UC method in the integration of the variable RE generator into the existing microgrid power system. The technique used in this UC implementation solution is a GA, which is intended to overcome the considerable variability of load changes and the uncertainty of the RES generator output fluctuation in the microgrid system. To improve GA implementation, PL techniques are used to help determine the initial population [START_REF] Fossati | Unit Commitment and Economic Dispatch in Micro grid[END_REF][START_REF] Sarjiya | Unit Commitment Solution Using Genetic Algorithm Based on Priority List Approach[END_REF]. Several parameters and constraints have been applied in this algorithm, allowing GA to avoid local optima traps and accelerate solution convergence. For this purpose, a microgrid system simulation model was built based on the Timor electric power system to test the new algorithm as a sample reference for the development of UC for islands microgrid power system. This paper is structured as follows. The Concept of unit commitment is described in Section 2, including renewable energy conditioning and optimization algorithms for unit commitment. Section 3 presents the application of the unit commitment in the microgrid system. The UC procedures based on the PL and GA algorithms are described in Section 4. The simulation studies results are described in Section 5, along with an analysis of computation time, reserve requirements, and operating costs. Conclusions are given in Section 6.

II. THE CONCEPT OF UNIT COMMITMENT (UC)

A. Unit commitment UC is a method of scheduling power plant units' daily operation to fulfill load requirements while meeting various operational, environmental, and technical limitations [START_REF] Kazarlis | A genetic algorithm solution to unit commitment problem[END_REF][START_REF] Anjos | Unit Commitment in Electric Energy System[END_REF]. Optimization of the UC implementation is subject to several global constraints such as power balance, spinning reserves, and several individual constraints such as unit output limits, maximum ramp rates, and minimum up/down hours. The UC scheme in power systems with RES variable generation is shown in Fig. 1 below. In microgrid power systems, the main objective of UC implementation is to minimize the total generation cost, denoted by TOC [START_REF] Jo | Improved Genetic Algorithm-Based Unit Commitment Considering Uncertainty Integration Method[END_REF]. On/off status of the i unit at time t, if the unit is on then Uit = 1, and if the unit is off then Uit = 0, Fit(Pit) fuel cost function of the i-th unit, with output of Pit generator, at the t-time, with i-th unit start-up cost at t-hour, where ai, bi, ci coefficient of unit fuel cost, Pit, generation output from the i-unit at t-hour,

𝑇𝑂𝐶 = ∑ ∑ ∑(𝛼

In UC, several objectives can be formulated as follows: Total operation cost (TOC):

𝑀𝑖𝑛 𝑇𝑂𝐶 = ∑ ∑(𝑈 𝑖𝑡 . 𝐹 𝑖𝑡 (𝑃 𝑖𝑡 ) + 𝑈 𝑖𝑡 . 𝑆 𝑖𝑡 ) 𝑁 𝑖=1 𝑇 𝑡=1 (2)
In equation ( 2) above, there are two functions used; first is the fuel cost per MW, which can be formulated as follows:

𝐹 𝑖𝑡 (𝑃 𝑖𝑡 ) = 𝑎 𝑖 𝑃 𝑖𝑡 2 + 𝑏 𝑖 𝑃 𝑖𝑡 + 𝑐 𝑖 (3)
And the second is start up cost of generator unit :

𝑆 𝑖𝑡 =  𝐻𝑆𝐶 𝑖 , 𝑖𝑓 𝑀𝐷𝑇 𝑖 < 𝑇 𝑜𝑓𝑓𝑖 ≤ 𝑀𝐷𝑇 𝑖 + 𝑇 𝑐𝑜 𝐶𝑆𝐶 𝑖 , 𝑖𝑓 𝑇 𝑑𝑜𝑤𝑛,𝑖 < 𝑀𝐷𝑇 𝑖 + 𝑇 𝑐𝑜𝑙𝑑,𝑖 (4)
In UC implementation, besides fulfilling the main objectives above, some constraints must still be fulfilled as follows:

A.1. Power balance In this constraint, the power of generation should always keep equal to the customer's electricity load. The mathematical equation can be seen in ( 5).

∑ ∑ 𝑃 𝑡,𝑛,𝑘 𝑘 𝑛 = 𝐷 𝑡 (5)

A.2. Spinning reserve (SR)

The total generation must be greater than the load plus the spinning reserve needed to guarantee the stability of the power system. The total SR must be larger than 10% of demand.

∑ ∑(𝑈 𝑛 . 𝑉 𝑡,𝑛 -𝑃 𝑡,𝑛 ) ≥ 10% 𝐷 𝑡 𝑛 𝑡 (6) ∑(𝑈 𝑖𝑡 . 𝑃 𝑚𝑎𝑘𝑠 𝑖 ) ≥ (𝑙𝑜𝑎𝑑 𝑡 + 𝑆𝑅 𝑡 ), 1 ≤ 𝑡 ≤ 𝑇 𝑁 𝑖=1 (7)
A. [START_REF] Zafeiratou | Sustainable island power system -Scenario analysis for Crete under the energy trilemma index[END_REF]

𝑇 𝑜𝑓𝑓 (𝑖, 𝑡) ≥ 𝑀𝐷𝑇(𝑖), 𝑖𝑓 𝑈(𝑖, +1) 𝑚𝑢𝑠𝑡 𝑏𝑒 𝑜𝑛

A.5. Ramp up/downtime This constraint is the generating unit's time limit to be ready for full operation (ramp-up) or time to go inoperative (ramp down).

𝑃 𝑖,𝑡 -𝑃 𝑖,𝑡-1 ≤ 𝑈𝑅 𝑖 , 𝑖𝑓 𝑝𝑜𝑤𝑒𝑟 𝑖𝑛𝑐𝑟𝑒𝑎𝑠𝑒 (11) 𝑃 𝑖,𝑡-1 -𝑃 𝑖,𝑡 ≤ 𝐷𝑅 𝑖 , 𝑖𝑓 𝑝𝑜𝑤𝑒𝑟 𝑑𝑒𝑐𝑟𝑒𝑎𝑠𝑒. (12)
A.6. Must run units According to the forecasted output, the RES power plants are treated as a must-run unit and directly calculated into the power supply mechanism for the power system.

𝑃 𝑅𝐸 (𝑡) = 𝑃 𝑜𝑢𝑡𝑝𝑢𝑡 (𝑡) 0 ≤ 𝑡 ≤ 24 (13)
A.7. System's Frequency and voltage Limits of system frequency and voltage are :

𝑉 𝑚𝑖𝑛 ≤ 𝑉 ≤ 𝑉 𝑚𝑎𝑥 (14) 𝐹 𝑚𝑖𝑛 ≤ 𝐹 ≤ 𝐹 𝑚𝑎𝑥 (15)
B. UC with intermittent RES power plant Electric power systems with renewable energy sources (RES), such as solar PV and wind power plants, the operating arrangements of each generating unit are formulated in the form of constraints as above, adjusted to its primary energy availability. The constraints are consist of generator power limit, minimum up/downtime, ramp up/downtime, and unit restriction status. Regarding the limitations and characteristics of the availability of its primary energy sources, RES power plants are classified as must-run units. Solar photovoltaic power plants are characterized by sharp temporal fluctuations, which can have considerable variability.

In preparation for the implementation of UC, the RES power plant must also be considered according to its availability. For intermittent RES, weather forecasting, historical data, and analysis are critical in preparing the generation schedule. With accurate weather forecasting, the UC of the RES generating unit can be performed with proper accuracy. In calculating the total operation cost, we will also calculate the solar PV cost according to the estimated selling price per kWh. The output of solar PV can be presented:

𝑃 𝑉 = 0, 𝐺 < 𝐺 𝑖 𝜇. 𝐺, 𝐺 𝑖𝑛 ≪ 𝐺 < 𝐺 𝑟𝑎𝑡𝑒 𝑃 𝑟𝑎𝑡𝑒 , 𝐺 ≫ 𝐺 𝑟𝑎𝑡𝑒 ( 16
)
Where Pv is the power output of solar PV, Gin and Grate denote the minimum and rated insolation, and μ is a coefficient of the solar panel.

The output from the solar PV can be included in the net load calculation concept. The net load concept simplified the formulation of the UC problem [START_REF] Jo | Improved Genetic Algorithm-Based Unit Commitment Considering Uncertainty Integration Method[END_REF].

𝑁𝑒𝑡 𝑙𝑜𝑎𝑑 = 𝑙𝑜𝑎𝑑 -𝑃 𝑉 [START_REF] Sarjiya | Unit Commitment Solution Using Genetic Algorithm Based on Priority List Approach[END_REF] Due to the energy sources' intermittency, it is necessary to predetermine a schedule that meets the netload. Based on the predicted values, hourly random loads and solar power are considered in each scenario.

C. Optimization Method for Unit Commitment

UC has two main interrelated parts: scheduling the order of operation of the generating unit and economic dispatch. Here, an optimization technique is needed to manage the UC in order to produce the best results. The objective function of the UC is to minimize the total operating cost while meeting the constraints: power balance, spinning reserve, power generation limit, maximum ramp rate, and minimum up/ downtime.

The optimization techniques used to solve the UC can be divided into classical, stochastic, and hybrid techniques. Moreover, several optimization techniques have been developed to find suitable solutions for implementing the unit commitment, such as priority list, dynamic programming, Lagrangian relaxation, genetic algorithm, etc. GA is a population-based search algorithm that follows the principle of "survival of the fittest" and is a stochastic technique that can be used to provide solutions to UC problems. UC is a non-convex, large-scale, non-linear, and mix-integer blend optimization problem, with several limitations. The non-convex condition is due to the binary on/off condition of the generating unit, and the non-linear characteristic comes from the non-linear generation cost curve. Therefore, the GA technique is believed to be a suitable solution for large non-linear engineering optimization problems such as UC. GA continues to be developed by combining it with other techniques for increased performance and accuracy.

Priority list (PL) techniques have been implemented in many existing optimization methods to solve the UC problem [START_REF] Delarue | Enhanced priority list unit commitment method for power systems with a high share of renewables[END_REF][START_REF] Tingfang | Methodological Priority List for Unit Commitment Problem[END_REF]. The method is simple and fast but produces less optimal solutions. The PL algorithm can be easily programmed and is suitable for handling UC problems involving multiple generators, and Delarue (2017), in his paper, proves that PL works well in electric power systems with intermittent RES generators under low load conditions [START_REF] Sarjiya | Unit Commitment Solution Using Genetic Algorithm Based on Priority List Approach[END_REF]. On the other hand, GAs are characterized by exhaustively explore the search space, thereby ensuring the possibility of finding a close to the optimum solution. As one of the stochastic methods, GA relies on a population of candidate solutions. Therefore, determining the candidate solution population as the initial population will make the GA process to achieve its convergence faster. Furthermore, we can accomplish this after the PL method's scheduling results are used as the initial population of the GA process [START_REF] Fossati | Unit Commitment and Economic Dispatch in Micro grid[END_REF][START_REF] Sarjiya | Unit Commitment Solution Using Genetic Algorithm Based on Priority List Approach[END_REF]. Therefore, the combination of the PL and GA methods could be a new option for UC optimization.

III. OVERVIEW OF TIMOR POWER SYSTEM

The Timor power system is an electrical system that supplies the electric power needs of the Timor island, in East Nusa Tenggara province, Indonesia [START_REF]Indonesia Electricity Supply Business Plan 2016-2025[END_REF]. This electrical system was developed based on a microgrid system with a diesel-powered power plant, strengthened by coal-fired steam power plants, and currently with RES-based power plants. The fuel mix ratio in the Timor electric power system is 65.09%, which means that most of the power plants still use fuel oil, so RES is still very potential. Since 2015, Timor power systems operate 5 MWp Solar PV plants and contribute to electric power production. The single line diagram of the Timor power system can be seen in Fig. 2, and the list of generating units is described in TABLE I. It is necessary to adjust the generating unit operation schedule in daily system operations to compensate for the intermittent RES generator output fluctuation. In this study, we built a simulation model of an electric power system with reference to the electric power system of Timor, East Nusa Tenggara, Indonesia [START_REF] Rendroyoko | Optimizing Unit Commitment Schemes for Variable RES Power Plant Integration in Microgrid Systems[END_REF]. The Timor electric power system was chosen as a case study because the electricity system is rapidly developing and has the most potential for utilizing solar power in Indonesia.

0 0 0 0 IS(jam) -1 -1 -1 -1 - 1
Information: U1, U2: steam generated power units, U3: Solar PV, U4 -U10: Diesel power units In addition to system data, the forecast data for the power output of solar PV plant is shown in PV (kW) 100 0 0 0 0 0 0 0 Fig. 3 Load profile of Timor electricity system [START_REF]Indonesia Electricity Supply Business Plan 2016-2025[END_REF] Timor electric power systems' load profile is shown in Fig. 3, which shows the typical characteristics of electricity use in developing areas with the domination of residential customers and few industrial and business customers with slow growth.

IV. PL & GA ALGORITHM FOR UC OPTIMIZATION IN MICROGRID SYSTEM

This section discusses applying the GA technique to unit commitment problems in power systems with RES. This optimization process for UC begins with the PL technique to get the population initials and continues with optimization using the GA technique. Following are the steps for implementing the PL and GA methods for UC solutions [START_REF] Fossati | Unit Commitment and Economic Dispatch in Micro grid[END_REF][START_REF] Sarjiya | Unit Commitment Solution Using Genetic Algorithm Based on Priority List Approach[END_REF]:

Initialization

The generation of the initial solution is essential, mainly for the UC problem. The initial solution is usually generated at random. The selection of the candidate solutions in this study was carried out using the priority list (PL) technique.

Evaluation

Stage of evaluating the fitness value of the initial population candidate solutions.

Selection

Selection allocates more copies of those solutions with higher fitness values and imposes the fittest mechanism's survival on the candidate solutions. 4. Recombination / crossover combine parts of two or more parent solutions to a new possible better solution. A properly designed recombination mechanism will determine success.

Mutation

The mutation is a genetic operator used to maintain genetic diversity from one generation of a GA chromosome population to the next. The flow chart of the optimization process for implementing UC in the microgrid system is shown in Fig. 4. The novelty of this algorithm is to use PL and GA techniques, with PL is done first for UC to get the initial population, then continue to dispatch with PL + GA, where one individual GA in the initial population stage uses the dispatch value from PL to accelerate the convergence of lowest operating costs.

Replacement Replacement of parent populations with offspring populations made by selection, recombination, and mutation

For solving the UC problem, the first step is to use the priority list (PL) technique to obtain the initial population. The priority order of the generating units in TABLE I can be obtained based on each unit parameter. Generators with a higher capacity and lower heat rate (HR) will receive higher priority. For machines with low capacity and higher HR will have less priority. The PL and calculated heat rate are shown in equations ( 18) and ( 19) [START_REF] Sarjiya | Unit Commitment Solution Using Genetic Algorithm Based on Priority List Approach[END_REF][START_REF] Tingfang | Methodological Priority List for Unit Commitment Problem[END_REF].

𝐻𝑅 𝑖 = 𝐹 𝑖 (𝑃 𝑖𝑚𝑎𝑥 (𝑡)) 𝑃 𝑖𝑚𝑎𝑥 (𝑡) (18)
𝐶𝑜𝑠𝑡 𝑖 = 𝐻𝑅 𝑖 × 𝐹 𝑖 𝑃 𝑖ℎ (19)
The PL method results are shown in Erreur ! Source du renvoi introuvable. (output power in MW) and will be used as the GA's initial population. This operating schedule is obtained from the PL calculation process based on the heat rate to meet the electrical load. 19:00

1 1 1 1 1 1 0 0 0 20:00 1 1 1 1 1 1 0 0 0 21:00 1 1 1 1 1 1 0 0 0 22:00 1 1 1 1 1 0 0 0 0 23:00 1 1 1 1 1 0 0 0 0 00:00 1 1 1 1 0 0 0 0 0
Once the initial population has been created, we must evaluate the fitness of each candidate solution. In the process, all constraints are checked. First, the spinning reserve (SR) limits will be reviewed hourly using the schedule from the initial population. If any SR is breached within one hour, the penalty count will be added to the penalty function. Second, the minimum up/downtime constraints must be met so that individuals can enter the load dispatch stage. If these constraints cannot be met, the individual will be given a calculated penalty.

After the load dispatch stage is completed, each individual will enter the process of calculating the fitness function. The penalty function is the multiplication of the penalty factor and penalty count (20. The penalty function, start-up costs, and fuel costs obtained from the load dispatch stage will be added to get each individual's fitness value in the population (21).

𝑃𝑒𝑛𝑎𝑙𝑡𝑦 𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛 = 𝑝𝑒𝑛𝑎𝑙𝑡𝑦 𝑓𝑎𝑐𝑡𝑜𝑟 × 𝑝𝑒𝑛𝑎𝑙𝑡𝑦 𝑐𝑜𝑢𝑛𝑡 (20)
𝐹𝑖𝑡𝑛𝑒𝑠𝑠 = 1 𝑠𝑡𝑎𝑟𝑡𝑢𝑝 𝑐𝑜𝑠𝑡 + 𝑓𝑢𝑒𝑙 𝑐𝑜𝑠𝑡 + 𝑝𝑒𝑛𝑎𝑙𝑡𝑦 𝑓𝑢𝑛𝑐 (21)
After evaluating all individual populations and calculating their fitness value, the selection stage can be carried out. The method used is the "roulette wheel" method due to the wheel's proportion being assigned to each of the possible selections based on their fitness value. Therefore, the best individual will be kept to be used in the next generation. The next stage is recombination, where several individuals are recombined to make new offspring. This process consists of selecting two individuals at random and recombining them with probability PR.

In the proposed Priority List -Genetic Algorithm method, the one-point cross-over is used between two large groups of parents from individual candidate solutions (Fig. 5). In this case, the cross-over site is randomly selected on the string's length, and then the parts of the two strings outside this point of intersection are exchanged. The individual chromosome for this recombination comes from two major groups: the father's chromosome and the mother's chromosome, cross-over to obtain new possible better solutions. While recombination operates on two or more parent chromosomes, the next stage of mutation is carried out to add diversity to the population and ensure that it is possible to explore the entire search space. The cross-over process shown by genes in both individual parents will be exchanged based on the cross-over point Fig. 5.

After recombination, the next stage is mutation. The mutation is a technique that moves randomly, approaches the optimal point, and possibly moves away. The purpose of mutation is to add diversity to the population and ensure that it is possible to explore the entire search space. The mutation performed is a variation of the well-known "bit flip" mutation where each bit in the string is converted from 0 to 1 or vice versa with probability Pm. The primary function of these mutations in GA is to avoid the local optima trap.

In the final stage, the offspring population made through selection, recombination, and mutation replaces the original parent population. The steps are repeated in an iterative process until stopping conditions are reached. During the evaluation process, the best individuals are saved from being reintroduced to the next generation. Thus, convergence is accelerated, and the best individuals' information is preserved to get the string with the lowest value of total operating costs. This computational process is coded and simulated in MATLAB R2016a. The results of the convergence check will produce the optimal value of the UC implementation.

V. SIMULATION RESULTS

A microgrid system simulation model has been prepared to test the proposed UC algorithm. The simulated microgrid consists of steam generating units, dispatchable diesel engines, and a solar PV plant. The implementation of the unit commitment is set based on the generating unit's capability and availability, and the dispatch arrangement is based on the lowest operating cost of the generating unit. The cost rate for each generator is modeled by the quadratic function for fuel costs and start-up costs. The energy cost of PV generation is also modeled by the multiplication of kWh of production and the selling price rate. The operating costs of all sources will be calculated using formula [START_REF]Indonesia Electricity Supply Business Plan 2016-2025[END_REF], and the characteristics of each element of the grid are listed in TABLE I. The time-based solar photovoltaic power output forecast is used to calculate the UC setting (see TABLE II). The optimal UC schedule together with dispatch is shown in TABLE IV and Fig. 7, and here it can be seen that the obtained UC schedule meets all constraints and produces a quality solution. 01:00 16,50 16,50 5,39 3,00 0,00 0,00 0,00 0,00 0,00 02:00 16,50 16,50 3,55 3,00 0,00 0,00 0,00 0,00 0,00 03:00 16,50 16,10 3,00 3,00 0,00 0,00 0,00 0,00 0,00 04:00 16,50 15,60 3,00 3,00 0,00 0,00 0,00 0,00 0,00 05:00 16,50 15,30 3,00 3,00 0,00 0,00 0,00 0,00 0,00 06:00 16,50 15,00 3,00 3,00 0,00 0,00 0,00 0,00 0,00 07:00 16,50 16,50 3,74 0,00 0,00 0,00 0,00 0,00 0,00 08:00 16,50 16,50 3,53 0,00 0,00 0,00 0,00 0,00 0,00 09:00 16,50 16,50 5,12 3,00 0,00 0,00 0,00 0,00 0,00 10:00 16,50 16,50 5,66 3,00 3,00 0,00 0,00 0,00 0,00 11:00 16,50 16,50 6,81 3,00 3,00 0,00 0,00 0,00 0,00 12:00 16,50 16,50 6,41 3,00 3,00 0,00 0,00 0,00 0,00 13:00 16,50 16,50 6,53 3,00 3,00 0,00 0,00 0,00 0,00 14:00 16,50 16,50 7,46 3,00 0,00 0,00 0,00 0,00 0,00 15:00 16,50 16,50 6,78 3,00 3,00 0,00 0,00 0,00 0,00 16:00 16,50 16,50 5,85 3,00 3,00 0,00 0,00 0,00 0,00 17:00 16,50 16,50 6,45 3,00 3,00 0,00 0,00 0,00 0,00 18:00 16,50 16,50 7,61 3,00 3,00 0,00 0,00 0,00 0,00 19:00 16,50 16,50 8.00 7,95 3,00 3,00 0,00 0,00 0,00 20:00 16,50 16,50 8,00 6,84 3,00 3,00 0,00 0,00 0,00 21:00 16,50 16,50 8,00 5,40 3,00 3,00 0,00 0,00 0,00 22:00 16,50 16,50 8,00 4,23 3,00 0,00 0,00 0,00 0,00 23:00 16,50 16,50 8,00 3,09 3,00 0,00 0,00 0,00 0,00 00:00 16,50 16,50 7,97 3,00 0,00 0,00 0,00 0,00 0,00

In this case, the spinning reserve is specified as 10% from the load, mutation probability 0.8, and maximum iteration 500 times. The iteration process will achieve convergence to get the least total operating costs of the microgrid system with the RES generator shown in Fig. 6. It has been shown here that the process of obtaining low operating costs is carried out properly using genetic algorithm techniques. As a comparison, the PL and GA algorithms' efficacy and efficiency are then compared with the PL algorithms. The results of calculations using the two algorithms are shown in TABLE V, which shows that using the PL and GA algorithms gives the least operating costs, even though the time is a little longer. These results can prove that to optimize unit commitment implementation on the island power system with the RES generator, PL and GA algorithms may be used for better and more accurate results in a relatively short time. In this paper, we propose a unit commitment solution using an optimization algorithm of the improved priority list and genetic algorithm, considering the contribution of renewable energy generation in the microgrid power system. In the microgrid system, to obtain appropriate and better optimization results for implementing the unit commitment, the priority list algorithm and genetic algorithm are applied. In order to simulate the effectiveness of the PL + GA algorithm, computation time and operating costs are compared with other optimization methods for deterministic cases. Here, the PL + GA algorithm shows better performance compared to other optimization methods.

For this reason, the use of a combination of this optimization algorithm can be suitable for application in microgrid systems with renewable energy generators. As the regulator, the government can also determine and recommend this method as a solution for supplying electricity in remote areas and islands.

In this study, it can also be seen that the penetration of renewable energy sources can help alleviate the cost of supplying electricity and stringent reserve requirements.

Fig. 1

 1 Fig. 1 Unit commitment in electric power system

Fig. 2

 2 Fig.2Timor electricity power system[START_REF]Indonesia Electricity Supply Business Plan 2016-2025[END_REF]

Fig. 4

 4 Fig. 4 Flow chart of PL and GA for UC in microgrid system

Fig. 5

 5 Fig. 5 Population chromosome in GA

Fig. 6

 6 Fig. 6 GA convergence for the interconnected configurationThe results obtained for implementing the PL + GA algorithm in managing the unit operating schedule to meet the electrical load requirements are shown in Fig.7.

Fig. 7

 7 Fig. 7 Generation schedule of Timor Microgrid system using PL+GA optimization

  𝑘,𝑛 𝑘 . 𝑃 𝑡,𝑛,𝑘 + 𝛽 𝑘,𝑛 . 𝑉 𝑡,𝑛,𝑘 )

		𝑡	𝑛	(1)
			+ ∑ ∑(𝑆 𝑡,𝑛 * 𝐶𝑆 𝑛 )
			𝑡	𝑛
	In the compilation of UC, several terms or notations are
	commonly used. Some of the notations used in this paper are
	as follows:
	𝑇𝑂𝐶	Total operation cost of power generation
	𝑃 𝑡,𝑛,𝑘	Power output
	𝑉 𝑡,𝑛,𝑘	on/off status of unit n, at time t, and at segment k
	𝛼 𝑘,𝑛	cost coefficients of segment k of unit n
	𝛽 𝑘,𝑛	cost coefficients of segment k of unit n
	𝑆 𝑡,𝑛	start-up status of unit n at time t
	𝐶𝑆 𝑛	start-up cost of unit n
	𝐷 𝑡	Total demand
	𝑈 𝑛	maximum output of unit n
	n	power generating units,
	t, T	time, total period of time
	k	segment of power generation,
	Uit

TABLE I

 I

		Power generations data of Timor power system [1]
		U1	U2	U3	U4	U5
	Pmax(MW)	16	16	5	6	6
	Pmin(MW)	10	10	0	3	3
	a($/Jam)	1000	1000	-	100	100
	b( $/kWh)	16,19	16,19	-	1500	1500
	c($/kW^2h)	0,00048	0,00048	-	0.022	0.022
	MUT(jam)	8	8	1	1	1
	MDT(jam)	8	8	1	1	1
	HSC($)	4500	4500	-	50	50
	CSC($)	9000	9000	-	100	100
	CS(jam)	0	0	0	0	0
	IS(jam)	-1	-1	-1	-1	-1
		U6	U7	U8	U9	U10
	Pmax(MW)	6	6	3	3	3
	Pmin(MW)	3	3	1	1	1
	a(Rp/Jam)	125	125	125	125	50
	b(Rp/kWh)	1500	1500	1500	1500	40
	c(Rp/kW^2h)	0.025	0.025	0.025	0.025	0.001
	MUT(jam)	1	1	1	1	1
	MDT(jam)	1	1	1	1	1
	HSC(1000 Rp)	50	50	50	50	20
	CSC(1000 Rp)	100	80	100	100	50
	CS(jam)	0

TABLE III

 III

			Generator scheduling from PL
	Hour/ Unit	P1	P2 P3 P4 P5 P6 P7 P8 P9
	01:00

TABLE IV Best

 IV dispatch power production of generator units

	Hour/ Unit	P1	P2	P3	P4	P5	P6	P7	P8	P9

TABLE V

 V

				VI. CONCLUSIONS
		Comparison of optimization techniques
	Method	Total fuel cost ($)	Total operational cost ($)	Elapsed time (s)
	PL	3.9885e+04	4.2049e+04	1.6566
	PL+GA	3.9447e+04	4.2028e+04	1.9581

ACKNOWLEDGMENT

We thank Mr. Muhammad Yasirroni, research assistant from University of Gadjah Mada, Yogyakarta, Indonesia, for assistance with technical coding and MATLAB programming of the simulation model in this paper.

His research interests include energy services industries, corporate finance, investment management, capital markets, statistics for business, and complex adaptive system.


Xiaoyu Mao
email: maoxiaoyu@smail.xtu.edu.cn

Jie Liu
email: liujie@smail.xtu.edu.cn

Jianxin Zhong
email: jxzhong@xtu.edu.cn

Rudolf A Römer

Disorder effects in the two-dimensional Lieb lattice and its extensions

We study the localization properties of the two-dimensional Lieb lattice and its extensions in the presence of disorder using the transfer matrix method and finite-size scaling. We find that all states in the Lieb lattice and its extensions are localized for W ≥ 1. Clear differences in the localization properties between disordered flat bands and disordered dispersive bands are identified. Our results complement previous experimental studies of clean photonic Lieb lattices and provide information about their stability with respect to disorder.

I. INTRODUCTION

Quantum descriptions of materials in the solid state can be carried out equally well in real space as in momentum space [START_REF] Ashcroft | Solid State Physics[END_REF][START_REF] Hook | Solid State Physics[END_REF]. In perfect crystals, the momentum space description offers the advantage that all spatial information is already encapsulated in the unit cell, hence allowing for a dramatic reduction in the complexity of characterising such materials. The structure of the resulting energy bands, their number, curvature, linearity, filling and behaviour at the Fermi energy are all important descriptors, essential for the determination of the material properties of said systems. In recent years, e.g., interest has broadened this approach towards a description of the topological properties of the bands [START_REF] Kane | [END_REF]4] with concepts such as Fermi arcs [5,6] and Weyl points [7] indicating highly-specialized band structure features.

In a similar context, so called flat bands have also received much attention [8]. The term refers to energy bands that are dispersionless in the whole of k-space [9][10][11][12]. This is equivalent to effectively having zero kinetic energy in these bands and hence to promoting the remaining terms in the Hamiltonian such as potential and interaction terms. Usually, such flat bands emerge only for specially "engineered" lattice structures such as quasi-1D lattices [13][14][15], diamond-type lattices [16], and Lieb lattices [17][18][19][20][21][22]. Lieb lattices were originally used to enhance magnetic effects in model studies [9,23,24] but have now been shown to be of experimental relevance in Wigner crystals [12], high-temperature superconductivity [10,18], photonic wave guide arrays [8,[25][26][27][28], Bose-Einstein condensates [29,30], ultra-cold atoms in optical lattices [31] and electronic systems [32].

The eigenstates in a flat band lattice are known to form so-called compact localized states [13], which due to the effectively zero group velocity, cannot participate in transport although they are in principle infinitely degen-erate [16]. Since the existence of flat bands relies on exact conditions, weak deviations from these conditions, such as potential disorder, could in principle fundamentally alter the nature of these localized states. Indeed, previous studies reveal a wide spread of tendencies depending on the nature of the perturbations [13,16,[33][34][35][36] as well as the environment of the bands [37][38][39][40]. Our focus here is slightly different: instead of concentrating on the localization properties of the flat band states, we will investigate how the localization properties in the neighboring dispersive bands are changed by the presence of the flat bands. In particular, we are interested in whether the presence of the flat bands changes the localization properties and the finite-size scaling behaviour in disordered two-dimensional (2D) Lieb model and its extensions.

II. MODELS AND METHOD

A. The Lieb lattices L(n)

The original Lieb lattice [23] is based on a square lattice with additional sites added along the horizontal and vertical directions as shown in Fig. 1. We indicate the original sites as A while B and C denote the added horizontal and vertical sites, respectively. In Ref. 41, extended Lieb lattices have been proposed with additional added sites as also shown in Fig. 1. We shall denote all these lattices symbolically as L(n) with n marking the number of sites added in both horizontal and vertical directions, i.e. here L(1) to L(4). Our numerical approach to study the localization properties of the Lieb lattices L(n) is based on (i) direct diagonalization of the Hamiltonian for small system sizes, and (ii) the transfer matrix method (TMM) [42,43]. TMM has been used extensively for various disorder problems [44,45] and is the currently accepted method-of-choice for high-precision estimates of parameter-dependent localization lengths [46].

We start with the standard Anderson tight-binding Hamiltonian where | r denotes the tight-binding state site r = (x, y), x and y indicate the exact position in the 2D plane and r = r indicates the nearest-neighbor connections as defined in Fig. 1. The onsite potential energy ε r at site r is randomly and uniformly distributed in [-W/2, W/2] such that W indicates the strength of the disorder. The t r, r are the hopping integrals for a particle going from site r to site r , and represent the kinetic energy part of the Hamiltonian. We only set nearest neighbours hopping in our model, so that t r, r = t ≡ 1 for nearest sites sets the energy scale and t r, r = 0 otherwise. We have chosen hard-wall boundary conditions, but also checked that results are similar for periodic boundaries.

H = r ε r | r r| - r = r t r, r | r r | , (1)
B. Transfer-matrix method for L(1)

For the L(1) lattice, in order to compute the localization length λ of the wave function from the TMM corresponding to the Schrödinger equation HΨ = EΨ, we consider a quasi-1D strip consisting of M chains of A and B sites with additional M sites, labelled C, in the vertical (y-) direction as shown in Fig. 1. We choose to propagate the TMM in the horizontal (x-) direction until convergence [? ]. This implies a TMM setup with two transfer matrices, T A→B and T B→A . From sites A to B in x-direction, we find

  ψ B x+1 ψ A x   = T A→B   ψ A x ψ B x-1   =    εx,y-E t - t εx,y-1-E - t εx,y+1-E 1 M - ty εx,y-1-E - t † y εx,y+1-E -1 M 1 M 0 M      ψ A x ψ B x-1   (2)
with M × M matrix

t y = t         0 1 0 • • • 0 0 0 0 1 • • • 0 0 . . . . . . . . . 0 0 0 • • • 1 0 0 0 0 • • • 0 1 (1) 0 0 • • • 0 0         (3)
1 M ≡ diag εx,1-E t , εx,2-E t , . . . , ε x,M -E t
while y -1 and y + 1 indicate the unrenormalized sites neighboring y for the hopping terms with t y and t † y in Eq. ( 2). From sites B to A, a more standard TMM form emerges as

  Ψ A x+1 Ψ B x   = T B→A   Ψ B x Ψ A x-1   =    εx,y-E t 1 M -1 M 1 M 0 M      Ψ B x Ψ A x-1   . (4)
The transfer matrices T B→A and T A→B transfer the wave vector amplitudes Ψ between atoms A and B. Along the transfer (x-) axis, we multiply these matrices, reorthogonalize the M Ψ states, i.e. the columns of the TMM, at least after every 10th multiplication, compute Lyapunov exponents γ i , i = 1, . . . , M and study their accumulated changes in variance until we reach a desired accuracy (error of typically 0.1%) for the smallest γ min [43,[48][49][START_REF] Beenakker | Random-matrix theory of quantum transport[END_REF]. The localization length is then found as λ(M, W, E) = 1/γ min for given energy E and disorder strength W .

C. TMM for L(2) and beyond

For the extended Lieb lattice L(2), we can derive the corresponding TMM recursion formula in a similar fash-ion. Starting again from the first sites in x-direction as given in Fig. 1 for L(2), we find for the transfer from sites

A → C,   Ψ C x+1 Ψ A x   = T A→C   Ψ A x Ψ B x-1   =    E1 M - tty (εx,y-2-E)(εx,y-1-E)-t 2 - tt † y (εx,y+2-E)(εx,y+1-E)-t 2 -1 M 1 M 0 M    ×   Ψ A x Ψ B x-1   , (5)
where

E = ε x,y -E t - t (ε x,y-2 -E) (ε x,y-2 -E) (ε x,y-1 -E) -t 2 - t (ε x,y+2 -E) (ε x,y+2 -E) (ε x,y+1 -E) -t 2 . ( 6
)
From sites C → B, we have

  Ψ B x+1 Ψ C x   = T C→B   Ψ C x Ψ A x-1   =   ( εx,y-E t )1 M -1 M 1 M 0 M     Ψ C x Ψ A x-1   (7a) and for B → A   Ψ A x+1 Ψ B x   = T B→A   Ψ B x Ψ C x-1   =   ( εx,y-E t )1 M -1 M 1 M 0 M     Ψ B x Ψ C x-1   . ( 7b
)
As before, the hopping in vertical (y-) direction has been included in renormalization of the onsite energies of the A sites as well as renormalized hopping between A sites, albeit more complicated in structure than for L(1). Still, following this procedure, it is not difficult -albeit time consuming -to construct the transfer matrices for the remaining L(3) and L(4) lattices and we will not include those lengthy terms here.

D. Finite-size scaling in 2D

According to renormalization group ideas [START_REF] Wilson | The renormalization group and the {lunate} expansion[END_REF], we can study the properties of large critical systems from smaller subsystems [START_REF] Wegner | [END_REF][53][54]. We define the dimensionless, reduced localization length Λ M (E, W ) = λ(E, W )/M . In order to get the infinite-size localization length ξ [55] from the Λ M (E, W ), we should apply the one-parameterscaling hypothesis, i.e.

Λ M (E, W ) = f (ξ(E, W )/M ). ( 8
)
As a numerical procedure, this means that the finitesize Λ M values can be scaled by a factor ξ, such that all values collapse onto a single scaling function f . We fit the parameter ξ by using a least squares procedure [42] either as a function of E for constant W , or, more standard, as a function of W for constant E. We then determine the absolute scale of ξ(W, E) by fitting the largest disorder with the form of a truly localised behaviour, i.e.

Λ M = ξ/M + b(ξ/M ) 2 , (9)
as expected for large W t. We emphasize that numerically, this procedure does not use the Taylor expansion methods nor does it assume any power-law diverging terms as is normally standard in studies of the 3D Anderson transition and its variants [46,56]. Rather, as there is no evident metal-insulator transition in the 2D Anderson model [43], there is no universal functional form to fit towards [42] and we hence have to revert back to the simpler methods of maximizing data overlap as done in most pre-1999 studies of Anderson localization [43].

III. RESULTS

A. Dispersion relations and disorder-broadened density of states for L(n)

In the clean lattice L(1), we can easily compute the dispersion relations from the tight-binding model (1) as

E 1 = 0, E 2,3 = ± 4 + 2 (cos k x + cos k y ) ( 10
)
where we have chosen the unit length to be 1, and the k x , k y are the wave vectors mapping onto the x and y axes. Fig. 2 plots the L(1) band structure and we can clearly see the flat band at energy E = 0, and a Dirac cone around the 'M ' point (k x , k y ) = (π, π). Similarly, we compute the band structures for L(n), n = 2, 3, 4 and show them in the last three plots in Fig. 1. We see that every L(n) exhibits n flat bands separated by n + 1 dispersive bands [41]. Furthermore, the odd n = 1, 3 lattices show Dirac cone at special points in the (k x , k y ) plane while the even n = 2, 4 lattices only seem to have parabolic-type dispersions. We now turn on the disorder, i.e. W > 0. Obviously, k x and k y are no longer good quantum numbers and the dispersion relations lose their meaning. Nevertheless, as discussed previously, the flat bands correspond to zero kinetic energy are dominated by the remaining term in the Hamiltonian, i.e. the random disorder term. Hence already small W values can quickly lift the degeneracy of the flat bands and the states can begin the overlap with states in the originally dispersive, neighboring bands. In order to study the interplay of flat bands and dispersive bands in the disordered situation, we have computed the disorder-dependence of the density of states (DOS) by exact diagonalization for small system sizes M 2 = 13 2 , 10 2 , 9 2 and 8 2 for L(1), L(2), L(3) and L(4), respectively. The DOS has been obtained after averaging over 300 samples. The results are given in Fig. 2. Clearly, the figure shows that in all cases the degeneracy of the flat bands is lifted very quickly and the originally degenerate states rapidly begin to move into energy regions occupied by the dispersive bands. Similarly, but somewhat less pronounced, the flat regions of the dispersive bands vanish upon increasing W . Overall, when W becomes of the order of 2-3, comparable to the band width of the largest subband in each case, the DOS loses it peak features.

B. Localization and finite-size scaling

In Fig. 3, we give the result of the TMM calculations for Λ M (E, W ) at E = 0. The transverse system widths were M = 10, 12, . . . , 20 with error 0.1% for L(1), and M = 10, 12, . . . , 22 with error 0.2% for L(2), L(3) and L(4). This results in about 10 8 transfer matrix multiplications for each pair (E, W ). We have chosen E = 0 since this corresponds to the flat bands for L(1) and L(3) while it is in the centre of the subbands for L(2) and L(4). Hence E = 0 allows us to study directly the influence of the disorder on originally flat and on originally dispersive bands. We find that in all cases, irrespective of the flatness of the original clean lattice bands, the values of Λ M (0, W ) reduce for all values of W when increasing M as shown in the inset of Fig. 3. This is the behaviour as expected for localized states where for sufficiently large

L(1)
L( 2) W and M , the values of λ M saturate such that Λ M decreases when M increases. Indeed, in the main panels of Fig. 3, we see that all Λ M (0, W ) data can be made to scale onto single scaling curves. These curves also decrease with increasing M , hence showing strong localization. Furthermore, for λ M M such that Λ M 1, we find that Λ M ∝ 1/M as expected for strong localization with λ M ∼ λ ∞ (where λ ∞ denotes the value of the localization length in the thermodynamic limit). We also note immediately that the localization lengths for the states originating from flat bands are nearly one order of magnitude smaller than the corresponding Λ M (0, W ) values for the originally dispersive bands. Hence, in perfect agreement with our aforementioned expectations, we see that the absence of kinetic energy in the flat bands for the clean systems translates into strong localization for the disordered systems.

L(3) L (4)
In Fig. 4 we show the resulting W dependence of the scaling lengths ξ(0, W ). According to Thouless [42,57], we expect ξ(W ) ∝ W -2 for a strictly 1D weakly disordered chain while in the 2D Anderson model the low W behaviour has been shown to follow a non-universal ξ(W ) = aW -α exp(βW -γ ) with a, β positive and α, γ of order unity [43]. These fits are also given in Fig. 4 and the resulting fit parameters are listed in Table I. As is clear from Fig. 4, the behaviour for L(1) and L(3), i.e. the flat bands, is clearly different from the behaviour of the dispersive bands for L(2) and L(4). This is already visible in the absolute magnitude of the ξ values where the L(1) and L(3) cases are about three order of magnitude below the values for L(2) and L(4) for low W . For L(1) and L(3), the simple power-law fit ξ(W ) ∝ W -2 seems to work reasonably well at 1 ≤ W ≤ 2 while at W > 2, in both cases the estimated ξ values drop below the fit. One might want to speculate that this apparent consistency with a 1D fit result might also indicate that the spatial localization itself is along 1D structures, but this is clearly beyond what TMM can reliably infer. We also note that the agreement with the α = 2 value in 1D localization is not perfect while the value of a should be 105 at E = 0 [58] for 1D and about 12M for quasi-1D stripes of width M [59]. Clearly, none of these values are obvious from the fit parameters given in Table I. On the other hand, for L(2) and L(4), none of our fits as given in Table I seem very convincing with the non-universal 2D behavior perhaps performing slightly better. Of course the lack of a simple fit to the low-W behaviour of ξ is not surprising and simply indicates that the disordered Lieb . The values for a, α, β and γ are chosen from Tab. I.

IV. CONCLUSIONS

We have studied the localization properties of the disordered 2D Lieb lattice and its extensions. After appropriately finite-size scaling, we find that all states are localized, down to disorder of W ∼ 1. For even lower dis-orders, we do not have sufficient computational resources for reach the same accuracies as presented here, but, as for the standard 2D Anderson model [59], we also do not find any indications that the qualitative behaviour would change. Our results indicate that the finite-size scaled localization lengths ξ for the energies corresponding to flat bands in the clean system (W = 0) show a behaviour reminiscent of perturbative result for 1D. This could potentially provide spatial information about the nature of these localized states. For energies corresponding to dispersive bands in the clean system there is a tendency towards much larger localization lengths as is expected from the 2D Anderson model of localization [55]. Hence, the clear differences observed between flat and dispersive bands in photonic Lieb lattices [41] are predicted to be largely robust against the presence of disorder.

where the ω = -1+ √ 3i 2 , ρ ± = 3 q ± q 2 -5 3 3 and q = cos k x + cos k y . For L(3), we have The remaining 5 dispersive bands emerge as solutions to a 5th order equation, i.e.

E 1 = 0, E 2,3 = ± √ 2, (A2a)
E 5 -7E 3 + 9E -2 (cos k x + cos k y ) = 0, (A3b) which we solve numerically.

LFIG. 1 .

 1 FIG.1. Schematic structure of lattices L(1), L(2), L(3) and L(4), respectively. Circles denote the lattice sites with blue/dark circles corresponding to the sites of a simple square lattice and the red/light circles denoting the added sites of the Lieb lattices. The thin black lines between the sites indicate hopping connections. The letters A, B, . . . , E give special labels used in the explanation of the TMM in the text. The x and y axes are given by arrows.

  denoting the hopping connections to neighboring A sites in y direction and (1) indicating a possible periodic boundary term in y-direction. Note that in this setup, the C sites have been effectively renormalized away into effective onsite energies and changed vertical hopping terms such that the TM has the standard 2M × 2M size with 1 M and 0 M denoting M × M identity and zero matrices. The M × M wave function amplitudes in the xth line are marked as Ψ A,B x , either A or B, with y = 1, . . . , M , labelling the position of the renormalized cubic sites [47]. In this notation the term εx,y-E t

FIG. 2 .

 2 FIG. 2. The band structure (top row) and the smoothed and normalized density of states (DOS, bottom row) for L(1) with width M = 13, L(2) with M = 10, L(3) with M = 9 and L(4) with M = 8. Colours in the top row vary from purple at low energies to red at high energies. In the bottom row, the colour indicates the numerical values of the DOS and the thin black lines are equal-DOS contours.

FIG. 3 .

 3 FIG. 3. Double logarithmic plot of scaled reduced localization length, ΛM , as a function of ξ/M for L(1), L(2), L(3) and L(4) as indicated. The chosen values for W are 1.0 (♦), 1.01, . . . , 1.05 ( ), . . . , 2.0, 2.1, 2.2, 2.4, 2.6 . . . , 10.0 ( ) and vary in color. Some data points are not shown with symbols for clarity. Error bars are within symbol size. Insets: ΛM as a function of 1/M for each L(n), respectively. Symbols and colors are as in the main panels, dashed lines are guides to the eye only.

FIG. 4 .

 4 FIG.4. The scaling parameter ξ(0, W ) (•) and fitted functional forms (lines) for L(1), L(2), L(3) and L(4), respectively. The blue solid line corresponds to aW -α , the purple dashed line is the fit function aW -2 e βW -1, while the green dotted line is the fit function aW -α e βW -γ. The values for a, α, β and γ are chosen from Tab. I.

E 4 , 5 , 6 , 7 = 1

 45671 ± 3 ± 5 + 2 (cos k x + cos k y ). (A2b)For L(4), the flat band solutions are given by E

TABLE I

 I

	. Fit functions for ξ(W ) and lattices L(1), L(2), L(3)
	and L(4) with a, α, β, γ are fit coefficients estimated from a
	Levenberg-Marquardt non-linear fitting procedure.
	models L(n) have indeed localization behaviour different
	from the standard 2D Anderson model.

ACKNOWLEDGMENTS

We wish to acknowledge the National Natural Science Foundation of China (Grant No. 11874316), the Program for Changjiang Scholars and Innovative Research Team in University (Grant No. IRT13093), and the Furong Scholar Program of Hunan Provincial Government (R.A.R.) for financial support. This work also received funding by the CY Initiative of Excellence (grant "Investissements d'Avenir" ANR-16-IDEX-0008) and developed during R.A.R.'s stay at the CY Advanced Studies, whose support is gratefully acknowledged. We thank Warwick's Scientific Computing Research Technology Platform for computing time and support. UK research data statement: Data accompanying this publication are available from the corresponding authors.

Appendix A: Expressions for the dispersion relations

For completeness, we give the dispersion relations for the Lieb lattices L(2), L(3) and L(4) here. For L(2), we find


Jie Liu
email: liujie@smail.xtu.edu.cn

Xiaoyu Mao
email: maoxiaoyu@smail.xtu.edu.cn

Jianxin Zhong
email: jxzhong@xtu.edu.cn

Rudolf A Römer

Localization, phases and transitions in the three-dimensional extended Lieb lattices

We study the localization properties and the Anderson transition in the 3D Lieb lattice L3(1) and its extensions L3(n) in the presence of disorder. We compute the positions of the flat bands, the disorder-broadened density of states and the energy-disorder phase diagrams for up to 4 different such Lieb lattices. Via finite-size scaling, we obtain the critical properties such as critical disorders and energies as well as the universal localization lengths exponent ν. We find that the critical disorder Wc decreases from ∼ 16.5 for the cubic lattice, to ∼ 8.6 for L3(1), ∼ 5.9 for L3(2) and ∼ 4.8 for L3(3). Nevertheless, the value of the critical exponent ν for all Lieb lattices studied here and across disorder and energy transitions agrees within error bars with the generally accepted universal value ν = 1.590(1.579, 1.602).

I. INTRODUCTION

Flat energy bands have recently received renewed attention due to much experimental progress in the last decade [1]. The hallmark of such flat bands is an absence of dispersion in the whole of k-space [2][3][4][5], implying an effectively zero kinetic energy. This leads to a whole host of effects in transport and optical response such as, e.g. localization of eigenstates without disorder [6] and enhanced optical absorption and radiation. Further studies explorations of flat-band physics have now been done in Wigner crystals [5], high-temperature superconductors [3,7], photonic wave guide arrays [1,[8][9][10][11], Bose-Einstein condensates [12,13], ultra-cold atoms in optical lattices [14] and electronic systems [15].

Systems that exhibit flat-band physics correspond usually to specially "engineered" lattice structures such as quasi-1D lattices [6,16,17], diamond-type lattices [18], and so-called Lieb lattices, [7,[19][20][21][22][23]. Indeed, the Lieb lattice, a two-dimensional (2D) extension of a simple cubic lattice, was the first where the flat band structure was recognized and used to enhance magnetic effects in model studies [2,24,25]. Most other flat-band systems cited above are also of the Lieb type and exists as either 2D, quasi-1D or 1D lattices [26]. Less attention has been given to 3D flat-band systems [18] or extended Lieb lattices [23,27]. Furthermore, until recently hardly any work has investigated the influence of disorder on flatband systems [16,28]. Recently, instead of concentrating on the properties of flat-band states, we investigated how the localization properties in the neighboring dispersive bands are changed by the disorder for 2D flat-band systems [27].

In the present work, we extend these studies to the class of 3D extended Lieb lattices. As is well known [29] the Anderson transition in a simple cubic lattice with uniform potential disorder x ∈ [-W/2, W/2] at each site x is characterized by a critical disorder W c = 16.0(5)t [30], with t denoting the nearest neighbor hopping strength. The full energy-disorder phase diagram is characterized by a simple-connected region of extended states ranging from ±6t at W = 0 and ending at W c = 16.530 (16.524, 16.536) for E = 0 [31]. The critical exponent of the transition has been determined with ever greater precision as close to, e.g., ν = 1.590(1.579, 1.602) [31] and 1.57(2) [32]. The 3D Lieb model, shown in Fig. 1 together with its extensions, is characterized by additional sites on the edges between the original site of the cubic lattice. As such, the transport along the edges should become more 1D-like and we expect that the phase diagram should have a smaller region of extended states.

II. MODELS AND METHOD

A. Transfer-matrix method for the 3D Lieb lattices and its extensions L3(n)

We denote the Lieb lattices as L d (n) if there are n equally-spaced atoms between two original nearest neighbors in a d-dimensional lattice. Here, we shall concentrate on L 3 (1), L 3 (2) and L 3 (3) as shown in Fig. 1. To explore the effects of disorder, we use the standard Anderson Hamiltonian

H = x x |x x| - x =y t xy |x y|. (1)
The orthonormal Wannier states |x describes electrons located at sites x = (x, y, z) of Lieb lattice with hard boundary condition (we have similar results for periodic boundary conditions as well). The hopping integrals t xy = t only for x, y being nearest neighbors as indicated by the lines in Fig. 1, otherwise t xy = 0. For L 3 (1), in order to calculate the localization length λ of the wave function by the transfer-matrix method (TMM), we consider a quasi-one-dimensional bar, with cross area M 2 and length L >> M . A unit length corresponds to original site-to-site distances as indicated by A sites in Fig. 1. Along the transfer axis in the z-direction, there are two different slices in L 3 (1). The first slice contains the original A sites, and the added B and C sites to form an A-B-C slice, the second (D-)slice only contains the added D sites as shown in Fig. 1. The TMM equation implementing HΨ = EΨ at energy E for the Hamiltonian (1) can be written as two parts. First, transferring from slice A-B-C to slice D slice, we have

  Ψ D z+1 Ψ A z   = T A→D   Ψ A z Ψ D z-1   =   E1 M 2 - 1 z,x-1,y -E t x-- 1 z,x+1,y -E t x+ - 1 z,x,y-1-E t y-- 1 z,x,y+1-E t y+ -1 M 2 1 M 2 0 M 2     Ψ A z Ψ D z-1   , (2)
where

E = z,x,y -E t - t z,x-1,y -E - t z,x+1,y -E - t z,x,y-1 -E - t z,x,y+1 -E , (3)
and 0 M 2 , 1 M 2 denote M 2 × M 2 zero and identity matrices, respectively. Similarly, t x+ , t x-, t y+ and t y-are M 2 × M 2 connectivity matrices in the positive/negative x/y directions. With this choice of TMM set-up, we effectively renormalize the added B, C (red) sites shown in Fig. 1(a). Taking M = 3 as an example, we can explicitely write the 3 2 × 3 2 matrices

t x-=t             
0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 (1) 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 (1) 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 (1) 0 0

             (4)
and

t x+ = t † x-. Similarly, t y-=t             
0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 (1) 0 0 0 0 0 0 0 0 0 (1) 0 0 0 0 0 0 0 0 0 (1) 0 0 0 0 0 0

             (5)
and t y+ = t † y-. In Eqs. ( 4) and ( 5), the matrix entries (1) can be chosen 0 for hard-wall boundaries and 1 for periodic boundaries. In this way, the effects of sites B and C have been renormalized into effective onsite energies E and hopping terms t x± , t y± keeping the transfer matrix T A→D in the standard 2M 2 × 2M 2 form. We emphasize that Ψ A,D z denotes a vector of length M 2 for wave function amplitudes in the zth slice [START_REF]and D are not needed to label Ψ A,D z , but we retain them here for the readers convenience[END_REF], either A or D, with x, y = 1, . . . , M , labelling the position of the original cubic sites in this slice. In this notation the term z,x,y -E t

1 M 2 ≡ diag z,1,1-E t , z,
  Ψ A z+1 Ψ D z   = T D→A   Ψ D z Ψ A z-1   =   z,x,y -E t 1 M 2 -1 M 2 1 M 2 0 M 2     Ψ D z Ψ A z-1   . (6)
in similar notation.

The TMM method proceeds by multiplying successively T A→D by T D→A along the bar in z-direction, using M 2 possible starting vector Ψ A z (1) = (1, 0, . . . , 0), Ψ A z (2) = (0, 1, . . . , 0), Ψ A z (M 2 ) = (0, 0, . . . , 1) to form a complete set. We regularly renorthogonalize these M 2 Ψ states, usually after every 10th multiplication. The Lyapunov exponents γ i , i = 1, 2, . . . , M 2 , and their accumulated changes are calculated until a preset precision is reached for the smallest γ min [29,[START_REF] Oseledets | [END_REF][35][START_REF] Beenakker | Random-matrix theory of quantum transport[END_REF]. The localization length λ(M, E, W ) = γ min > 0, the dimensionless reduced localization length is Λ M (E, W ) = λ(M, E, W )/M . These considerations set out the TMM for L 3 (1). For the extended Lieb lattices, we follow a similar strategy, leading to an even more involved renormalization scheme which we refrain to review in the interest of brevity.

B. Finite-size scaling

The metal-insulator transition (MIT) in the Anderson model of localization is expected to be a second-order phase transition, characterized by a divergence in a cor-

relation length ξ(W ) ∝ |W -W c | -ν at fixed energy E, and ξ(E) ∝ |E -E c | -ν at fixed disorder W [37]
, where E c is the critical energy and ν, W c as before.

We determine the reduced correlation length ξ/M in the thermodynamic limit assuming the single parameter scaling i.e. Λ M (M, E, W ) = f (ξ/M ) [30]. For a system with an MIT this scaling function consists of two branches corresponding to localized and extended phases. Using finite-size scaling (FSS) [38], we can obtain estimates of the critical exponent. Here, we use a method [32,[START_REF] Eilmes | [END_REF] that models two kinds of corrections to scaling: (i) the presence of irrelevant scaling variables and (ii) non-linearity of the scaling variables. Hence one writes Λ = F (χ r M 1/ν , χ i M y )., where χ r the relevant scaling variable and χ i the irrelevant scaling variable. We next Taylor-expand Λ and F up to order n i and n r such that

Λ= ni n=0 χ n i M ny F n (χ r M 1/υ ), F n = nr k=0 a nk χ k r M k/ν .(7)
Furthermore, we also expand χ i and χ r by ω

= (W c - W )/W c (or (E c -E)/E c ) to consider the importance of the nonlinearities, χ r (ω) = mr m=1 b m ω m , χ i (ω) = mi m=0 c m ω m . (8)
In order to fix the absolute scales of Λ in (7) we set b 1 = c 0 = 1. We then perform the FSS procedure for various values of n i , n r , m i , m r , in order to obtain the best stable and robust fit by minimizing the χ 2 statistic. We quote goodness of fit p values to allow the reader to judge the quality of our results.

III. RESULTS

A. Dispersion and disorder-broadened density of states for L3(n)

For a clean L 3 (1) system, the dispersion relation can be derived from (1) as

E 1,2 = 0, E 3,4 = ± 6 + 2 (cos k x + cos k y + cos k z ), (9)
where the k x , k y , k z are the reciprocal vectors corresponding to the x, y and z axes, respectively. Fig. 2(a) shows the energy structure of L 3 (1), where we can see two dispersive bands which meet linearly at the R point (k x , k y , k z ) = (π, π, π) at E = 0. This coincides in energy with the doubly-degenerate flat band. Analogously, we calculate the energy structures for L 3 (n), n = 2, 3, 4 and plots them in Figs. 2(b),(c) and (d), respectively. We can see that each L 3 (n) lattice has n doubly degenerate flat bands separating n + 1 dispersive bands. Furthermore, the two dispersive bands at high and low energies are separated by energy gaps for these models. We also note that for L 3 (3) two dispersive bands again meet linearly, as for L 3 (1), but in this instance at the Γ point (k x , k y , k z ) = (0, 0, 0) at E = 0. No such linear behaviour can be found for L 3 (n) with n even.

We now include the disorder, i.e. W > 0, and we calculate the disorder-dependent density of states (DOS) by direct diagonalization for small system sizes M 3 = 5 3 , 5 3 , 4 3 and 4 3 for L 3 (n), n = 1, 2, 3, 4, respectively. The DOS is generated from W = 0 to W = 5.2 in step of 0.05 with 300 samples for L 3 (n), n = 1, 2, 3, while we have 100 samples for L 3 (4). We also apply a Gaussian broadening of the energy levels to obtain a smoother DOS. The results are shown in Fig. 2. For weak disorders we can clearly identify the large peaks in the DOS with the flat bands for all L 3 (n) models. From W ∼ 3 onward, the various peaks have merged into one broad DOS. Also, the energy gaps for L 3 (n), n = 1, 2, 3, 4, vanish quickly with increasing W .

B. Phase Diagrams

Fig. 3 shows the energy-disorder phase diagram for L 3 (1). The phase diagram was determined from the scaling behaviour of the Λ(E, W ) for small system sizes M = 6, M = 8 and M = 10 with TMM error ≤ 0.1% [START_REF] Eilmes | [END_REF]. Data for W < 1 fluctuates too much to give useful results and hence has been omitted from the figure. Clearly, the phase diagram is qualitatively similar to the phase diagram of the standard 3D Anderson model, although the band width and the critical disorder at E = 0 are different. In particular, the critical disorder is reduced by about 50% compared to the Anderson model. This is in agreement with the discussion in section I. Close to the band edges for small W ≤ 4 we also see a small reentrant region as is also found in the 3D Anderson model. However, the shoulders that develop at E ∼ ±2.75 and W = 6 are a novel feature. The DOS at such strong disorder does not exhibit clearly any similar signatures.

For L 3 (2) and L 3 (3), we show the phase diagrams in Fig. 4, determined with TMM errors of ≤ 0.2% and with the same system sizes as for L 3 (1). As before, small disorder results have to be excluded. Our numerical results support, as for L 3 (1), a mirror symmetry at E = 0 and the results as shown in Fig. 4 have been explicitly symmetrized. For both L 3 (2) and L 3 (3), the phase boundaries of the central dispersive band support a reentrant behaviour, although this is less so for L 3 (3).

The obvious difference between the phase diagrams of L 3 (1), L 3 (2) and L 3 (3) is that the extended region for L 3 (1) lattice is simply connected, while for L 3 (2) and L 3 (3) it is disjoint. This difference can be attributed to the presence of the energy gaps in L 3 (2) and L 3 (3) as in Fig. 2. Let us denote, as in the cubic Anderson model, a critical disorder W c as the disorder value at the transition point from extended to localized behaviour at energy E = 0. Then we see that the critical disorders are W c ∼ 16.530 for the cubic lattice [31], ∼ 8.6 for L 3 (1), ∼ 5.9 for L 3 (2) and ∼ 4.8 for L 3 (3). Hence as expected, in the Lieb lattices the last extended states vanish already at much weaker disorders and the trend becomes stronger with increasing n in each successive L 3 (n).

C. High-precision determination of critical

properties for the Lieb models

1. Model L3(1)
In order to determine the critical properties at the phase boundaries for the Lieb models, we have to go to larger system size for a reliable FSS. In all cases, the results are collected up to M = 20 and with TMM convergence errors ≤ 0.1%. Using the phase diagram as in Fig. 3 as a rough guide, we pick out 4 points of special interest, namely, two transitions as a function of W at the band centre at constant E = 0 and outside the band centre at E = 1. Furthermore, we also study two transitions as function of E corresponding to the point marking the reentrant behaviour at constant W = 3 and the kink in the phase boundary at constant W = 6. In Fig. 5, we show the Λ M (E, W ) data, the resulting scaling curves and the variation of the scaling parameter ξ for typical examples of FSS results.

In Table I we present fits for all 4 cases shown in Fig. 5 with higher expansion coefficients n r and m r that show that our results are stable with respect to an increase in an expansion parameter. We have also checked that they are stable with respect to slight changes in the choice of parameter intervals δW and δE for fixed energy and fixed disorder transitions, respectively. The reader will have noticed, however, that the accuracy of the data is not good enough to reliably fit irrelevant scaling contributions and hence the results in Table I are all for n i = m i = 0 although we have indeed performed our FSS allowing for these additional parameters. Furthermore, one can see in Fig. 5 that the accuracy of the TMM data becomes worse for the fixed disorder transitions at W = 3 and especially W = 6. The reason for this behaviour is in principle well understood since at the points, the DOS has an appreciable variation which leads to extra corrections not well captured in the FSS [39]. Usually, larger system sizes M can reduce these variations but this is not possible here due to computational limitations.

Models L3(2) and L3(3)

We follow a similar strategy as in the previous section in order to finite-size scale the localization lengths for L 3 (2) and L 3 (3). The TMM convergence errors were chosen as ≤ 0.1% up to M = 16 and, due to the increased complexity of these models, as ≤ 0.2% for the largest system size with M = 18. Fig. 6(a) shows Λ M (E = 0, W ) and the scaling curve for L 3 (2) at energy E = 0 with n r = 3, m r = 3. From the panel with the Λ M (E = 0, W ) data, it is very hard to observe a clear crossing at W c . The situation improves for Λ M (E, W = 4) in Fig. 6(b) which exhibits a clear crossing of Λ M around E c ∼ 1.70. For L 3 (3) shown in Fig. 6(c) the crossing for Λ M (E = 0, W ) is again somewhat less clear. Nevertheless, in all three cases, the FSS results produce stable and robust fits with estimates for W c , E c and ν as shown in Table I. As for L 3 (1), the FSS fits L 3 (2) and L 3 (3) do not resolve potential irrelevant scaling corrections.

IV. CONCLUSIONS

There are two ways to understand the Lieb lattices as originating from the normal simple cubic lattices: (i) as shown in Fig. 1, one can view the L 3 (n) lattices as a cubic lattice with additionally added sites between the vertices of the cube, effectively allowing for additional backscattering and interference along the original site-to-site connections and hence potentially leading to more localization. On the other hand, one might argue that (ii) the L 3 (n) lattices can be constructed by deleting sites from a cubic lattice, for example a central site in Fig. 1(a) and the 6 face-centered sites. In this view, the decrease of possible transport channels should give rise to stronger effective localization. Both constructions lead to the same predictions and agree with what we find here, namely, the localization properties in all L 3 (n) lattices show an increased localization with respect to the cubic Anderson lattice and become stronger when n increases. This is, for example, clear from looking at the behaviour of W c (n) in Table I. It is instructive to study the behaviour as n → ∞. From Fig. 2, we see that the overall band width decreases as n increases. At the same time, the number of flat bands increases and the extremal energy of these bands extends as well towards |E| = 2. Hence I.

for very large n, where L 3 (n) is simply a 2 3 renormalized lattice, but n renormalized sites apart, with proliferating flat bands. Our results for the critical exponent then suggest that as n increases and the dispersive bands become smaller, the critical properties in each band still retain the universality of the 3D Anderson transition -at least up to n = 3 that we have been able to compute (cp. Fig. 7. This is in good agreement with previous results in loosely coupled planes of Anderson models in which the universal 3D behaviour was also retained [40]. However, for loosely coupled planes, the MIT was retained even for small interplane coupling -a truly 2D localization behavior only emerged when the interplace coupling was zero. The point of view of this work is different, i.e. the change from 3D dispersive bands with an MIT to a solely 1D system without MIT is not a continuous change, but rather an eventual replacement and shrinking of dispersive bands by a proliferation of flat bands as n grows.

L3 E), expansion orders nr, mr are listed as well as resulting critical disorders Wc (or energies Ec), their 95% confidence intervals (CI), the critical exponent ν, its CI, and the goodness of fit probability p. The averages contain the mean of the three preceding Wc (or Ec) and ν values, with standard error of the mean in parentheses. The bold Wc, Ec and ν values highlight the fits used as examples in Figs. 5 and6.

Appendix A: Dispersions

For completeness, we here include the dispersion relations shown in Fig. 2. For L 3 (2), we have

E 1,2 = 1, E 3,4 = -1, E 5 = ρ + + ρ -, (A1a)
E 6 = ωρ + + ω 2 ρ -, E 7 = ωρ -+ ω 2 ρ + , (A1b)
where ω = -1+i

√ 3 2 , ρ ± = 3 -q(k) 2 ± q(k) 2 2 -7 3
3
and q(k) = 2 (cos k x + cos k y + cos k z ). For L 3 (3), we find

E 1,2 = √ 2,

FIG. 1 .

 1 FIG. 1. (a) The Lieb lattice L3(1) and its extensions (b) L3(2) and (c) L3(3). The blue spheres denote the original nearest-neighbor sites in the underlying cubic lattice while the red spheres show the added sites. The solid lines indicate the cubic structure. The coordinate system is to help identify the TMM setup used in our study as are the labels A, B, C, and D.

FIG. 2 .FIG. 3 .√ 3 and Emax = 2 √ 3 .

 2323 FIG. 2. (a)-(d) Dispersion relations for clean systems and (e)-(h) dependence of the normalized DOS on W for L3(1) to L3(4). In all cases, the flat bands are doubly degenerate. Different colors in the dispersion relations denote different bands while the colors in the DOS indicate different DOS values as also emphasized by the contour lines.

0 FIG. 4 .

 04 FIG. 4. Phase diagrams for (a) L3(2) and (b) L3(3) lattices. The symbols, lines and colors are as in Fig. 3, i.e. representing small M estimates with M = 6, 8 and 10. The solid squares ( ) denote high-precision FSS results from ΛM with an TMM error ≤ 0.1% for width M ≤ 16 and ≤ 0.2% for width M = 18. The diamonds ( ) denote the maximal band edges from W = 0 at ±3 for L3(2) and ±2 √ 2 for L3(3).

FIG. 5 .

 5 FIG. 5. (a) FSS of the localization lengths for L3(1) with E = 0, (b) E = 1, (c) W = 3, and (d) W = 6. System sizes M are 14 ( ), 16 ( ), 18 (♦), 20 ( ). The left half in each panel denotes a plot of ΛM disorder W or energy E, the solid lines are fits to the data acquired by Eqs. (7)-(8) with (a+b) nr = 3, mr = 1, (c) nr = 2, mr = 1 and (d) nr = 1, mr = 1. The right half in each panel shows the scaling function F (solid line) and the scaled data points with the same nr and mr as in the corresponding left half while each inset gives the scaling parameter ξ as a function of disorder strength W , in (a) and (b), or energy E in (c) and (d). The parameters of the fits are shown in detail in TableI.

FIG. 6 .FIG. 7 .

 67 FIG. 6. FSS of the localization lengths for (a) L3(1) at E = 0 and (b) W = 4 as well as for (c) L3(3) at E = 0. System sizes M are 10 ( ), 12 ( ), 14 ( ), 16 ( ), 18 (♦). The arrangement in each panel is as in Fig. 5, i.e. scaling curves (solid lines) and scaled ΛM data (symbols) in the left half of each panel, scaling curve F (lines) with scaled data (symbols) in the right half and ξ in the inset. The chosen expansion coefficients are (a) nr = 2, mr = 2, (b) nr = 2, mr = 1 and (c) nr = 2, mr = 1 as highlighted in TableI.

TABLE I .

 I Critical parameters of the MIT for L3(n), n = 1, 2, 3 and 4. The columns denoting system width M , fixed E (or W ), range of W (or

						1)
	∆M	E	δW	nr	mr	Wc	CI(Wc)	ν	CI(ν)	p
	16-20	0	8.25-8.9	3	1	8.594	[8.585, 8.604]	1.57	[1.49, 1.65]	0.15
	16-20	0	8.25-8.9	2	2	8.598	[8.586, 8.610]	1.55	[1.46, 1.63]	0.08
	16-20	0	8.25-8.9	3	2	8.595	[8.582, 8.607]	1.57	[1.48, 1.66]	0.13
	Averages:					8.596(4)		1.56(3)
	∆M	E	δW	nr	mr	Wc	CI(Wc)	ν	CI(ν)	p
	14-20	1	8.0-8.8	3	1	8.435	[8.429, 8.441]	1.60	[1.54, 1.65]	0.18
	14-20	1	8.0-8.8	2	2	8.439	[8.432, 8.447]	1.57	[1.53, 1.62]	0.19
	14-20	1	8.0-8.8	2	3	8.438	[8.431, 8.446]	1.57	[1.53, 1.62]	0.21
	Averages:					8.437(3)		1.58(2)
	∆M	W	δE	nr	mr	Ec	CI(Ec)	ν	CI(ν)	p
	16-20	3	3.725-3.785	2	1	3.748	[3.747, 3.749]	1.75	[1.68, 1.82]	0.88
	16-20	3	3.725-3.785	2	2	3.748	[3.747, 3.749]	1.76	[1.67, 1.84]	0.86
	16-20	3	3.725-3.785	3	1	3.748	[3.747, 3.749]	1.75	[1.68, 1.82]	0.86
	Averages:					3.748(1)		1.75(3)
	∆M	W	δE	nr	mr	Ec	CI(Ec)	ν	CI(ν)	p
	16-20	6	3.04-3.11	1	1	3.077	[3.070, 3.083]	1.54	[1.08, 2.01]	0.14
	16-20	6	3.04-3.11	2	1	3.076	[3.069, 3.082]	1.54	[1.09, 1.99]	0.24
	16-20	6	3.04-3.11	2	2	3.077	[3.069, 3.084]	1.54	[1.07, 2.00]	0.21
	Averages:					3.077(3)		1.54(14)
						L3(2)
	∆M	E	δW	nr	mr	Wc	CI(Wc)	ν	CI(ν)	p
	12,14,18	0	5.85-6.05	2	2	5.964	[5.958, 5.969]	1.75	[1.57, 1.92]	0.08
	12,14,18	0	5.85-6.05	2	3	5.965	[5.959, 5.970]	1.70	[1.51, 1.89]	0.08
	12,14,18	0	5.85-6.05	3	2	5.963	[5.956, 5.971]	1.75	[1.57, 1.92]	0.07
	Averages:					5.964(3)		1.73(6)
	∆M	W	δE	nr	mr	Ec	CI(Wc)	ν	CI(ν)	p
	10,12,14	4	1.6-1.8	2	1	1.704	[1.701, 1.708]	1.55	[1.43, 1.68]	0.18
	10,12,14	4	1.6-1.8	1	3	1.705	[1.701, 1.709]	1.56	[1.43, 1.70]	0.1
	10,12,14	4	1.6-1.8	2	2	1.703	[1.700, 1.707]	1.53	[1.40, 1.66]	0.2
	Averages:					1.704(2)		1.55(5)
						L3(3)
	∆M	E	δW	nr	mr	Wc	CI(Wc)	ν	CI(ν)	p
	12-18	0	4.7-4.875	2	1	4.79	[4.786, 4.794]	1.63	[1.48, 1.78]	0.49
	12-18	0	4.7-4.875	1	2	4.791	[4.786, 4.795]	1.63	[1.48, 1.78]	0.47
	12-18	0	4.7-4.875	2	2	4.791	[4.786, 4.795]	1.63	[1.48, 1.78]	0.47
	Averages:					4.790(2)		1.63(5)

ACKNOWLEDGMENTS

We wish to acknowledge the National Natural Science Foundation of China (Grant No. 11874316), the Program for Changjiang Scholars and Innovative Research Team in University (Grant No. IRT13093), and the Furong Scholar Program of Hunan Provincial Government (R.A.R.) for financial support. This work also received funding by the CY Initiative of Excellence (grant "Investissements d'Avenir" ANR-16-IDEX-0008) and developed during R.A.R.'s stay at the CY Advanced Studies, whose support is gratefully acknowledged. We thank Warwick's Scientific Computing Research Technology Platform for computing time and support. UK research data statement: Data accompanying this publication are available from the corresponding authors.

Last, for L 3 (4), the four doubly degenerate flat bands are given as

and the remaining five dispersive bands are the solutions of the 5th order equation