Z Ren

Thomas W Cornelius

C Leclere

A Davydok

J.-S Micha

O Robach

G Richter

O Thomas

Three-point bending behavior of a Au nanowire studied by insitu Laue micro-diffraction

Keywords: plasticity, geometrically necessary dislocations (GNDs), nanowire, in-situ threepoint bending test, Laue micro-diffraction

The elastic and plastic deformation of a gold nanowire tested in three-point bending configuration using the custom-built scanning force microscope SFINX was studied in-situ by Laue micro-diffraction. A new data treatment method, which bases on the integration of diffraction patterns recorded along the deformed nanostructure, is introduced visualizing both movement and shape of the diffraction peaks as a function of the measurement position.

Besides bending, torsion is evidenced during the elastic deformation originating from a misalignment of the SFINX-tip of the order of 60 nm with respect to the nanowire center. As demonstrated by post-mortem Laue micro-diffraction maps, the plastic deformation is governed by the storage of geometrically necessary dislocations. Analyzing the shape of the diffraction peaks, the activation of two unexpected slip systems is found which does not coincide with the slip systems with the highest resolved shear stress. These unexpected slip systems are probably related to the dislocation nucleation process at the clamping point, which is influenced by the local curvature.

Introduction

Since the first mechanical tests on microwhiskers in the 1950s [1,2], studies have demonstrated that the mechanical behavior of micro-and nanostructures differ from their bulk counterpart e.g. showing higher yield strengths with values near the theoretical limit of the materials for some of the defect-free macro-pillars and nanowires [3,4]. Most of these microand nano-mechanical tests were conducted in uniaxial testing configuration applying either a tensile or a compressive stress. When plasticity sets in, dislocations tend to glide through the structures before they can interact with each other [5], at variance with what happens at macroscopic scales. In tension tested single-crystalline Au nanowires, for example, dissociated dislocations are generated as leading partials with higher Schmid factors. After nucleation, these partial dislocations rapidly glide through the nanowire leaving stacking faults behind, which accumulate and eventually lead to Σ3 twin formation [6][7][8].

Bending is evidently a more complex mechanical test, which is less well studied especially at the nanoscale. However, the understanding of bending behavior is crucial for many practical situations such as for flexible electronics or MEMS devices. As pointed out by Nye in the 1950s [9], in bending tests, dislocations can be non-uniformly stored to minimize the strain energy during mechanical loading, and are thus required for the compatible deformation of various parts of the crystal. Recently, Konijnenberg et al. demonstrated that the bending deformation for FIB fabricated cantilever micro-beams is mainly accommodated by the storage of geometrically necessary dislocations (GNDs) [10]. Contrary to these micron-scale cases, some recent molecular dynamics (MD) simulations suggest that at the nanoscale the bending is accommodated by the formation of wedge-shaped twins similar to the ones observed in uniaxial tension tests [11,12]. While most studies reported in literature focused on the cantilever bending of micrometer sized beams prepared via FIB milling [10,[13][14][15][16][17][18], reports on plastic bending experiments for defect-scarce nanostructures are limited. To further explore the deformation mechanism for the bending at the nanoscale, in-situ three-point bending test is performed on an as-grown Au nanowire in combination with Laue micro-diffraction in the present work. Based on the rotation of the Au crystal along the nanowire, the deformation profile is extracted, which is further compared with finite element method (FEM) simulations.

Combined with post-mortem Laue micro-diffraction maps and SEM characterization, the activated slip systems are discussed shedding additional light on mechanisms of plastic bending in defect-scarce nanomaterials.

Experimental

In-situ three-point bending tests are performed using the custom designed scanning force microscope SFINX [19], which is installed at the French CRG-IF beamline BM32 at the European Synchrotron (ESRF) in Grenoble (France) as schematically illustrated in the Fig. 1(a). The incident polychromatic X-ray beam with an effective energy range between 5 and 25 keV is focused down to 400 nm (V) x 500 nm (H) on the sample surface using a pair of Kirkpatrick-Baez (KB) mirrors. The fluorescence yield is monitored using a Röntec XFlash 1001 energy-resolved point detector and the diffracted X-rays are recorded by a MarCCD165 detector with 2048 x 2048 pixels of 80 µm in size, which is mounted at a distance of 70 mm and 90° with respect to the incident beam (angular resolution of the detector around 0.065˚).

Single-crystalline and nominally defect-free Au nanowires were grown by physical vapor deposition on carbon coated tungsten or molybdenum substrates [3]. For three-point bending tests, the growth substrate was rubbed against a Si substrate patterned with 10 µm wide and 1.5 µm deep Si micro-trenches, transferring the as-grown Au nanowire. Some of the randomly distributed Au nanowires cross Si micro-trenches forming self-suspended nano-bridges. To avoid any sliding during mechanical testing, these Au nanowires are clamped to the Si supports by depositing carbon from the residual gas in a SEM chamber. As illustrated in Fig. 1(b), the growth direction of these Au nanowires is along 011 [7] and the Au nanowires lie down on the Si substrate with one of their <111> side facets [20]. The length of the nanowire is around 19 µm and the width of the nanowire is about 360nm, which is similar to the size of the X-ray beam.

Bragg coherent X-ray diffraction imaging as well as Laue micro-diffraction on similar nanowires revealed that the suspended nanowire parts are free of both strain and defects However, blurred diffraction patterns recorded for nanowire parts on top of the Si supports suggest that the transfer process might affect the strain state in these nanowire sections [20].

For nano-mechanical testing, firstly, pre-selected nanowires are located by AFM imaging using SFINX (noise around 10 nm) and mapping of the Au L III fluorescence with the X-ray beam.

Subsequently, the SFINX-tip is positioned above the center of the self-suspended nanowire and then lowered at a constant speed of 2 nm/s using the piezoelectric z-stage installed in SFINX.

The loading process is interrupted after each 50 nm displacement of the piezoelectric stage and Laue micro-diffraction patterns are recorded along the deflected nanowire in steps of 770 nm by rotating the downstream KB mirror in steps of 0.004° (calibration factor = 192 µm/degree).

Each KB scan takes around 3 min, where 13 diffraction patterns are recorded along the suspended part of the nanowire. In this way, the complete deformation profile of the deflected nanowire is measured via the recently developed KB scan approach, which avoids potentially destructive vibrations to the nano-mechanical testing equipment [21].

After the complete retraction of the SFINX-tip, a Laue micro-diffraction map is recorded along the deformed Au nanowire in steps of 500 nm employing the sample translation stage.

The finer step size yields a diffraction map with higher intensity and spatial resolution. However, the measurement positions are also slightly difference between the in situ and post mortem diffraction map, which is compensated by comparing both of them to the SEM image. A background Laue micro-diffraction pattern is recorded on the bare Si surface a few micrometers distant to the Au nanowire. The nanowire deformation is additionally simulated within continuum elasticity by the finite element method (FEM) using COMSOL Multiphysics taking into account the exact nanowire geometry, the elastic constants of bulk Au, and considering the Au nanowire being thoroughly clamped.

Laue microdiffraction data analysis

For data treatment, the background diffraction pattern is first subtracted from the nanowire diffraction patterns, removing the Si diffraction peaks and enhancing the visibility of the Au Laue spots (see Fig. S1 in supplementary material). The Au diffraction patterns are then indexed employing the LaueTools software [START_REF] Micha | LaueTools: Open Source Python Packages for X-ray MicroLaue Diffraction Analysis[END_REF], and the Au crystalline orientation is extracted as an orientation matrix (called UB matrix in the following) considering at least eight Laue spots (The X direction of the lab coordinates is along the beam direction, the Z direction perpendicular to the 2D detector). Employing this indexing process, one-dimensional orientation maps along the nanowire are generated, which are presented as the deviation angles θ for three distinct Au crystalline directions (Au [111], Au 011 , Au 2 11 ) (further details can be found in the supplementary material). These three crystalline directions correspond to the top surface direction, the nanowire growth axis, and the direction perpendicular to the loading plane, respectively.

One-dimensional orientation plots present the rotations of the crystal lattice which are caused either by elastic deformation or by the presence of unpaired dislocations. However, for weak diffracting objects such as the nanowires, extracting the absolute orientation by indexation is sometimes hindered by some eventual blurring of the diffraction peaks.

Furthermore, as pointed out by Nye [9], three-dimensional measurements of the crystal orientation are needed to calculate the curvature of the atomic lattice and to deduce the density of stored GNDs. Extracting the GND configuration from one-dimensional orientation maps requires additional assumptions, which compensate for components missing in Nye's tensor [START_REF] Hardin | [END_REF]. Thus, in the present work, an alternative experimental strategy is developed to further track the lattice rotation along the nanowire under mechanical testing. In this approach, the background subtracted diffraction patterns recorded at the same deformation state are summed up generating integrated diffraction patterns. Thus, each integrated diffraction pattern reflects a complete orientation profile of the deflected Au nanowire. Analogous to diffraction patterns measured using a pencil beam oriented parallel to the nanowire long axis, the displacements and the shape of the Au Laue spots along the nanowire are clearly visualized as continuous paths on the detector (see Mov. S1 in the supplementary material). Compared to a pencil beam measurement, the experimental approach described in the present work, however, provides additional information about both the measurement position and the shape of individual Laue spots. Similar to ref. [24], quantitative analysis of each individual diffraction pattern is then achieved by comparing these paths to the ones produced by the relative rotation of the Au crystal along the nanowire around defined crystallographic axes.The orientation matrix after a certain rotation is described with the following equation

UB UB • C • cos sin 0 sin cos 0 0 0 1 • ( 1
)
where ω is the rotation angle and C is the transformation matrix between the crystal coordinates (x, y, z axis along the three principal crystalline directions) and rotation coordinates (z' axis is the rotation axis, x', y' axes are orthogonal to the z'). Considering the geometrical parameters of the detector, these different orientation matrices can be expressed as different positions of the diffraction intensity on the detector, which eventually predicts the movement and elongation direction of the diffraction peaks after a certain rotation.

The movement of the diffraction peaks gives access to the rotation along the nanowire similar to 1D orientation maps. For a bending in the vertical plane, the Au crystal along the nanowire rotates around the 2 11 axis in the present work. The calculated deviation angles for the top surface direction (Au 111 ) and the growth direction (Au 011 ) of the nanowire coincide with the bending angles of the nanowire, while the deviation angle for the direction perpendicular to the loading plane (Au 2 11 ) is zero. The transformation matrix C can be written as

0 1/√3 2/√6 1/√2 1/√3 1/√6 1/√2 1/√3 1/√6 (2)
Combined with equation 1, the expected orientation matrix after this rotation is calculated, which predicts the expected path of the diffraction peaks for a pure vertical bending.

Considering edge dislocations, the 12 slip systems in FCC crystals induce rotations around the 12 different <112> axes, thus generating 12 independent movements and streaking directions of the Laue spots (as illustrated by Fig. S3 in the supplementary material). In Laue diffraction formalism [25,26], the principal elongation direction of the Laue spots is usually referred as ξ direction, which is perpendicular to both the dislocation line direction and the momentum transfer vector . While the activation of a single slip system causes movement and streaking of the diffraction peaks in a well-defined crystallographic direction, the simultaneous activation of multiple slip systems results in a movement that is defined by a linear combination of the corresponding crystalline rotations. Hofmann et al. [27] illustrated experimentally that already a single dislocation in a nanostructure leads to large bending angles as compared with what happens in the bulk. The induced pronounced rotation field may significantly afffect the peak shape as it will be demonstrated for a plastically deformed Au nanowire in the present work. This combination of diffraction peak shape analysis with a simple model provides additional information about the local rotation field within the probed volume giving access both to the activated slip systems and to the spatial distribution of the stored dislocations (see section 4.2).

Results

A sequence of integrated background subtracted Laue micro-diffraction patterns at different deformation states during a three-point bending test is presented in Fig. 2 On the other hand, the deviation angle along the Au 2 11 , which is the direction perpendicular to the loading plane, does not exceed 4° all along the nanowire. Since the diffraction map is recorded after the complete unloading, the elastic strain within the nanowire is neglected, in a first approximation, and the large variation of the crystalline orientation observed is thus attributed to the storage of geometrically necessary dislocations (GNDs).

Activated slip systems

In the post-mortem Laue micro-diffraction map (Fig. 4 axes with ratios of 9.5 and 5.8, respectively (see Fig. 6). While this linear combination describes well the shape of the diffraction peaks between the left Si support and the kink (positions i and ii) as well as for some of the diffraction patterns close to the right Si support (position iv), the measured Laue spots increasingly deviate from this prediction when approaching the kink from the right-hand side. This presumably is caused either by elastic strain or by additionally activated slip systems.

Number of stored GNDs

According to Nye [9], the plastic rotation of a crystal by an angle Δθ is directly related to the density of stored GNDs:

ρ ∆ ( 3
)
where b is the Burger's vector, and L is the probe length along b. Considering a single edge dislocation and t being the effective thickness along the slip plane direction, the lattice rotation induced by a single dislocation is described by the following equation:

∆θ b/t (4)
The effective thickness in nanomaterials is limited by the specimen size. As illustrated experimentally by Hofmann et al. [27], a single dislocation line within a nanostructure induces a sharp change of the lattice rotation similar to a tilt boundary in a bulk specimen (see Fig. given by the slip system and the angle are thus fixed. In the case of dislocations with screw components, however, the line direction of the slip system is no longer fixed, which may generate the variation of the line directions. Taking this into account would, however, require three-dimensional calculations which are out of the scope of this work.

Discussion

The traditional indexation process measures the crystalline orientation on a 3D real space grid allowing for decomposing the curvature tensor into contributions from different slip systems [9]. The peak shape analysis developed in this work can be considered as a projected 3D measurement of the crystalline orientation. Since there are only 9 independent components in Nye's tensor but 18 distinct dislocation types (edge and screw), no unique solution exists for the slip systems. However, considering the finite dimension of our specimen and the continuity of the dislocation lines within the specimen, the activation of the slip systems described above is the most reasonable solution.

Contrary to the findings by Laue micro-diffraction, the resolved shear stress (RSS) distribution, calculated for the 12 slip systems using FEM in the framework of continuum anisotropic elasticity considering the complete experimental configuration (angle between nanowire and Si micro-trench, misalignment of SFINX tip) for Δz = 1100 nm is maximal for the 1 01 11 1 , 110 11 1 , 101 1 1 1 , and 1 10 1 1 1 slip systems (see Fig. 7). It is thus similar to a perfect plane strain bending (see Table 1 in the supplementary material) despite a significant amount of torsion being present during the in-situ three-point bending test and tilt of the nanowire with respect to the Si support. Therefore, in spite of a very precise account of the imperfections of the nanomechanical test, they cannot account for the observed plastic behavior. The discrepancy of our experimental findings with respect to the simulations might be related to the dislocation nucleation process within the nanowire. Actual three-point bending tests contain points of stress concentrations at the loading position and at the clamping points (as illustrated by maximal RSS at these positions in Fig. 7). respectively. With increasing load, the shape as well as the movement of the diffraction peak gradually approaches the one after unloading, which is indicated by yellow lines (position i to ii in the post-mortem diffraction map). As pointed out by Kirchlechner et al [18], in the nanomaterial this broadening may be caused either by the storage of dislocations or by an inhomogeneous strain distribution. While the broadening of the diffraction peak for Δz ≥ 1200 nm is mainly due to the storage of GNDs, the slight deviation of the peak shape and the increased peak width at highest load (Δz =1300 nm) are attributed to a remaining torsion and an additional inhomogeneous strain within the loaded nanowire, respectively. The less pronounced streaking of the diffraction peaks from the Si support towards position ii implies that the GND density diminishes (see also Fig. 5

Conclusions

In conclusion, the elastic and plastic properties of three-point bent Au nanowires are investigated by in-situ Laue micro-diffraction. A new data treatment method employing the integration of Laue micro-diffraction patterns is described, which allows for visualizing the movement of the diffraction peaks along the nanowire and extract the actual elastic and plastic deformation. For the elastic deformation, both bending and torsion are present. While the former is evidently expected, the latter is caused by a misalignment of the SFINX-tip with respect to the center of the nanowire width of the order of 60 nm. For the plastic deformation, detailed peak shape analysis further allows for counting the number of stored dislocations and points to the activation of the 011 111 and the 101 1 11 slip system which are unexpected and in contradiction with recent MD simulations predicting the formation of wedgeshaped twins [11] as well as with the calculated critical resolved shear stresses. The experimental findings suggest that for defect-scarce nanomaterials, the final GNDs distributions might be influenced not only by the complete deformation profile, but also by the local rotation field at the dislocation nucleation position. Further insight in the plasticity of three-point bent nanostructures may be obtained by future dedicated experiments concentrating on the nucleation of the very first dislocations as well as by three-dimensional atomistic simulations.

  with the inset displaying the enlarged area around the Au 222 Laue spot. A movie showing the complete evolution of the integrated diffraction patterns can be found in the supplementary material (see Mov. S2). In the initial state (see Fig.2(a)), when the SFINX-tip does not touch the nanowire, a slight spread of the diffraction peak positions is apparent, which corresponds to a small variation of the Au crystalline orientation on the order of 2°. This small deformation might be related to an unintentional folding of the cantilever part of the nanowire during the alignment process (see post-mortem SEM image in Fig.4(a)). At a movement of the piezoelectric z-stage of 200 to 300 nm the SFINX-tip gets into contact with the nanowire and the integrated diffraction patterns start evolving continuously. With increasing load, the Laue diffraction peaks describe larger and larger elliptical paths on the detector (see Fig.2 (b-d)). At the piezo movement of 1200 nm, the paths of the diffraction peaks split into two semi-circles, which are separated around the loading point. In a first approximation, the continuous change observed at small movements of the piezoelectric stage is attributed to the elastic bending of the Au nanowire, while the discontinuous change at z = 1200 nm is ascribed to plasticity.

4. 1

 1 Fig. 3(c) and (d), respectively. With increasing movement of the piezoelectric stage the peak-to-

4. 2 Plasticity 4 . 2 . 1

 2421 Fig. 4(b)), the nanowire is divided into three sections, which are marked with four points. In the

  (b)), the movement of the Laue spots does not coincide with any of the expected directions for the 12 slip systems generated by edge dislocations in FCC crystals (see Fig. S3 in the supplementary material). Instead, the change of the crystalline orientation shown in Fig. 4(a) can well be explained by the combination of two simultaneously activated slip systems. Inverse pole figures calculated for the <112> poles (as presented in Fig. S4 in the supplementary material) reveal that the measured crystalline orientation follows very accurately the expected rotation induced by the 011 111 slip system for the 1 2 1 pole and the 2 1 1 pole, whereas it does not match any of the expected rotations generated by other slip systems for the other <112> poles. This leaves only two possibilities for a second activated slip system: either 101 1 11 or 011 11 1 . The shape of the diffraction peaks give access to the local rotation field induced by GNDs stored in the material, and thus provide additional information on the activated slip systems. Comparison of the streaking directions induced by the two aforementioned slip systems with the shape of the diffraction peaks (as illustrated by Fig. S5 in the supplementary material) yields a good accordance between the streaking caused by 101 1 11 dislocations and the vertical spread of the Laue spots, whereas the streaking direction induced by 011 11 1 dislocations does not show any correspondence with the shape of the Laue spots. Thus, the GNDs stored are inferred as combination of 011 111 and 101 1 11 slip systems, which induces rotations around the Au 2 11 and the Au 121 axes. The upper and lower borders of the trapezoidal shaped diffraction peaks can be well reproduced by a linear combination of rotations around these two

6

 6

  (a) and (b)). When two edge dislocations of two different slip systems are present in a nano-object (illustrated by the red and blue lines in Fig. 6(c)), each of them induces a separate rotation around its corresponding <112> axes. With different line directions, the two dislocations divide the crystal into four parts exhibiting four different combinations of rotations. This distribution of lattice rotations spreads the diffraction peak into two directions on the detector. For the 011 111 and the 101 1 11 slip systems observed in the nanowire, multiple dislocations are stored within the nanowire. As indicated by the blue and red lines in Fig. 6(d), the dislocation line for edge dislocations of 011 111 type is perpendicular to the nanowire axis while it is tilted by 30° for the 101 1 11 type dislocation. Considering a nanowire thickness of 90 nm, each 011 111 and 101 1 11 type dislocation induces a rotation of 0.18° around Au 2 11 and 0.06˚ around Au 121 , respectively. Within the probed volume, for one rotation angle around the Au 2 11 axis, a multitude of rotation angles around 121 is present leading to the vertical lines separating the diffraction peaks in Fig. 5. The ratio of these two rotations differs across the nanowire width corresponding to the upper and lower border of the diffraction peak. Connecting the positions with the same rotation angles on the two sides of the nanowire width (see Fig. S6(e) in the supplementary material), the distribution of the stored GNDs is extracted as presented by Fig. 6(e). A maximum of 15 011 111 and 6 101 1 11 type GNDs is found close to the left clamping point (position i) and the number of stored dislocations steadily decreases towards position ii. The angle between the line directions of the two slip systems varies between 27˚ and 58°. Considering only edge dislocation, the line direction is

  At the loading position (between position ii and iii in the post-mortem diffraction map), a kink is observed after plastic deformation. While the post-mortem SEM image may suggest a glide band where most of the dislocations escaped at the nanowire free surfaces, the nondetectable low diffraction intensity at this position indicates a strongly perturbed crystalline lattice caused by the storage of a large number of stored GNDs at this position. Comparing the integrated Laue micro-diffraction patterns recorded in-situ with the post-mortem pattern provides further information on the starting point and the evolution of plasticity. A sequence of the integrated Au 331 diffraction peak extracted from the left-hand side of the loading point for piezo movements of 1100, 1200, and 1300 nm are presented in Fig. 8(a), (b), and (c),

  (a)). Hence, it is reasonable to assume that dislocations observed after unloading are nucleated at the clamping position. As demonstrated in section 4, the streaking of the Laue spots mainly corresponds to the rotation of the Au crystalline lattice.In the vicinity of the clamping point, the peak broadening remains the same during and after loading indicating that the local lattice rotation at this position does not change. Thus, at this position the GNDs are stored according to the curvature tensor as described by Nye[9]. Along the nanowire, the broadening gradually spreads from the clamping point. The final GND distribution might thus not only be influenced by the stress state determined by the complete wire profile, but also by the local rotation at the nucleation position. The local rotation in the vicinity of the clamping point is affected by the orientation the nanowire with respect to the Si ridges which is about 114˚ as illustrated in Fig.1. Since the Si ridges are more rigid than the Au nanowire, their edges impose a certain shape to the bent nanowire which differs from a perfect three-point bending test, thus leading to localized rotations. As shown by the FEM simulations in Fig.3(b), slight variations of the torsion exist close to the two clamping points, which deviate from the linear increase observed in the rest of the nanowire. The unexpected slip systems are energetically the most favorable ones accommodating the best the deformation geometry in the vicinity of the clamping. They might be the L 2 minimized solutions for Nye's tensor at the clamping point[9].Before the discontinuity in the evolution of the integrated diffraction patterns at Δz = 1200 nm (Fig.8(a)), broadened diffraction peaks are already observed near the Si support, which might be related to the nucleation of dislocations. While the in-situ study presented here shows the overall behavior (elastic and plastic) of the Au nanowire during a three-point bending test, the loading steps are too coarse to detect the very first dislocations nucleated and stored within the nanowire. In addition, the inhomogeneous elastic strain within the bent nanowire broadens the diffraction peaks and, thus hinders the detection of the first defects. Therefore, dedicated in-situ investigations are needed with smaller loading steps and intermittent unloading to verify the elastic deformation and determine the nature of the GNDs at the clamping points.

Fig. 1 :

 1 Fig. 1: (a) Schematic diagram of experimental configuration at BM32. (b) Scanning electron

Fig. 2 :

 2 Fig. 2: Integrated background-subtracted Laue micro-diffraction patterns of KB scans recorded

Fig. 3 :

 3 Fig. 3: (a) Enlarged view of the area around the Au 222 diffraction peak at the piezo movement

Fig. 4 :

 4 Fig. 4: (a) 1D orientation map of the plastically deformed Au nanowire superimposed on the

Fig. 5 :

 5 Fig. 5: Integrated diffraction patterns from position i to ii in Fig. 4 of (a) the Au 222, and (b) the

Fig. 6 :

 6 Fig. 6: Schematic illustration of the lattice rotation induced by a single edge dislocation: (a) side

Fig. 7 :

 7 Fig. 7: The resolved shear stress distribution along the nanowire for the 12 slip systems in the

Fig. 8 :

 8 Fig. 8: The Au 331 diffraction peak integrated from the left Si support to the loading position

Acknowledgements

The authors gratefully acknowledge the financial support from the French National Research Agency through the project ANR-11-BS10-01401 MecaniX. They further thank the French CRG program committee for the allocated beamtime at the BM32 beamline at ESRF.

Supplementary Material

The supplementary online material demonstrates the background subtraction of Laue microdiffraction patterns, thus removing the Si Laue spots from the substrate and enhancing the visibility of the Laue diffraction peaks originating from the Au nanowire. It further describes the calculation of one-dimensional orientation maps from Laue diffraction patterns and gives additional details on the finite element method simulations. The analysis of the shape of the diffraction peaks that is introduced in this paper is further described and additional information on the determination of the activated slip systems is given. Moreover, a table with the Schmid factors for the tensile side of a bent nanowire is presented.


Martine Labbé
email: mlabbe@ulb.ac.be[

Miguel A Pozo
email: miguelpozo@us.es.

Justo Puerto
email: puerto@us.es

Computational comparisons of different formulations for the Bilevel Minimum Spanning Tree Problem

Keywords: Minimum Spanning Tree, Bilevel optimization, Stackelberg game

Let be given a graph G = (V, E) whose edge set is partitioned into a set R of red edges and a set B of blue edges, and assume that red edges are weighted and contain a spanning tree of G. Then, the Bilevel Minimum Spanning Tree Problem (BMSTP) is that of pricing (i.e., weighting) the blue edges in such a way that the total weight of the blue edges selected in a minimum spanning tree of the resulting graph is maximized. In this paper we present different mathematical formulations for the BMSTP based on the properties of the Minimum Spanning Tree Problem and the bilevel optimization. We establish a theoretical and empirical comparison between these new formulations and we also provide reinforcements that together with a proper formulation are able to solve medium to big size instances random instances. We also test our models in instances already existing in the literature.

Let G be a given a graph whose edge set is partitioned into a set of red edges and a set of blue edges, and assume that red edges are weighted and contain a spanning subgraph of G. Then, the Bilevel Minimum Spanning Tree Problem (BMSTP) consists in pricing (i.e., weighting) the blue edges in such a way that the total weight of the blue edges selected in a Minimum Spanning Tree (MST) of the resulting graph is maximized.

F o r R e v i e w O n l y
An example of the BMSTP is the following. Suppose a telecommunications company (TC) owns several connections (blue edges) between different nodes of a network. A new provider wants to enter into the market building a network that connects all nodes at minimum cost. It may use the connections of TC or those of its competitors (red edges). The target is to maximize the profit of the connections that TC could sell to the new provider.

The BMSTP can be seen as a bilevel optimization problem where the second level is a MST and where the objective functions are bilinear at both levels. Optimization problems related to spanning trees, or simply Spanning Tree Problems (STP), are among the core problems in combinatorial optimization. On the one hand, the combinatorial object that represents spanning trees has important structural properties. On the other hand, from a practitioner point of view, spanning trees are found in a wide range of applications in many fields (e.g. computer networks design, telecommunications networks, transportation, etc). Furthermore, they often appear as subproblems of other more complex optimization problems.

In game theory, a bilevel optimization problem is known under the name of Stackelberg game [START_REF] Von Stackelberg | Marktform und Gleichgewicht (Market and Equilibrium)[END_REF] and it consists in a game between a leader and a follower who play sequentially. Those players compete with each other: the leader makes the first move, and then the follower reacts optimally to the leader's action. This kind of hierarchical game is asymmetric in nature, where the leader and the follower cannot be interchanged. The leader knows ex ante that the follower observes its actions before responding in an optimal manner. Therefore, to optimize its objective, the leader anticipates the optimal response of the follower. In this setting, the leader's optimization problem contains a nested optimization task that corresponds to the follower's optimization problem.

Several papers have been published covering the problem in which the second level consists in choosing shortest paths between pairs of origins and destinations. The problem was first introduced by [START_REF] Labbé | A bilevel model of taxation and its application to optimal highway pricing[END_REF]. [START_REF] Roch | An approximation algorithm for stackelberg network pricing[END_REF] show that the problem is strongly NP-hard even when the second level consists in one single shortest path. More references regarding that bilevel optimization problem can be found in the surveys of [START_REF] Van Hoesel | An overview of stackelberg pricing in networks[END_REF] and [START_REF] Labbé | Bilevel programming and price setting problems[END_REF]. [START_REF] Gassner | Maximal spannende Baumprobleme mit einer Hierarchie von zwei Entscheidungstrgern[END_REF] studied a discrete variant of the bilevel minimum spanning tree problem where a partition of the set of edges into leader-and follower-edges is given. The leader's action is to choose a subset of his edges while the follower's reaction is to build up a spanning tree that includes the edges chosen by the leader. Hence, the leader's and follower's decision vectors are discrete. [START_REF] Cardinal | The stackelberg minimum spanning tree game[END_REF] proved the APX-hardness of the BMSTP even when the number of red edge costs is 2, and gave an approximation algorithm with guaranteed worst case performance. They also give an integer programming formulation for the problem and study its linear programming relaxation. Further [START_REF] Cardinal | The stackelberg minimum spanning tree game on planar and bounded treewidth graphs[END_REF] proved that the problem remains NP-hard even if G is planar, while it can be solved in polynomial time once that G has bounded treewidth. [START_REF] Bilò | Specializations and generalizations of the stackelberg minimum spanning tree game[END_REF] point out that the hardness in finding an optimal solution for the BMSTP lies in the selection of the optimal set of blue edges that will be purchased by the follower. Since once a set of blue edges is part of the final MST, their best possible pricing can be computed in polynomial time, as shown in [START_REF] Cardinal | The stackelberg minimum spanning tree game[END_REF]. [START_REF] Morais | A branch-and-cut-and-price algorithm for the stackelberg minimum spanning tree game[END_REF] introduce a reformulation and a Branch-and-cut-and-price algorithm for BM-STP. The reformulation is obtained after applying KKT optimality conditions to a BMSTP non-compact Bilevel Linear Programming formulation and is strengthened with a partial rank-1 RLT and with valid inequalities from the literature. They also implemented a Branch-and-cut algorithm for an extended formulation derived from another one in the literature and a preliminary computational study comparing

F o r R e v i e w O n l y both methods is presented.
In this paper we will present different mathematical formulations for the BMSTP based on the properties of the MSTP and the bilevel optimization paradigm. We establish a theoretical and empirical comparison between these new formulations that are able to solve medium to big size random instances. We also test our models in instances already existing in the literature taken from [START_REF] Morais | A branch-and-cut-and-price algorithm for the stackelberg minimum spanning tree game[END_REF].

The remainder of the paper is organized as follows. In Section 2 we formally define BMSTP and provide a new heuristic algorithm. Sections 3 and 4 present the catalogue of formulations that we study for the BMSTP including those corresponding to the MST subproblem. The empirical performance of the resulting BMSTP formulations is analyzed in Section 5, where we present extensive numerical results and a comparison with existing ones. Finally, some conclusions are summarized in Section 6.

Problem description and preliminary results

The BMSTP can be formally defined as follows. Let be given a graph G = (V, E) whose edge set E is partitioned into a set B of blue edges (controlled by the leader) and a set R of red edges, and assume that red edges are weighted and contain at least one spanning tree of G, thus, |R| ≥ |V | -1. A positive cost c e is associated to each red edge e ∈ R and a positive price T e (T = [T 1 , ..., T |E| ]) has to be determined for each blue edge e ∈ B. We denote by x = [x 1 , ..., x |E| ] the design variables used to describe the spanning tree polytope T .

Then, the Bilevel Minimum Spanning Tree Problem (BMSTP) consists in determining T e for each e ∈ B in such a way that the total weight of the blue edges selected in a minimum spanning tree of the resulting graph is maximized. A very general non-linear model for the BMSTP is then the following: BMSTP : max

Objective function (1a) maximizes the total weight of the blue edges selected in the solution. Constraints (1b) return the MST in the graph for a given cost vector T .

Example 1. Let G be a graph as depicted in Figure 1 left where red edges provide a spanning tree of total cost 20. Blue edges can be priced as in Figure 1 center or Figure 1 right in order to provide a BMSTP solution of value 15. We first observe (see [START_REF] Cardinal | The stackelberg minimum spanning tree game on planar and bounded treewidth graphs[END_REF] that in every optimal BMSTP solution the prices T e , e ∈ B take values in the set of red edges cost.

Property 1. In every optimal BMSTP solution T e ∈ C R = {c e : e ∈ R} for each e ∈ B.

Second, we observe (see [START_REF] Cardinal | The stackelberg minimum spanning tree game on planar and bounded treewidth graphs[END_REF] that the cost of each blue edge belonging to an optimal solution is bounded from above by the minimum among the maximum of the red costs of each cycle that contain the blue edge. Let C(e, S) be the set of cycles of G that include edge e and edges of the set S ⊆ E\{e}. Then, [START_REF] Cardinal | The stackelberg minimum spanning tree game[END_REF] provides a simple approximation algorithm, called Best-Out-Of-k algorithm: Let c 1 < ... < c |K| be the |K| different edge costs that appear in the initial red set of edges, where k ∈ K is the index of the k -th cost and C R = {c 1 < ... < c |K| } is the set of costs. The Best-Out-Of-k algorithm consists in choosing the best cost c k to be assigned to all T values. We can also rewrite this algorithm as a property in the following way: Property 5. (BMSTP feasibility) A feasible BMSTP solution (lower bound) is given by: In this subsection we generalize the Best-Out-Of-k algorithm by means of a local search algorithm that we denote BMSTP-H. Basically, this algorithm starts fixing the prices of the blue edges with a given value (current best solution) and iteratively, a subset of blue edges is chosen and their associated prices are increased/decreased. The BMSTP solution is then evaluated fixing the new prices of the blue edges. If an improvement is achieved, the current best known solution is updated. Otherwise, we reset the modified

F o r R e v i e w O n l y
T = (c k ) 1×|B| / k = arg
F o r R e v i e w O n l y
blue edges prices and we iterate until a stop condition is fulfilled. For a particular value of the search parameters, this algorithm has as a particular case the Best-Out-Of-k algorithm. However, BMSTP-H allows one to intensify the search by varying (1) the number of blue edges prices to modify, (2) increasing or decreasing the prices, (3) choosing or not blue edges already chosen in previous iterations and (4) modifying the stop condition. Additionally the algorithm can be run several times in a row for different values of the search parameters.

Note that the Best-

Out-Of-k algorithm is equivalent to BMSTP-H(b = |B|, p1 = 1, p2 = 0, sol best = ∅, ST OP c = "repeat |K| times").
In our experiments we apply three runs of BMSTP-H sequentially. The current best solution achieved on each module is used as an input in the following:

1. BMSTP-H(b = |B|; p1 = 1; p2 = 0; ST OP c = "repeat |K| times") → sol 1 . 2. BMSTP-H(b = 1; p1 = 0; p2 = 0; sol best = sol 1 ; ST OP c = "repeat until the set of chosen edges have size |B|") → sol 2 . 3. BMSTP-H(b = 1, p1 = 0, p2 = 1, sol best = sol 2 ; ST OP c = "repeat until the set of chosen edges have size |B|") → sol 3 .
We report on the performance of solutions found by the heuristic in Section 5. As an indicator Table 4 in Section 5 shows the number of times (in %) that the BMSTP-H algorithm reached the best lower bound for each BMSTP formulation.

Primal-dual BMSTP formulations

Let min x≥0 {cx : x ∈ T } be a continuous linear MSTP formulation and max µ≥0 {dµ : µ ∈ T D (c)} its dual form. We also denote by T D (c, T ) the polyhedron resulting by replacing c e by T e for each e ∈ B of T D . We therefore can express the BMSTP as

F 0 : max T ≥0 e∈B T e x e (3a) s.t.: x ∈ T (3b) µ ∈ T D (c, T ) (3c)
e∈B

T e x e + e∈R c e x e = dµ. (3d)
Since T and T D (c, T ) define linear domains, the strong duality theorem holds and the primal and dual solutions can be forced to take the same value as in (3d). Therefore, the revenue (3a) can be optimized ensuring that an optimal minimum spanning tree is chosen by the follower. In addition, if c e are integers, at least one optimal solution has integers variables T e as well (as previously mentioned) and, therefore, integrality constraints are not needed for variables x.

Note that the product of T and x variables makes F 0 non-linear. In this section, we firstly propose two linearizations of F 0 . Second we develop a polyhedral description of the Spanning Trees (ST) of G (coming from the one given by [START_REF] Martin | Using separation algorithms to generate mixed integer model reformulations[END_REF] and its dual form. Finally we propose different alternatives to describe the spanning tree polyhedron T within F 0 .

F o r R e v i e w O n l y 3.1. Linear BMSTP formulations
The product of variables T x can be linearized defining a new set of variables p e = x e T e , e ∈ B as the profit of edge e. In the following we denote by F p the linearization of (3a)-(3d) using variables p, that is [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part i -convex underestimating problems[END_REF].

F p : max T ≥0 e∈B p e (4a) s.t. x ∈ T (4b) µ ∈ T D (c, T ) (4c
Note that once T x is linearized, the integrality property of the problem is lost and integrality conditions (4i) are required.

We observe that variables T e and p e can be discretized, thus obtaining alternative formulations. Indeed, let {0, c 2 , ..., c |K| } be the set made up of the zero value and the |K -1| different edge costs that appear in the initial red tree, where k ∈ K is the index of the k -th cost and K is a sorted set in non-decreasing order. More precisely, we can also define the set K for each edge as

K e = {k ∈ K : m e ≤ c k ≤ M e }.
Let z k e be a binary variable equal to one if and only if the price of edge e is equal to the k-th cost, that is,

T e = k∈Ke z k e c k e ∈ B. (5)
Analogously, the values of p can be also discretized by using the binary variable zk e equal to one if and only if edge e is priced with the k -th cost, that is,

p e = k∈Ke zk e c k , e ∈ B. (6)
Therefore F p can be rewritten as follows:

F o r R e v i e w O n l y F z : max e∈B k∈Ke zk e c k (7a) s.t. x ∈ T (7b) µ ∈ T D (c, z) (7c) e∈B k∈Ke zk e c k + e∈R c e x e = dµ (7d) k∈Ke z k e = 1 e ∈ B ( 7e
)
k∈Ke

zk e = x e , e ∈ B ( 7f
)
zk e ≤ z k e , e ∈ B, k ∈ K e (7g) z k e ≤ zk e + (1 -x e ), e ∈ B, k ∈ K e (7h) z k e , zk e ∈ {0, 1} e ∈ E, k ∈ K e (7i) x e ∈ {0, 1} e ∈ E (7j)
Constraint (7e) ensures that each edge is priced with one of the costs. Constraint (7f) ensures that each edge leaves a profit coming from one of the costs if the edge is chosen by the follower and a null profit otherwise. Constraint (7g) ensures that the benefit of a blue edge is lower or equal than the profit of the edge. Constraint (7h) together with (7g) ensures that if a blue edge is chosen, its profit is equal to its price.

In addition, we note in constraints (7c) by T D (c, z) the polyhedron coming from T D (c, T ) when variables T e , e ∈ B are replaced by the values given in (5).

Therefore, note that (4e)-(4h) are implied by (7e)-(7j) as we show in the following property.

Property 6. Let Ω p,T ;x LR be the projection of the polyhedron defined by constraints (4b)-(4g) over the x variables and Ω z,z;x LR the projection of the polyhedron given by (7b)-(7h) over the x variables. Then

Ω z,z;x LR ⊆ Ω p,T ;x LR . Proof. ⊆) Let (x, z, z) ∈ Ω z,z;x
LR , we prove that (x, z, z) verifies (4e)-(4h). First, by replacing relations (5) and (6) in (7c) and (7d) we obtain (4c) and (4d). Now, we prove that (x, z, z) verifies (4e), ( 4f) and (4g):

• For each e ∈ B, (6) ⇒ p e = k∈Ke c k zk e ≤ M e k∈Ke zk e (7e)
= M e x e ⇒ (4e).

• For each e ∈ B, k ∈ K e , (7g) ⇒ zk e ≤ z k e ⇒ c k zk e ≤ c k z k e ⇒ k∈Ke c k zk e ≤ k∈Ke c k z k e ⇒ p e ≤ T e ⇒ (4f). • For each e ∈ B, k ∈ K e , (7g) ⇒ 0 ≤ z k e -zk e ≤ k∈Ke c k (z k e -zk e ) ≤ M e k∈Ke (z k e -zk e ) = M e (1 -x e ) ⇒ (4g). F o r R e v i e w O n l y 3.2.

KM-MST formulation and dual form

We recall that the general schemes proposed in the previous section require a linear STP formulation, T , and its dual form, T D . For that aim, we present next the Kipp Martim MSTP formulation and its dual form. [START_REF] Martin | Using separation algorithms to generate mixed integer model reformulations[END_REF] proposes a TDI formulation for the MST problem with a number of variables and constraints which are polynomial in the input size. For this, an arborescence rooted at each vertex k ∈ V is modeled. The arcs of such arborescences are then related to the design variables x defined above. For k ∈ V, (u, v) ∈ E, let us denote by q kuv and q kvu to the decision variables that respectively indicate whether or not arc (u, v) and (v, u) belong to the arborescence rooted at k, where A stands for the set or arcs. Then, the formulation is as follows:

min e∈E c e x e (8a)
s.t. e∈E x e = n -1 (8b)
s∈V :(k,s)∈A

q kks ≤ 0 k ∈ V (8c) v∈V :(u,v)∈A q kuv ≤ 1 k, u ∈ V : u = k (8d) q kuv + q kvu = x uv k ∈ V, (u, v) ∈ E (8e) x uv ≥ 0 (u, v) ∈ E (8f) q kuv ≥ 0 k ∈ V, (u, v) ∈ A (8g)
Constraint (8a) ensures that the tree has n -1 edges. Constraints (8b)-(8d) break cycles that could be generated by the n -1 edges. Note that if there is a cycle of undirected edges containing vertex k, then by (8b) there is a corresponding cycle of directed edges defined by q kij which also contains vertex k (refer to this set of directed arcs as the k-arcs). However, there cannot be a cycle of k-arcs which contains vertex k. This is impossible by (8d) and the cycle cannot be directed. If the cycle is not directed, then there is at least one vertex i with two k-arcs directed out of it. This is impossible by (8c). Thus, there are no cycles in the solution and we have by (8a) a spanning tree.

Formulation (8a)-(8f) can be simplified removing variables q kkv = 0, k ∈ V, (k, v) ∈ A as follows:

F o r R e v i e w O n l y min e∈E c e x e (9a)
e∈E

x e = n -1 (9b) (u ,v)∈E: (u =k∧v=u)∨ (u =u∧v=k) x u v + (u,v)∈A:v =k q kuv ≤ 1 k, u ∈ V : u = k (9c) q kuv + q kvu = x uv k ∈ V, (u, v) ∈ E : u, v = k (9d) x uv ≥ 0 (u, v) ∈ E (9e) q kuv ≥ 0 k ∈ V, (u, v) ∈ A : v = k (9f)
In the following, we denote by T km to the polyhedron associated to the KM-MSTP formulation, that is, equations ( 9b)-(9f). In addition, the dual of ( 9a)-( 9f) can be written as:

max α(n -1) - k∈V u∈V :v =k β ku (10a) α -β uv -β vu - k ∈V :k =u,v γ k uv ≤ c uv (u, v) ∈ E (10b)
-

β ku + (u ,v )∈E: (u =u∧v =v)∨ (u =v∧v =u) γ k u v ≤ 0 k ∈ V, (u, v) ∈ A : u, v = k (10c)
β ku ≥ 0 k, u ∈ V : u = k (10d) γ k uv k ∈ V, (u, v) ∈ E : u, v = k (10e)
Finally, we denote by T D to the polyhedron associated to the KM-MSTP dual formulation, that is, equations (10b)-(10e). Note that now, both T km and T D can be implemented in previous BMSTP formulations to set effective valid representations.

Other primal-dual BMSTP formulations

Previous BMSTP formulations assume constraint x ∈ T km . Note that this constraint can be replaced by any other STP polyhedron, namely subtour elimination, Miller-Tucker-Zemlim (MTZ), flow, etc (see [START_REF] Magnanti | Optimal trees[END_REF][START_REF] Fernández | Ordered weighted average optimization in multiobjective spanning tree problems[END_REF]. For example, the Miller-Tucker-Zemlim formulation [START_REF] Miller | Integer programming formulation of traveling salesman problems[END_REF] uses variables y uv which take the value 1 if and only if arc (u, v) belong to the arborescence and continuous variables l u , denoting the position that node u occupies in the arborescence with respect to the root node.

It is well-known that STP formulations exist with the integrality property. Unfortunately, when they are embedded within the BMSTP framework the integrality property is lost. Then, alternative STP formulations without such property may now be superior (in computational terms) and this explains why

F o r R e v i e w O n l y
some of the formulations we have used lack the integrality property. In Table 0 we resume the main properties of the STP formulations that we have considered. The criteria that have guided the selection of the formulations are either their good theoretical properties or some characteristic that seemed useful as, for instance, a small number of variables or constraints.

T sub e∈E(S) xe ≤ |S| -1, ∅ = S ⊂ V O(|E|) Exp(n) Yes Kipp Martin
Martin (1991)

T km (u,v)∈δ + (u) q kuv ≤ 1, k ∈ V, u ∈ V : u = k 0, k ∈ V, u = k ∀k O(n|E|) O(n|E|) Yes
Miller-Tucker-Zemlim

Miller et al. (1960)

T mtz lv ≥ lu + 1 -n(1 -yuv), (u, v) ∈ A r O(|E|) O(|E|) No Flow Gavish (1983) T f low (u,v)∈δ + (u) ϕuv - (v,u)∈δ -(u) ϕvu =      n -1, u = r -1, u ∈ V \ {r} r O(|E|) O(|E|) No KM extended Fernández et al. (2017) T km2 (u,v)∈δ + (u) quv ≤ 1, u ∈ V : u = r 0, u = r r O(|E|) Exp(n) Yes
In particular, the subtour and Kipp Martin formulations present the advantage of the integrality property (as you see in the last column) but they exhibit the inconvenience of an exponential number of constraints in the case of the subtour, or a cubic number of variables and constraints as it is the case of Kipp Martin. As an alternative, the MTZ and FLow formulations present lower dimensions and are, therefore, easier to handle.

We propose as an alternative STP formulation a relaxation of Kipp Martin that instead of building an arborescence at each node, it builds only one of them (see [START_REF] Fernández | Ordered weighted average optimization in multiobjective spanning tree problems[END_REF]. Therefore cut-set inequalities will be required to be included dynamically in a Branch&Cut algorithm. The separation of this inequalities can be carried out in polynomial time by finding the cut Gomory-Hu tree.

A path-based BMSTP formulation

In this section we present an alternative BMSTP formulation that does not require the strong duality property. Instead, we impose minimum cost optimality constraints to a path-based STP formulation. In this way, we can impose in the objective function a maximization of the benefits of the blue edges which in turns ensures the validity of this approach.

First, we propose the announced path-based formulation for the STP. Let P denote the set of pairs of nodes such that i < j. We define now ϕ ij uv as the flow through edge (u, v) going from origin i to

F o r R e v i e w O n l y
destination j with (i, j) ∈ P . The following set of constraints define a polyhedral description of the spanning trees of G.

T path :

v∈V :(i,v)∈A

ϕ ij iv = 1 (i, j) ∈ P (11a) (u,v)∈A ϕ ij uv - (v,u)∈A ϕ ij vu = 0 (i, j) ∈ P, v ∈ V : v = i, j (11b)
(u,j)∈A

ϕ ij uj = 1 (i, j) ∈ P (11c) ϕ ij uv + ϕ ij vu ≤ x uv (i, j) ∈ P, (i, j ) ∈ P, (u, v) ∈ E : u, v = i, j (11d)
(u,v)∈E

x uv = n -1 (11e) ϕ ij uv ≥ 0 (i, j) ∈ P, (u, v) ∈ A : v = i, u = j (11f) 0 ≤ x e ≤ 1 e ∈ E (11g)
Formulation ( 11a)-( 11g) defines a tree on the graph G. Constraints (11a)-( 11c) guarantee that the flow subnetwork is connected. Constraints (11d) ensure that an edge is selected if there is a positive flow traversing any of its arcs.

Note that we do not define variables ϕ ij ui and ϕ ij jv in T path since flow that is sent from i to j does not arrive to i or depart from j. In addition, constraints (11d) can alternatively be defined as

ϕ ij uv + ϕ ij vu ≤ x uv (i, j) ∈ P, (u, v) ∈ E : u, v = i, j. (12)
Replacing ( 11d) by ( 12) reduces significantly the number of variables and constraints in T path . However, with this enhancement the integrality property of variables x in T path is lost. Next, we can add an additional constraint to T path in order to guarantee this tree to be optimal (minimal cost):

(ϕ ij uv + ϕ ij vu )c uv ≤ c ij (1 -x ij ) (i, j) ∈ P, (u, v) ∈ E : (u, v) = (i, j) (13a)
Constraint (13a) impose that if there is flow sent from i to j along arcs

(u, v) or (v, u) (that is ϕ ij uv + ϕ ij vu = 0, what implies x ij = 0) then c uv ≤ c ij .
This let us formulate the BMSTP as follows:

F o r R e v i e w O n l y F path : max e∈B T e x e (14a) s.t. (x, ϕ) ∈ T path (14b) (ϕ ij uv + ϕ ij vu )(c uv + T uv ) ≤ (T ij + c ij )(1 -x ij ) (i, j) ∈ E, (u, v) ∈ E : (u, v) = (i, j) (14c)
vars:

x e ∈ {0, 1} e ∈ E (14d) ϕ ij uv ≥ 0 (i, j) ∈ P, (u, v) ∈ A (14e) T e ≥ 0 e ∈ B ( 14f
)
Note that for each (u, v) ∈ E constraints (14c) are well defined because we assume c e = 0, e ∈ B and T e only defined in B.

The linear formulation that removes the non-linearity of the product of p and T variables and the product of ϕ and T variables can be obtained by standard techniques. However, adding the new required inequalities destroys the integrality property of the tree polytope. Therefore, to have a valid formulation, constraint (15m) must be augmented. e ∈ B (15e)

t ij uv ≤ (ϕ ij uv + ϕ ij uv )M e (i, j) ∈ P, (u, v) ∈ B (15f) t ij uv ≤ T uv (i, j) ∈ P, (u, v) ∈ B (15g) T uv ≤ t ij uv + M e (1 -ϕ ij uv -ϕ ij vu ) (i, j) ∈ P, (u, v) ∈ B (15h) (ϕ ij uv + ϕ ij vu )c uv ≤ c ij (1 -x ij ) (i, j) ∈ R, (u, v) ∈ R : (u, v) = (i, j) (15i) t ij uv ≤ c ij (1 -x ij ) (i, j) ∈ R, (u, v) ∈ B : (u, v) = (i, j) (15j)
(ϕ ij uv + ϕ ij vu )c uv ≤ T ij -p ij (i, j) ∈ B, (u, v) ∈ R : (u, v) = (i, j) (15k)
t ij uv ≤ T ij -p ij (i, j) ∈ B, (u, v) ∈ B : (u, v) = (i, j) (15l) x e ∈ {0, 1} e ∈ E (15m)
ϕ ij uv ≥ 0 (i, j) ∈ P, (u, v) ∈ A (15n) T e ≥ 0 e ∈ B (15o)
t ij uv ≥ 0 (i, j) ∈ P, (u, v) ∈ A (15p) p e ≥ 0 e ∈ B (15q)
Constraints (15c)-(15e) provide a linearization of the product p e = x e T e , e ∈ B by means of the standard McCormick linearization. Constraints (15f)-(15h) provide a linearization of the product ϕT by means of variable t defined as

t ij uv = (ϕ ij uv + ϕ ij vu )T uv (i, j) ∈ P, (u, v) ∈ B (16)
and by means of the standard McCormick linearization. Constraints (15i)-(15l) provide a linearization of (14c) by means of variable t and distinguishing the different cases where

(i, j) ∈ R or B and (u, v) ∈ R or B.
There exists a similar formulation translating variables T and p into variables z and z. Given that such formulation has not provided good results in computational terms, we skip its description for the sake of readability.

Computational results

Next, we report on the results of some computational experiments that we have run, in order to compare empirically the proposed formulations. We have studied the BMSTP combining the different formulations proposed for the STP.

Instances G = (V, E) are generated as in [START_REF] Morais | A branch-and-cut-and-price algorithm for the stackelberg minimum spanning tree game[END_REF]. We first generate a complete graph G = (V, E c ) according to different values of |V | and the components of the cost vectors are randomly chosen from a set K of |K| random integers drawn from a uniform distribution on [1, c max ]. We then compute a MST solution (V, Ê) and we initialize E ← Ê, R ← Ê, B ← ∅. Additional edges are then randomly picked from E c \ Ê, until a desired graph density d is obtained. If a given edge e ∈ E c \ Ê is added to E, we choose with a probability p if R ← R ∪ {e} otherwise, B ← B ∪ {e}. If e is added to R, the cost assigned to e (c e ) is randomly chosen from a set K of |K| random integers. In particular, we choose |V | ∈ {20, 30, 50, 70}, c max = 150, d ∈ {10%, 20%, 30%, 50%}, p = 0.5, and |K| = {3, 5, 7}. Note that in this case the mean value of

|B| is |V |(|V |-2) 2 .
In all tables, reported results correspond to groups of 10 instances with the same triplet (|V |, d, |C|). We present average results (and some maximum values) for each group. This way, in total, we have a set of 240 benchmark instances. All instances were solved with the MIP Xpress 7.7 optimizer, under a Windows 10 environment in an Intel(R) Core(TM)i7 CPU 2.93 GHz processor and 16 GB RAM. Default values were initially used for all parameters of Xpress solver and a CPU time limit of 1800 seconds was set. We have also tested different combinations of parameters for the solver cut strategy and intensity of heuristics but, unless it is specified, the best results were obtained with the parameters of the solver set to the default values. An initial solution was given to the problem by the three modules of BMSTP-H described in Section 2.2. The separation of the cutset inequalities in formulation T km2 was implemented using a max-flow based algorithm [START_REF] Gusfield | Very simple methods for all pairs network flow analysis[END_REF].

Tables are grouped in blocks. The first block contains three columns with the values of the instances parameters. Then, we give a block of 7 columns for each tested formulation. The columns of each block are the following:

1. Columns gRL give the percentage relative gap, computed as 100 objR-obj L objR , where obj R denotes the optimal value of the linear relaxation at the root node and obj L denotes the best known lower bound obtained in all our experiments. 2. Columns gUL give the percentage relative gap, computed as 100 objU -obj L objU , where now obj U denotes the upper bound at termination. termination, respectively), and gUL provides a quality measure of the lower bounds. In addition, gUL and gUL * provide measures of both upper and lower bounds for average and worst case performance respectively. Entries with the symbol "-" indicates that the average/maximum gaps are 0, or in other words, that all instances were solved to optimality.

The caption just below each block gives the formulation the block refers to. Throughout the section, F mor denotes the formulation with the best results reported in [START_REF] Morais | A branch-and-cut-and-price algorithm for the stackelberg minimum spanning tree game[END_REF] 3 shows the number of times (in %) that the BMSTP-H algorithm reached the best lower bound for each formulation. 4. Table 4 choses the best blocks so far, and extends the time limit of some rows of these blocks from 1800 seconds to 5 hours. 5. Table 5 displays a comparison between the results reported in [START_REF] Morais | A branch-and-cut-and-price algorithm for the stackelberg minimum spanning tree game[END_REF] and the results obtained with our best formulation, with the same set of instances.

In order to facilitate the comparison among tables, best results in each table are marked in bold. In this sense, Tables 1 and2 In Table 1, the results of block F f low p exhibit the worst values of gRL in the group. However, most of these instances have similar values in terms of gUL than the other formulations of the group. As colummn nod shows, the number or required nodes to reach the optimal values in instances (20, 30, |C|) is the biggest in the group and consequently, gaps are bigger than in other blocks. On the contrary, the results of block F km p exhibit the best values of gRL. Howerver, in many cases these gaps only improve slightly (or not improved at all) other gRL values of the group. In addition, gaps gUL and gUL * remain far from the best values of the group. Note that according to the low average number of explored nodes in the B&B tree, in particular for sizes (|V |, d) = (70, 20), solving the LP relaxation of the problem becomes quite hard so the corresponding gaps at termination remain quite large in comparison with other formulations. We recall that F km p uses the largest number of variables and constraints, which can be too high in larger graphs. Block F mtz p shows good average gaps gRL and gUL, and it also provides some best values of the group for gaps gUL * in the largest instances. Block F km2 p shows a similar performance than F mtz p in gaps gRL and gUL, and it also provides the best values of the group for gaps gUL, gUL and gUL * in instances |V | < 70. Since many of the constraints in this formulation are added on the fly within the B&B search tree, we can observe that consecuently, the number of explored nodes is the largest of the group. From this table, we conclude that the most promising formulation is F km2 p for medium to large size instances. In Table 2, the results of block F f low z exhibit the worst values of gRL in the group. However, most of these instances do not remain far, in terms of gUL, from the other formulations of the group. As in the F shows good averge gaps gRL, gUL and gUL * but it is outperformed in general by block F km2 z that shows the best performace of this group. From this table, we conclude that the most promising formulation is F km2 z for medium to small size instances. Table 3 shows the number of times (in %) that the B&B search returned the same lower bound as the one given by the BMSTP-H algorithm. We observe that the percentages are smaller for the F (.) p formulations compared to the F (.) z formulations. This means that F (.) z formulations are harder to tackle in order to find feasible solutions (lower bounds) because the number of binary variables is highly superior to the one of the F (.) p formulations. Besides, we conclude from this table the BMSTP-H algorithm provides a reasonably good feasible initial lower bound that in many cases is no able to be outperformed by the solver after 1800 seconds of running time.

Table 4 shows a comparison of the best blocks so far, namely F km2 p and F km2 z , with a time limit of 1800 seconds and a time limit of 5 hours. In this case, we display results for only some rows of instances that correspond to the same combinations of (|V |, d, |C|) that were studied in [START_REF] Morais | A branch-and-cut-and-price algorithm for the stackelberg minimum spanning tree game[END_REF]. Obviously, mostly all gaps are improved when the time limit is extended to 5 hours but we observe a bigger improvement for F km2 z in small-medium instances and a bigger improvement for F km2 p in medium-big instances.

Table 5 shows a comparison between the results provided in [START_REF] Morais | A branch-and-cut-and-price algorithm for the stackelberg minimum spanning tree game[END_REF] and the results with our best formulations F km2 p and F km2 z . Note that in this case, we have used the same set of instances as in [START_REF] Morais | A branch-and-cut-and-price algorithm for the stackelberg minimum spanning tree game[END_REF] and only results for a single instance per row are provided. Block F mor shows the best results reported by [START_REF] Morais | A branch-and-cut-and-price algorithm for the stackelberg minimum spanning tree game[END_REF] implemented in C++ and tested with a 2.4GHz Intel XEON E5645 machine, with 32 GB of RAM, under Linux Operating System. From this table, we observe first that both F km2 p and F km2 z are able to provide better lower bounds objL and upper bounds objU than those in F mor . Consequently, gaps gU L are smaller in mostly all cases showing also that 11 out of 16 instances were able to be solved to optimality. In addition, running times (displayed in column t) show that 6 out of 16 instances were solved to optimality for F km2 p in less than 10 minutes. In general terms, we conclude that blocks F km2 p and F km2 z outperform significantly block F mor .

Conclusions

In this paper we have presented different mathematical programming formulations for the BMSTP based on the properties of the MSTP and the bilevel optimization paradigm. We have established theoretical and empirical comparisons between these new formulations that have shown to be effective for efficiently solving medium to big size instances. In addition we are able to outperform previous existing computational results in the literature.

Fig. 1 .

 1 Fig. 1. Graph with edge costs (left) and two BMST optimal solutions for |B| = |E| -(|V | -1)

  . (x, ϕ) ∈ T path (15b) p e ≤ M e x e e ∈ B (15c) p e ≤ T e e ∈ B (15d) T e ≤ p e + M e (1 -x e )

  3. Columns gUL give the percentage relative gap, computed as 100 obj U -objL obj U , where now obj U denotes the best known upper bound obtained in all our experiments and obj L denotes the lower bound at termination. 4. Columns gUL give the percentage relative gap, computed as 100 objU -objL objU . 5. Columns gUL * give the maximum of the gUL values among the 10 instances of the row. 6. Columns |#| indicate the number of instances in the group that could be solved to optimality within the CPU time limit. 7. Columns nod indicate the average number of nodes explored in the branch-and-bound tree.Note that, while gRL and gUL provide quality measures of the upper bounds (at the root node and at

  are treated as a sigle one, so as to highlight the best values among the eight proposed formulations, namely F f low p

p

  case, we observe that F f low z shows nthe worst gap values in general terms. The results of block F km z exhibit the best values of gRL and in some cases, these gaps are improved giving rise to best values of gaps gUL and gUL * in medium size instances. As in F km p solving the LP relaxation in sizes (|V |, d) = (70, 10) becomes quite hard giving rise to no more than three explored nodes in the B&B tree for these cases. Block F mtz p

  max According to p 1 , a subset B S of b blue edges is chosen from B ∪ S or from B\S. Edges e ∈ B S verifying T e > M e or T e < m e are removed from B S . According to p 2 , for each e ∈ B S increase/decrease by one unit k in c k e .

		k∈K	e∈B	x∈T c k x e : x = argmin { e∈B	c k x e +	e∈R	c e x e }
	2.2. A general framework for providing BMSTP feasible solutions
	Algorithm 1: BMSTP-H algorithm
		input :
		• sol best : Current best solution (by default T e = c |K| , ∀e ∈ B).
		• b: Number of blue edges to modify (by default b = |B|)
		• S: Set of edges that has been chosen in previous iterations (by default S = ∅)
		• p 1 : probability of choosing edges from B or from B\S (by default p 1 = 0)
		• p 2 : probability of choosing a direction of movement where "moving up" is chosen with probability
		p 2 and "moving down" with probability 1 -p 2 (by default p 2 = 0)
		• ST OP c : STOP condition (by default "repeat |K| times")
		output: sol best : Current best solution.
	1 while ST OP c = f alse do
	2
	6	Evaluate the BMSTP solution updating c k e for all e ∈ B S .
	7	if sol best is outperformed then
	8	Update sol best
	9	else
	10	reset values c k e for all e ∈ B S

3

S ← S ∪ B S . 4 5

  for the BM-STP. Otherwise, we denote by F

				(.) p the combination of the BMSTP F p formulation together with
	a spanning tree polyhedron T (.) (idem with F z ). In this way, we report results of formulations (.)
	F f low p	, F km p , F mtz p	, F km2 p	and F f low z	, F km z , F mtz z	, F km2 z	that have shown the best performance
	in preliminary tests. Note that we do not report results of the path formulations F path p	since those results
	were clearly outperformed by F p and F (.) z in preliminary tests. Then, we have summarized the results (.)
	in five tables.
	1. Table 1 shows results for F p formulations, namely F f low (.) p	, F km p , F mtz p	and F km2 p	.
	2. Table 2 shows results for F z formulations, namely F f low (.) z	, F km z , F mtz z	and F km2 z	.
	3. Table

Table 1

 1 BMSTP results for the F

	(.) p formulations.

  formulations: % of times BMSTP-H algorithm returned the lower bound of the B&B search Table4BMSTP results comparison for the best formulations with time limits of 0.5h and 5h.

	Table 3
	BMSTP-H results for F	(.) p and F	(.)
		|V | d F f low p	F km p	F mtz p	F km2 p	F	f low z	F km z	F mtz z	F km2 z
		20 30 33.3	33.3 33.3	33.3	33.3	33.3 33.3	33.3
		20 50 23.3	23.3 23.3	23.3	26.7	26.7 26.7	26.7
		30 30 30.0	30.0 33.3	26.7	66.7	50.0 50.0	33.3
		30 50 33.3	30.0 33.3	30.0	73.3	60.0 50.0	43.3
		50 10 30.0	20.0 16.7	16.7	46.7	16.7 20.0	30.0
		50 20 43.3	16.7 20.0	16.7	80.0	70.0 70.0	70.0
		70 10 60.0	50.0 50.0	43.3	93.3	86.7 90.0	90.0
		70 20 63.3	66.7 43.3	50.0	100.0 96.7 76.7	83.3
		total 39.6	33.8 31.7	30.0	65.0	55.0 52.1	51.3

z

Acknowledgements

The research of the first author has been partially supported by the Fonds de la Recherche Scientifique -FNRS under Grant(s) no PDR T0098.18. The second author was partially supported by the V Plan Propio de Investigación (Universidad de Sevilla). The second and third authors were partially supported by project MTM2016-74983-C02-01 (MINECO/FEDER). This support is gratefully acknowledged.

We also would like to acknowledge Vinicius Morais, Alexandre Salles da Cunha and Philippe Mahey for providing us their set of instances presented in Morais et al. (2016) that we have used in the present paper.

 1e4 20 50 3 4.8 1.6 -1.6 8.8 7 2e5 3.4 1.4 -1.4 8.6 8 5e4 3.4 1.8 -1.8 8.1 6 1e5 3.4 0.8 -0.8 7.8 9 1e5 20 50 5 7.1 5.8 -5.8 12 0 4e5 7.1 5.2 -5.2 10.6 1 2e5 7.1 5.2 -5.2 11.1 1 3e5 7.1 5.2 -5.2 11.4 1 4e5 20 50 7 8.9 7.8 1.7 7.9 16 1 3e5 8.9 7.9 1.6 7.9 16 1 2e5 8.9 8 1.7 8 16 1 3e5 8.9 8 1.7 8 16 1 4e5 30 30 3 6.9 3.4 0.3 3.4 14.5 6 1e4 4.4 3.3 0.3 3.3 14.8 5 1e4 4.4 3.2 0.3 3.2 15 5 2e4 4.4 2.8 0.3 2.8 13.3 6 2e4 30 30 5 9.8 6 3.1 6.1 14.1 2 3e4 7.2 6 2.9 6 14 2 2e4 7.2 6 2.9 6 14 1 5e4 7.2 6.4 2.9 6.4 14 1 6e4 30 30 7 13.4 9.4 6.8 9.5 16.2 0 6e4 10.1 9.4 6.8 9.4 15.2 0 2e4 10.1 9.4 6.9 9.5 15.2 0 5e4 10.1 9 6.8 9 13.3 0 7e4 30 50 3 0.4 0.2 -0.2 1.6 9 3e3 0.2 0.2 -0.2 1.6 9 1e3 0.2 0.2 -0.2 1.6 9 4e3 0.2 0.2 -0.2 1.6 9 4e3 30 50 5 4.1 3.8 0.8 4.1 10.3 1 3e4 3.8 3.8 0.5 3.8 7.4 1 1e4 3.8 3.8 0.5 3.8 7.6 1 3e4 3.8 3.8 0.5 3.8 7.4 1 5e4 30 50 7 5.9 5.7 4.4 6.5 21 0 3e4 5.7 5.7 4.3 6.3 19.5 0 1e4 5.7 5.7 4.5 6.5 21 0 3e4 5.7 5.7 4.2 6.2 18.5 0 6e4 50 10 3 16.2 4.1 1 4.4 13 2 1e4 3.7 1.9 0.6 1.9 10.9 7 2e3 4.4 2 0.6 2 10.4 7 2e3 4.8 1.8 0.6 1.8 9.6 6 5e3 50 10 5 19.3 5.3 1.5 6.3 9.7 1 1e4 5.5 2.4 0.5 2.5 7 5 3e3 6.3 2.9 0.5 2.9 7.2 4 7e3 6.8 2.3 0.6 2.4 7.3 5 1e4 50 10 7 22.5 6.3 2.8 7.4 18.2 2 1e4 8 3.7 1.9 4 14.9 5 2e3 8.6 4.7 1.7 4.8 13.2 4 7e3 8.6 3.4 1.6 3.4 10.5 5 9e3 50 20 3 4.1 2.6 0.6 2.6 6.6 1 5e3 2.8 2.4 0.6 2.4 6.6 2 1e3 2.8 2.5 0.6 2.5 6.6 1 4e3 3 2.7 0.6 2.7 6.6 0 8e3 50 20 5 9.4 8.2 14.4 14.8 25.1 0 5e3 8.2 8.2 9.7 10.2 17.3 0 1e3 8.2 8.2 9.9 10.3 16.6 0 5e3 8.2 8.2 8.7 9.1 14.5 0 8e3 50 20 7 13.8 12.1 20.1 20.4 30.8 0 5e3 12.1 12.1 15.7 16 25.9 0 1e3 12.1 12.1 15.5 15.8 24.1 0 5e3 12.1 12.1 14.3 14.7 26.6 0 9e3 70 10 3 16.3 7.6 7.3 11.2 21.3 0 1e3 7.3 6.8 6 9.2 19.2 1 5e2 7.5 6.5 6.4 9.3 19.1 1 9e2 7.5 6.7 6 9.1 19.5 1 1e3 70 10 5 18.7 13.1 15.7 16.4 26.6 0 9e2 13 12.4 15.8 16 26.7 0 1e2 13.1 12.8 15.4 15.9 27.3 0 8e2 13.1 12.5 14.4 14.6 23.7 0 9e2 70 10 7 21 14.8 22.7 23.4 33.2 0 9e2 15.1 14.4 22 22.3 33.2 0 1e2 15.2 14.5 22 22.3 30.4 0 7e2 15.2 14.6 19.9 20.4 32.7 0 1e3 70 20 3 1.1 1.1 4.7 5.1 18.6 0 1e3 1.1 1.1 5.7 6.1 24.6 0 25 1.1 1.1 1 1.4 3.8 0 7e2 1.1 1.1 0.9 1.3 3.8 0 5e2 70 20 5 4.7 4.6 7.7 7.7 23.9 0 2e2 4.6 4.6 7.7 7.8 24.3 0 3 4.6 4.6 5.4 5.4 13.2 0 4e2 4.6 4.6 5.2 5.3 13.8 0 6e2 70 20 7 8.1 8 11.9 11.9 23.7 0 2e2 8 8 11.9 11.9 23.7 0 0 8 8 10.1 10.1 21.8 0 4e2 8 8 10.8 10.8 21.8 0 8e2

Table 2 BMSTP results for the F -10 1e5 20 50 3 4.8 0.7 -0.7 6.8 9 5e4 3.4 0.8 -0.8 8.1 9 1e4 3.4 0.7 -0.7 6.7 9 3e4 3.4 ----10 2e4 20 50 5 7.1 2.5 0.1 2.6 8 1 5e5 7.1 0.8 0.6 1.4 8.1 5 7e4 7.1 0.7 -0.7 6.2 7 2e5 7.1 0.2 0.4 0.6 4.2 6 2e5 20 50 7 8.9 4.2 2.6 5.1 10.5 0 5e5 8.9 3.9 2 4.2 10.1 3 1e5 8.9 2.3 1.7 2.4 5.8 3 2e5 8.9 2.4 1.8 2.5 5.5 3 3e5 30 30 3 6.9 2.1 1.1 2.9 12 6 5e4 4.4 1 0.3 1 7.2 8 3e3 4.4 1.1 0.3 1.1 5.7 7 5e4 4.4 0.7 0.4 0.8 4.8 8 2e4 30 30 5 9.8 5 4.2 6.3 21.8 1 8e4 7.2 5.7 4 6.7 22.1 2 1e4 7.2 4.7 3.5 5.3 14 2 4e4 7.2 3.6 3.1 3.8 11.8 3 6e4 30 30 7 13.4 8.2 11.9 13.3 24.9 0 9e4 10.1 9.4 10.3 12.9 24.9 0 1e4 10.1 8.5 9.1 10.8 19 0 6e4 10.1 7.5 10.5 11.2 23.9 0 1e5 30 50 3 0.4 0.2 0.2 0.3 3.1 9 1e4 0.2 ----10 1e3 0.2 ----10 6e2 0.2 ----10 2e2 30 50 5 4.1 3.8 6.1 9.1 29.5 1 7e4 3.8 3.6 5.4 8.3 29.5 2 1e4 3.8 2.4 2.3 4.1 18.5 5 2e4 3.8 1.5 0.5 1.6 4.6 4 8e4 30 50 7 5.9 5.7 9.1 11 21 0 8e4 5.7 5.7 6.4 8.4 21 0 1e4 5.7 5.7 7.2 9 21 0 6e4 5.7 4.7 7.3 8.3 20 0 1e5 50 10 3 16.2 2.7 1.1 3.2 13.5 3 1e4 3.7 1.9 0.6 1.9 10 6 1e3 4.4 1 0.8 1.2 12.2 9 1e3 4.8 1.2 0.8 1.4 8.6 6 6e3 50 10 5 19.3 5.7 2.4 7.5 10.5 0 2e4 5.5 2.6 0.6 2.7 6.6 4 2e3 6.3 2.5 1.1 3.1 12.8 3 9e3 6.8 3 1.2 3.7 12.2 2 2e4 50 10 7 22.5 7.1 3.7 9 24.2 1 1e4 8 4.8 2 5.2 13.3 3 3e3 8.6 5.2 2.4 5.9 13.5 3 8e3 8.6 4.5 2.3 5.2 14 3 1e4 50 20 3 4.1 2.4 6.5 8 25.2 2 1e4 2.8 2.3 4 5.5 25.2 3 7e2 2.8 2.4 2.3 4 13.9 2 4e3 3 1.6 1.3 2.3 6.6 2 2e4 50 20 5 9.4 8.2 16.1 16.5 26.4 0 1e4 8.2 8.2 16 16.4 26.4 0 1e3 8.2 8.2 15.1 15.6 26.4 0 8e3 8.2 8.2 15.3 15.7 26.4 0 1e4 50 20 7 13.8 12.1 20.8 21.1 30.8 0 1e4 12.1 12.1 20.6 20.9 30.8 0 9e2 12.1 12.1 20.7 21 30.8 0 8e3 12.1 12.1 21 21.3 30.8 0 1e4 70 10 3 16.3 5.8 7.7 9.9 21.2 1 1e3 7.3 6.1 7.8 10.3 19.7 1 1e2 7.5 6.4 7.8 10.5 19.7 1 6e2 7.5 6.4 7.1 10 19.3 1 8e2 70 10 5 18.7 13.2 16.6 17.5 28.1 0 1e3 13 12.6 16.3 16.6 25.5 0 1e2 13.1 12.5 16.5 16.6 26 0 6e2 13.1 12.7 16.6 17 28.1 0 1e3 70 10 7 21 14.9 22.8 23.6 33.2 0 1e3 15.1 14.2 22.7 22.8 33.2 0 1e2 15.2 14.7 22.8 23.3 33.2 0 7e2 15.2 14.6 22.8 23.3 33.2 0 1e3 70 20 3 1.1 1.1 9.3 9.7 28.4 0 2e3 1.1 1.1 9.2 9.6 28.4 0 3 1.1 1.1 4.3 4.6 28.4 0 1e3 1.1 0.9 8.3 8.5 28.4 2 2e3 70 20 5 4.7 4.6 7.7 7.8 24.3 0 3e3 4.6 4.6 7.7 7.8 24.3 0 3 4.6 4.6 7.7 7.8 24.3 0 2e3 4.6 4.6 7.7 7.8 24.3 0 5e3 70 20 7 8.1 8 12.1 12.1 24.8 0 2e3 8 8 12.1 12.1 24.8 0 3 8 8 12.1 12.1 24.8 0 1e3 8 8 12.1 12.1 24.8 0 4e3

2e2 30 50 5 3.8 3.8 0.6 3.8 7.4 1 5e4 3.8 1.5 0.7 1.6 4.6 4 8e4 3.8 3.5 0.6 3.5 7.4 2 5e5 3.8 0.6 0.7 0.7 3.4 7 7e5 30 50 7 5.7 5.7 4.2 6.2 18.5 0 6e4 5.7 4.7 7.3 8.3 20 0 1e5 5.7 5.7 3.7 5.7 16.1 0 6e5 5.7 3.7 6.2 6.2 16.8 0 1e6 50 10 5 6.8 2.3 1.1 2.4 7.3 5 1e4 6.8 3 1.7 3.7 12.2 2 2e4 6.8 1.3 1 1.3 5.3 7 6e4 6.8 2.1 1.6 2.7 12.2 6 1e5 50 10 7 8.6 3.4 2.2 3.4 10.5 5 9e3 8.6 4.5 2.9 5.2 14 3 1e4 8.6 2.2 2.2 2.2 9 6 8e4 8.6 3.2 2.3 3.3 11 4 1e5 50 20 3 3 2.7 0.8 2.7 6.6 0 8e3 3 1.6 1.5 2.3 6.6 2 2e4 3 2.5 0.8 2.5 6.6 2 6e4 3 0.9 1 1.1 3.4 4 2e5 50 20 7 12.1 12.1 14.3 14.7 26.6 0 9e3 12.1 12.1 21 21.3 30.8 0 1e4 12.1 12.1 11.9 12.2 21.3 0 9e4 12.1 11.8 18.1 18.1 30.3 0 9e4 70 10 3 7.5 6.7 7.9 9.1 19.5 1 1e3 7.5 6.4 9.1 10 19.3 1 8e2 7.5 6.6 6 7.1 19.3 1 1e4 7.5 5.5 7.1 7.1 18.8 2 7e3 70 10 7 15.2 14.6 20.1 20.4 32.7 0 1e3 15.2 14.6 23 23.3 33.2 0 1e3 15.2 14.5 14.4 14.5 26.9 0 1e4 15.2 14.5 21.2 21.3 33.2 0 1e4 70 20 3 1.1 1.1 0.9 1.3 3.8 0 5e2 1.1 0.9 8.3 8.5 28.4 2 2e3 1.1 1 0.7 1 2.6 1 5e3 1.1 0.7 6.3 6.3 28.4 4 3e4 70 20 5 4.6 4.6 5.2 5.3 13.8 0 6e2 4.6 4.6 7.7 7.8 24.3 0 5e3 4.6 4.6 4.6 4.6 10.7 0 7e3 4.6 4.6 6.5 6.5 24.3 0 5e4