S Le Roux

F Rézaï-Aria

Topological and metric properties of microscopic crack patterns: application to thermal fatigue of high temperature tool steels

Topological and metric properties of microscopic crack patterns: application to thermal fatigue of high temperature tool steels

Sabine Le Roux, Farhad Rezai-

Introduction

When a brittle material (like an oxide film on a metallic substrate) is subjected to transient thermal gradients under thermal fatigue (TF), in-plane tensile stresses are generated during the cooling phase, due to thermal expansion mismatch between the layer and the substrate. According to their amplitude, these tensile stresses may be relieved by the formation of a multi-directional interconnected network of cracks on the surface [1]. Similar 'cooling-induced' or 'shrinkage' cracks patterns are produced in ceramics under thermal shock experiments [2,3], in the glaze of ceramics [4], in solidified basaltic lava flows [5], or in thin layers of desiccating mud or gel [4,6,7]. These networks partition the free surface of the material in a two-dimensional random tessellation of cells, comparable to that formed by grains in planar sections of polycrystalline materials, or by bubbles in soap froths [8].

The topological properties of natural cellular patterns have been studied for a long time, in many scientific fields (e.g. biology, metallurgy, geology, geography) and at widely different size scales (from m 2 to µm 2 ) [8]. Most often, the universal Lantuejoul's procedure [9], developed for polycrystalline ceramics and based on mathematical morphology, is used to determine the distribution of the number of first neighbours (or sides) of the cells. This method is applicable to all kind of networks, whatever the shape of polygons. However, the computation time can be very long for analysing a structure containing a large number of cells, due to the iterative nature of processing. Bandeira et al [10,11] proposed a multi-layer approach to increase the processing speed of the algorithm, since the spatial resolution of their remotely acquired images was very high. This improved algorithm reduces the computation time significantly, but its application is restricted to tri-and tetravalent networks. Bohn et al [12] investigated the dynamic evolution of the topological structure of 'hierarchical' crack patterns, using the 'genealogical tree' of cell division. However, such an approach requires that the cracking phenomenon is observed in situ.

Some semi-empirical statistical rules are commonly used to relate topological and geometrical properties of the cells. The Aboav-Weaire's law [13,14] states a linear relationship between the number of neighbours (or sides) of a cell, n, and the average number of neighbours (or sides) of its neighbouring cells, m(n):

m(n) = ⟨n⟩ -a + (⟨n⟩a + µ 2 )/n, (1)
which can otherwise be expressed as

nm(n) = (⟨n⟩ -a)n + µ 2 + ⟨n⟩a, (2)
where ⟨n⟩ is the average number of sides (or neighbours) of the cells (⟨n⟩ = 6 in infinite cellular structures), a is a structural constant and µ 2 is the second moment (or variance) of the probability distribution P (n) of cell sides (µ 2 = P (n)(n-⟨n⟩) 2 ). This structural equation of state (called 'nearestneighbours correlation') is a quantitative expression of the empirical fact that in planar space-filling cellular structures, many-sided cells tend to be surrounded by few-sided cells, and vice versa. It describes a statistical equilibrium, assuming that there is no correlation in cell shape beyond the nearest neighbours [8]. Since then, some authors [15][16][17] have shown that Aboav-Weaire's law can be generalized to more distant neighbours.

Studying the cellular structure of vegetal tissues, Lewis [18] found that the average area of a n-sided cell, ⟨A(n)⟩, was linearly related to n as

⟨A(n)⟩ = α(n -n 0 ), (3)
where both α (slope) and n 0 (intercept) depend on the properties of the cellular structure. This law reflects the tendency for many-sided cells to be larger, on average, than fewer-sided cells. Using the maximum-entropy principle, Rivier [8,20,21] showed that agreement to this equation of state means that the shaping of the network is governed only by mathematical constraints imposed by the need to completely fill the space without gaps. Furthermore, he defined λ (Lagrange multiplier related to the time) as an important structural parameter with a physical meaning, since it measures the topological ageing in a coarsening structure (like a soap froth or a polycrystalline aggregate). He demonstrated that λ is related both to the slope (α = λ⟨A⟩) and the intercept (n 0 = 6 -1/λ) of Lewis' law [20][21][22]. Equation (3) thus becomes:

⟨A(n)⟩ = ⟨A⟩ [1 + λ(n -6)] , (4)
where ⟨A⟩ is the average area of all cells in the network. The above-mentioned author has also suggested in [21] that the average area of a cell is linearly dependent on the time, in agreement to Von Neumann's law [START_REF] Neumann | Discussion-shape of metal grains Metal Interfaces[END_REF].

Differentiating Lewis' law (3) with respect to time, the rate of change of the average area of the cells, d⟨A(n)⟩/dt, can be expressed as a function of their local topology by the following relationship:

d⟨A(n)⟩/dt = γ (n -6), ( 5
)
where γ is the diffusion constant of the system (γ > 0). Concretely, this equation implies that, in average, the cells with less than six sides grow, those which have more than six sides shrink, whereas the area of 6-sided cells does not evolve over time. This 'topological growth' correlation, established by Von Neumann for soap bubble froths [START_REF] Neumann | Discussion-shape of metal grains Metal Interfaces[END_REF], is in fact applicable to any evolving cellular mosaic, such as biological tissues, soap froths and metallurgical grain structures [21].

Both Aboav-Weaire's and Lewis' laws have been empirically validated for a wide variety of random cellular systems [8,[START_REF] Chiu | [END_REF], as biological and vegetable tissues [16,18,22,[25][26][27], soap froths [13,28], polygonal convective structures [29], colloidal aggregations [30] and polygonal networks on the surface of stars [31,32]. However, experiments reported in soap foams [26,33] and polycrystalline materials [21,[34][35][36] have shown that Lewis' law is not always obeyed. For these structures, it is the cell perimeter (or radius), rather than the area, which depends linearly on n (as expressed by Desch's law [21,37]). The failure of Lewis' law has been attributed to particular constraints acting on the cellular network, other than strictly the mathematical and universal constraints of space-filling. These additional specific forces, which may be physical, biological or chemical, are linked to the energy carried by the cell interfaces (surface tension in soap froths, or grain boundary in metallurgical aggregates [22]), and associated with the mechanism driving the evolution towards statistical equilibrium [21,30].

In this paper, the topological properties of microscopic crack patterns, produced by TF laboratory experiments on the oxidized surface of a X38CrMoV5 tool steel, are investigated using image analysis and Lantuejoul's algorithm [1]. Based on many experimental data, topological and metric correlations (such as Aboav-Weaire's and Lewis' laws) are verified and compared with other natural cellular structures. The effects of the number of cycles and of the conditions of the TF test are shown and discussed.

Experimental procedure

TF tests

TF experiments were carried out on X38CrMoV5 tool steel, heat-treated to achieve either 42 or 47 HRC hardness. The tubular cylindrical-shaped specimens present an external diameter of 30 mm with a polished surface, and a wall thickness of 10 mm in the central gauge area [38]. During the test, the external surface of the specimen was alternatively fast-heated using a high-frequency induction system, and air-cooled by natural convection. The inner wall was continuously cooled by a flow of water circulating through the central hole of the specimen, thus generating a high thermal gradient.

Several thermal cycles were investigated, with a minimal temperature (T min ) always fixed to 100 • C, while the maximum

( • C) ht (s) ( • C s -1 ) max (MW m -2 ) ct (s) ( • C s -1
) of stops of cycles) temperature (T max ) was varied from 550 to 685 • C. For the tests performed between 100 and 650 • C, various heating periods (1.2 < ht < 6.5 s) were applied, leading to maximum heatflux densities ( max ) imposed on the specimen surface ranging from 2.13 to 4.95 MW m -2 (table 1). The experiments were regularly stopped in order to investigate the crack pattern evolution, and determine the number of cycles to crack initiation and to achieve a stabilized crack density (see details in [38]). In table 1, for each TF test, the first interruption of cycling (corresponding to the number of cycles when the micro-cracks were detected first), and the last interruption are given (end of the test, when the crack density was fully stabilized).

Crack pattern analysis

After each cycling interruption, the surface of the TF specimen was observed with a scanning electron microscope and a back-scattered electron detector. By selecting a magnification of 250, an interconnected microscopic crack pattern (named 'heat-checking' [1]) was revealed on the superficial oxide layer covering the steel. A matrix of (3 × 3) images was acquired on a randomly chosen area of 0.825 mm 2 in the centre of the specimen. The nine SEM images were assembled to form a cartography sizing about (1740 × 1210) pixels 2 , with a spatial resolution of 0.5 µm/pixel (figure 1(a)).

Using an image analysis algorithm [1] and the Aphelion ® software, the micro-cracks were segmented and skeletonized (i.e. thinned to one pixel). An example of the resulting two-dimensional cellular network is illustrated in figure 1(b).

To analyse the topological features of the cellular structure, the cells touching the field border were eliminated, since their neighbourhood was incomplete. Within the 'measuring mask' consisting of the remaining cells, between 320 and 4750 cells (N c ) were counted in the 70 heat-checking maps investigated for this study (table 2). Despite their high number, an exhaustive analysis of all cells was systematically carried out, leading to large computation time (typically several hours for most of cartographies). On the other hand, a statistical reliable description of the cellular networks was achieved. A previous study [1] has revealed that the crack-junctions were predominantly tri-, tetra-and pentavalent, excluding the possibility to apply the algorithm developed by Bandeira et al [11]. Therefore, the analysis was carried out following the procedure of Lantuéjoul [9], considering that two cells are neighbours if they share at least one side. This method consists of five steps, as illustrated in figure 2: individual extraction of each cell in the measuring mask (a) ; dilatation of size 2 of the selected cell using a square structuring element, which provides markers both in the selected cell and its first neighbours (b); geodesic reconstruction of the cell and its neighbourhood, using image b as a 'marker' in the initial cellular mosaic (c); logical subtraction of the central cell (d); labelling and counting of the first neighbours of the central cell (e). These operations were iteratively applied to each cell of the network, leading to a labelled image wherein the cells are classified with a grey level related to their topological class n (figure 2(f )). The counting of the cells of each topological class n allowed determining the discrete probability function P (n), as shown in figure 2(g). To compute the average number of neighbours of second order m(n) (i.e. the average number of neighbours of the cells that are adjacent to nneighboured cells), the previous procedure was repeated for each adjacent cell of the selected n-sided cells. According to Lantuejoul's algorithm, in the particular case of two different cracks crossing at right angle, thus determining four adjacent cells, all these cells are considered as neighbours since they share a common point on their border (i.e. the intersection of the two cracks).

It should be emphasized that detailed investigations [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF] have revealed that under the test conditions examined in this study, spalliation of the oxide layer did not occur in the central useful area of the specimen (although some spalled areas may appear at geometrical singularities such as corners or free edges, due to growth effects leading to high strains of the oxide on highly curved surfaces). However, for maximum temperatures (T max ) of the thermal cycle lower than 550 • C, oxide scale spalling becomes a non-negligible damage mechanism. The oxide layer then breaks preferentially by a buckling mechanism, and the heat-checking network is absent or very limited [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. In such cases, the analysis procedure proposed here would not be appropriate.

Topological properties of the heat-checking networks

Distribution of number of cell sides

Table 2 gives some topological descriptors, such as the density of cells in the heat-checking network (d c , i.e. the number of cells N c per mm 2 ) and the statistical distribution of the number of cell sides n (minimum n min , maximum n max , average ⟨n⟩, variance µ 2 , fraction of 6-sided cells P (6) and cumulative fraction of 5-6-7-sided cells 7 n=5 P (n)). These results show that the investigated cellular patterns are composed of polygons whose number of sides lies in the range from 3 to between 10 and 14. The average number ⟨n⟩ varies between 5.91 and 6.03, depending on the TF test. This is very close to the ergodic theoretical value of 6 expressed by Euler's theorem [8,9,14], indicating that the polygonal cells are predominantly 6-sided (as in trivalent networks). In all heat-checking networks, the cells having 5, 6 or 7 sides represent 50-80% of the population. The contribution of 4-sided cells and 8-sided cells is minor, whereas cells with 3 sides and 9 sides or more are rarely observed. Rather comparable distributions of cell sides, mainly composed of a pentagonal-hexagonal-heptagonal mix, have also been reported for other cooling-induced crack networks (in thermal shocked ceramics [2] and basaltic lava flows [5]) as well as for many cellular systems (such as soap froths [8,13], polycrystalline materials [9,14,34], Bénard-Marangoni convective structures [29] or polygonal terrains on Mars [32]). On the other hand, shrinkage-crack patterns generated by desiccation (such as in thin layers of dried slurries [6,7] or in a ceramic glaze [4]) form an arrays of predominantly 4-sided polygons. The variance µ 2 , which measures the spread in the probability distribution of cell sides P (n), is a key indicator of the topological disorder which characterizes the deviation of the cellular structure with respect to the perfect hexagonal pattern [29]. The higher the µ 2 , the more the cells are randomly distributed in the network and the probability of finding 4-sided cell and/or 8-sided cells is high. Delannay and Le Caër [START_REF] Delannay | [END_REF] reported that µ 2 never exceeds 5 in natural cellular structures. The values determined on the heat-checking patterns fall in the interval 1.27-4.30 (table 2), with an average value close to 2.2. A rather comparable degree of disorder was measured in thermal shocked ceramics (where µ 2 = 2.75 ± 0.07 [2]), Al metal grains (µ 2 = 2.90 ± 0.8 [34]) and polygonal terrains on Mars (1.24 < µ 2 < 3.46 [32]), despite the very different size scales of these cellular networks (from µm 2 to m 2 ). However, lower values (µ 2 1.5) were reported, for example in biological cells, cork patterns and Bénard-Marangoni convective structures (see [41] and references therein).

Aboav-Weaire's law

As indicated in table 3, all heat-checking networks satisfy the Aboav-Weaire's law remarkably well. The slope of the regression line is close to 5, and the coefficient of determination R AW is better than 0.99 in 97% of cases (the minimum R AW is 0.97). The system-constant a (calculated from the slope (⟨n⟩a) of equation ( 2)) ranges from 0.44 to 1.47. These values are in agreement with those reported in [41] for most cellular networks (0.5 < a < 1.5). Different interpretations of this parameter are proposed in the literature. According to Zsoldos and Szasz [41], a is a numerical expression of the collectivity of the cellular pattern: the more the increase in a, the stronger the collectivity (with a maximum value of 2). A strong collectivity corresponds to the state of minimum energy for the system. It reflects the presence in the hexagonal structure of dipoles (or chains) of pentagonal and heptagonal cells, which can be regarded as 'topological dislocations' [8] (like the crystallographic defects or irregularities within the atomic structure of a crystalline material). These local 'defects' result from topological transformations (for example successive cell divisions) that change the number of sides of the cells involved. By considering the structure as an elastic solid, the topological dislocations constitute a local manifestation of plasticity [22]. In addition, Cerisier et al [29] have stated that a should be 1 when the cellular mosaic is evolutive, meaning that the cells can coalesce, divide or disappear (such as in polycrystallization processes [13], soap froths [28,42] or biological tissues [25]). In our experiments, the mean value of a is slightly less than 1. We can therefore consider that our heatchecking networks are non-equilibrium structures in which the collectivity is not very strong (the results presented next in section 4.2 will confirm that the cell organization evolves effectively upon cycling).

Lewis' law

The linear Lewis' law (3), relating the average area of the cells ⟨A(n)⟩ to their number of sides n, is broadly verified by all heat-checking networks (table 4). Determination coefficients higher than 0.85 are found, except in four cases (0.635 < R L < 0.998). The slope α of the linear regression shows a wide variation between 32 and 559, indicating that the mean size of the heat-checking cells may vary considerably from one test to another and during cycling. Parallelly, the intercept n 0 ranges from about 0.5 to 4. In most cellular systems, n 0 generally stand between 2 and 2.5 [START_REF] Chiu | [END_REF]29], which is close to the average value found for the 70 heat-checking maps investigated here (n 0 = 2.06 +2.05 -1.58 ). The parameter λ (calculated from (4), as the slope of the normalized curve ⟨A(n)⟩/A versus n) lies in the interval from 0.16 to 0.62, with a mean value of 0.27. For comparison, values of the same order of magnitude were reported for Bénard-Marangoni convective structures (where λ = 0.23 in the steady regime [29]), and polygonal networks on Mars (0.26 < λ < 0.51 [32]).

Discussion

Relations between topological quantities

For the stabilized crack networks, the average number of cell sides ⟨n⟩ seems to be strongly linked to the cell density d c (figure 3), according to a power law with a high coefficient of determination (R = 0.936):

⟨n⟩ = 5.76d 0.0052 c . (6)
This relationship can be attributed to the fact that in finite cellular structure, ⟨n⟩ increases as the number of cells increases (as stated by Euler's theorem: ⟨n⟩ 6 -12/N c [27,30]). Such observation has also been reported in other cellular structures (two-dimensional soap froths [42], inflorescences of plants [27] and polygonal networks on Mars [32]). In our heat-checking networks, N c varies from 323 to 4749, meaning that ⟨n⟩ should theorically vary between 5.96 and 6.00. In fact, the range of variation found experimentally is larger than expected (5.91 ⟨n⟩ 6.03), and values greater than 6.0 are even achieved, as also noted by Saraiva et al. Presumably, this could be explained by the large number of cells (much higher than in other studies) analysed in both cases.

Both the fraction of 6-sided cells P (6) and the topological variance µ 2 reflect the inhomogeneity of the cellular structure, the first inversely and the second directly since it measures the topological disorder [30]. Figure 4, plotting the variation of the topological variance µ 2 as a function of P (6), confirms that these two thermodynamic variables are linked, despite the diversity of their origin. As can be seen, our data and those collected in the literature for others 2D cellular structures [2,18,29,31,32,34,35] lie on a single curve. This nonlinear relationship between µ 2 and P (6) was first discovered empirically by Lemaître et al in two-dimensional foams [43,44]. Lemaître's law has been expressed by Le Caër and Delannay [45], and later by Rivier (based on maximum entropy analysis [START_REF] Rivier | Maximum entropy for random cellular structures From Statistical Mechanics to Statistical Inference and Back[END_REF][START_REF] Rivier | Statistical thermodynamics of foam Foams and Emulsions[END_REF]), as µ 2 •P (6) 2 = 1/(2π) = 0.159 (for 0.3 < P (6) < 0.7), (7.1) and µ 2 = 1 -P (6) (for P (6) > 0.7).

(7.2)

This structural equation of state, which is the equivalent of the virial equation of state in liquids and gases, is very robust for both experimental and simulated cellular systems [22,30,45]. The second expression (7.2) reflects an asymptotic behaviour for P (6) > 0.7. It corresponds to the case of a 'topological gas' (TG), which is an ideal model of uncorrelated arrangement of cells introduced by Fradkov [START_REF] Fradkov | [END_REF]. Our experimental data obey (7.1) with good accuracy, confirming the universality of Lemaître's law for cellular mosaics as a consequence of maximum entropy inference [START_REF] Rivier | Maximum entropy for random cellular structures From Statistical Mechanics to Statistical Inference and Back[END_REF]. Figure 4 further shows the reliability of this equation in the interval 0.18 < P(6) < 0.3 (hitherto not reported in the literature), indicating that Lemaître's law is also valid for highly disordered cellular networks with µ 2 3.

It seems that the parameters of Aboav-Weaire's law (a and µ 2 ) are interdependent [START_REF] Delannay | [END_REF]49]. Figure 5 shows that there is a tendency for the decrease of the system-constant a when the variance µ 2 becomes larger:

a = -0.168µ 2 + 1.3482. ( 8
)
Even if the scattering is high (R = 0.335), it can be noted that the slope of the linear regression curve calculated from our data (s = -0.168) is very close to that reported by Vincze et al (s ′ = -1/6, with a coefficient of determination even lower, R ′ = 0.21 [49]).

In figure 6 the variation of ρ (quotient a/µ 2 ) as a function of µ 2 (with µ 2 ranging from 0.14 to 4.3) is plotted, for various natural structures including our crack networks. As already established by Delannay and Le Caër [START_REF] Delannay | [END_REF], this figure reveals that ρ decreases quasi-universally according to a curve near 1/µ 2 , and approaches an asymptotic value equivalent to that of a TG (ρ TG = -1/6) for cellular structures with a high disorder µ 2 .

Effect of the number of TF cycles

For most TF tests the topological properties of the heatchecking patterns were investigated in the early stages of cracking (for the first 500 or 1000 cycles, as indicated in table 1). Figure 7 shows that the frequency distributions of cell sides, P (n), are quasi-normal and centred on 6. However, it should be emphasized that for a low number of cycles (i.e. when the heat-checking density is not fully stabilized or saturated [1]), P (n) is slightly asymmetric and flattened, with a smaller proportion of 6-sided cells (figure 7(a)) compared with larger numbers of cycles (figures 7(b) and (c)). Whatever the number of cycles N , P (6) is always higher than P (5), and P (5) higher than P (7), as also reported for crack patterns formed in quenched ceramics [2] and basaltic lava flows [5], or for polygonal networks on Mars [32]. In addition, one can observe that the number of cells with 9 or more sides decreases or disappears with increasing N (figures 7(a)-(c)).

With rare exceptions, the variance µ 2 is less than 3 (table 2). According to Stavans and Glazier [8,28,42], such value means that the spatial disorder is 'long-lived in time' and characteristic of a non-asymptotic behaviour. This statement is confirmed by figure 8, plotting the evolution of µ 2 versus number of cycles N , which shows a decreasing trend and tends to a minimal value (µ 2 < 2) when N increases, regardless of T max . This figure suggests that the heat-checking pattern evolves towards a more ordered geometry in the stabilized regime, without however achieving a markedly asymptotic behaviour. In soap froths, a somewhat different evolution was reported for µ 2 versus time [42]: the initially ordered cellular structure, composed mostly of hexagonal bubbles (with µ 2 < 0.5), widens to a maximum (µ 2 = 2.65) before narrowing to a steady intermediate value (µ 2 = 1.5 ± 0.2) at longer times. In their experiments, Stavans and Glazier observed the evolution of the soap bubbles in situ, while our TF tests must be periodically interrupted to check the crack initiation. Therefore, we are not sure that the very early stages of fragmentation were observed, and a behaviour similar to soap froths cannot be completely excluded.

Both Aboav-Weaire's and Lewis' laws are satisfied throughout the whole TF cycling, as shown for test 5 in figures 9 and 10(a), respectively. It should be emphasized that the coefficient of determination of both laws slightly increases with N. A higher deviation from the linear regression is, however, observed for many-sided cells (with n > 9), that can be explained by a lower number of data for this population. This is particularly true for Lewis' law at low number of cycles, since the heat-checking cells have not yet completely reached their stabilized size [1]. Figure 10(b), showing the evolution of the slope of Lewis' law (α) as a function of the number of cycles N , indicates that α decreases to an asymptotic minimal value. It was previously shown [1,38] that the heat-checking density (which is inversely proportional to the mean cell size expressed here by α) follows a sigmoidal trend versus N, whatever the thermal cycling conditions. After an incubation period during which no crack is observed (since the oxide layer is too thin and has not achieved a critical thickness to rupture), the crack network initiates and stabilizes rapidly to a saturated density. Figure 10(b), which shows a bi-linear dependence of α on N , reflects the same phenomenon. Thus, α decreases strongly at the beginning of the cracking process (as shown by the first highly sloped linear portion of the curve), indicating that the heat-checking cells undergo rapidly a strong shrinking after successive cell divisions. Then, beyond a critical number of cycles (N sat ) which may change depending on the test conditions [38], α remains quasi-constant when N increases (corresponding to the second linear part with a nearzero slope), meaning that the cells have reached their stabilized (or saturated) size. Thus, the heat-checking networks are non-equilibrium systems, whose topological and metric properties evolve with time. In figure 11(a), the rate of change of the average area of the cells, d⟨A(n)⟩/dN , is plotted versus their topological class n, for different intervals of cycling (test 5). By excluding the cells with n > 9 (which are very few in the heat-checking network), a linear evolution is observed as a function of n, for each interval between two successive cycling stops:

d⟨A(n)⟩/dN = k(N ) • n + b, (9)
where k(N ) is called the 'rate constant', which depends on the TF cycling interval and indicates either a growth (k > 0) or a shrinkage (k < 0) of the cells, and b is a constant. This expression can be considered as a 'modified' Von Neumann's law, by substituting the time t by the number of cycles N in (5) (it must be recalled that N is proportional to t for a given thermal cycle). At the beginning of the test (i.e. for low number of cycles, up to 2500 cycles in the case of test 5), a negative rate constant is observed (k < -2.10 -2 µm 2 /cycle), meaning that the heat-checking cells shrink (due to their fragmentation). This behaviour is analogous to that of biological tissues (such as epidermis), where successive divisions occur under a microscopic process called 'mitosis' [22]. It can be also noted that the linear regression curve predicts negative values of d⟨A(n)⟩/dN even for values of n higher than 6, in contradiction with Von Neumann's law (which states a linear dependence of d⟨A(n)⟩/dt with (n -6)). This is also in contradiction with experimental observations reported for bubbles in soap froths [28] and grains in polycrystalline aggregates [34,35], which may shrink or expand in agreement with Von Neumann's law. However, even if the sign of d⟨A(n)⟩/dN does not change according to the topological class of the cells, manysided cells have a higher 'shrinkage rate' than few-sided cells (figure 11(a)). In addition, the cells with n > 10 tend to disappear during the fragmentation process, as previously shown in the frequency distributions P (n) (figure 7(d)). Figure 11(b) indicates that the rate constant k, which expresses the rate constant of the cell size variation (slope of equation ( 8)), varies according to the number of cycles N . It may be noted that this curve follows the opposite trend to that of α versus N (figure 10(b)). As explained above, the heat-checking cells globally shrink at the beginning of the fragmentation process, whatever their number of neighbours. Figure 11(b) better shows that the rate of shrinkage k decreases rapidly between the beginning of the fragmentation process and the number of cycles to saturation N sat (as k varies from -0.041 to -0.0009 µm 2 /cycle). When the saturated regime is achieved (i.e. before 5000 cycles), k stabilizes to a constant value near zero, indicating that the cell size no longer evolves under further cycles. Thus, this curve reflects the fact that the fragmentation process is time dependent. Even if the stabilization of the crack network is rapidly achieved, this is not an instantaneous phenomenon since it takes a few thousand cycles before the cells reach their saturated size. This time dependence of k (which is equivalent to γ in original Von Neumann's law) was not observed in soap froths, where γ remains constant with time [28].

Effect of TF test conditions

Focusing on the saturated heat-checking networks (when the crack density is stabilized at a high number of cycles), the effects of TF test conditions (maximal temperature T max and heating period ht of the thermal cycle) on the topological properties are investigated. Whatever T max , a Gaussian probability distribution of the number of cell sides P (n) is obtained, with a slightly more flattened shape when T max < 600 • C (figure 12(a)). In addition, figure 12(b) shows that the kurtosis of the distribution is lower when the heating rate hr decreases (i.e. ht increases), leading to a positively skewed distribution for the slowest thermal cycle (ht = 6.5 s). The average number of cell sides does not drastically change with neither T max (⟨n⟩ = 6.01 ± 0.2), nor ht (⟨n⟩ = 5.99 +0.03 -0.05 ). However, it seems that a relationship exists between ⟨n⟩ and the maximum heat-flux density ( max ) applied to the specimen during the heating period of the TF test (figure 13). It must be recalled that max depends on the test conditions and especially on the heating rate, as indicated in table 1 (2.13 < max < 4.95 MW m -2 when 85 < hr < 458 • C s -1 for the 100-650 • C thermal cycle). In figure 13, data points are best-fitted by a power law (with R = 0.82): ⟨n⟩ = 5.915 0.0107 max .

This relationship can be related to the tendency of ⟨n⟩ to increase with the cell density d c , as previously shown in figure 3. Indeed, the thermal cycles performed with a low max produce heat-checking networks with fewer cells than those performed with a high max , and vice versa (table 2). No very significant effect of the steel hardness and the TF test conditions is observed on Aboav-Weaire's correlation, although the system-constant a appears to be lower when the test is performed at low T max (<600 • C) or at high ht (table 3). The topological disorder µ 2 exhibits greater variations. More disordered cellular networks (with µ 2 > 2) are obtained at low T max , while the most ordered organization seems to be reached at 650 • C (figure 14(a)). On the other hand, figure 14(b) indicates that µ 2 increases by increasing ht, meaning that low heating rates tend to produce a greater deviation from geometric equilibrium, compared with fastest cycles.

Figure 15(a) shows the effect of T max on Lewis' linear correlation, whose slope α does not vary much in the saturated regime. However, it can be noted that when T max < 600 • C (tests 1 and 2), the dependence of ⟨A(n)⟩ upon n seems rather quadratic than linear (figure 15(b)), thus improving the average correlation coefficient from 0.91 to 0.97. The failure of Lewis' law under these conditions would suggest that additional constraints to inescapable mathematical spacefilling requirements are involved in shaping the heat-checking network [20][21][22]. In addition, figure 16 indicates that the slope of the linear regression, α, is drastically increased by increasing ht (i.e. decreasing the mean heating rate hr). The scattering from the linear regression is also larger by increasing ht, because the number of cells becomes lower (especially when ht 3.7 s). In fact, figure 17 reveals that α, which is linked to the average cell area, decreases with max following a power law:

α = 2414 -2.362 max . (11)
This curve, which fits our experimental points with a remarkably high coefficient of determination (R = 0.95), clearly shows that α reaches an asymptotic minimum when max is greater than 4 MW m -2 . It means that beyond this value, the saturated heat-checking cells would have achieved their minimum size, and could not be more fragmented. A previous study [38] revealed that the saturated crack density ρ sat is dependent on the input energy (applied at the specimen surface during the heating period of the thermal cycle), which is itself linked to the maximum heat-flux density max . Thus, ρ sat increases by increasing max , and beyond a critical limit of the transmitted energy per unit-area, ρ sat remains quasiunchanged. Therefore, the stabilized size of the cells may vary in a wide range depending on the test conditions (and especially the heating rate of the thermal cycle). The higher the max (and therefore the energy received by the specimen surface, i.e. the oxide layer), the smaller the resulting cell size. With our test facility, the maximum heat-flux density that can be achieved (by adjusting the power of the generator to 100%) does not exceed 5 MW m -2 , leading to heat-checking networks with a minimum cell size of about 15-20 µm [38]. Once the heat-checking network is initiated and stabilized, increasing the number of thermal cycles leads to propagation of the micro-cracks of the oxide layer in the steel substrate. Under such conditions, the thermally induced mechanical energy is dissipated by the crack propagation inwards of the steel, rather than by subdividing the heat-checking cells [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. The interaction between the thermal stresses and oxidation ahead of the crack tip modifies the oxidation diffusion paths. In fact, due to in-depth crack propagation, localized pit oxidation occurs preferentially to the formation of oxide on the surface of the specimen. Therefore, fewer driving forces are available to subdivide the existing cells.

By plotting the normalized area of the cells, ⟨A(n)⟩/⟨A⟩, versus their topological class n (normalized Lewis' law) for all test conditions (figure 18), all data points tend to cluster on a single curve, namely the linear law:

⟨A(n)⟩/⟨A⟩ = 0.267n -0.574. ( 12
)
Despite a higher scattering for the many-sided cells (n > 8), a quite good coefficient of determination is obtained (R = 0.88). This result indicates that irrespective of the TF test conditions, the heat-checking networks are similar apart from a size scaling factor, suggesting that the same rules govern the shaping of the cellular structure in the stabilized state. However, as previously shown, the cell organization could be affected, under some test conditions, by additional forces to the inescapable space-filling constraints. Actually, during our TF experiment, the maximum heat-flux density imposed on the specimen generates thermomechanical loading by self-constraining [38]. This thermo-physical constraint contributes to the shaping of the heatchecking network, under certain conditions (low T max and low hr).

The steel grade investigated in this study is commonly used as a tool material in hot-forming processes, such as diecasting or forging. In service, the tools suffer from thermal and mechanical solicitations, which can be reproduced at a lower scale by our TF test [1,38]. In the run of industrial application, this work shows that by changing the test parameters (in particular the heat-flux density by modifying the heating rate of the thermal cycle), heat-checking networks with various cell sizes can be generated. It is well known that the TF life is considerably reduced when the superficial heat-checking network is denser, because these micro-cracks contribute to weakening the steel by accelerating the oxidation mechanisms and the initiation of macro-cracks that propagate through the tool material as a result of accumulated thermo-mechanical solicitations [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. Therefore, better controlling the heat exchanges between the tool and the processed parts could be an interesting way to improve the life of hot-working tools in service [38].

Conclusions

This paper presents an experimental investigation on the topological properties of microscopic heat-checking patterns formed under thermal fatigue on the oxidized surface of a X38CrMoV5 steel specimen. 300 and 6000 cells were analysed in each polygonal crack network, leading to statistically relevant results. An average number of sides per cell very close to the ideal value expressed by Euler's theorem is found (⟨n⟩ = 6.0 +0.03 -0.09 ), indicating that the cellular networks are pseudo-hexagonal. The applicability of general topological laws (such as Aboav-Weaire's and Lewis') is verified, meaning that geometric correlations do not extend beyond adjacent cells in the polygonal network, and that the shaping of the network is governed by space-filling constraints [8]. The topological properties of the microscopic heat-checking networks are broadly consistent with those reported in the literature for crack patterns in thermal shocked ceramics or for polygonal terrains on Mars, which suggests that these random cellular systems are governed by the same rules despite their very different scale. The universality of Lemaître's law is confirmed, and its domain of validity extended to highly disordered structures.

It is shown that the heat-checking networks are nonequilibrium systems, as their topological and geometrical properties change with the number of cycles. The topological disorder, estimated by the variance of P (n), is rather high (µ 2 = 2.22 +2.08 -0.95 ). It decreases when increasing the number of cycles, and evolves towards a geometric equilibrium with µ 2 < 2. A greater disorder is achieved for lower maximal temperatures (T max < 600 • C) and lower heating rates. It is shown that the heat-checking cells undergo a strong shrinkage at the beginning of the thermal cycling. A modified version of Von Neumann's law is established, relating the rate of change of the average area of the cells to their topological class.

Whatever the number of cycles and the TF test conditions, the linear Aboav-Weaire's law is satisfied with a correlation coefficient always greater than 0.97. The system-constant is near 1 on average (a = 0.96 ± 0.52). The Lewis' law, which states a linear relationship between the average area of the cells and their number of sides, is generally obeyed with a slope (32 < α < 559) mainly depending on the thermal cycling conditions (and especially the heating rate). However, better correlation coefficients are achieved when the heat-checking pattern is stabilized (for a high number of cycles), and at higher maximum temperatures (T max > 600 • C). Below this T max , the best-fit correlation is quadratic rather than linear, suggesting that physical forces (related to the maximum heatflux density applied to the specimen surface) are involved in the fragmentation process. Neither the intercept (n 0 = 2.06 +2.05 -1.58 ) nor the 'ageing' parameter (λ = 0.27 +0. 35 -0.11 ) of Lewis' law appears related to the test conditions. It is shown that for stabilized cracks networks, the maximum heat-flux density ( max ) influences both the average number of sides (⟨n⟩) and the average size of the heat-checking cells (i.e. the slope α of Lewis' law). Whatever the test conditions, a unique normalized Lewis' law is found, indicating that the stabilized cellular networks are indistinguishable apart from a scale factor (the average cell size).

Figure 1 .

 1 Figure 1. (a) Typical SEM cartography of the microscopic crack pattern observed on the surface of the specimen; (b) enlarged binary image of the skeletonised crack network (cracks in black and polygonal cells in white).

  6.01-6.02 1.75-2.22 0.19-0.22 0.56-0.60 a Test performed on a 42 HRC heat-treated specimen.

Figure 2 .

 2 Figure 2. Procedure to determine the number of first neighbours (n) of the heat-checking cells [9]: (a) selection of a cell in the measuring mask; (b) dilatation of size 2 of the selected cell; (c) reconstruction of the cell and its neighbourhood; (d) subtraction of the central cell; (e) labelling and counting of its first neighbours; (f ) cellular mosaic in which the cells are classified with a grey level related to their number of neighbours; (g) histogram of frequency distribution P (n) (probability of finding a cell with n neighbours or sides).

Figure 3 .

 3 Figure 3. Relation between the average number of sides ⟨n⟩ of the cells and their density dc in the stabilized heat-checking networks.

Figure 4 .

 4 Figure 4. Universality of Lemaître's law: correlation between the variance µ 2 and P (6) for the heat-checking networks investigated in this work and other cellular structures reported in the literature.

Figure 5 .

 5 Figure 5. Linear relation between a (expressing the collectivity in Aboav-Weaire's law) and the variance µ 2 , for the heat-checking networks.

Figure 6 .

 6 Figure 6. Variation of ρ(= a/µ 2 ) versus µ 2 for the heat-checking networks investigated in this work and other cellular structures reported in the literature.

Figure 7 .

 7 Figure 7. Frequency distributions P (n) of the number of cell sides at increasing number of cycles N (test 5): (a) N = 1000 cycles; (b) N = 5000 cycles; (c) N = 20 000 cycles; (d) evolution of P (n) versus N, for n = 5 to 7. (⟨n⟩ is the average number of sides of the cells, and µ 2 is the variance of P (n)).

Figure 8 .

 8 Figure 8. Evolution of the variance of µ 2 versus number of cycles, indicating a decreasing topological disorder in the heat-checking pattern.

Figure 9 .

 9 Figure 9. Applicability of the Aboav-Weaire's law (linear relation relating the average number of sides of cells surrounding a n-sided cell, m(n), to (n), for test 5. The data are the mean values of n • m(n) calculated between the first (1000 cycles) and last (20 000 cycles) interruption of the test, and the error bars represent their minimum and maximum values.

Figure 10 .Figure 11 .

 1011 Figure 10. Applicability of Lewis' law (test 5): (a) linear relation relating the average area ⟨A(n)⟩ of the cells to n, at different numbers of cycles; (b) evolution of the slope α versus number of cycles.

Figure 12 .

 12 Figure 12. Frequency distributions of the number of cell sidesP (n) in the stabilized regime, for various maximum temperatures T max (a) and heating periods ht of the thermal cycle (b).

Figure 13 .

 13 Figure 13. Effect of the maximum heat-flux density max on the average number of cell sides ⟨n⟩, for the saturated heat-checking networks.

Figure 14 .Figure 15 .

 1415 Figure 14. Effect of the maximal temperature T max (a) and of the heating period ht of the thermal cycle (b) on the variance µ 2 in the stabilized regime (µ 2 is an indicator of the topological disorder in the cellular network).

Figure 16 .

 16 Figure 16. Effect of the heating period ht of the thermal cycle on Lewis' law in the stabilized regime (note that the slope α increases with increasing ht).

Figure 17 .

 17 Figure 17. Effect of the maximum heat-flux density max on the slope of the Lewis' law α, for the saturated heat-checking networks.

Figure 18 .

 18 Figure 18. 'Universal' normalized Lewis' law (⟨A(n)⟩/⟨A⟩ versus n) applied to all the heat-checking networks produced with various test conditions (maximum temperature T max and heating period ht).

Table 1 .

 1 Thermal fatigue test conditions performed on 42 and/or 47 HRC treated X38CrMoV5 steels.

		Thermal	Heating period	Maximum	Cooling period	Cycling interruptions
	TF	cycle	Mean	heat-flux	Mean	First and last
	test T min -T max Duration rate hr	density	Duration rate cr	Number stops (number
	N •

  Tests non-interrupted in the early stages of crack initiation.

	1	100-550	0.90	500	4.75	14.5	31	2	30 000-80 000 a
	2	100-575	0.95	500	4.84	14.8	32	7	5000-50 000
	3 4	100-600 100-625	1.0 1.1	500 477	4.85 4.93	15.0 16.0	33 33	9 5	2500-60 000 10 000-50 000 a
	5	100-650	1.2	458	4.95	17.0	32	9	1000-20 000
	6	100-650	2.0	275	3.93	20.0	28	9	1000-30 000
	7	100-650	3.7	149	2.78	22.0	25	8	2000-30 000
	8	100-650	6.5	85	2.13	30.0	18	7	3000-50 000
	9	100-685	1.6	366	4.75	22.0	27	7	500-3000

a

Table 2 .

 2 Range of variation of topological properties (minimum and maximum values between the first and last cycle interruptions) calculated on the microscopic heat-checking cells produced under various TF test conditions on 47 HRC heat-treated specimens.

	TF	Measuring mask		Distribution of the number of cell neighbours (or sides) n
	test Total number Cell density
	N •	of cells N c	d c (mm -2 )	n min n max	⟨n⟩	µ 2	P (6)	7 n=5 P (n)
	1	2812-3002	3705-4015	3	12

-13 6.01-6.02 2.37-2.63 0.19-0.25 0.57-0.66 2 2495-3439 3311-4523 3 11-13 6.01-6.03 2.11-2.94 0.18-0.26 0.52-0.69 3 2256-3381 4109-4606 3 11-12 6.01-6.02 1.73-2.

Table 3 .

 3 Least-squares fit parameters of Aboav-Weaire's law (minimum and maximum values between the first and last cycle stops) calculated on the heat-checking networks produced under various TF test conditions on 47 HRC heat-treated specimens.

		TF test conditions	Aboav-Weaire's law (2) parameters
	N •	T max ht	max	Slope of linear System	Determination
	test ( • C) (s)	(MW m -2 ) regression	constant a coefficient R AW
	1	550	0.90 4.75	5.1-5.5	0.51-0.96 0.992-0.995
	2	575	0.95 4.84	5.1-5.3	0.68-0.93 0.995-1.000
	3	600	1.0	4.85	4.9-5.1	0.93-1.09 0.991-1.000
	3 a	600	1.0	4.85	4.6-5.2	0.78-1.39 0.993-1.000
	4	625	1.1	4.93	4.9-5.2	0.85-1.11 0.997-1.000
	5 5 a	650 650	1.2 1.2	4.95 4.93	4.7-5.2 4.7-5.2	0.79-1.26 0.994-1.000 0.79-1.25 0.996-1.000
	6	650	2.0	3.93	4.7-5.1	0.85-1.47 0.980-1.000
	7	650	3.7	2.78	4.9-5.1	0.88-1.15 0.992-0.999
	8	650	6.5	2.13	4.5-5.5	0.44-1.44 0.973-0.999
	9 9 a	685 685	1.6 1.6	4.75 4.75	4.8-5.2 4.8-5.2	0.81-1.20 0.995-1.000 0.84-1.16 0.991-1.000

a Test performed on a 42 HRC heat-treated specimen.

Table 4 .

 4 Least-squares fit parameters of Lewis' law (minimum and maximum values between the first and last cycle stops) calculated on the heat-checking networks produced under various TF test conditions on 47 HRC heat-treated specimens.

		TF test conditions		Lewis' law parameters in (3) and (4)
	N •	T max ht	max	Slope α	Intercept n 0 Slope λ	Determination
	test ( • C) (s)	(MW m -2 ) (3)	(3)	(4)	coefficient R L
	1	550	0.90 4.75	69-78	2.56-2.63	0.30-0.32 0.856-0.956
	2	575	0.95 4.84	67-170	2.61-4.11	0.31-0.62 0.790-0.972
	3 3 a	600 600	1.0 1.0	4.85 4.85	52-60 32-73	1.57-2.46 0.48-2.50	0.24-0.34 0.940-0.981 0.16-0.30 0.854-0.994
	4	625	1.1	4.93	39-52	1.38-2.19	0.22-0.28 0.760-0.992
	5 5 a	650 650	1.2 1.2	4.95 4.93	47-96 45-95	0.71-2.78 1.02-2.25	0.20-0.32 0.920-0.997 0.20-0.28 0.911-0.998
	6	650	2.0	3.93	63-102	0.96-1.88	0.20-0.34 0.923-0.990
	7	650	3.7	2.78	196-259 1.34-2.02	0.20-0.25 0.899-0.993
	8	650	6.5	2.13	303-559 1.33-2.43	0.21-0.30 0.635-0.993
	9 9 a	685 685	1.6 1.6	4.75 4.75	39-77 41-62	1.61-2.59 1.50-3.00	0.23-0.32 0.924-0.995 0.23-0.36 0.692-0.994

a Test performed on a 42 HRC heat-treated specimen.

Acknowledgments

The authors wish to thank Professor C Levaillant, former head of the research centre 'ICA-Albi', for the fruitful discussions and encouragement to undertake this investigation. Dr F Medjedoub is thanked for his contribution to thermal fatigue experimentations.


Napoléon Gutierrez

Arnaud Fernandez

Olivier Llopis

Stéphane Calvez

Stéphane Balac

Numerical studies on Kerr comb generation in Si3N4 resonators with frequency dependent access coupler properties

Keywords: Racetrack Resonator, Frequency-dependent coupling, Kerr Comb

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Kerr Frequency Comb (KFC) generation is the process by which an optical spectrum consisting of evenly spaced lines are generated by taking advantage of the Kerr effect. Optical resonators are attractive platforms for KFC generation because of their resonant nature i.e. because of their mode selectivity that and their power enhancement properties. In continuous wave (CW) pumped micro resonators, KFC generation arises in the anomalous dispersion regime when the parametric gain engendered by light circulating inside the cavity overcomes losses and the first modulation sidebands are generated by degenerate Four Wave Mixing (FWM). As more power is accumulated inside the cavity, the generated combs lines interact with each other by nondegenerate FWM forming an evenly space comb with a multiple-FSR spacing [START_REF] Herr | Universal formation dynamics and noise of Kerr-frequency combs in microresonators[END_REF][START_REF] Herr | Temporal solitons in optical microresonators[END_REF] . As these comb lines become strong enough to create parametric sidebands strong enough to overcome total cavity losses, they broaden occupying the resonator modes between them and merge together creating a FSR-spaced frequency comb [START_REF] Herr | Universal formation dynamics and noise of Kerr-frequency combs in microresonators[END_REF] .

The power threshold for KFC generation in optical resonators can be optimized by maintaining adequate coupling properties for optimum power enhancement over a large bandwidth. Furthermore, total roundtrip losses can also be optimized so that the losses that have to be overcome by parametric gain in order for FWM to occur are further reduced. Hence, this is why the work of Arlotti et al. on vertically coupled whispering-gallery-mode (WGM) Si 3 N 4 resonators design for wideband critical coupling, which has demonstrated coupling geometries that allow critical bandwidths larger than 330 nm [START_REF] Arlotti | Achromatic critically coupled racetrack resonators[END_REF] , is so important. These results prove very promising in the context of Kerr frequency comb generation in these platforms as access coupling properties are a crucial factor in a resonator's power enhancement properties and finesse [START_REF] Yariv | Universal relations for coupling of optical power between microresonators and dielectric waveguides[END_REF] . Furthermore, Si 3 N 4 is a material that proves interesting in the context of KFC generation due to its high nonlinear refractive index (n 2 =2.5 10 -15 cm 2 .W -1 ), its compatibility with CMOS technologies and the absence of two-photon absorption in telecommunication band [START_REF] Levy | CMOScompatible multiple-wavelength oscillator for on-chip optical interconnects[END_REF] . However, as optimized as the access coupling can be, these structures can present spectral regions where the resonator's quality (Q) factor (and thus finesse) is strongly degraded. Moreover, the resonator's access coupler dispersion has a strong influence on the overall resonator dispersion. All this creates problems in numerically simulating Kerr comb generation with the Lugiato-Lefever Equation (LLE) [START_REF] Coen | Modeling of octave-spanning Kerr frequency combs using a generalized mean-field Lugiato-Lefever model[END_REF] .

This issue can be circumvented by adopting a model that separates intra-cavity dispersion from dispersion in the access coupler such as the Ikeda map 8 , already used in references [START_REF] Hansson | Frequency comb generation beyond the Lugiato-Lefever equation: multistability and super cavity solitons[END_REF][START_REF] Gutierrez | An iterative method for the dynamic modeling of ultra-short pulse generation in nonlinear optical ring resonator[END_REF] . While this approach is certainly more computationally cost-intensive due to the fact that the step-size has to be inferior to the cavity length, we have put in place a RK4IP method 2 to reduce computation time and keep computation accuracy high. In addition to this, Si 3 N 4 having a higher absorption coefficient than materials such as fluoride glasses (CaF 2 and Mg 3 F 2 11, 12 ), the quality factor of such resonators doesn't usually exceed 10 7 making numerical simulations of these cavities less time-intensive.

In this paper, we investigate on Kerr comb generation for the three resonators modelled by Arlotti et al. in reference [START_REF] Arlotti | Achromatic critically coupled racetrack resonators[END_REF] . We start by presenting the modelling approach and parameters used to analyze the frequencydependent cavity transfer properties and Kerr frequency comb generation simulations in section 2. Section 3 describes the transfer properties of the three optical resonators studied in this paper. Section 4 presents simulation results of Kerr comb generation in resonators with frequency-dependent and complex-valued access coupler coefficients.

NUMERICAL MODEL AND PARAMETERS

Our numerical approach consists in simulating the round-trip to roundtrip evolution of the intra-cavity field of a pumped optical resonator by separating intra-cavity propagation from propagation through the coupler. This is a two time-scale approach that associates a slow time-scale to the variable t=mt R to the m-th roundtrip of light circulating inside the resonator with t R the cavity roundtrip time, and a fast time-scale to the variable τ comprised between -t R /2 and + t R /2 that describes the temporal distribution of the optical signal circulating inside the optical cavity. A schematic of the device modelled is presented in figure 1, along with the absolute value of the coupler's transmission and coupling coefficients at each simulated resonator mode In figure 1, A in represents the pump signal at the input of the optical resonator. In our simulation approach, the pump laser is a continuous wave (CW) signal that is superposed to a white noise floor 150 dB weaker than the pump power. A m represents the intra-cavity signal after propagation through the cavity at the end of roundtrip m. A m-1 represents the intra-cavity field at round-trip m before propagation through the cavity. Coupler coefficients are represented as a function of mode number which normalizes frequency with respect to the pumped mode and the resonator's FSR. Intra-cavity propagation is described the Non Linear Schrödinger Equation (NLSE) which, limiting dispersion up to the second order, is given by

A z, τ 2 A z, τ 2 ∂ A z, τ A z, τ |A z, τ | (1)
In equation ( 1), parameter β 2 models dispersion during propagation, parameter α models field attenuation and γ is the Kerr coefficient that models the nonlinear Kerr effect during propagation. Material and geometrical properties (e.g. bending radius) of the racetrack resonator have a direct impact on the dispersive properties and propagation losses of the ring cavity [START_REF] Moss | New CMOS-compatible platforms based on silicon nitride and Hydex for nonlinear optics[END_REF] . Spectral dependence of and depicted in figure 2 shows respectively an attenuation of 29.69 and an anomalous chromatic dispersion of 18.88 at 1.55 µm (pump laser wavelength). The non-linear index ! combined with mode effective area " #$$ lead to a Kerr coefficient 1200 & . Access coupler coefficients are complex valued and frequency dependent. They are used to model not only linear propagation through a coupler of non-zero length but also dispersive effects. In fact, the magnitude values of ρ 1 , ρ 2 , κ 1 , κ 2 model the amplitude of the transfer between input waveguide and ring cavity while the phase of these coefficients models changes in the optical path of the different spectral components propagating through the coupler, i.e. dispersion during propagation through the coupler. The passage of light through the access coupler is modelled by a transfer matrix that describes signal transfer between the input/output of the resonator and the cavity in the spectral domain.

' " ( ) " ( *+, -λ / 0 1 1 0 2 ' " ( 34 " ( ) -λ (2)
Propagation through the access coupler is performed numerically at the beginning of every iteration in the spectral domain with the help of fast Fourier transforms (FFT), respecting relation (2). Then, intra-cavity propagation is performed using a RK4IP solver of the NLSE through the length of the cavity taking into account fluctuations of the dispersion and absorption parameters.

This approach has a higher computation cost than ones based in coupled modes equation or split-step simulations of the LLE because it forces the user to adopt a step-size that is lower than the cavity physical length. However, by separating the treatment of intra-cavity propagation and propagation through the access coupler, this two time-scale approach allows for the implementation of non-reciprocal access coupler properties. Furthermore, this approach accurately models Kerr Frequency Comb generation in optical resonator for any optical cavity regardless of its finesse, quality factor or dispersion.

MODELLING OF THE CAVITY TRANSFER PROPERTIES

Transfer function

A resonator's transfer properties give useful insight about its resonance linewidth, cavity enhancement properties and resonance mode structure. To determine the transfer properties of an optical resonator, it is necessary to consider the ring resonator behavior in linear regime i.e. with very low optical power circulating in the cavity such that the non-linear term doesn't manifest itself in the NLSE (see eq. 1). A very straightforward solution of the propagation equation under this circumstances allows us to determine the frequency dependent Q factor, free spectral range and power enhancement factor (PEF) of our optical resonator. Indeed, having access to the frequency-dependence of the quantities listed above provides substantial information concerning frequency dependent resonator losses as well as resonance mismatch that subsequently affect KFC generation and reveals the impact on the coupler design. By assuming the absence of Kerr nonlinearity in equation ( 1), the intra-cavity field propagation admits an analytical solution given by equation (3).

" ( ) 5 6 " ( ) 5 6 7 8 9 :; <= > ? @ ? AB

Here, 5 6 corresponds to the frequency at mode in an evenly spaced frequency grid 5 6

. CDE ∈ ℤ . By combining eq.

(3) with the coupler's transfer function (2), on can easily obtain the resonator's transfer function H I in the steady state regime at a given frequency which is given by the relation

H I " ( *+, " ( 34 1 0 0 1 1 J K L M 7 ;N 1 0 J K L M 7 ;N (4)
In equation ( 4), the ring attenuation and dispersion per roundtrip are represented respectively through operators J K 7 O ?

B and L M 7 ; <= > ? P ? ? B . The variable I represents the linear phase accumulated by the intra-cavity field over one round-trip with respect to the FSR spaced frequency grid. This linear phase is linked to a shift in frequency Q ∈ R CDE/2; CDE/2U which represents the frequency detuning in the vicinity of the k th optical mode: 2V I 2V = > :$ WXY .

Quality factor and Power enhancement factor

The resonator's quality factor Q is linked to the photon storage time of the optical resonator. In the spectral domain the Q factor is inversely proportional to the full-width at half maximum ΔQ [ of the cold cavity transfer function squared modulus |H I | at resonance. It is expressed as following: \ 5 ] /ΔQ [ , in which 5 ] represents the absolute frequency at resonance. Experimental spectro-RF metrology technique invoke this Q factor expression for measurement [START_REF] Abdallah | Radio frequency spectral characterization and model parameters extraction of high Q optical resonators[END_REF] . In this numerical work, we compute the local Q factor at each resonator mode by referring to equation ( 4) for each of the three simulated optical resonators. Increasing resonator's Q factor reduces the total cavity losses that have to be overcome by parametric gain in order to start Kerr comb generation and thus reduces the pump power threshold. Another way to decrease the comb generation threshold is to optimize the accumulation of optical power inside an optical cavity at resonance. Pump storage efficiency is quantified by a resonator's power enhancement factor (PEF). It is determined, again, using equations ( 2) and ( 3) to determine the ratio of intra-cavity power with respect to pump power. The analytical description of the cavity PEF is given by equation ( 5)

^_C `" ( ) ` `" ( 34 ` |κ | |1 J K 1 | (5)
For a passive ring cavity (which is the case in our study), total cavity losses `J K 1 `b 1 and critical coupling only occurs for the singular condition: J K 1 corresponding to a maximum achievable PEF with a given propagation medium. Since the ring resonator intensity losses J K and the coupling properties are both frequency dependent, this inevitably degrades the PEF and affects coupling regime. The frequency-dependent quality factor and PEF of the three resonators studied in this letter are depicted in figure 3. We observe a similar Q factor c 1.4 10 e at the pumped mode for all three cavities; however the coupling design affects drastically the Q factor value as mode number increases. It is easy to distinguish the case corresponding to Coupler 2 as it corresponds to a high quality factor over a considerably narrower bandwidth compared to that of the two other cases.

Cold resonance mismatch: Intra-cavity and coupler dispersion

Another critical aspect in Kerr comb generation is total roundtrip dispersion. As stated in section 1, KFC generation in optical resonators arises from the compensation of dispersive effects by nonlinear effects. The fact that the optical path of light is frequency dependent produces changes in the spectral position of each resonance of an optical resonator. If dispersion is adequate, this allows the Kerr effect to pull the resonator modes into an evenly spaced frequency-grid as pump power is accumulated inside the cavity. In order for modulation instability to occur and give rise to an optical frequency comb, the mode structure of the cavity has to present a concave shape at the pump wavelength. This means having a second order dispersion parameter that is negative at the pump wavelength. There are two sources of dispersion in the model presented in this paper and those are intra-cavity dispersion modeled by coefficient in the propagation equation and the phase of the complex-valued coupler coefficients. The phase modulation term 2V5 6 /2 in equation ( 3) models the resonance frequency mismatch due to intra-cavity dispersion of a resonator and can be used to obtain its contribution to the cold cavity resonance mismatch of a resonator with the equation

∆5 g h ! !! 2V5 n !k1 CDE 2V (6)
Since the contribution of the access coupler to the cold cavity resonance mismatch is modelled by the coupler's transfer matrix coefficients, it is necessary to use equation ( 4) to obtain the global cold cavity resonance mismatch. Numerically, we find the resonance position of each resonator mode using equation ( 4) and then compute the difference between resonance position and the evenly spaced frequency grid 5 6

. CDE ∈ ℤ obtaining the cold resonance frequency mismatch of each resonator. This effect is mutual to all the simulated optical resonators and is obtained by using equation ( 6). Figure 4(b) depicts the cold cavity resonance mismatch computed using equation ( 4) to determine the exact spectral position of each resonator mode with respect to an equally, FSR-spaced frequency grid. Figure 4 shows how in the resonators modelled by Arlotti et al., total round-trip dispersion is strongly affected by coupler geometry. This highlights the importance of taking into account dispersion effects during propagation through the access coupler.

Furthermore, the simple mathematical analysis that was put in place in this section allows us to determine which of the three resonators presents better properties for Kerr comb generation. In terms of frequency-dependent losses, a high quality factor over a large bandwidth is desired so that parametric gain engendered by the pump signal circulating inside the cavity has lower losses to overcome in order for degenerate FWM to occur and give rise to an optical frequency comb. In this category, resonators 1 and 3 are clearly better options than resonator 2. With respect to the cavity PEF, it is desirable to have the strongest PEF possible at the pumped mode to best take advantage of the power enhancement properties of the cavity and reduce the pump power threshold necessary for Kerr comb generation even further. Concerning resonator dispersion, as long as dispersion is anomalous at the pumped mode, a flatter resonance mismatch curve lends itself to the formation of broader frequency combs [START_REF] Coen | Universal scaling laws of Kerr frequency combs[END_REF] . Figure 4 b) shows that resonator 2 possesses the best adapted mismatch curve even suggesting the possibility of dispersive-wave formation around modes -140 and 140; [START_REF] Bao | High-order dispersion in Kerr comb oscillators[END_REF] .

SIMULATION RESULTS

Kerr comb simulations rely on the model presented in section 2. Because of Kerr nonlinearity, the phase of light propagating in the optical medium is modulated by its own intensity. The term |" ) | " ) in the nonlinear Schrödinger equation accounts for the Kerr non-linear phase shift accumulated by light during propagation. This phase shift is due to the modification of the medium's refractive index by the Kerr effect effectively lengthening the optical path of light circulating in the optical resonator as power is increased.

In the context of nonlinear optical cavities, as optical power is accumulated inside the resonator, the pumped mode resonance shifts towards longer wavelengths (lower frequencies) due to self-phase modulation. Because of this, controlling the pump laser frequency detuning with respect to the pumped mode resonance is crucial in practice and in Kerr comb simulations. In order to efficiently accumulate power inside a non-linear optical cavity, the frequency of the pump laser has to be controlled in such a manner that it follows the Kerr-induced resonance shift that takes place as optical power is accumulated inside the optical cavity. This means that the pump laser has to be tuned from higher to lower frequencies across the pumped mode resonance in order to follow the resonance drift caused by the Kerr effect as light is accumulated in the cavity. In simulations, the phase detuning of the pump laser with respect to the cold cavity resonance is controlled via a detuning parameter l ] such that " 34 " 7 3m n

With " corresponding to the envelope of the pump optical field at round-trip . In simulations we use a CW pump with a white noise floor 150dB lower than the pump that serves as a seeding signal for degenerate FWM interactions. The quantity l ] is the linear phase detuning between the pump signal and the closest cold cavity resonance.

In section 3, the Power enhancement factor of each resonator was established. As this quantity strongly differs at the pumped mode (1550 nm), pump power used in simulations has been adapted for each resonator so that the quantity ^_C 1550 ! ^p+)p stays constant for all three cases in order to be able to compare the simulation results in a compelling way. In simulations, pump power is adjusted to keep the quantity ^_C 1550 ! ^p+)p 210&. This comes to be pumping powers in simulations of 1.7 W for resonator 1; a simulated pump power of 4.2W for resonator 2 which has the lowest PEF at the pumped mode and a simulated pump power of 1W for resonator 3. Figure 5 shows the average intra-cavity power during laser scan simulations performed over the cavity resonance for each case of study with normalized pump power as well as the intra-cavity field temporal distribution and the intra-cavity field spectra for the 3 resonators. . These patterns correspond to a phase-locked frequency comb with a comb spacing equal to an integer number of cavity FSRs and are not observed in resonator 2. As the detuning of the pump laser continues, these generated sidebands widen and overlap filling the resonator modes between them [START_REF] Herr | Universal formation dynamics and noise of Kerr-frequency combs in microresonators[END_REF] . This leads to a broad but incoherent spectrum with a comb spacing equal to the resonator's free spectral range. The incoherent phase relationship between the spectral components of the comb translates to a chaotic temporal signal as seen for resonators 3 and 1. Finally, when the stable multiple pulse regime is attained, a spectrum with a well-defined phase relationship between comb teeth is generated inside the resonator [START_REF] Guo | Universal dynamics and deterministic switching of dissipative Kerr solitons in optical microresonators[END_REF][START_REF] Karpov | Universal dynamics and controlled switching of dissipative Kerr solitons in optical microresonators[END_REF] . The envelopes of the spectra generated are not smooth due to the fact that multiple pulses circulate inside the resonator at the same time. However, the temporal spacing between pulses remains constant and so the stable spectrum in the pulsed regime results from the superposition of each spectrum associated to each pulse circulating in the cavities.

Resonator 1 has a wide band high quality factor with a peak around mode -100 and a dip around mode 150. Low total losses allow for the formation of a large comb with circulating CW power of around 30W inside the cavity. Because losses remain of the same order of magnitude over the whole frequency span in simulations, the shape of the intra-cavity spectra generated in this resonator is shaped mainly by total resonator dispersion. Resonator 2 has a high quality factor over a narrow bandwidth close to the pumped mode that quickly degrades as mode number increases reaching lows of the order of magnitude 10 3 around modes 140 and -140 (fig 3.). A substantially higher CW power circulating inside the cavity (around 95W) is necessary in order for parametric gain to overcome total cavity losses as seen in figure 5. Furthermore, in the case of this resonator, the generated KFCs are shaped by total resonator losses as there is a very strong contrast of a factor of around 100 between the quality factor at the pumped mode and the lowest quality factor. The cold cavity mode structure of this resonator (fig 4 .) even suggest the potentiality of the generation of dispersive waves not far from modes 140 and -140 [START_REF] Bao | High-order dispersion in Kerr comb oscillators[END_REF][START_REF] Wang | Broadband Kerr frequency combs and intracavity soliton dynamics influenced by high-order cavity dispersion[END_REF] but with losses being most important in these particular regions, this effect does not take place. This resonator is an example of a case where frequency comb bandwidth is limited by the degradation with frequency of total resonator losses.

Resonator 3 presents what seem to be the best KFC generation properties of the 3 resonators with a frequency-dependent Q factor profile that remains highest on average of all three couplers and a flatter cold cavity mismatch curve than resonator 1. The first modulation sidebands are generated with circulating powers of around 30W as in resonator 1 due to the fact that both resonators have quality factor spectral profiles that maintain and order of magnitude of at least 10 5 . As is the case with resonator number 1, the spectrum generated with resonator 3 parameters is shaped by the cold cavity mode structure (i.e. total resonator dispersion).

CONCLUSIONS

In summary, we have presented in this paper the impact that access coupler geometry can have in the wideband transfer function of an optical resonator and we have presented by the means of numerical simulations, the consequences this entails for the generation of Kerr frequency combs. The modelling approach used in this work predicts that access coupler geometry will have a strong impact in the spectral distribution of the quality factor of the resonator as well as the resonator's mode structure (even more so than intra-cavity dispersion in this particular study). Given that both, losses and cold resonator mode structure, are affected by the coupler geometry, we have run numerical simulations based the Ikeda Map of the resonators in order to study Kerr Frequency comb generation in these structures. Simulations show that the degradation of a resonator's quality factor can strongly hinder the comb generation process and in turn have allowed us to determine that losses can become the limiting factor in comb bandwidth even when dispersion lends itself well for Kerr comb generation. These results highlight the importance of access coupler design particularly with vertically coupled Si 3 N 4 waveguide resonators. Furthermore, simulations predict promising Kerr comb generation results with access couplers 1 and 3 as designed by Arlotti et al an applied to 220 GHz high confinement Si 3 N 4 race-track resonators. Resonator 3 particularly predicts the formation of large frequency combs with as low as 1W of pump power with a structure that favors CMOSintegration (Si 3 N 4 ) and Kerr nonlinearity due to high confinement.

Figure 1 .

 1 Figure 1. Schematic of the simulated resonator and absolute value of the coupling and transmission coefficients of the 3 access couplers used in simulations. ρ 1 , ρ 2 : transmission coefficients. κ 1 , κ 2 : outward and inward coupling coefficients respectively. Coupler 1: Bold black lines. Coupler 2: Bold red lines. Coupler 3: Dashed black lines

Figure 2 .

 2 Figure 2. a) Representation of the cavity's second order dispersion parameter and wave vector as functions of mode number. b) Material absorption coefficient and power attenuation per roundtrip as functions of mode number.

Figure 3 .

 3 Figure 3. a) Representation of the Quality factor and Power enhancement factor as functions of modes number for the three simulated resonators. Black line: Resonator 1; Red line: Resonator 2; Gray line: Resonator 3.

Figure 4 .

 4 Figure 4. a) Cold cavity resonance mismatch induced by intra-cavity dispersion only (obtained using equation (6)). b) Cold cavity resonance mismatch induced by total roundtrip dispersion (intra-cavity and coupler) of each simulated resonator.

Figure 4

 4 Figure 4 depicts the cold cavity resonance mismatch computed without taking coupler dispersion into account in figure (a).This effect is mutual to all the simulated optical resonators and is obtained by using equation(6). Figure4(b) depicts the cold cavity resonance mismatch computed using equation (4) to determine the exact spectral position of each resonator mode with respect to an equally, FSR-spaced frequency grid. Figure4shows how in the resonators modelled by Arlotti et al., total round-trip dispersion is strongly affected by coupler geometry. This highlights the importance of taking into account dispersion effects during propagation through the access coupler.

Figure 5 .

 5 Figure 5. Evolution of intra-cavity field temporal distribution and intra-cavity spectrum during laser scan simulations for each simulated coupler. For resonators 1 and 3: I. Apparition of the first modulation sidebands. II. Formation of temporal patterns inside the cavity. III. Chaotic comb regime. IV. Multiple pulse regime. The different functioning regimes in frequency scan simulations can be identified in figure 5. After the apparition of the first modulation sidebands in fig 5.I) comes the formation of temporal patterns (II) as the detuning of the pump laser is performed

17

 17


Juan Arroyo
email: arroyo@us.es

Sometimes, sex is in the head

A recommendation of

Torices R, Afonso A, Anderberg AA, Gómez JM, and Méndez M. Architectural traits constrain the evolution of unisexual flowers and sexual segregation within inflorescences: an interspecific approach. However, attention has focused mostly on cases related to sex differentiation (dioecy and associated conditions [START_REF] Geber | Gender and sexual dimorphism in flowering plants[END_REF]). Separate unisexual flower types on the same plant (monoecy and related cases, rendering the plant functionally hermaphroditic) have been much less studied, apart from their possible role in the evolution of dioecy [START_REF] Charlesworth | Theories of the evolution of dioecy[END_REF] or their association with particular modes of pollination [START_REF] Friedman | A phylogenetic analysis of the evolution of wind pollination in the angiosperms[END_REF].

Two specific non-mutually exclusive hypotheses on the evolution of separate sexes in flowers (dicliny) have been proposed, both anchored in Lloyd's views and Darwin's legacy, with selfing avoidance and optimal limited resource allocation. Intermediate sex separation, in which sex morphs have different combinations of unisexual and hermaphrodite flowers, has been crucial for testing these hypotheses through comparative analyses of optimal conditions in PEER COMMUNITY IN EVOLUTIONARY BIOLOGY | DOI: 10.24072/pci.evolbiol.100069 suggested transitions. Again, cases in which floral unisexuality does not lead to sex separation have been studied much less than dioecious plants, at both the microevolutionary and macroevolutionary levels. It is surprising that the increasing availability of plant phylogenies and powerful methods for testing evolutionary transitions and correlations have not led to more studies, even though the frequency of monoecy is probably highest among diclinous species (those with unisexual flowers in any distribution among plants within a population [START_REF] Renner | The relative and absolute frequencies of angiosperm sexual systems: dioecy, monoecy, gynodioecy, and an updated online database[END_REF]).

The study by Torices et al. [START_REF] Torices | Architectural traits constrain the evolution of unisexual flowers and sexual segregation within inflorescences: an interspecific approach. bioRxiv 356147, ver 3 peer-reviewed and recommended by[END_REF] aims to fill this gap, offering a different perspective to that provided by Diggle and Miller [START_REF] Diggle | Developmental plasticity, genetic assimilation, and the evolutionary diversification of sexual expression in Solanum[END_REF] on the evolution of monoecious conditions. The authors use heads of a number of species of the sunflower family (Asteraceae) to test specifically the effect of resource limitation on the expression of sexual morphs within the head. They make use of the very particular and constant architecture of inflorescences in these species (the flower head or "capitulum") and the diversity of sexual conditions (hermaphrodite, gynomonoecious, monoecious) and their spatial pattern within the flower head in this plant family to develop an elegant means of testing this hypothesis. Their results are consistent with their expectations on the effect of resource limitation on the head, as determined by patterns of fruit size within the head, assuming that female fecundity is more strongly limited by resource availability than male function.

The authors took on a huge challenge in choosing to study the largest plant family (about 25 thousand species). Their sample was limited to only about a hundred species, but species selection was very careful, to ensure that the range of sex conditions and the available phylogenetic information were adequately represented. The analytical methods are robust and cast no doubt on the reported results. However, I can't help but wonder what would happen if the antiselfing hypothesis was tested simultaneously. This would require self-incompatibility (SI) data for the species sample, as the presence of SI is usually invoked as a powerful antiselfing mechanism, rendering the unisexuality of flowers unnecessary. However, SI is variable and frequently lost in the sunflower family [START_REF] Ferrer | Macrophylogenetic analyses of the gain and loss of selfincompatibility in the Asteraceae[END_REF]. I also wonder to what extent the very specific architecture of flower heads imposes an idiosyncratic resource distribution that may have fixed these sexual systems in species and lineages of the family. Although not approached in this study, intraspecific variation seems to be low. It would be very interesting to use similar approaches in other plant groups in which inflorescence architecture is lax and resource distribution may differ. A whole-plant approach might be required, rather than investigations of single inflorescences as in this study. This study has no flaws, but instead paves the way for further testing of a long-standing dual hypothesis, probably with different outcomes in different ecological and evolutionary settings. In the end, sex is not only in the head.

PEER COMMUNITY IN EVOLUTIONARY BIOLOGY | DOI: 10.24072/pci.evolbiol.100069

  bioRxiv 356147, ver 3 peer-reviewed and recommended by PCI Evol Biol (2019). DOI: 10.1101/356147 Plants display an amazing diversity of reproductive strategies with and without sex. This diversity is particularly remarkable in flowering plants, as highlighted by Charles Darwin, who wrote several botanical books scrutinizing plant reproduction. One particularly influential work concerned floral variation [3]. Darwin recognized that flowers may present different forms within a single population, with or without sex specialization. The number of species concerned is small, but they display recurrent patterns, which made it possible for Darwin to invoke natural and sexual selection to explain them. Most of early evolutionary theory on the evolution of reproductive strategies was developed in the first half of the 20th century and was based on animals. However, the pioneering work by David Lloyd from the 1970s onwards excited interest in the diversity of plant sexual strategies as models for testing adaptive hypotheses and predicting reproductive outcomes [1]. The sex specialization of individual flowers and plants has since become one of the favorite topics of evolutionary biologists.

Appendix