Gia-Diem Pham

Khanh-Duy Cong Do

Duy-Hung Chung

Sung-Goon Park

Frédéric Plourde

Kwang-Yong Kim

Cong-Truong Dinh

HEAT TRANSFER CHARACTERISTICS OF AN INTERNAL COOLING CHANNEL WITH PIN-FINS USING RECESSED-EXTRUDED ENDWALL COUPLING IN GAS TURBINE BLADE

This work numerically explores the vortices creation and preservation in a channel with pin-fins for two endwallturbulation systems: the flat endwall and the recessed-extruded endwall, while studying the flow structure therein. The heat transfer characteristics, including Nusselt number, friction factor, and HTEI, are analyzed and compared with the smooth endwall case at Reynolds numbers ranging from 7400 to 36000. The simulations are carried out with the RANS method. It is reported in the results that, with the new endwall configuration, the high heat transfer areas near the pin-fins are noticeably enlarged compared to the flat endwall at all the tested Reynolds numbers. Besides, the new endwall design increases the heat transfer level near the pin-fins, represented by an enhancement of 64.3% of the Nusselt number compared to the flat endwall at Re = 36000. However, the growth in friction factor due to the flow obstruction at all the Reynolds numbers is trivial. It has been discovered that the enhancement in heat transfer capacity is induced mainly by the combination of the indentations and protrusions that the geometry creates. When considering the heat transfer performance, the HTEI of the recessed-extruded endwall is increased by approximately 62.4% at Re = 21500 compared to the flat endwall, with a height of 2.5% channel height. These findings demonstrate the tremendous potential for enhancing the pin-fin's heat transfer by altering the endwall design.

μ Fluid dynamic viscosity, Kgm -1 s -1 ρ Fluid density, kgm -3 Abbreviations RANS Reynolds-Averaged Navier-Stokes HV Horseshoe Vortex/Horseshoe Vortices HTEI Heat Transfer Efficiency Index 1. INTRODUCTION As the need for a gas turbine with a high thermal efficiency grows, gas turbines can achieve improved thermal efficiency by raising the temperature of the gas entering the turbine component. In this context, the scientists are facing enormous material-related difficulties since the operating conditions for turbine blades are particularly harsh when the temperature exceeds 2000K. This temperature is far higher than the material's melting point in modern gas turbines. Therefore, to maintain the temperature within an acceptable range for strength and service life, elements of the turbine that come into touch with the hot gas need to be cooled or protected. Turbines are frequently cooled internally, which involves passing colder air from the compressor through ducts into the vanes and blades of the turbine. The heat transfer capacity of the cooling ducts can be improved by using ribs, pin-fins, and dimples as turbulators.

This study focuses on pin-fins, a common turbulator located in the trailing edge domain of turbine vanes and blades. This type of turbulator is used to enhance heat transfer in this domain and, in the meantime, promoting structural support for the turbine blade. Conventionally, it is understood that the presence of features like pin-fins in a cooling route increases the wetted flow area and intensifies turbulence, improving the convective heat transfer efficiency. According to APLab's studies [START_REF] Dinh | Heat Transfer Performance of Internal Cooling Turbine Blades Using Cutted-Root Rib Design[END_REF][START_REF] Dinh | Numerical investigation of truncated-root rib on heat transfer performance of internal cooling turbine blades[END_REF][START_REF] Hoang | Large eddy simulation of the turbulence heat and mass transfer of pulsating flow in a V-sharp corrugated channel[END_REF] on the heat transfer characteristics of channels with modified turbulators, there are many benefits to use turbulators like ribs and pin-fins in order to increase heat transfer of the turbine blade's cooling channel.

Many researchers have studied the flow characteristics and heat transfer properties in cooling ducts with pin-fins using experimental and computational methods. Most of those studies investigated the effects of geometric variables on heat transfer, e.g., pin-fin arrangement, cross-section shape, height, and spacings in the streamwise and spanwise directions. The earliest studies about pin-fins appeared in the late eighties with the investigation of Metzger et al. [START_REF] Metzger | Developing Heat Transfer in Rectangular Ducts With Staggered Arrays of Short Pin Fins[END_REF]. They came to the conclusion that the heat transfer coefficient in a channel with pin-fins is much greater than that on the endwalls or in the smooth channels. Siw et al. [START_REF] Siw | The Effects of Different Pin-Fin Arrays on Heat Transfer and Pressure Loss in a Narrow Channel[END_REF] studied pin-fin arrangements. They achieved that staggered pin-fin arrays are more likely to cause increased heat transfer capacity and a more potent friction factor than inline configuration. Ames and Dvorak [START_REF] Ames | Turbulent Transport in Pin Fin Arrays: Experimental Data and Predictions[END_REF] discovered that the first few rows show a progressive rise in heat transfer coefficient for both inline and staggered pin-fin arrays. As the flow enters the arrays, turbulence gradually intensifies because of the array's congestion, which speeds up the flow, increases flow separation, and enhances HV.

Numerous studies of pin-fin designs on the cross-sectional shape, including circle, cube [START_REF] Sahiti | Performance comparison of pin fin in-duct flow arrays with various pin cross-sections[END_REF], ellipse [START_REF] Li | Multiobjective Optimization of Staggered Elliptical Pin-Fin Arrays[END_REF], and fan-shape [START_REF] Moon | Analysis and optimization of fan-shaped pin-fin in a rectangular cooling channel[END_REF][START_REF] Moon | Heat Transfer Performance of a New Fan-Shaped Pin-Fin in Internal Cooling Channel[END_REF], attempted to determine the ideal cross-sectional shape. The pin-fin's turbulence-related mechanism alters significantly depending on the pin-fin's cross-section [START_REF] Ghosh | Shape Optimization of Pin Fin Arrays Using Gaussian Process Surrogate Models Under Design Constraints[END_REF]. It is found that the pin-fins with sharp edges, e.g., cubic and diamond forms, highly raise the pressure loss of the channel and also the heat transfer capability [START_REF] Sahiti | Performance comparison of pin fin in-duct flow arrays with various pin cross-sections[END_REF]. Moreover, the pin-fins with aerodynamic shapes seem attractive for increasing the total heat transfer efficiency index (HTEI) by lowering the friction factor [START_REF] Li | Multiobjective Optimization of Staggered Elliptical Pin-Fin Arrays[END_REF]. Over these cross-sections, elliptical pin-fins show the highest HTEI by lengthening the pin-fin's streamwise diameter, which results in engendering the friction factor and concurrently expanding the contact region of the pin-fin with the fluid.

On the pin-fin height-to-diameter ratio (H/D), short pin-fins generally induce better heat transfer capacity than long pin-fins. Park et al. [START_REF] Park | Heat Transfer on Rotating Channel With Various Heights of Pin-Fin[END_REF] investigated the heat transfer of pin-fins by varying H/D in a range of 2-4 experimentally, and found that the level of the heat transfer increases with the increase of H/D. On the other hand, Brigham and VanFossen [START_REF] Brigham | Length to Diameter Ratio and Row Number Effects in Short Pin Fin Heat Transfer[END_REF] studied short pinfins and found that the heat transfer capacity only varies with the channel's Reynolds number for H/D < 2. The Nusselt number is influenced by both H/D and Re for longer pin-fins (H/D > 2). On the pin-fin spacings, it was concluded that the heat transfer influences of each pin-fin on the side walls it is attached to, and its neighbors change with the changes in these spacings. The modification of the pin-fin arrangement (staggered and inline), as well as the pin-fin's spanwise and streamwise spacing-todiameter ratio (Sx/D and Sy/D), were studied by Jeng and Tzeng [START_REF] Jeng | Pressure drop and heat transfer of square pin-fin arrays in inline and staggered arrangements[END_REF]. The overall findings demonstrated that the staggered disposition often results in increased heat transfer capacity and friction factor penalty yielded by the flow obstructing features.

Of experimental studies, most results are obtained by timeaveraged techniques. Only Ames et al. [START_REF] Ames | Turbulent Transport in Pin Fin Arrays: Experimental Data and Predictions[END_REF][START_REF] Ames | Turbulent Augmentation of Internal Convection Over Pins in Staggered-Pin Fin Arrays[END_REF], Ostanek and Thole [START_REF] Ostanek | Wake development in staggered short cylinder arrays within a channel[END_REF], and Anderson and Lynch [START_REF] Anderson | Time-resolved stereo PIV measurements of the horseshoe vortex system at multiple locations in a low-aspect-ratio pin-fin array[END_REF] were able to resolve the flow's unsteadiness through experiments. Park et al. [START_REF] Park | Heat Transfer on Rotating Channel With Various Heights of Pin-Fin[END_REF] used the naphthalene sublimation technique on time-averaged outcomes and observed the conceptualization of secondary flow and stagnant vortices induced by pin-fins. Chyu and Goldstein [START_REF] Chyu | Influence of an array of wall-mounted cylinders on the mass transfer from a flat surface[END_REF] pointed out the Horseshoe Vortices (HV) that contribute to the high heat transfer of the endwalls. On time-resolved studies using particle image velocimetry, Ostanek and Thole [START_REF] Ostanek | Wake development in staggered short cylinder arrays within a channel[END_REF] showed that the shear layer eddies and Karman vortices dominated the momentum transport in the near wake downstream of pin fins. Anderson and Lynch [START_REF] Anderson | Time-resolved stereo PIV measurements of the horseshoe vortex system at multiple locations in a low-aspect-ratio pin-fin array[END_REF] investigated the HV at Reynolds numbers between 10000 and 50000 using stereoscopic particle image velocimetry. They discovered that the HV system oscillates between two quasi-periodic modes.

Of the computational methods, most of the works are based on the steady Reynolds-averaged Navier-Stokes (RANS) analysis due to its cost compared to its usefulness. In a straight channel, Li et al. [START_REF] Li | Assessment of Six Turbulence Models for Modeling and Predicting Narrow Passage Flows, Part 2: Pin Fin Arrays[END_REF] evaluated six turbulence models: k-ω, shear-stress transport, explicit algebraic Reynolds stress, realizable k-ε, renormalization group k-ε, and v 2 -f. Their research revealed that all the models tested accurately predicted the qualitative characteristics of flow and heat transfer, but their quantitative predictions showed 10%-30% discrepancies from experimental data. Jin et al. [START_REF] Jin | Effect of shape and distribution of pin-fins on the flow and heat transfer characteristics in the rectangular cooling channel[END_REF] investigated the effect of different variants of turbulence models, which include standard k-ε, realizable k-ε, renormalization group k-ε, and standard k-ω models, on heat transfer of a channel with short pin-fin arrays. The group of k-ε models has enhanced wall treatment enabled. It is discovered that the error between numerical results and experimental data of Nu when using the standard k-ε model is within ±20%. Besides RANS analysis, large-eddy simulation (LES) or hybrid RANS-LES has also been utilized in pin-fin heat transfer prediction. To forecast flow and heat transfer in a pin-fin array, Lee and Shih [START_REF] Lee | Effects of heat loads on flow and heat transfer in the entrance region of a cooling duct with a staggered array of pin fins[END_REF] investigated the heat transfer of pin-fin arrays using a scale-adaptive model (SAS), a form of hybrid method that utilizes the von Karman length scale in the turbulent structure. Their conclusions include a maximum of 5% uncertainty in the results of Nu compared to experimental data.

On the interaction between flow and a pin-fin, the HV pattern that forms around the cylindrical pin-fin fixed on a plate was investigated by Choi et al. [START_REF] Choi | Endwall heat transfer and fluid flow around an inclined short cylinder[END_REF] and Sircar et al. [START_REF] Sircar | Turbulent flow and heat flux analysis from validated large eddy simulations of flow past a heated cylinder in the near wake region[END_REF]. They observed the form of a HV around the pin-fin and the vortex shedding caused by shear layer instability. Moreover, the wake in the rear of the pin-fin is the leading region of turbulence production through shear in the mean velocity field. Without a doubt, the channel's heat transfer properties depend on the robustness of these vortices and their interaction. Uzol and Camci [START_REF] Uzol | Elliptical Pin Fins as an Alternative to Circular Pin Fins for Gas Turbine Blade Cooling Applications: Part 2 -Wake Flow Field Measurements and Visualization Using Particle Image Velocimetry[END_REF] used particle image velocimetry to monitor and measure the wake flow for elliptic and circular pin-fin arrays. The findings demonstrate that the enormous low-momentum wake zone, created by a separation of the flow via the pin-fin, causes the loss downstream of the circular pin.

A new approach has been opened by Do et al. [START_REF] Do | Numerical Investigation of Heat Transfer Characteristics of Pin-Fins with Roughed Endwalls in Gas Turbine Blade Internal Cooling Channels[END_REF]. They implied that almost all studies of pin-fins heretofore investigated the effects of pin-fins geometry and their spacing on heat transfer characteristics, where the impacts of the endwall are left behind. They suggested that modification of the endwall can be used to spread and preserve the vortices around the pin-fins. In their study, the changing endwall designs modified the flow behavior and vortex system around the pin-fins. They also showed that the recessed endwall could enlarge HV while deducting the wake size. Their study suggested that changing the endwall shape increased the Nusselt number by 78.0% and the HTEI by 77.7%. Hence, the investigation of the effect of endwall configuration on the channel's heat transfer performance still looks compelling and promising to all researchers.

By simulating and performing experiments, Bai et al. [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF][START_REF] Bai | Numerical investigation on heat transfer and pressure drop of pin-fin array under the influence of rib turbulators induced vortices[END_REF] induced pin-fins in a ribs-roughed channel and investigated the heat transfer of the combination of pin-fins and ribs. They discovered that different shapes of ribs induced different secondary flows leading to distinct pin-fin domain's inlet flow phenomena. Chang et al. [START_REF] Chang | Endwall heat transfer and pressure drop in rectangular channels with attached and detached circular pin-fin array[END_REF] researched the heat transfer of channel with attached and detached pin-fins by experiments. They figured out that the channel with detached pin-fins shows a considerably lower friction factor when compared with attached pin-fins but a moderately lower Nusselt number.

In this work, the heat transfer properties of a pin-fin channel with modified endwalls are investigated numerically. Here, a new design of a recessed-extruded endwall is suggested. The heat transfer characteristics of this model, i.e., Nusselt number, friction factor, and HTEI, are differentiated from those of the flat and recessed endwalls.

NUMERICAL METHODS

Description of the problem 2.1.1 Configuration of the calculation model

An experimental model developed by Bai et al. [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF] served as the basic configuration of the channel for building the present computational domain. This paper proposes the recessed-extruded endwall design of the heated endwalls to modify the vortex system surrounding the pin-fins. The geometric illustration of the flat endwall, recessed endwall from previous work [START_REF] Do | Numerical Investigation of Heat Transfer Characteristics of Pin-Fins with Roughed Endwalls in Gas Turbine Blade Internal Cooling Channels[END_REF], and the recessedextruded endwall designs are shown in Fig. 2. For the recessed endwall, the indentations are placed at the 1 st , 3 rd , and 5 th rows and represented by the red frame in Fig. 2b. In Fig. 2c, the protrusions (blue frame) are positioned between the indentations (red frame), which are situated in the pin-fin rows. The protrusions and indentations are arranged in both the upper and Figure 1. Schematic diagram of Bai et el. [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF] lower heated walls. The protrusions and indentations here have an equal height of 0.5mm, i.e., 2.5% of the height of the channel. Besides, the values of inclined angles for the two objects are equal.

Boundary conditions

In this work, the incompressible dry air with constant characteristics that enters the cooling channel has a uniform temperature of Tin = 318.15K and a 5% turbulence intensity. The range of Reynolds numbers used to compute fluid velocity at the intake is from 7400 to 36000. The gauge pressure at the channel's outlet is 0Pa. A temperature of 298.15K and the no-slip wall condition are given to the pin-fins, the endwalls of the test section and the plain heated section. The temperature conditions at the inlet and heated wall are set differently to investigate the ability to transfer heat from the fluid to the walls. All the other walls are adiabatic and non-slip. Thus, even though the momentum boundary layer starts at the duct's inlet (x/D = -12.5), the thermal boundary layer does not start until the inlet of the plain heated section (x/D = 0). The reference pressure of the computational domain is 1 atm, and the reference temperature is 287K. These hypotheses are acceptable since the channel is relatively short and the heat load is rather modest, indicating a minimal pressure loss. Boundary conditions are established according to the experimental and simulation conditions set out by Bai et al. [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF]. The results of the present simulations will be compared with their experimental and simulation results.

Numerical analysis

The commercial software Ansys CFX ® 19.1 [29] was used in the present study to numerically solve the governing differential equations, which comprise continuity, RANS, and energy equations, in order to study the flow characteristics and heat transfer capacity in the channel with pin-fins. The advection terms are discretized using a high-resolution scheme and the finite volume method. The k-ω turbulence model was used throughout the study.

Figure 3. The mesh structure of the channel

The convergence requirement is satisfied if the roots-meansquare values of the continuity, momentum, energy, and turbulent kinetic energy equations are less than 10 -6 . Additional monitoring metrics are developed to confirm convergence, such as the mass flow at the outlet, the changes in gas temperature at the inlet and outlet surfaces, and the heat flux at the heated wall surfaces. The structured grid division in the ANSYS-ICEM ® program was utilized in this work. The elements near the wall and the pin-fins are specified with a height of 1x10 -6 m to comply with the demands of the turbulence model and have y + values less than 1 over the whole Reynolds number range. The structural grid in Fig. 3 is constructed using hexahedral elements, and an O-grid is placed around the pin-fins to enhance mesh quality and depict their round shape.

Data reduction

The Reynolds number (Re) of the channel is defined as:

Re = ρuD h μ ( 1
)
where u is the inlet velocity of the channel; ρ is the air density; Dh and μ are the channel hydraulic diameter and dynamic viscosity of the coolant, respectively. The local heat transfer coefficient (h) of the heated surfaces is determined by:

h = q T w -T b ( 2
)
where q is the heat flux of the heated wall, Tw is wall temperature, and Tb is the bulk temperature of the fluid, which is interpolated based on the x location along the streamwise direction according to the formula of Bai et al. [START_REF] Bai | Numerical investigation on heat transfer and pressure drop of pin-fin array under the influence of rib turbulators induced vortices[END_REF] The local Nusselt number (Nu) is used to represent the local heat transfer capacity, which can be written as:

Nu = hD h λ ( 3
) where λ is the thermal conductivity of coolant fluid.

According to the correlation formula of Dittus-Boelter [START_REF] Dittus | Heat transfer in automobile radiators of the tubular type[END_REF], the Nusselt number (Nu0) for fully developed turbulent flow in a smooth channel with a rectangular cross-section is:

Nu 0 = 0.023.Re 0.8 Pr 0.3 (4) The friction factor, which represents the pressure loss of the channel, is defined as:

f = ∆PD h 2ρ∆Lu 2
(5)

The friction factor obtained for fully developed flow in a smooth channel [START_REF] Dittus | Heat transfer in automobile radiators of the tubular type[END_REF] is expressed as: f 0 = 0.316Re -0.25 (6) The HTEI (η) is defined as:

η = ( Nu Nu 0 ) / ( f f o ) 1/3 (7)

RESULTS AND DISCUSSION

Verification and Validation

The dimensionless wall distance y + on all the heated wall surfaces is evaluated, with Re = 36000. The conditions of the kω turbulence model are satisfied by the y + value area-averaged on these walls, which is 0.07. The maximum y + value on all the walls is 0.39, which is completely insured to be smaller than 1. The fact that y + is proportional to the velocity at the entry, which is proportional to the Reynolds number, also indicates that the y + value is fulfilled for all other Reynolds numbers.

A grid sensitivity study was performed to determine the optimum grid system. In this process, the tested mesh size varies from 1.9 to 5.6 million. Figure 4 shows how changing the mesh size affects the averaged Nusselt number of all heated surfaces. The results show that the averaged Nusselt number barely varied by 0.09% when the mesh size changes from 3.3 million to 4 million elements. Hence, for all of the numerical simulations performed in this work, the numerical grid size of the 3.3 million elements was chosen to achieve a compromise between processing resources and analytical precision. To seize the mesh size dependency in detail, pressure and heat flux results obtained from three different vital meshes are compared. Selected meshes are indicated in Fig. 4. Figure 5 illustrates how mesh refinement affects the convergence of the channel's characteristics. The pressure and wall heat flux at y = D/2 along the heated channel are analyzed with three different meshes. Because the wall temperature is lower than the inlet flow temperature, the actual heat flux is negative. The heat flux is thus assumed to be the absolute value in this instance for the sake of consistency and simplicity. Mesh 1 is chosen with 1.9 million meshes. The baseline mesh chosen for the project is Mesh 2, which has a 3.3 million element count. Mesh 3 has a 5.6 million element count. Figure 5 shows that the presented lines of the baseline mesh (Mesh 2) and the finest mesh (Mesh 3) nearly coincide for both the pressure and the heat flux. It was further discovered that the highest relative error between these two meshes is less than 0.5%. Furthermore, the outcomes are within an asymptotic range of convergence, as seen by the magnified photos and the extracted data. As a consequence, the dependability of the baseline mesh is reasonably validated.

For validation, the present numerical results are compared with previous numerical and experimental results [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF][START_REF] Chang | Endwall heat transfer and pressure drop in rectangular channels with attached and detached circular pin-fin array[END_REF] for the averaged Nusselt number on heated walls. The configurations of the pin-fin array of the reference literatures are the same as that used in this study. Furthermore, all the endwalls in the reference cases are flat. Note that the no ribs configuration in the study of Bai et al. [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF] is the same as the flat endwall design in this study. Figure 6 compares the numerical results using the k-ω model with those obtained by Bai et al. [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF], as well as with two sets of experimental data from Bai et al. [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF] and Chang et al. [START_REF] Chang | Endwall heat transfer and pressure drop in rectangular channels with attached and detached circular pin-fin array[END_REF]. The main factor that differs the numerical method between the reference [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF] and this work is the mesh type. Which, reference [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF] uses auto-generated tetrahedral and prism elements near the walls. Moreover, to ensure y + around 1, the first layer thickness may be very low, even though the mesh size is limited. Thus, they used a mesh with high volume transition; this corresponds to high aspect ratio elements that may affect the results.

Figure 6. Comparison of numerical and experimental results at

Re=21500 It is found that, in every case, the averaged Nusselt number is proportional to the Reynolds number. Compared to the numerical results from Bai et al. [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF], the findings from the current investigation are more in line with both the experimental data sets. The mean Nu calculated on the heated walls using the k-ω turbulence model is lower than the experimental data by 24.9%, 9.6%, and 1.9% at Re = 7400, 14000, and 21500, respectively. The numerical values are successively higher than the experimental results at Re = 29000 and 36000 by 0.2% and 4.4%, respectively. This is partially due to the experimental uncertainty of 1% of the mean velocity error of the flowmeter, 0.75C of the thermocouple, and 2% of the Nusselt number measurement error [START_REF] Bai | Investigation of ribs disturbed entrance effect of heat transfer and pressure drop in pin-fin array[END_REF]. Thus, the validity of the numerical technique has been rationally confirmed with a maximum acceptable deviation of about 25% and a credible range of Reynolds number deviation of less than 5%. Those errors are acceptable for the prediction of heat transfer when taking into account the error of the experimental method and the simplification of the numerical processes.

Flow structure in a rectangular channel with circular pin-fins

Figure 7. Distribution of Nusselt number on the lower heated wall at Re = 21500

The introduction thoroughly discusses the primary characteristics and properties of the pin-fin array. As mentioned, the heat transfer characteristics induced by the pin-fins can be grouped into two main regimes: the HV and the wake vortex.

Figure 7 shows the distribution of the Nusselt number on the lower heated wall. Figure 8a shows the vortex structure around the pin-fins. It can be observed that the vortex system around the pin-fins includes some complex vortices, e.g., HV and wake vortex. It is remarked that before the flow reaches the first row of the pin fin arrays, the streamlines are entirely uniform. When the flow reaches the first row, a stagnation point forms. Then, the flow accelerates from this point and separates on the rear part of each pin-fin due to the adverse pressure gradient, forming the wake vortices. Here, the poor heat transfer zone is mainly caused by the wake vortices. Although these vortices constantly come into contact with the heated surfaces of the pin-fins and endwalls, they do not interfere with the flow of cooling fluids. They are the primary reason for the poor heat transfer zones located downstream of each pin-fin.

(a) Streamlines and velocity distribution near the endwall (b) Contours of Nu on the lower endwall and profiles of vorticity X near pin-fin

Figure 8. Vortex structures induced by pin-fins at Re = 21500

Besides, Fig. 8b shows several vorticity contours in the xdirection at different sections (parallel with the yz plane) around the middle pin-fin in the first row, which also describes the 3D shape of the horseshoe vortex and its operational regions near the endwall. This vortex appears in a horseshoe shape, hence called a HV. This vortex shape fits the shape studied by other investigations.

At the junctions between pin-fins and endwalls, HV forms upstream of each pin-fin and fades downstream. Under the effects of wake vortices, HV, and high-speed regions caused by the narrow spaces between pin-fins, the flow gradually becomes more turbulent, and the pressure drop increases when moving to the rear rows. Turbulence enhances the heat transfer coefficient significantly. Observing Fig. 7 more carefully, the HV and the wakes dominate the heat transfer in the first two rows due to lack of turbulence. However, the HV gradually lose its indistinctness in the following rows and gradually becomes more apparent. Besides, by the cause of the low heat transfer regions, the wakes from the first two rows are in direct contact with the pin-fins of the third row and so on. This affects the flow characteristics of the rear rows, which makes the flow more turbulent and also enhances the Nusselt number downstream. Thus, turbulence is the leading cause of the enhancement of heat transfer. In this and below subsection, the results of the reference recessed endwall from previous literature [START_REF] Do | Numerical Investigation of Heat Transfer Characteristics of Pin-Fins with Roughed Endwalls in Gas Turbine Blade Internal Cooling Channels[END_REF] is added to compare the flow and heat transfer properties of the three designs. Be aware that the height of the recessed endwall in the study is 1 mm, which is double the height of the recessed endwall here. When considering Fig. 9, which shows the velocity contours and streamlines for the three cases, it is achieved that the wake vortices behind the pin-fins of the flat endwall case are distinct. On the other hand, the recessed endwall may partially inhibit the growth of these vortices, leading to a discernible cutoff in both their breadth and intensity. However, the local vortices in this case are still presented, evidenced by the low heat transfer zones in Fig. 14b. In the case of the recessed-extruded endwall, the wake vortices of the rows with five pin-fins (the 1 st , 3 rd , and 5 th rows) are totally excluded, while the wake vortices in the other two rows are pushed away from the endwall. As a consequence, low heat transfer zones in Fig. 14c are reduced as a result of this phenomenon because the impact of local vortices on the endwall is lessened. Local vortices of the downhill endwall, which play a role in the development of HV downstream, are the persistent vortices in this situation. Moreover, the highest velocity zones of all three cases are found in the 2 nd and 4 th rows of the recessed endwall case due to the more narrowed channel. According to the velocity distribution at the mid-channel cross-section (Fig. 10), the maximum velocity is slightly greater in the flat case than in the recessed and recessed-extruded instances. This is explained by the fact that, in accordance with the continuity equation, the recessed-extruded and recessed exhibit concave geometry lower than the initial flat endwall, increasing the cross-sectional area and decreasing the flow velocity. The low-velocity vortices downstream the pins in the flat and recessed cases are more equivalent, although being more prominent in the recessed-extruded case. These are the wake vortices that are pushed away from the endwalls. Additionally, towards the endwalls, the effects of the roughed endwall cases are also more pronounced.

Effects of the roughed endwall on the flow field

The velocity distribution near the endwall (at the section that is offset 0.1mm from the lower endwall) at Re = 7400 is shown in Fig. 11. It is observed that the velocity distributions near the pin-fins of the cases with altered endwalls are enhanced compared to the flat endwall case. In the figures, the highvelocity zones of the three cases coincide with the high-Nu zones in Fig. 14 despite the different Re. It is then achieved that the phenomenon of increased heat transfer is the same for all Reynolds numbers.

Furthermore, it appears that the upward profile of the protrusion causes the high-velocity zones in the case of roughed endwalls to drift toward the pin-fins. This is evidenced by the high-velocity zones tend to stick with the pin-fins more in the indented rows in Figures 11b and11c. This is particularly clear in the recessed-extruded case, where the flat endwall design induces an adverse pressure gradient, the high-velocity zone is less protracted at the stagnation point and splits swiftly as it advances downstream. This leads to the disappearance of lowvelocity zones behind pin-fins in Fig. 11c compared to Fig. 11a. Besides, the velocity distribution of the recessed endwall case near the pin-fin is more similar to the flat endwall instance because the 2 nd row of the recessed case is the row where the profile is the most stable and least affected by the front rows. Moreover, the jet-like beam caused by the remained HV of the recessed endwall case is observable in the 2 nd and 4 th rows of Fig. 11b. However, it almost disappeared in the recessed-extruded endwall case. Figure 12 demonstrates the velocity distributions near the endwall at Re = 36000. In this figure, due to the altered concave wall, the vortex textures downstream of each pin-fin differ between the two cases of roughed endwall. Common features of this case include the widening of the high-velocity zone around each pin-fin and the shrinking of the rear wake vortex compared to the flat endwall in all rows. This is associated with the diminishment of the wake vortices discussed above.

In addition, the flow through it behaves similarly to the interaction with a rib turbulator. The continuous upward and downward endwall profile are likewise regarded as a type of rib, or so named in this study as protrusion. When the flow strikes the top of each protrusion, the flow splits apart, rolls downhill, and collides with the pin-fin at the stagnation point, which helps to widen the HV. In addition, the protrusion also produces a wake region in its back, as can be observed in Fig. 9c. The vortices in this region interact directly with the rear HV. The heat transfer characteristics of this protrusion can be observed in Fig. 12c. Due to the local vortices presented, the protrusion exhibits high heat transfer areas upstream and low heat transfer regions downstream. This phenomenon passively strengthens the HV, as previously considered, while intensifying the turbulence in the channel. A quick peek at the figures revealed that the Nusselt number on the endwalls increases for both the modified endwall designs. The Nusselt number of the channel is increased by roughly 54-64% in the recessed-extruded endwall design compared to the flat endwall, or about twice as much as in the recessed endwall, which is increased by just 30-32%. Moreover, the friction factor of the recessed-extruded endwall drops slightly at Re of 7400 compared to the flat endwall. However, due to the flow restriction caused by the protrusions, all the other Re values show rises in f. At Re = 7400, the friction factor decreases by 2.6%; at all other Re values, i.e., Re = 14000, 21500, 29000, 36000, it increases by 1.2%, 3.1%, 4.0%, and 4.7%, respectively compared to the flat endwall. Overall, the enormous growth in the Nusselt number outweighs the growth in the friction factor, which leads to an overall rise in the HTEI at all the Re values. Accordingly, the overall increases in η range from 55.6% to 62.4%. While it is cut off halfway from the changed height, these statistics are around 1.7 times more than the gains in HTEI of the recessed endwall for all the Re numbers, which is about 33.9-36.7%.

For three cases of flat endwall and roughed endwall at Re = 21500, contours of the Nusselt number on the lower heated wall are shown in Fig. 14. Previous research [START_REF] Do | Numerical Investigation of Heat Transfer Characteristics of Pin-Fins with Roughed Endwalls in Gas Turbine Blade Internal Cooling Channels[END_REF] found that the indentations increase the amount of the flow-endwall interaction upstream of each pin-fin. Heat transfer is intensified by this action, providing momentum for the flow to advance backward. Expansion of the HV of the rows with the indentation, as seen in the 1 st , 3 rd , and 5 th rows of Fig. 14b, proves that this further boosts heat transfer in the rear rows of pin-fin arrays. As a result, the three-indentation recessed endwall design is replaced with a five-indentation recessed-extruded endwall design. As anticipated, the channel's heat transfer is increased in all rows, with the most significant differences in the 2 nd and 4 th rows.

In terms of the primary cause of the most heat transfer, the HV. Rows with indentations of recessed endwall, i.e., the 1 st , 3 rd , and 5 th rows, experience increased heat transfer due to the enlargement of HV regions, as was previously recognized [START_REF] Do | Numerical Investigation of Heat Transfer Characteristics of Pin-Fins with Roughed Endwalls in Gas Turbine Blade Internal Cooling Channels[END_REF]. On the other hand, the recessed-extruded endwall case has an indentation in each row, which increases the ability of this case to transfer heat across all rows. The HV of the recessed-extruded endwall are intensified in strength and widened in HV regions surrounding the pin-fins when comparing Figs. 14a and14c. In addition, as shown in Fig. 14c, the HV of the recessed-extruded endwalls are further increased in both the flow directions as well as to both sides of the pin-fins in the 3 rd and 5 th rows. This is because the small protrusions act like a rib cooling and lifting the flow above it, raising turbulent kinetic energy and passively enlarging the HV at back rows. When it comes to the wake vortices, in contrast to Fig. 14a, where the local vortices behind each pin-fin are observed at ease, it is found that the local vortices behind them have almost vanished in Figs. 14b and14c. This is because the HV continue to accelerate due to the endwall's upward narrowing geometry after passing through the pin-fins (Figs. 9b and9c). This action of the roughened endwalls prevents local vortices from forming behind each pin-fin. It lessens turbulence downstream of the pin-fins in the center to alleviate the loss of airflow that develops to the next rows.

Moreover, the upward endwall of the indentations engages in flow interaction and promotes heat transfer there. Additionally, the steeper endwall of the recessed-extruded endwall design produces increased heat transfer in this region compared to the recessed endwall, as evidenced in the 1 st row of Figs. 14b and14c. The higher Nusselt number in the 1 st row of Fig. 15a is primarily the result of this phenomenon. The jet-like beams between the pin-fins of Fig. 14c are hardly visible when the HV are magnified and are only discernible close to the side walls of the 2 nd and 4 th rows. However, in the case of recessedextruded endwall, the heat transfer of these beams is more significant than the recessed endwall case.

Figure 15 indicates the Nusselt number and pressure drop development for different rows. The Nusselt number by row is achieved by averaging the local Nusselt number by the area at both the upper and lower heated endwall of the current row. It is observed in Fig. 15a that the increment of averaged Nusselt number for the case of the recessed-extruded endwall is due to the enlargement of Nu in all rows. The highest Nusselt number of all rows appears in the 3 rd row, where the most turbulence is perceived with five pins. On the other hand, only the 1 st , 3 rd , and 5 th rows, where the indentations are, benefit considerably from the recessed endwall design. The Nusselt number increases by around 2.6 times for the recessed-extruded endwall compared to the flat endwall case in the 5 th row, which is the most prominent change. The minimum increase in Nu of the recessed-extruded endwall case is 2.17 times the Nusselt number of the flat endwall at the 4 th row. The maximum deviation for the recessed endwall is just 2.3 times in the 5 th row. The most significant difference between the two modified endwall is in the 4 th row, with a 44% increment for the recessed-extruded endwall design. The minimum change between these two new models in Nu is in the 5 th row, with a 12.8% increase for the recessed-extruded endwall. Furthermore, when considering the pressure drop variation, the recessed-extruded endwall induces more pressure drop in the odd rows than the recessed endwall case, along with increased heat transfer capacity. This phenomenon is due to the enlargement of HV in these rows of the recessed-extruded endwall design compared with the two other models. However, the even row's lower pressure drop of the recessed-extruded endwall is caused by the absence of protrusions in these rows. The pressure drop in the 2 nd row of the recessed-extruded endwall case is 18.6% higher than that of the flat endwall case, which is 15.8% lower than the case of the recessed endwall. This row has the highest pressure drop increase of the two altered endwall designs compared with the flat endwall design. Besides, the recessed-extruded endwall induces 8.5% and 8.3% reductions in pressure drop compared to the flat endwall at the 1 st and 5 th rows, correspondingly.

CONCLUSION

A computational study based on RANS analysis was performed to study the heat transfer of the pin-fin channel with roughed endwalls. The heat transfer characteristics, including Nusselt number, friction factor, and HTEI, were analyzed and compared with different cases.

It is achieved that the HV around the pin-fins are relatively enlarged with the recessed-extruded endwall design compared to the flat endwall at all tested Re numbers. The indentations enlarge the area of the HV around the pin-fins and protect them when they are on the verge of collapse. Meanwhile, the protrusions in the recessed-extruded endwall design inhibit the development of wake vortices which show low heat transfer. On the other hand, these protrusions function as rib turbulators, creating local vortices both upstream and downstream. These vortices raise the turbulence level in the channel and aid in enlarging the HV magnitude around the pin-fins. As a consequence, the channel with reference recessed-extruded endwalls has a better capacity for transferring heat, as seen by increases of 54.0%, 59.2%, 63.9%, 64.1%, and 64.3% in Nu when compared to the flat endwall at Re = 7400, 14000, 21500, 29000, and 36000, respectively. On the other hand, the growth in friction factor brought by the flow obstructions is insignificant. In general, the HTEI of the channel with this design is raised by 55.6%, 58.7%, 62.4%, 62.1%, and 61.9% compared to the flat endwall at Re = 7400, 14000, 21500, 29000, and 36000, respectively. In general, the increase in the HTEI of the recessedextruded endwall at Re = 21500 is approximately 1.7 times the increase in the HTEI of the recessed endwall design with double the recessed height.

In general, the recessed endwall design is more suited for the turbine, which requires less friction while yet having the ability to endure deeper cuts into the surface, with an increase of 36.7% in HTEI, higher than the flat endwall. In contrast, it is concluded that the recessed-extruded endwall elevates the HTEI to a peak value of 62.4% increase, higher than the flat endwall, but also raises the friction factor, even though it is minor. However, the design has a reduced recessed height, which is half the height ratio of the reference recessed endwall case.

  Fig. 1 shows a schematic diagram of the tested channel, where the pin-fin diameter (D = 8mm = 0.315inches) is used to indicate all dimensions. The length of the tested channel is Lt = 87.5D, the height is H = 2.5D, and the width is W = 10D. The tested channel comprises four sections: three smooth sections, which include a development section at the entrance, a smooth heated section, an output section at the outlet, and a test section measuring L = 12.5D in length. The development section with a length of Ld = 12.5D and the plain heated channel with a length of Lp = 37.5D assure homogeneity of the entry flow for the pin-fin domain. The output channel is designed to reduce backflow productively. Five rows of pin-fins are positioned in a staggered pattern in the test section. Each pinfin is positioned with the neighbors with the spacings, Sx/D = 2 and Sy/D = 2. The cross-section of the computed domain is a rectangle with an aspect ratio AR = W/H = 4. The upper and lower walls correspond to z = 0 and z = H, and the two side walls are located at y = -W/2 and y = W/2. The equivalent diameter of the channel, Dh = 2HW/(H+W), is 32mm. The direction of the flow corresponds to the direction of the x-axis.

Figure 2 .

 2 Figure 2. Demonstration of all designs and the design parameters of the recessed-extruded endwall

Figure 4 . 5 .

 45 Figure 4. Influence of mesh size on the result of Nusselt number at Re = 21500

Figure 9 .

 9 Figure 9. Velocity contours at the middle section (y = 0 mm) and streamlines near the endwall of the channel at Re = 21500

Figure 10 .

 10 Figure 10. Velocity contours at the middle section (z = 10mm) of the channel at Re = 21500

Figure 11 .

 11 Figure 11. Velocity contours near lower endwall at Re = 7400

Figure 12 . 3 . 4

 1234 Figure 12. Velocity contours near lower endwall at Re = 36000

Figure 13 .

 13 Figure 13. Comparison of the new designs with the flat endwall

Figure 14 .Figure 15 .

 1415 Figure 14. Nusselt number contours of flat endwall and roughed endwall cases on the lower heated wall at Re = 21500

ACKNOWLEDGEMENTS

This research is funded by Hanoi University of Science and Technology (HUST) under project number T2022-PC-018, and the cooperation research between HUST and Viettel Aerospace Institute (VTX).


Jianhao Zhang

Carlos Alonso-Ramos

Laurent Vivien

Sailing He

Eric Cassan

Self-adaptive waveguide boundary for wideband multi-mode four-wave mixing

Keywords: Nonlinear optics, phase matching, index manipulation, subwavelength

We propose a new approach to provide wideband multi-mode four-wave mixing, independent of the intrinsic waveguide dispersion. We adopt concepts from quantum mechanics and sub-wavelength engineering to design an effective photon well, with a graded potential along the waveguide cross section, that provides flexible control over the mode confinement. The selfadaptive nature of the waveguide boundary allows different spatial modes with equi-spaced frequencies and shared propagation wavevector, automatically fulfilling both, energy conservation and wavevector phase matching conditions. Capitalizing on this concept, we show

INTRODUCTION

Nonlinear processes like four-wave mixing (FWM) have garnered a great interest due to their unique capabilities for on-chip light generation, with an immense potential for the implementation of wideband sources for silicon photonics [START_REF] Cassan | Nonlinear Fiber Optics[END_REF][START_REF] Hansryd | Fiber-based optical parametric amplifiers and their applications[END_REF][START_REF] Leuthold | Nonlinear silicon photonics[END_REF][START_REF] Borghi | Nonlinear silicon photonics[END_REF]. Harnessing 3rd nonlinear Kerr effect in silicon already allowed the demonstration of promising frequency combs [START_REF] Xue | Second-harmonic-assisted four-wave mixing in chip-based microresonator frequency comb generation[END_REF][START_REF] Singh | Octave-spanning coherent supercontinuum generation in silicon on insulator from 1.06μm to beyond 2.4μm[END_REF][START_REF] Stern | Battery-operated integrated frequency comb Generator[END_REF][START_REF] Guo | Mid-infrared frequency comb via coherent dispersive wave generation in silicon nitride nanophotonic waveguides[END_REF], optical parametric amplification [START_REF] Salem | Signal regeneration using low-power four-wave mixing on silicon chip[END_REF][START_REF] Foster | Broad-band optical parametric gain on a silicon photonic chip[END_REF][START_REF] Foster | Broad-band continuous-wave parametric wavelength conversion in silicon nanowaveguides[END_REF][START_REF] Turner-Foster | Frequency conversion over two-thirds of an octave in silicon nanowaveguides[END_REF] and mid-infrared light sources in silicon [START_REF] Liu | Mid-infrared optical parametric amplifier using silicon nanophotonic waveguides[END_REF][START_REF] Zlatanovic | Mid-infrared wavelength conversion in silicon waveguides using ultracompact telecom-band-derived pump source[END_REF][START_REF] Liu | Bridging the mid-infrared-to-telecom gap with silicon nanophotonic spectral translation[END_REF]. A great effort has been devoted to compensate both the intrinsic material dispersion and the nonlinearityinduced dispersion, that hamper phase matching which is key to maximize the efficiency of FWM processes [START_REF] Hansryd | Fiber-based optical parametric amplifiers and their applications[END_REF]. Nevertheless, achieving broadband phase-matching in simple-to-fabricate and fabrication-tolerant silicon waveguides remains an open question. Optimization of transversal dimensions of conventional strip waveguides, provided phase-matching in a relatively narrow wavelength ranges, thus compromising the bandwidth of the nonlinear wavelength conversion processes [START_REF] Salem | Signal regeneration using low-power four-wave mixing on silicon chip[END_REF][START_REF] Foster | Broad-band optical parametric gain on a silicon photonic chip[END_REF][START_REF] Foster | Broad-band continuous-wave parametric wavelength conversion in silicon nanowaveguides[END_REF][START_REF] Turner-Foster | Frequency conversion over two-thirds of an octave in silicon nanowaveguides[END_REF][START_REF] Liu | Mid-infrared optical parametric amplifier using silicon nanophotonic waveguides[END_REF][START_REF] Zlatanovic | Mid-infrared wavelength conversion in silicon waveguides using ultracompact telecom-band-derived pump source[END_REF][START_REF] Liu | Bridging the mid-infrared-to-telecom gap with silicon nanophotonic spectral translation[END_REF]. Broadband phase-matching has been shown based on optimization of high order diffraction terms, e.g. 4th order waveguide dispersion [START_REF] Foster | Broad-band continuous-wave parametric wavelength conversion in silicon nanowaveguides[END_REF], or by implementing rib geometries [START_REF] Turner-Foster | Frequency conversion over two-thirds of an octave in silicon nanowaveguides[END_REF]. Yet, the proposed solutions require complex fabrication processes, with deposition of different materials, or tight control of rib and slab thicknesses. On the other hand, photonic crystal (PhC) waveguides with flexible dispersion properties, high confinement and large group velocities were also considered as a FWM photonic platform [START_REF] Mcmillan | Observation of four-wave mixing in slow-light silicon photonic crystal waveguide[END_REF][START_REF] Monat | Four-wave mixing in slow light engineered silicon photonic crystal waveguides[END_REF]. However, the bandgap guiding in PhC waveguides results in a limited bandwidth and a high sensitivity to fabrication imperfections that substantially increase propagation loss and distort dispersion properties. An alternative approach to dispersion-engineered waveguides is the use of photonic cavities to effectively enhance the light-matter nonlinear interactions by resonant enhancement and phase-matching among different cavity modes modes. This method was demonstrated in various configurations, including single rings [START_REF] Xue | Second-harmonic-assisted four-wave mixing in chip-based microresonator frequency comb generation[END_REF][START_REF] Singh | Octave-spanning coherent supercontinuum generation in silicon on insulator from 1.06μm to beyond 2.4μm[END_REF][START_REF] Stern | Battery-operated integrated frequency comb Generator[END_REF][START_REF] Guo | Mid-infrared frequency comb via coherent dispersive wave generation in silicon nitride nanophotonic waveguides[END_REF][START_REF] Azzini | Ultra-low power generation of twin photons in a compact silicon ring resonator[END_REF], coupled rings [START_REF] Morichetti | Travelling-wave resonant four-wave mixing breaks the limits of cavity-enhanced alloptical wavelength conversion[END_REF][START_REF] Zeng | Design of triply-resonant microphotonic parametric oscillators based on Kerr nonlinearity[END_REF][START_REF] Zeng | Four-wave mixing in silicon coupled-cavity resonators with port-selective, orthogonal supermode excitation[END_REF],

and coupled nanobeam cavities [START_REF] Azzini | Stimulated and spontaneous four-wave mixing in silicon-on-insulator coupled photonic wire nano-cavities[END_REF][START_REF] Lin | High-efficiency degenerate four-wave mixing in triply resonant nanobeam cavities[END_REF]. Concurrently, PhC cavities with parabolic mirror designs were also proposed for nonlinear frequency comb generation [START_REF] Sumetsky | Surface nanoscale axial photonics[END_REF][START_REF] Oreshnikov | Multiple nonlinear resonances and frequency combs in bottle microresonators[END_REF][START_REF] Suchkov | Frequency comb generation in SNAP bottle resonators[END_REF][START_REF] Combrie | Comb of high-Q Resonances in a Compact Photonic Cavity[END_REF]. Though these demonstrations showed promising performance, their broadband operation is hampered by tradeoffs in phase matching and spatial overlaps. Fig. 1. Schematic of degenerate four wave mixing operating in a single-mode waveguide with anomalous dispersion (a) and operating within an inter-modal scheme regardless the absolute dispersion provided that all dispersion curves are obtained by translating the same curve with a constant frequency step (e.g. 21=10 here).

Here, we propose a new strategy to satisfy both phase matching and energy conservation conditions in an ultra-broad wavelength range. Rather than tuning the waveguide dimensions to yield anomalous dispersion, we shape the index profile of the waveguide to support different spatial modes with the same propagation constant and equal frequency spacing, ensuring phase matching and energy conservation, respectively. In the conventional single-mode FWM approach, depicted in Fig. 1(a), energy is transferred from the pump into signal and idler propagating in the fundamental waveguide mode with different propagation constants. Thus, precise control of dispersion is required to fulfill phase matching condition (2kzP=kzS+kzI). In the proposed multi-mode FWM approach, see Fig. 1(b), energy is transferred from the pump into signal and idler propagating in different waveguide modes with the same propagation constant (kzP=kzS=kzI). Therefore, the phase matching condition is automatically satisfied. Concurrently, energy conservation requires equal frequency spacing. Then, the bandwidth of the proposed scheme does not depend on the exact dispersion of the waveguide, but on the relative slopes of the dispersion curves of the modes, determining the wavelength range where energy conservation is fulfilled. We adopt concepts from quantum mechanics and sub-wavelength engineering to design an effective photon well with a graded potential along the waveguide cross section, yielding a different spatial confinement for each waveguide mode. This self-adaptive boundary simultaneously provides equal frequency spacing, parallel dispersion curves and large spatial overlap among different waveguide modes, overcoming the major bandwidth and efficiency limitations of conventional approaches.

Compared to previous demonstrations using inter-modal FWM [START_REF] Cheng | Intermodal four-wave mixing in a higher-order-mode fiber[END_REF][START_REF] Friis | Inter-modal four-wave mixing study in a two-mode fiber[END_REF][START_REF] Signorini | Intermodal four-wave mixing in silicon waveguides[END_REF][START_REF] Kittlaus | On-chip inter-modal Brillouin scattering[END_REF], our work presents a more systematic and controllable strategy with clear analytical explanation. Unlike previous works that relied on parabolic photonic wells along the propagation direction in photonic cavities [START_REF] Combrie | Comb of high-Q Resonances in a Compact Photonic Cavity[END_REF], we exploit a self-adaptive boundary effect in waveguides with a photonic well in the transverse direction to light propagation. This method facilitates broadband wavelength operation compared to the cavity-based four-wave mixing. We also demonstrate that a parabolic guide profile is not essential to obtain a quasi-universal phase matching condition in multimode waveguides with propagative modes. Last, we also show the remarkable flexibility of the proposed approach, allowing large conversion span, with a material-independence design method compatible with sub-wavelength and/or suspended membrane waveguides.

FREQUENCY SPACINGS OF A STEP-INDEX TWO-DIMENSION (SLAB)

WAVEGUIDE

In order to present the step-by-step concepts, we start here with the simplified case of a slab step-index waveguide, while realistic 3D waveguides corresponding to photonic integration standards are discussed in sections 3 and 4. Let us consider the two-dimension (infinite depth along the x axis) step-index slab waveguide presented in the inset of Fig. 2(a). The eigen equation for modes propagating along the 𝑧 axis with an electric field polarized along 𝑦 axis reads as [START_REF] Okamoto | Fundamentals of Optical Waveguides[END_REF]:

ℎ𝑎 = 𝑚𝜋 2 + 𝑎𝑟𝑐𝑡𝑎𝑛( 𝛾𝑛 𝑤 2 ℎ𝑛 𝑐 2 ) (1)
where 𝑛 𝑤 and 𝑛 𝑐 are the material index of waveguide core and cladding, while ℎ = √(𝑘 0 2 𝑛 𝑤 2 -𝑘 z 2 ) and 𝛾 = √(𝑘 z 2 -𝑘 0 2 𝑛 𝑐 2 ) are the wavevector along the 𝑦 axis, inside and outside the waveguide core, respectively. Silicon and silicon dioxide are chosen here as the core and cladding materials, considering material dispersion. The waveguide width 2𝑎 was set at 700nm which is enough to support 4 modes with effective index values higher than 2. From Equ. (1) we obtain the dispersion curves for the first four order modes (with corresponding mode order 𝑚=0, 1, 2, 3), shown Fig. 2(a). The frequency spacings change in a nonlinear manner with the mode order 𝑚, resulting in an uneven spacing that precludes the satisfaction of energy conservation required for the proposed multi-mode FWM approach. Schematic of a photonic well described by the cut-off frequency for photon propagating along 𝑧 axis with wavevector 𝑘 𝑧 .

To better understand the evolution of frequency spacing ∆𝜔 𝑚 , we consider the derivative of the mode frequency with respect to the mode order. which gives a good explanation on the monotonous and nonlinear evolution of ∆𝜔 𝑚 observed in Fig. 2(b). This confirms that the analytical frequency spacing 𝑑𝜔 𝑑𝑚 is a useful and simple tool to investigate how the frequency spacing evolves with the waveguide dimension and the index profile.

The step-index slab waveguide can be understood from the point of view of quantum wells, just by considering the cut-off frequencies for the core (𝑘 𝑧 𝑐/𝑛 𝑤 ) and cladding (𝑘 𝑧 𝑐/𝑛 𝑐 ), for a given wavevector 𝑘 𝑧 [START_REF] Eichenfield | Optomechanical crystals[END_REF][START_REF] Alpeggiani | Effective bichromatic potential for ultrahigh Q-factor photonic crystal slab cavities[END_REF][START_REF] Simbula | Realization of high-Q/V photonic crystal cavities defined by an effective Aubry-André -Harper bichromatic potential[END_REF]. Then, the behavior of frequency spacings is similar to the solutions of harmonic oscillator in a finite-depth potential well. The square potential well formed by the stepindex waveguide, depicted in Fig. 2(d), results in unevenly spaced frequencies. Conversely, it is well known that parabolic potential wells yield equi-spaced frequencies [START_REF] Combrie | Comb of high-Q Resonances in a Compact Photonic Cavity[END_REF]. Therefore in the following section we study the influence of the index profile in the frequency separation of the modes.

FREQUENCY SPACINGS OF A GRADED-INDEX SLAB WAVEGUIDE WITH SELF-ADAPTIVE BOUNDARY

In this section we study the frequency spacing ∆𝜔 𝑚 for a 2D graded-index slab waveguide, with a nonuniform index profile described by 𝑛(y) = (𝐴 + 𝐵𝑦) 𝑝 . By changing the value of 𝑝, the index profile can be made linear (𝑝 = 1) or nonlinear (𝑝 ≠ 1). The waveguide width is 2𝑎, and the indices for the waveguide core, waveguide edge and cladding are 𝑛(0) = 𝑛 𝑐𝑒𝑛𝑡 , 𝑛(𝑎) = 𝑛 𝑏 , and 𝑛(|𝑦| > 𝑎) = 𝑛 𝑐 , as displayed in Fig. 3(a). From the analytical solution of the wave equation for this graded-index waveguide [START_REF] Okamoto | Fundamentals of Optical Waveguides[END_REF] (see Supplementary Information), it follows that the field distribution inside the waveguide core follows a cosinus-like profile while the guide index is larger than the mode effective index (𝑛(𝑦) > 𝑛 𝑒𝑓𝑓 ) and an exponentially decaying profile otherwise (𝑛(𝑦) < 𝑛 𝑒𝑓𝑓 ). An effective mode width, [-𝐿 𝑚 , 𝐿 𝑚 ] can be defined as the region where the 𝑛(𝑦) > 𝑛 𝑒𝑓𝑓 , i.e. where mode profile exhibits a cosinus-like shape. Then, as each mode has a different effective index, each mode also has a different effective width. This different confinement effect is illustrated in Fig

3(a)
, that shows the analytical field profile for the first three waveguide modes. Such a wave confining method, based on the index profile with condition 𝑛 𝑒𝑓𝑓 > 𝑛 𝑏 , is then defined as "selfadaptive boundary (SAB)" hereafter. I propose to move the explanations of eq. 3-11 to the supplementary information. First, we consider a linear index profile, i.e. 𝑝 = 1 and 𝑛(y) = 𝐴𝑦 + 𝐵 . With this profile, the frequency spacings can be described as (see Supplementary Information S2): In the following section we present several realistic implementations of this type of multimode graded waveguides, chosen to illustrate the method. However, we would like to highlight the very general nature of the proposed approach, which can be seamlessly adapted to all types of geometries, photonic platforms (Si, III/V, polymer guides, etc.), and spectral ranges (near infrared, medium infrared).

𝑑𝜔 𝑑𝑚 = -𝐴𝜋𝑐/ [𝑛 𝑒𝑓𝑓 2 𝑙𝑜𝑔 (√ 𝑛 𝑐𝑒𝑛𝑡 2 𝑛 𝑒𝑓𝑓,𝑚 2 -1 + 𝑛 𝑐𝑒𝑛𝑡 𝑛 𝑒𝑓𝑓,𝑚 ) + 𝑛 𝑐𝑒𝑛𝑡 √𝑛 𝑐𝑒𝑛𝑡 2 -𝑛 𝑒𝑓𝑓,𝑚 2 ] (3)

GRADED-PROFILE THREE-DIMENSION WAVEGUIDE WITH SELF-ADAPTIVE BOUNDARY

The self-adaptive boundary presented in previous section, which conceptually inherits from the ability of non-uniform potential well, provides equi-spaced frequency modes sharing the same propagation constant, as shown in Fig. 4. One further step is to verify this method in a more practical situation, i.e. with three-dimensional waveguides. As a large index change along the waveguide cross-section is not always feasible from the material point of view, an alternative approach is to implement this graded index profile through nanostructured subwavelength grating (SWG) engineering. By periodically combining high-index and low-index sections with a pitch shorter than half of the wavelength, subwavelength gratings allow the implementation of a material with a synthetic refractive index that can be tuned at will between those of the high and low index materials [START_REF] Ortega | Analysis of "Quasi-Modes" in Periodic Segmented Waveguides[END_REF][START_REF] Halir | Waveguide sub-wavelength structures: a review of principles and applications[END_REF][START_REF] Cheben | Subwavelength integrated photonics[END_REF]. In Fig. 4(a), the unit cell of two types of periodically-structured subwavelength waveguides are presented for illustration. Here, the gradual index variation is implemented by apodization of the waveguide length or size of the engraved holes. Based on this index-equivalent effective material method and the Marcatili's waveguide approximation, we are again able to semi-analytically investigate the three-dimensional waveguide (Detailed in Supplementary Information S3). ). We have to mention that this result came from the operation of first 3 modes as a demonstration. However, this design strategy can be scaled up to higher order modes, e.g. 5 modes. As mentioned earlier, these capabilities arise from the self-adaptive boundary (SAB) condition 𝑛 𝑒𝑓𝑓,𝑚 > 𝑛 𝑏 , ∀𝑚 which is more than a simple graded-profile condition. In the Supplementary Material S4, we analytically show that a multimode gradedindex waveguide with a classical condition 𝑛 𝑒𝑓𝑓 < 𝑛 𝑏 cannot equalize frequency spacings between its modes since it adopts the same effective boundary for different modes as in a stepindex waveguide.

In addition to the demonstrated capability and flexibility of the proposed approach for degenerate four-wave mixing, the nonlinearities-induced phase mismatch can be considered here as well.

The phase matching considering nonlinearities is described [START_REF] Hansryd | Fiber-based optical parametric amplifiers and their applications[END_REF][START_REF] Foster | Broad-band continuous-wave parametric wavelength conversion in silicon nanowaveguides[END_REF] as ∆𝑘 = 2𝛾𝑃 𝑃 -(2𝑘 𝑧𝑃 -𝑘 𝑧𝑆 -𝑘 𝑧𝐼 ), which is governed by the nonlinear part 2𝛾𝑃 𝑃 and linear part ∆𝑘 𝐿 = 2𝑘 𝑧𝑃 -𝑘 𝑧𝑆 -𝑘 𝑧𝐼 . Since 2𝛾𝑃 𝑃 is normally positive in silicon, the linear dispersion ∆𝑘 𝐿 needs to be a bit larger than zero to fulfill the global phase matching condition, which is classically addressed by tuning the dispersion to its anomalous regime in a classical waveguide. In contrast, in our case, the condition is translated to the fact that: the frequency spacing ∆ω SP should be slightly different from ∆ω PI , which can be easily achieved just by slightly moving the operating point of pump wave, as illustrated in Fig. 5(a).

If the condition that ∆ω SP = ∆ω 21 - the effective nonlinear Kerr nonlinearity γ can be described [START_REF] Serna | Suspended SOI waveguide with sub-wavelength grating cladding for mid-infrared Experimental GVD engineering in slow light slot photonic crystal waveguides[END_REF][START_REF] Santagiustina | Theory of slow light enhanced four-wave mixing in photonic crystal waveguides[END_REF]

∭ 𝐸 𝑦𝑆 •𝐸 𝑦𝑃 •𝐸 𝑦𝑃 * •𝐸 𝑦𝐼 * 𝑑𝑟 3 ∭ |𝐸 𝑦𝑃 | 2 𝑑𝑟 3 √∭ |𝐸 𝑦𝑆 | 2 𝑑𝑟 3 •√∭ |𝐸 𝑦𝐼 | 2 𝑑𝑟 3
in a single unit cell [START_REF] Signorini | Intermodal four-wave mixing in silicon waveguides[END_REF], which gives a value of 0.12.,𝐸 𝑦𝑆 , 𝐸 𝑦𝑃 and 𝐸 𝑦𝐼 corresponding to the field of signal wave, pump wave and idler waves. ) obtained from Fig. 5(b), we can predict an 3dB tunable bandwidth of around 1THz (~10nm) for a 1mm long waveguide. 𝑉 𝑔2 and 𝑉 𝑔0 are corresponding to the 2nd and fundamental modes, respectively By generalizing this approach, the 3dB tunable bandwidth as a function of waveguide length can be plotted in Fig. 5(d). It can therefore be seen that, even in a non-optimized waveguide, the spectral operating band of the FWM process is as wide as few tens of nanometers for sub-millimeter long waveguide. With all these results, we have successfully demonstrated the feasibility of using a multimode SWG waveguide with self-adaptive boundary for degenerate FWM, based on the interaction among the first three modes. However, nonlinear processes are not strictly limited to this scenario. Using higher order modes (e. g. fundamental mode, 2nd mode and 4th mode), FWM process could be further scaled to even larger conversion spans by using the same method, for example to directly convert near-infrared waves to nearly mid-infrared waves. Most importantly, since the proposed approach no longer requires specific normal/anomalous dispersion, we can easily adopt this approach to any other platforms, no matter what kinds of dispersion the waveguide modes exhibit.

As an illustrative example of the highly adaptable nature of the method to different platforms with different thicknesses, targeting different frequency ranges, etc, we investigated its application to different situations. Using, 3D FDTD method and considering material dispersion, we calculated the dispersion curves for Si sub-wavelength waveguides with Si thicknesses ranging between 220nm and 600nm. For the 220nm-thick Si waveguide with parameters 𝑎=700nm, 𝑏=300nm, the optimized position for ∆ω 21 =∆ω 10 is located at wavevector around 𝑘 𝑧 = 0.95 × 10 7 m-1, leading to a conversion span from 220 THz to 180 THz, as shown in Fig. 6(a). By raising 𝑎 and reducing 𝑏 to 750nm and 250 nm, we have been able to shift the optimized position continuously to 𝑘 𝑧 = 1.53 × 10 7 m-1,with starting and stop frequencies of 280 THz to 245 THz, respectively. In order to push the operation point to a lower frequency, with the purpose of generating light wavelengths up to 2𝜇𝑚, larger waveguide cross-sections can be considered. Very interestingly, by simply enlarging the thickness to 600nm, the working frequency for the fundamental mode can be shifted to around 155THz, with conversion span over 40THz (from 1.53μm to 1.93μm, i.e. 400 nm), as shown in Fig. 6(c), with almost no displacement on 𝑘 𝑧 . Similar result can also be observed in another configuration (𝑎=750nm, 𝑏=250nm), with a slightly shifted working point 𝑘 𝑧 .These gathered results unambiguously show the simplicity and flexibility of our approach for degenerate four-wave mixing, i. e., 1) for each waveguide thickness the waveguide is capable to offer adjustable working conditions in a wide frequency range, within a varying index-profile; 2) for a certain optimal lineshape, the strategy for shifting the working wavelength, is to simply adjust the thickness.

CONCLUSION

With thorough analytical explanation and numerical confirmation, we show that waveguides with material-engineered graded index profile (unambiguous condition that 𝑛 𝑒𝑓𝑓 > 𝑛 𝑏 ) can be designed to support the modes that adapt themselves to different effective boundaries ("selfadaptive boundary"), according to their effective index values. This self-adaptive behavior, that relies on similar concepts as those found for potential wells, provides new degrees of freedom to achieve simultaneous satisfaction of both the energy conservation and phase matching, regardless the intrinsic dispersion of the considered optical waveguide modes. This strategy, that can be adapted to different wavelength ranges and material platforms, opens a new design space for degenerate FWM. Capitalizing on this concept, we show that phase-matching condition can be satisfied over 400nm (bridging from telecom wavelength to almost 2μm) by employing three waveguide modes. This range could be further expanded by using higher order modes enabled by flexible index profile optimization. We foresee that the self-adaptive boundary concept will expedite the development of a new generation of nonlinear circuits with an immense potential for light generation in mid-infrared wavelength, but also for nonlinear process beyond FWM and applications in which dispersion manipulation is of major relevance.

SUPPLEMENTARY MATERIAL

S1 Calculation for a step-index two-dimensional waveguide; S2 Calculation for a gradedindex two-dimension waveguide with self-adaptive boundary (𝑛 𝑒𝑓𝑓 > 𝑛 𝑏 ); S3 Semi-analytical calculation for a graded-profile three-dimensional waveguide with self-adaptive boundary (𝑛 𝑒𝑓𝑓 > 𝑛 𝑏 ); S4 Discussion on a graded-index waveguide with classical condition 𝑛 𝑒𝑓𝑓 < 𝑛 𝑏 .

ACKNOWLEDGMENT

As introduced in the text, for a two-dimension (infinite depth along 𝑥 axis) step-index waveguide, the eigen equation for modes propagating along 𝑧 axis with electric field polarizing along 𝑦 axis can be described as ℎ𝑎 =

𝑛 𝑏 2 √(𝑛 𝑒𝑓𝑓 2 -𝑛 𝑐 2 ) 𝑛 𝑐 2 √(𝑛 𝑏 2 -𝑛 𝑒𝑓𝑓 2 )
) accounts for the abrupt index change on the boundary. The related first derivatives can be written as:

𝑑𝑓 𝑚1 𝑑𝑛 𝑒𝑓𝑓 = 2𝑘 𝑧 𝜋 𝑑 𝑑𝑛 𝑒𝑓𝑓 (𝑎√ 𝑛 𝑤 2 𝑛 𝑒𝑓𝑓 2 -1) = -2𝑎𝑘 𝑧 𝑛 𝑤 2 𝜋𝑛 𝑒𝑓𝑓 2 √𝑛 𝑤 2 -𝑛 𝑒𝑓𝑓 2 (S3) 𝑑𝑓 𝑚2 𝑑𝑛 𝑒𝑓𝑓 = - 2 𝜋 • 𝑑 𝑑𝑛 𝑒𝑓𝑓 𝑎𝑟𝑐𝑡𝑎𝑛( 𝛾𝑛 𝑤 2 ℎ𝑛 𝑐 2 ) = - 2 𝜋 • 1 1+ 𝑛 𝑏 4 (𝑛 𝑒𝑓𝑓 2 -𝑛 𝑐 2 ) 𝑛 𝑐 4 (𝑛 𝑏 2 -𝑛 𝑒𝑓𝑓 2 ) • 𝑑 𝑑𝑛 𝑒𝑓𝑓 ( 𝑛 𝑏 2 √(𝑛 𝑒𝑓𝑓 2 -𝑛 𝑐 2 ) 𝑛 𝑐 2 √(𝑛 𝑏 2 -𝑛 𝑒𝑓𝑓 2 ) ) = - 2 𝜋 • 𝑛 𝑒𝑓𝑓 -1+𝑛 𝑒𝑓𝑓 2 (1/𝑛 𝑐 2 +1/𝑛 𝑤 2 ) • 1 √(𝑛 𝑒𝑓𝑓 2 -𝑛 𝑐 2 ) • 1 √(𝑛 𝑤 2 -𝑛 𝑒𝑓𝑓 2 ) (S4)

S2. CALCULATING THE 𝒅𝝎/𝒅𝒎 OF A GRADED-INDEX TWO-DIMENSION WAVEGUIDE WITH SELF-ADAPTIVE BOUNDARY (𝒏 𝒆𝒇𝒇 > 𝒏 𝒃 )

We consider a 2D graded-index slab waveguide (then with graded cut-off potential) with a waveguide width of 2𝑎. The indices for the waveguide core, waveguide edge and cladding are 𝑛(0) = 𝑛 𝑐𝑒𝑛𝑡 , 𝑛(𝑎) = 𝑛 𝑏 and 𝑛 𝑐 , respectively, as displayed in Fig. S1. we propose to systematically rely on a nonuniform index profile 𝑛(𝑦). When we introduce a condition that the effective indices of the modes are larger than the physical boundary index (i.e. 𝑛 𝑒𝑓𝑓 > 𝑛(𝑎) = 𝑛 𝑏 for all the guided modes) into the index profile, the waveguide can be considered as splitting

𝐻 𝑥 = 𝐴 1 -𝑐𝑜𝑠ℎ[𝐷(𝑦)] + 𝐵 1 -𝑠𝑖𝑛ℎ[𝐷(𝑦)], -𝑎 < 𝑦 ≤ 𝐿 (S9) 𝐻 𝑥 = 𝐴 2 -𝑒𝑥𝑝 [𝑘 0 √𝑛 𝑒𝑓𝑓 2 -𝑛 𝑐 2 • (𝑦 + 𝑎)] , -𝑎 ≥ 𝑦 (S10)
in which the 𝐴 1,2 + and 𝐵 1,2 + are the amplitudes of the decaying components for the positive direction, while 𝐴 1,2 -and 𝐵 1,2 -stand for the negative direction. 𝜑 is the biased phase that is related to the mode order (𝜑 = 0 or 𝜑 = 𝜋 2 ⁄ , with respect to modes of symmetric and anti-symmetric parities). The pattern of zone 0 and zones, ∓2 jointly determine the hyperbolic form of part 1.

Due to the varying index profile, currently the total phases of the wave propagation along the 𝑦 axis inside the waveguide can be transformed to the integrals of position-related wavevector √(𝑘 0 2 𝑛 2 (𝑦) -𝑘 𝑧 2 ) to 𝑦 coordinate (𝑁 𝑚 and 𝐷 𝑚 stand with respect to zone 0 and zones ∓1, respectively): As an example of the possible design of three-dimensional graded-profile waveguides with a self-adaptive boundary, a subwavelength grating (SWG) waveguide with a period of 150nm is chosen, which the geometry is depicted in Fig. S2. The material considered here are silicon and air, for the high-index and low-index regions, respectively. As the SWG is able to be considered as being made of equivalent materials, the strategy here to achieve the graded-profile (gradedindex) waveguide is tapering the component width of silicon, from the center to the edge, while keeping the period fixed. In this case, the width of silicon components is changed from 150nm to 40nm, i. e. adjusting the filling factor from 1 to 0.267. According to the well-known equivalent properties of waveguide, the index of SWG waveguide can be described as:

𝑁 𝑚 (𝑦) = ∫ √(𝑘 0 2 𝑛 2 (𝑦) -𝑘 𝑧 2 ) 𝐿 𝑚 0 𝑑𝑦 (
𝑛 𝑆𝑊𝐺 (𝑦) = √𝜂𝑛 𝑐𝑒𝑛𝑡 2 + (1 -𝜂)𝑛 𝑐 2 = √𝜂(𝑛 𝑐𝑒𝑛𝑡 2 -𝑛 𝑐 2 ) + 𝑛 𝑐 2 (S19)
which leads to an index ranging from around 3.48 to 1.8 for the C-band wavelengths. The 𝜂(𝑦) = (𝐴𝑦 + 𝐵) 𝑜𝑟𝑑𝑒𝑟 is the filling factor of silicon of the SWG waveguide at a transverse position 𝑦. The schematic of potential well that describes e cut-off frequency, as introduced before, is presented in the inset of Fig. S2. In order to approximately investigate the frequency spacings of this three-dimension (3D) SWG waveguide, the Marcatili's method [START_REF] Cassan | Nonlinear Fiber Optics[END_REF] is adopted to find an equivalent simpler 2D waveguide with eight cladding regions, then the phase item of the dispersion modes Equ. (S14) can be rewritten as:

𝑁(𝐿 𝑚 ) = ∫ √(𝑘 0 2 𝑛 𝑆𝑊𝐺 2 (𝑦) -𝑘 𝑥 2 -𝑘 𝑧 2 ) 𝐿 𝑚 0 𝑑𝑦 = 𝑚𝜋 2 + 𝜋 4 (S20) 𝑘 𝑥 𝑡 = 𝑎𝑟𝑐𝑡𝑎𝑛 ( 𝛾 𝑥 𝑘 𝑥 ) (S21) = √(1 -𝐷)[𝜂(𝑛 𝑐𝑒𝑛𝑡 2 -𝑛 𝑐 2 ) + 𝑛 𝑐 2 ] -𝐶 (S23)
With this equivalent index, (S12) can be simplified to the general form:

𝑁(𝐿 𝑚 ) = ∫ √(𝑘 0 2 𝑛 2 (𝑦) -𝑘 𝑧 2 ) 𝐿 0 𝑑𝑦 = 𝑚𝜋 2 + 𝜋 4 (S24)
Now the phase item shares the same expression as (S14). Thus, we can adopt the method discussed before to optimize the waveguide. Very importantly, the effective index 𝑛 𝑒𝑓𝑓 is shifted to a lower level due to the limited thickness, which gives a maximum equivalent index of 𝑛(0)=√(1 -𝐷)[(𝑛 𝑐𝑒𝑛𝑡 2 -𝑛 𝑐 2 ) + 𝑛 𝑐 2 ] -𝐶.

In order to ease the effect coming from 𝑛 𝑥 and simultaneously promote the effective index values of all the modes, we introduce additionally a shared section in which the index is fixed as the central section shown in Fig. S2, with a width of 2𝑏. By setting this, a minimum phase can be ensured and the total phase can be rewritten as: The only thing that matters is the absolute value frequency spacing. This result unambiguously confirms us that, without the assistance of a variable "effective width" supported by the selfadaptive boundary, there is no way to reshuffle the frequency spacings in the multi-mode scheme, since the integral range has a significant effect on the left part of eigen dispersion equation.

(1) K. Okamoto, "Fundamentals of Optical Waveguides," Elsevier 2006.

Figure 2 (

 2 b) depicts the evolution of the frequency spacing ∆𝜔 𝑚 = 𝜔 𝑚+1 -𝜔 𝑚 and effective index as a function of the mode order. With the increasing mode order m (thus decreasing 𝑛 𝑒𝑓𝑓 ), the frequency spacing ∆𝜔 𝑚 monotonously and rapidly increases.

Fig. 2 .

 2 Fig. 2. (a) Dispersion curves of first 4 modes of a two-dimension silicon waveguide with silica

Fig. 3 .

 3 Fig. 3. (a) Sketch of a non-uniform index profile and the mode distribution of the first three

Figure 3 (

 3 Figure 3(b) shows d𝜔/d𝑚 as a function of 𝑛 𝑒𝑓𝑓 , calculated from Equ. (3), considering

Fig. 4 .

 4 Fig. 4. (a) Sketches of a photonic well described by the cut-off frequency for photon propagating

2 =

 2 ∆ω SI can be satisfied at the shifted position characterized by 2𝑘 𝑧𝑃 > (𝑘 𝑧𝑆 + 𝑘 𝑧𝐼 ), then the nonlinearity-induced phase mismatch can be compensated, as in Fig.5(a). This condition requires 𝑑(∆ω 21 -∆ω 10 ) 𝑑𝑘 𝑧 > 0, which can be easily satisfied by adjusting the structure-profile. For example, in the in Fig.5(b) we plot the 𝛿𝜔 for the optimized waveguide shown in Fig.4. rightward the working point, we see that ∆ω 21 > ∆ω 10 with increasing 𝑘 𝑧 (grey line). The negative-to-positive trend well validates the possibility of compensating the nonlinear mismatch. For silicon waveguide working at telecom wavelengths,

Fig. 5 .

 5 Fig. 5. (a) Schematic of the proposed SWG waveguide for the simultaneous energy conservation

Fig. 6 .

 6 Fig. 6. (a), (b), (c) and (d) are the dispersion curves of first 3 modes in SWG waveguide with

From 2 ⁄ 2 - 1

 221 in which 𝑛 𝑤 and 𝑛 𝑐 is the material index of waveguide core and cladding while ℎ = √(𝑘 0 2 𝑛 𝑤 2 -𝑘 𝑧 2 ) and 𝛾 = √(𝑘 𝑧 2 -𝑘 0 2 𝑛 𝑐 2 ) are the wavevector along the 𝑦 axis, inside and outside the waveguide, respectively. Assuming that the 𝑘 𝑧 is a constant and 𝜔 is a function of 𝑛 𝑒𝑓𝑓 and 𝑚, the frequency spacing between modes can therefore be expressed as: Equ. (S1) we can infer that if 𝑑𝑚 𝑑𝑛 𝑒𝑓𝑓 ⁄ is proportional to -1 𝑛 𝑒𝑓𝑓 , then the frequency spacing is fixed. Though the mode order number 𝑚 is a discrete integer indicating the phase solution, mathematically it still can be represented as a function of 𝜔 and 𝑛 𝑒𝑓𝑓 , from Equ. describing the phase from standing waves, while 𝑓 𝑚2

2 . 2 .

 22 shown in Fig.S1, modes with different mode orders m are confined in different spatial spans, which sizes increase with m. Each spatial span 𝐿 𝑚 (for m=1,2,3, …) which can be called "effective length" for mode 'm' is defined through 𝑛(𝐿 𝑚 ) = 𝑛 𝑒𝑓𝑓𝑚 . The condition 𝒏 𝒆𝒇𝒇 > 𝒏 𝒃 is defined as "self-adaptive boundary (SAB)" in the manuscript. By solving the equation (S6)-(S12), the eigen equation can be written as: Eigen equation for modes propagating along 𝑧 axis with an electric field polarized along the 𝑦 axis can be rewritten as: the step-index waveguide, to respond to the approximately linear phase increas on the right part of Equ. (1), the changes of 𝜔 and 𝑛 𝑒𝑓𝑓 were correlated thorugh the relationship When the proposed self-adaptive boundary 𝑛 𝑒𝑓𝑓,𝑚 > 𝑛 𝑏 ∀𝑚 is introduced, the spatial integral range is automatically selected which exactly gives the room to trim the frequency spacing through the following condition : To consider the improvement from this new condition, a similar analysis for frequency spacings is carried out as previously with a general index profile described by 𝑛(𝑦) = (𝐴 + 𝐵𝑦) 𝑝 . First considering a linear index profile, i.e 𝑛(𝑦) = 𝐴𝑦 + 𝐵 the phase 𝑁(𝐿 𝑚 ) can be recalculated as: of the self-adaptive behavior that the endpoint 𝑛(𝐿 𝑚 ) is equal to 𝑛 𝑒𝑓𝑓𝑚 , then the (S15) integral can be greatly simplified in a common form for different modes, with 𝑛(0) = 𝐵 = 𝑛 𝑐𝑒𝑛𝑡 , to: 𝑁(𝐿 𝑚 ) = -𝑘 𝑧 𝑛 𝑒𝑓𝑓,𝑚 2𝐴 𝑙𝑜𝑔[𝑘 𝑧 ] -(

Fig. S2 .

 S2 Fig. S2. Schematic of a graded-profile subwavelength waveguide. The equivalent cut-off

𝑛 2 2 - 1 +

 221 𝑛(𝐿 𝑚 ) = 𝑛 𝑒𝑓𝑓,𝑚 , if order = 1, that is to say 𝜂(𝑦) = (𝐴𝑦 + 𝐵), then 𝑁(𝐿 𝑚 ) = 𝑏𝑘 𝑧 √ 1st derivative of 𝑁(𝐿) is modified as: 𝑛 𝑒𝑓𝑓,𝑚 𝑙𝑜𝑔 [

  𝑎)], which we found, is not much different in the index (𝑛 𝑒𝑓𝑓,𝑚 ) range we adresse (i.e. 𝑛 𝑏 >3). Based on Equ. (S34), frequency spacings in different configurations are investigated and presented in Fig.S3. Compared to a step-index waveguide with 𝑛 𝑏 =3 and 𝑎 =350nm, the frequency spacings of linear-shape graded-index waveguide presents a faster change but in a basically similar line-shape. Even though using a wider graded-index waveguide with 𝑎=500nm or smaller boundary index 𝑛 𝑏 =2.6, the lineshapes in the concerned range of effective index (𝑛 𝑒𝑓𝑓 < 𝑛 𝑏 , i.e.[1.6, 3] and [1.6, 2.6], respectively) still behave monotonously.

Fig. S4 .

 S4 Fig. S4. Analytical Frequency spacings as a function of effective index 𝑛 𝑒𝑓𝑓 , with different

  Confirmed by the numerical methods, for 𝑎=350nm, 𝑛 𝑤 =3.48, 𝑛 𝑐 =1.445 and 𝑛 𝑒𝑓𝑓 from 1.6 to 3,

	we have	𝑑𝑓 𝑚2 𝑑𝑛 𝑒𝑓𝑓	∈ (-1.7, -0.19) while	𝑑𝑓 𝑚1 𝑑𝑛 𝑒𝑓𝑓	∈ (-5, -3), from which we may assume that the
	𝑑𝑚 𝑑𝑛 𝑒𝑓𝑓	is almost governed by	𝑑𝑓 𝑚1 𝑑𝑛 𝑒𝑓𝑓	for simplifying the analysis. Therefore,	𝑑𝑚 𝑑𝑛 𝑒𝑓𝑓	can be
	approximately expressed as:
	𝑑𝜔 𝑑𝑚	=	𝑘 𝑧 𝑐 -𝑛 𝑒𝑓𝑓 2	/(	𝑑𝑓 𝑚1 𝑑𝑛 𝑒𝑓𝑓	+	𝑑𝑓 𝑚2 𝑑𝑛 𝑒𝑓𝑓	) ≈	𝑘 𝑧 𝑐 -𝑛 𝑒𝑓𝑓 2	/(	𝑑𝑓 𝑚1 𝑑𝑛 𝑒𝑓𝑓	) =	𝑐𝜋√𝑛 𝑤 2 -𝑛 𝑒𝑓𝑓 2 2𝑎𝑛 𝑤 2	(S5)

  Comparing Equ. (S34) to Equ. (S5), the only change we can clearly see is the replacing from

											√𝑛 𝑐𝑒𝑛𝑡 2 -𝑛 𝑒𝑓𝑓,𝑚 2 √𝑛 𝑏 2 -𝑛 𝑒𝑓𝑓,𝑚 2	+𝑛 𝑐𝑒𝑛𝑡 +𝑛 𝑏	])
	(S32)
	Similarly,
	𝑑𝜔 𝑑𝑚	=	k z 𝑐 -𝑛 𝑒𝑓𝑓,𝑚 2	/	𝑑 𝑑𝑛 𝑒𝑓𝑓,𝑚	[	2 𝜋	𝑁(𝑎) + 𝑓 𝑚2 ]	(S33)
	=	𝜋k z 𝑐 -2𝑛 𝑒𝑓𝑓,𝑚 2	/ { 𝛽k z 2𝐴	[	𝑛 𝑐𝑒𝑛𝑡 √𝑛 𝑐𝑒𝑛𝑡 2 𝑛 𝑒𝑓𝑓,𝑚 -𝑛 𝑒𝑓𝑓,𝑚 2 2	-	𝑛 𝑏 √𝑛 𝑏 2 -𝑛 𝑒𝑓𝑓,𝑚 2 𝑛 𝑒𝑓𝑓,𝑚 2	+ 𝑙𝑜𝑔 ( √𝑛 𝑐𝑒𝑛𝑡 2 √𝑛 𝑏 2 -𝑛 𝑒𝑓𝑓,𝑚 -𝑛 𝑒𝑓𝑓,𝑚 2 2	+𝑛 𝑐𝑒𝑛𝑡 +𝑛 𝑏	)] +
	-𝑛 𝑒𝑓𝑓,𝑚 -1+𝑛 𝑒𝑓𝑓,𝑚 2 (1/𝑛 𝑐 2 +1/𝑛 𝑏 2 )	•	1 √(𝑛 𝑒𝑓𝑓,𝑚 2	-𝑛 𝑐 2 )	•	1 2 -𝑛 𝑒𝑓𝑓,𝑚 √(𝑛 𝑏 2	) }	(S34)
	𝑑 𝑑𝑛 𝑒𝑓𝑓	𝑓 𝑚1 to 𝑑 𝑑𝑛 𝑒𝑓𝑓

The French ANR agency is acknowledged for the its support through the SITQOM project.

The European Research Council (ERC) is acknowledged for its support through the POPSTAR project. National Natural Science Foundation of China (NSFC) (61774131) is also acknowledged.

Supplementary materials

S1. CALCULATING THE 𝒅𝝎/𝒅𝒎 OF A STEP-INDEX TWO-DIMENSION

WAVEGUIDE

into five zones (zones 0, ∓1, ∓2 in the Fig. S1). The central part (part 0) is within the [-𝐿 𝑚 , 𝐿 𝑚 ] range in which 𝑛(𝑦) is larger than 𝑛 𝑒𝑓𝑓 and can be expressed by a 𝑐𝑜𝑠 function as usually, instead of the physical boundaries (i. e. 𝑦 = ∓𝑎). The other four zones (zones ∓1, ∓2) are described by a decaying form because 𝑛(𝑦) is there smaller than 𝑛 𝑒𝑓𝑓 . To satisfy the Maxwell's equations and the corresponding boundary condition for an electric field polarizing along 𝑦, the field can be expanded as follows we rewrite the wave equation inside the waveguide [START_REF] Cassan | Nonlinear Fiber Optics[END_REF]: Based on this approximation, this case can be simplified as a two-dimension waveguide with an equivalent "material index" of:

We then obtain similarly the frequency spacings as:

Currently, the frequency spacing 𝑑𝜔 𝑑𝑚 is no longer simply dominated by the effective index 𝑛 𝑒𝑓𝑓,𝑚 .

Instead, also being affected by the new introduced item with the boundary-induced phase item 𝑎𝑟𝑐𝑡𝑎𝑛(

), as: If 𝑛(y) = 𝐴𝑦 + 𝐵, that is to say, for a linear-index profile waveguide, then