Valentina Bonino

Angelo Agostino

Carmelo Prestipino

Olivier Hernandez

Matteo Fretto

Lorenzo Mino

Marco Truccato

Structural and functional modifications induced by X-ray nanopatterning in Bi-2212 single crystals

Keywords: synchrotron, nanobeam, direct-write patterning, high-temperature superconductors, Bi-2212, oxygen, grain boundary, X-ray nanopatterning

Introduction

Lithographic methods lie at the heart of modern-day nanotechnology. Their development in the reduction of pattern motif in functional materials has enabled the fabrication of high-performing and cost-effective integrated circuits, storage devices, displays, biosensors, DNA microarrays and advanced materials 1 . The main improvements in the past few decades were mainly based on the decrease of the radiation wavelength, which has enabled a corresponding decrease of the minimum feature size, so that nowadays 193 nm ArF lasers represent the standard for high volume production.

However, further decrease towards shorter wavelengths at the moment is limited only to research and development activities, mainly due to the lack of suitable laser sources or photoresists compliant with industrial standards. For instance, this is still the case for Extreme UltraViolet (EUV) lithography exploiting λ = 13.5 nm, which is considered a good candidate for the sub-20 nm technological nodes and indeed has already shown the capability to produce patterns with resolutions down to 8 nm halfpitch in its interference version 2 .

A very different approach could be represented by the use of intense, highly-focused, hard X-ray beams to tailor the material growth and properties at the nanoscale. Indeed, X-rays produced by synchrotrons and free electron lasers have already proved to be able to stimulate the organized growth of nanostructures 3 , to manipulate defects [START_REF] Agostini | Characterization of semiconductor heterostructures and nanostructures[END_REF] and to modify the properties of materials [START_REF] Poccia | [END_REF] . From this point of view, the viability of an approach based on maskless X-ray nanopatterning has been recently demonstrated by Truccato et al. by fabricating a proof-of-concept electrical device out of the superconducting oxide Bi2Sr2CaCu2O8+δ (Bi-2212) 6 . This patterning method represents a conceptual novelty because it does not imply any material removal from the substrate and is just based on the spatial modulation of the lattice electrical properties, with great potential advantages in terms of heat dissipation, chemical contamination and mechanical stability of the devices.

However, the full development of this patterning technique requires a good understanding of the underlying chemical-physical mechanisms that connect the X-ray exposure with the changes in the electric properties, in order to finely tune the desired modifications. The processes involved are expected to depend on the material absorption properties and on dis-excitation mechanisms, that in turns depend on the material composition and structure, the X-ray wavelength, irradiation dose and possibly from the X-ray flux and exposure time. Indeed, several mechanisms have been evoked so far in literature for material modifications induced by X-rays, including for instance sample heating, formation of electronic defects or excited states, photoreduction of metallic ions due to photoelectron generation, production of highly reactive species from solutions as a consequence of radiolysis, and modifications of the liquid surface tension due to surface charge accumulation 7 . On the other hand, A c c e p t e d m a n u s c r i p t in the specific case of Bi-2212, the creation of oxygen vacancies and the increase of crystallographic mosaicity have been considered as possibly responsible for the increase of the superconducting critical temperature after irradiation 6 .

However, a detailed study of these modifications has not been performed yet, to the best of our knowledge. The present paper is devoted to fill this gap by studying the effect of successive X-ray exposures on Bi-2212 in order to obtain more insight on the structural phenomena underlying X-ray nanopatterning.

Experimental

The experimental strategy we have followed consists in repetitively exposing high-quality Bi-2212 microcrystals to X-rays, and in monitoring step by step the corresponding effects on their crystal structure by means of X-ray micro-diffraction, whereas changes in the superconducting and normal state properties have been monitored by means of electrical measurements. Due to time and experimental constrains at the synchrotron, a parallel sample strategy has been adopted where the samples have been produced with very similar initial features and underwent the same irradiation procedure. However, only one specifically prepared sample also underwent electrical measurements, which required the X-ray diffraction acquisition procedure to be adapted to its special features.

Sample preparation.

A set of straight whisker-like Bi-2212 microcrystals showing mirror surfaces and regular shapes has been selected under an optical microscope (100× magnification) from a batch produced following the procedure described elsewhere 8 . From this set two similar whisker crystals have been selected, the first one was used to study the irradiation effect in the absence of any physical constrain (the so-called freestanding crystal, i.e. with just one of its extremities fixed on a glass rod by glue), while the second one was integrated in an electric chip. More precisely, the second whisker was mounted onto a sapphire substrate with its c-axis normal to the substrate plane and then four electrical contacts have been fabricated via Ag physical vapour deposition. A picture of this sample is shown in Figure 1. The so-obtained chip has undergone an annealing at 450°C for 5 minutes (heating and cooling rates equal to 5 K/min) in 1 bar pure oxygen atmosphere in order to secure good electrical contacts minimizing possible changes in the crystal oxygen content. Indeed, the total oxygen diffusion length in the ab-plane during this thermal cycle is 1.7 µm 9 , which in principle could induce some variations in the oxygen content, but the fact that superconducting critical temperature Tc for this crystal in pristine conditions is Tc =78.3 K testifies that these changes, if any, are small enough to keep the crystal in its native overdoping regime 10,11 (see the following for more details). Finally, in order to create a reference point for the irradiation procedure, two Pt pillars about 1 µm high were deposited on the sapphire substrate in close proximity of the crystal by Focused Ion Beam (FIB)-assisted vapor deposition (FEI Quanta TM 3D FIB, with Ga-ion current of 10 pA). Scanning electron microscope (SEM) analysis has determined the sizes of the region of this crystal between the voltage contacts as 168.7×4.69×0.33 µm 3 (along the a-, b-, and c-axis directions, respectively), whereas the thickness of Ag contacts is about 2.1 μm. Figure S1 of the Supplementary Information displays the appearance of the central region of the sample at the end of the preparation procedures.

A c c e p t e d m a n u s c r i p t

Electrical measurements. The sample was electrically characterized along its ab crystallographic plane by the standard four-probe method. Measurements of the electrical resistance R were performed in a continuous flow Janis ST-100 helium cryostat in the temperature range 40 K ≤ T ≤ 290 K. The R A c c e p t e d m a n u s c r i p t vs T behaviors were acquired on warming by means of temperature ramps at the constant rate of 0.5 K/min and with a current value of 1 µA. The sample was stabilized at low temperature for at least one hour before every measurement. The substrate was accurately cleaned from any residuals and a thin layer of cryogenic high vacuum grease was applied to provide good thermal contact to the sample holder. Tc has been calculated as the inflection point of the R vs T curves, however, in order to analyze the presence of secondary superconducting phases, also a zero-resistance superconducting temperature Tc0 has been defined as the first non-zero resistance point above noise level.

X-ray set-up. The X-ray irradiation and characterization was performed at the microdiffraction ID13 beamline at the European Synchrotron Radiation Facility (ESRF) in Grenoble, France. The microbeam generated by an 18 mm period in-vacuum undulator was monochromatised to E 0 =13.0 keV and focused to a spotsize Sb = 2×2.5 µm 2 (horizontal × vertical). The beam flux Φ0 was estimated in 4×10 11 ph/s and fully used for crystal irradiation, whereas during the characterization the flux has been lowered by 2 orders of magnitude by an opportune absorber in order to prevent further significant exposure. The experimental station was equipped with a Frelon CCD camera 12 placed about 50 mm away from the sample. Tilt and distance of the camera were aligned using a LaB6 standard in a rotating capillary, finer adjustment has been done during data analysis.

The chip and the freestanding whisker-like crystal were alternately mounted on a motor stage composed by a three dimensional piezo translation stage (x,y,z) installed over a rotation stage (ϕ) with the rotation axis perpendicular to the storage ring plane, and intersecting the beam in the focus point, with the crystallographic a-axis indicatively parallel to the ϕ-axis. However, due to the sphere of confusion of the setup, the small size of the whisker and the limited allowed range of rotation for the mounted crystal, the alignment of the center of rotation, the whisker and the beam was not precise enough to keep the beam on the crystal along all the precession. For this reason, a horizontal translation motor perpendicular to the beam was scanned during irradiation or diffraction measurements with a step equal to the half width half maximum (HWHM) of the beam over a range greater than the alignment error estimated during sample alignment. This approach has allowed to deliver a constant energy profile to the crystal during the irradiation and to obtain the most significant sampling during diffraction.

Diffraction patterns have been acquired as a function of ϕ position every 0.2°. Due to the fact that the sapphire substrate is not transparent to the X-rays, the diffraction acquisition of the mounted-on-chip sample was limited to a ϕ range of 20° from the position with the b axis aligned with the beam, whereas the freestanding sample has been investigated over the maximum allowed range of 270°.

A c c e p t e d m a n u s c r i p t

The region of interest of the mounted-on-chip sample was located by taking as a reference the diffraction of the two Pt pillars.

X-ray data treatment.

As previously mentioned, for each angular position ϕ a set of frames (from 8 to 10, depending on sample alignment) has been collected, but due to the strong anisotropy and small sizes of the crystal along both the b-and the c-axis, only in one or two images the beam was actually on the sample. For this reason, a software has been developed to select the most significant frame for each angular position. For the freestanding sample the algorithm is very simple and selects for each set only the frame with the higher average pixel value, since when the sample is out of the beam only the small signal coming from air scattering is detected, whereas, when the sample is within the beam the signal generated by Bi-2212 is much stronger. For the mounted-on-chip sample the algorithm cannot use this scattering contrast anymore, because the average pixel intensity is dominated by the scattering of the substrate. In order to select the correct frames, we have taken advantage of the alignment of the BSCCO a-axis with the ϕ-axis: in this case all of the (00lm) reflections fall in a small region along the detector equatorial axis, so that the average of the pixel values has been calculated in this region only. Moreover, in order to avoid selection errors induced by streaks due to strongly saturated substrate peaks, the translation position of the selected frames has been fitted by a polynomial and the outlier positions have been revised manually. The subset of images obtained in this way has been analysed either by the CrysAlis software made available from Rigaku or subjected to an image treatment software developed ad hoc. The latter software allows performing pixel-wise maximum projections over a set of data (implementing a maximum intensity filter to reduce the contribution to the images stemming from saturating substrate reflections), and mosaic figures showing a small range of the detector along the ϕ rotation.

The evaluation of the c-axis length has been performed over the radially integrated average image by using the Le Bail fit algorithm 13 implemented on the FullProf code 14 . Radial integration has been performed with the help of the Fit2D code 15 with calibration parameters refined from LaB6 and corundum standards. Before integration, the averaged images were totally masked, with the exception of a stripe corresponding to the (00lm) reflections. This preliminary operation was mandatory in order to avoid the interference from the substrate reflections that are several orders of magnitude higher. It is worth noticing that the fits used to calibrate the Frelon camera were highly unstable due to the spottiness of the images collected for the standards, for this reason the sample-detector distance has been recalibrated by about 2% in order to obtain the expected value of the c-axis for the pristine sample.

A c c e p t e d m a n u s c r i p t

Irradiation sessions and XRD analysis

The aim of the present paper is to monitor the changes induced in Bi-2212 by multiple and prolonged X-ray irradiation sessions via both scanning X-ray diffraction and electrical measurements. The irradiation procedure was performed by scanning during each session the same region, which was placed 20 m above the Pt reference and extended over 7.5 µm along the a-axis and over the whole crystal thickness (0.33 µm) along the c-axis. The experiment followed a multistage procedure, performing in-series: i) a preliminary µ-XRD analysis, ii) an irradiation session, iii) a post-irradiation µ-XRD analysis, and iv) an off-line electrical characterization, and then repeating the sequence. The electrical measurement was carried out only for the mounted-on-chip sample. In order to assess the optimum irradiation time per point, the first irradiation session was performed in three stages of 480 s per point, at the end of which the first electrical measurement was performed. Subsequently, all of the three successive irradiations were carried out using 1440 s per point.

Irradiation # ∆t (s/point) F (10 11 J/m 2 ) CF (10 11 J/m 2 ) D (10 12 Gy) CD (10 Table 1 reports the irradiation details with the delivered fluence F and the absorbed dose D for each single irradiation, along with the corresponding cumulative values. The absorbed dose D is calculated from the Beer-Lambert law as:

𝐷 = Φ 0 Δt E 0 η m Φ = Φ 0 Δt E 0 (1-𝑒 - t 𝜆 𝑎 ) V Φ ρ ,
where the Δt is the irradiation time per point and mΦ the mass of the irradiated volume, which corresponds to the product between the crystal density ρ and the sample volume VΦ struck by the beam. The factor η = 1-exp(-t/λa) takes into account the decaying absorption of the incident radiation A c c e p t e d m a n u s c r i p t within a crystal with length t and attenuation length λa. Given E 0 =13.0 keV, λa for Bi-2212 is equal to 26.91 μm 16 .

Results and discussion

Bi2Sr2CaCu2O8+δ was discovered as a high-Tc superconductor by Maeda et a1. 17 Indeed, it is well-known that δ modulates the normal state electrical resistivity ρ (i.e. ρ decreases with increasing δ) and also affects Tc, even if in a non-monotonic way, because a single δ value exists corresponding to the highest Tc (the so-called optimal doping regime, OP), with lower Tc values corresponding both to lower δ values (the so-called underdoped regime, UD) and to higher δ values (the so-called overdoped regime, OD). Moreover, δ directly affects also the c-axis length, with higher δ values corresponding to shorter c-axis lengths. J/m 2 (c, see Table 1). Inspection of the pattern of the pristine state (Figure 2(a)) confirms that the sample remains crystalline after the electrical contact process, and that the frame is composed of three different features: i) a few isolated very strong peaks, due to reflections of the sapphire substrate, ii)

four continuous rings ascribed to the (111), ( 200), (220), and (311) reflections of Ag used for contact deposition, that are present only at lower angles, and finally iii) seven rows of peaks that can be ascribed to BSSCO reflections. In agreement with our experimental setup in which the exploration of the reciprocal space is strongly limited by the Ewald sphere curvature and by a rotation limit of 20°, the series of peaks can be identified with the (00lm), (10lm), (20l1), and (31lm) reflection series of Bi-2212, from the equatorial axis of the detector to the top, and with the corresponding series featuring opposite values for the h index in the lower part of the detector, with a small intensity asymmetry between the two regions of the detector due to the misalignment between the a and ϕ axes. It is possible to notice that the Ag contribution disappears after the irradiation cycles, probably because of the interaction between the Ag particles and the beam. Also a few changes in the substrate reflections are visible, probably due to a different alignment between the sapphire substrate and the crystal, which could have been originated from the thermal contractions and dilations induced by the cryogenic electrical measurements taken between the diffraction measurements; however, a nonconstant efficiency of the algorithm used to reduce the contribution of the substrate reflections to the image cannot be excluded. From Figure 2 it is also visible that the main effect of the radiation damage is on the Bi-2212 crystal mosaicity: indeed, the sample remains a single crystal but the reflections that look like circular spots in the pristine sample decrease in intensity and assume an elongated shape with increasing the cumulated fluence.

A c c e p t e d m a n u s c r i p t 1).

More detailed information can be obtained from Figure 3, which focuses on the (0 0 14) reflection and shows the evolution of the shape of the peak as a function of the fluence and of the ϕ position. In the pristine sample the peak can be considered as gaussian along the two directions of the detector axes and also along the ϕ rotation. Upon irradiation, the peak deforms irregularly in the vertical direction forming a small arc and spreads over a wider range in the ϕ rotation. The fact that the 2 direction less affected and the irregular nature of the other deformations indicate that such behavior should be ascribed to an increase of the crystal mosaicity and not to an increase of internal strain of the whisker. Indeed, the deformation of the gaussian spots of the pristine material into arc-like features with increasing the cumulative fluence is compatible with the interpretation that the irradiated portion of the crystal is fragmented in multiple nano-domains slightly misaligned and connected to each other by [010]-tilted grain boundaries. Moreover, the increase of the amplitude of the ϕ range where the reflection can be observed from about 2° for the pristine material up to about 4.4° for CF=30.708 × 10 11 J/m 2 seems to indicate that the irradiation also promotes the appearance of [100]tilted grain boundaries. We cannot exclude that this fragmentation process could be somehow affected The c-axis lattice parameters evaluated by Le Bail fit of the radially integrated intensities are listed in Table 2 and shown in Figure 4. The plots showing the integrated intensity and the quality of the fit are available in the Supplementary Information (see Figure S2 and Tables ST1 andST2). It can be noticed that the length of the c-axis increases with increasing the radiation fluence. This behavior confirms the general trend already observed in Bi-2212 by our group 18 , although the variation of the c-axis length shows some discrepancies with respect to our previous study. The key feature of the present experiment is represented by the acquisition of diffraction patterns via sample precession. This approach, compared to the previous nano-diffraction setup 18 , presents the advantage of allowing a more certain indexing of the pattern and to fully explore the spread of the reflections in the reciprocal space so that smaller uncertainties can be obtained and appearance of sub-domains can be detected.

In any case, these data confirm that, under some circumstances, part of the radiation damage induced by X-ray nano-beams in Bi-2212 is compatible with the mechanism of oxygen depletion, as recently modelled by our group 19 . The electrical resistance R of the mounted-on-chip sample is reported in Figure 5 irradiation the crystal was still continuous, even if clearly deformed (see Fig. S3). These results provide valuable information about both the superconducting and the normal state.

Concerning the superconducting behavior, the corresponding transition is highlighted in the inset of Figure 5(a), where the appearance of a tail with increasing the cumulative fluence is visible. This phenomenon can be properly described by comparing the critical temperature Tc to the zero-resistance temperature Tc0. Their behavior as a function of both CF and CD is shown in Figure 5(b). It can be noticed that no clear variation can be observed for Tc beyond the experimental error. Considering that the irradiated and non-irradiated portions of the sample are electrically in series and that the corresponding lengths are 7.5 µm and 161.2 µm, respectively, it is expected that the main contribution to the sample resistance comes from the non-irradiated part, unless the resistivity ρ of the irradiated part increases above a level about 161.2/7.5≈21 times greater than the pristine one. However, this is not the case for the curves of Figure 5(a), as it is discussed in more details in the following. Therefore, it is reasonable to associate the inflection point of the whole R vs T signal to the Tc of the nonirradiated part, which turns out to be practically constant, as expected. On the contrary, Tc0 is clearly affected by the irradiation process and its difference from Tc increases with increasing the cumulative fluence, up to a maximum difference of about 24 K after a cumulative fluence CF=15.354 x 10 11 J/m 2 . These modifications of the superconducting transition with increasing the irradiation dose are in agreement with the results obtained from previous studies, especially about the appearance of a tail even when only a portion of the sample with the shape of a trench is irradiated 6,18,20 . This experimental observation is compatible with a picture where the 7.5 µm long irradiated portion actually consists of a large number of very short (e.g. 100-200 nm) crystal domains, each of them with different Tc, critical current density Jc and ρ values that are functions of the local cumulative damage of the crystal structure and of the kind of interconnection with the other crystallites. Indeed, if we assume that some interaction between different irradiation points has taken place from the point of view of the damage, which somehow extends beyond the HWHM limits of the beam 19 , a spatial distribution of Tc's can be expected with the lowest values in the middle of the irradiated portion of the sample and increasingly higher values towards its boundaries. A similar behavior can also be expected for the local normal state resistivity ρ(T), with higher values in the center or the irradiated region and decreasing values approaching the pristine one towards its limits. With increasing the temperature, each damaged crystal domain switches to the corresponding normal state curve ρ(T), starting from the middle of the irradiated region and progressively adding its contribution to the sample resistance. However, each individual domain contribution to the sample resistance R is expected to be small, even if the local ρ(T) values can be considerably greater than the pristine value, A c c e p t e d m a n u s c r i p t because of the very short length of the domain. The superposition of all of these behaviors is expected to result in a slowly increasing R value of the sample as a function of T (i.e. the tail), until the Tc value of the large non-irradiated part is reached and a steep transition is observed.

Of course, this picture is in agreement with the indications coming from the diffraction measurements reported in Figure 3 about the progressive transformation of the single crystal towards a polycrystalline structure. From the point of view of the relationship between the changes of the crystal structure and the local values of Tc and ρ(T), in principle the possible interpretation is twofold. On one hand, if different crystal domains appear, then grain boundaries originate as well. Their effect on the properties of high temperature superconductors has been extensively studied 21 , and it has been proved that they deeply affect transport properties, mainly depending on the magnitude of the misalignement angle between adjacent grains. As an example, in YBa2Cu3O7-x a steep decrease of the intergrain critical current was found for misorientation angles above 10° 22 . Also in Bi-2212 Li Q. et al. proved that a 7° grain boundary was sufficient to suppress by about 9 K the Tc of a bicrystal film 23 . A similar effect was also shown in Bi-2223 for grain boundaries misorientations of the order of 30° 24 . In all of these studies, the modifications of the transport properties have been associated to the introduction of local structural disorder at the grain boundaries, which induces the presence of local inhomogeneity. On the other hand, irradiation could directly affect the structure of the individual crystal nano-domains. Indeed, concerning for instance YBa2Cu3O7-x , it has already been shown that electron irradiation can induce oxygen displacement and hole depletion 25 . About Bi-2212, it has been proved that ion irradiation can directly modify the non-stoichiometric oxygen content δ, inducing an increase of the c-axis length up to about 0.4 Å, corresponding to a transition from superconducting to insulating electrical behavior 26,27 . Moreover, oxygen removal from single crystals has also been induced by means of ordinary thermal annealing under controlled atmosphere: in the case of whisker-like Bi-2212 single crystals very similar to the ones used in the present experiment, a change in the c-axis length by about 0.15 Å has been observed to correspond to a change in the doping status from OD to UD 11 .

Actually, the data shown in Figure 3 and Figure 4 testify that both effects are present in our irradiated sample. Although each of them cannot individually explain the observed effects on the electrical behavior we have observed, they seem to act in a synergistic way to enhance each other. Indeed, the maximum crystal domain misalignment of about 4.4° that can be deduced from Figure 3(c) for CF=30.708 × 10 11 J/m 2 , if compared to previous results [22][23][24] should not significantly affect Tc, which would imply a practically constant value of Tc0 that is not confirmed by our measurements.

In the same way, an elongation of the c-axis length of about 0.09 Å like the one observed for CF=23.031 × 10 11 J/m 2 can be expected at most to result into a shift from the OD to the UD regime, This synergistic action of crystal fragmentation and oxygen depletion represents the key result of the present experiment and introduces a more complex mechanism with respect to the picture emerging from our previous work on nano-diffraction 15 . This mechanism also clarifies why oxygen removal by photogenerated electrons is not able to explain by itself all of the experimental data available for Bi-2212 19 , suggesting that the influence of grain boundaries is expected to be larger for dense irradiation meshes. Moreover, this could also imply that much care has to be taken when comparing experimental results originated by thermodynamic equilibrium processes, like for instance thermal annealing, to experimental results originated by processes quite far from thermodynamic equilibrium, like for instance ion and X-ray irradiation.

Concerning the normal state resistivity ρ of Bi-2212 material, its behavior as a function of T is displayed in Figure 6(a) for the pristine sample and for two different levels of irradiation. It is worth noticing that the ρ values of the irradiated material have been calculated under the simplifying assumption that the irradiated portion of the crystal features uniform structural damage and corresponding uniformly changed electrical properties. Consequently, the sample has been schematized as an electrical resistance representing the irradiated portion of the crystal (whose uniform resistivity is a function of irradiation) that is electrically in series with the remaining pristine part. The corresponding mathematical model allows for the determination of ρ(T) for each irradiation level, apart from the region very close to Tc where some artifacts appear and that has not been reported for this reason. It is possible to see from Figure 6(a) that the resistivity at room temperature ρ(RT) increases from about 2 × 10 -6 Ωm for the pristine material up to about 10 × 10 -6 Ωm after a cumulative fluence CF= 7.677 × 10 11 J/m 2 and to 16.6 × 10 -6 Ωm after CF= 15.354 × 10 11 J/m 2 . The fact that X-ray irradiation has increased the material resistivity by a maximum factor of about 8 confirms that most of the resistance signal comes from unirradiated part of the sample, as previously hypothesized for the interpretation of Figure 5(b), therefore making our description of the system self-consistent. The monotonous increase of the resistivity with increasing the fluence shown in Figure 6(a) is in general agreement with what already observed for the radiation damage originated by both particles and photon sources 18,20,27,[29][30][START_REF] Aldica | [END_REF] . More in details, it is possible to notice that the ρ(RT) values obtained from our analysis lie approximately in the same range as the ones reported by Watanabe et al. 32 and by Chen et al. 28 for Bi-2212 single crystals with different oxygen content δ induced by means of If compared to the results by Watanabe et al. 32 following the procedure recently exploited by Torsello et al. 19 , the ρ(RT) values of Figure 6(a) would correspond to the non-stoichiometric oxygen contents δ= 0.279, δ= 0.209 and δ= 0.199 for the pristine, the CF= 7.677 × 10 11 J/m 2 and the CF=15.354 × 10 11 J/m 2 curves, respectively. However, in the case of the irradiated material, these values certainly represent an underestimation of the real oxygen content δ because they ascribe all of the resistivity increase induced by irradiation only to the oxygen loss, neglecting the effect of grain boundaries. In this sense, these estimations have to be considered only as lower limits for δ of irradiated Bi-2212 in our experiment.

On the other hand, concerning the contribution of grain boundaries to the resistivity of the irradiated material, we can say that, from a qualitative point of view, our data seem not to support the brickwall model for the current flow across grain boundaries of high temperature superconductors 33,34 , mainly because no contribution can be seen in the ρ(T) curves from the c-axis resistivity, which should show a semiconducting behavior with a typical peak above Tc. Conversely, the railway-switch model 35, 36 seems to be closer to our experimental situation, even if our observation of a substantial increase in resistivity and of the presence of [100]-and [010]-tilted grain boundaries indicates that some adjustment is necessary. In any case, detailed measurements of Jc and of its dependence on the magnetic field is required to tackle this problem on a quantitative level.

Figure 6(a) also shows the fit of the high temperature portion of the ρ(T) curves to the linear relation ρ(T)= ρ0+αT. Indeed, it is well-known that in high-temperature superconductors the deviation from linearity of ρ(T) depends on the doping level of the material [37][38][39][40][41] . More specifically, the temperature where the ρ(T) curve departs from linear behavior is considered a good estimation of T*, i.e. of the temperature where the so-called pseudogap opens in the electronic DOS of the normal state 42 , as firstly observed both by ARPES 43 and STS 44 measurements. Figure 6(b) displays the same information in the form of the rescaled dimensionless behavior ρ'=(ρ(T)-ρ0)/αT of the resistivity, according to the representation proposed by Ito et al. 37 , where T* corresponds to the departure from ρ'=1 and is indicated by the arrows. It is possible to observe that the pristine curve shows the typical behavior of OD materials, with a smooth upturn below T* ≈172 K before a steep decrease near Tc.

On the contrary, the curve corresponding to CF= 7.677 × 10 11 J/m 2 shows a deviation from linearity at T* ≈ 247 K, whereas the one corresponding to CF= 15.354 × 10 11 J/m 2 deviates at T* ≈261 K. The general features of these latter two behaviors correspond to the gradual decrease of ρ' typical of UD materials, but a small oscillation is also noticeable on top of these curves just below T*, whose origin is not clear. Nevertheless, the increase of T* with the fluence confirms that X-ray irradiation has A c c e p t e d m a n u s c r i p t promoted some shift of the material doping from the OD to the UD regime, in agreement with the indications coming from the length of the c-axis and from the qualitative estimation of δ.

Conclusions

In the present work we have investigated the mechanism of X-ray nanopatterning in a Bi-2212 single crystal via combined XRD and electrical resistance measurements. It has been shown that cumulative fluence values of the order of 10 12 J/m 2 (corresponding to cumulative doses of the order of 10 13 Gy) induce both oxygen loss from the material and an increase of crystal mosaicity, with the appearance of multiple nano-domains and corresponding grain boundaries. Both effects contribute in a synergistic way to change the superconducting and normal state electrical properties, which suggests the possibility of non-linear effects in these changes as a function of the cumulated fluence. This possibility could be very important to finely tune the structural defects as desired, and highlights how a careful assessment of the irradiation time is crucial to achieve a good control of X-ray nanopatterning. However, the understanding of its microscopic mechanisms is still far from being complete and further investigation is necessary.

Figure 1 :

 1 Figure 1: Micrograph of the mounted-on-chip sample (WBLP04). The Bi-2212 crystal is the thin horizontal line at the centre of the image. The four Ag strips vertically crossing the crystal represent the electrical contacts and are labelled according to their connection to the measurement circuit.

  and features a complex, incommensurably modulated, layered structure composed of perovskite and rock salt subunits with superspace group Bbmb (0β1) and a pseudo-tetragonal cell (a ≈ b ≈ 5.4 Å, c ≈ 30.7 Å) with modulation vector q*=[0, 0.210, 1]. Its transport properties are strongly dependent on the nonstoichiometric oxygen content δ, which refers to the amount of oxygen atoms hosted in the interstitial sites between the BiO and SrO layers and determines the amount of hole-doping of the material.

Figure 2 :

 2 Figure 2: Maximum projection maps of the micro-diffraction patterns collected during a 20° scanning around the a-axis (ϕ-axis). The panels refer to the following irradiation conditions of the sample: pristine

Figure 2

 2 Figure 2 shows the maximum projections of a 20° ϕ-scan of the contacted whisker, with the a-axis aligned along the rotation axis itself in the vertical direction. This figure compares the sample in the pristine state (a), after a cumulative fluence CF= 7.677 × 10 11 J/m 2 (b), and after CF=30.708 × 10 11

Figure 3 :

 3 Figure 3: Mosaic image of the (0 0 14) reflection as a function of the ϕ rotation in steps of 0.2°, selected from the maximum projections of Figure 2. Panel (a) corresponds to the pristine sample, panel (b) to cumulative dose CD=3.909 x 10 12 Gy and cumulative fluence CF =7.677 x 10 11 J/m 2 (irradiation #1.c) and panel (c) to CD=15.630 x 10 12 Gy and CF=30.708 × 10 11 J/m 2 (irradiation #4, see Table1).

A c c e p t e d m a n u s c r i p t also by individual

  features of the crystals, including for instance the release of some strain imposed to the whisker during the procedure of chip fabrication. Indeed, Figure3(a)shows for the mountedon-chip pristine sample a significant intensity of the (0 0 14) peak for about 2° during the ϕ-rotation, therefore revealing a more broadened rocking curve compared to the one of the free-standing sample, which, under the same instrumental conditions, shows the presence of this peak only in a single diffraction image, with a corresponding rocking curve width of about 0.2° during the ϕ-rotation.

Figure 4 :Table 2 :

 42 Figure 4: c-axis length versus cumulative dose CD and cumulative fluence CF for the mounted-on-chip crystal at different stages of the experiment.

A c c e p t e d m a n u s c r i p t

Figure 5 .

 5 Figure 5. Panel (a): Four-probe R versus T behavior for the mounted-on-chip sample: before irradiation (blue), after a cumulative dose CD=3.909 x 10 12 Gy and cumulative fluence CF =7.677 x 10 11 J/m 2 (red), and after CD=7.816 x 10 12 Gy and CF =15.354 x 10 11 J/m 2 (green). The inset shows a magnification of the region around the superconducting transition. Panel (b): Critical temperature Tc (circles) and zero-resistance critical temperature Tc0 (squares) versus cumulative dose CD and cumulative fluence CF.

  (a) as a function of temperature for the pristine condition, after a cumulative fluence CF=7.677 × 10 11 J/m 2 (irradiation #1.c) and after CF=15.354 × 10 11 J/m 2 (irradiation #2). A subsequent irradiation up to CF=23.031 × 10 11 J/m 2 (irradiation #3) has induced a transition to an insulating state (R> 20 MΩ), preventing the measurement of the corresponding electrical behavior. A SEM inspection has revealed that after this

A

  c c e p t e d m a n u s c r i p t with Tc still in excess of 70 K 11 , whereas in other cases it has resulted in a change of Tc by 5-7 K only28 . On the contrary, our measurements have shown that values of CF ≥ 23.031 × 10 11 J/m 2 have turned the material into an insulating state, which is not expected on the basis of each individual effect.

Figure 6 .

 6 Figure 6. Panel (a): Temperature dependence of electrical resistivity ρ along with the straight lines corresponding to the extrapolation of their linear fit to the equation ρ(T)= ρ0+αT in the range 260-295 K. Blue lines refer to the sample before irradiation, red lines to the sample after a cumulative fluence CF= 7.677 × 10 11 J/m 2 and a cumulated dose CD = 3.909 × 10 12 Gy, and green lines after CF= 15.354 × 10 11 J/m 2 and CD = 7.816 × 10 12 Gy. Panel (b): Corresponding rescaled curves (ρ(T)-ρ0)/αT. Arrows indicate the temperature T* at which the pseudogap opens. OD and UD refer to the overdoped and to the underdoped regime, respectively.

  A c c e p t e d m a n u s c r i p t equilibrium thermal treatments, but are lower by a factor of 2-3 with respect to the values corresponding to -irradiated single crystals 27 .

Table 1 :

 1

	12 Gy)

Irradiation time ∆t per point of the rastered mesh, fluence F, cumulative fluence CF delivered to the exposed surface, dose D and cumulative dose CD estimated for the irradiated volume of the Bi-2212 samples.

Acknowledgments

The authors thank the ESRF for the beamtime allocated and the ID13 beamline staff for their valuable help. They also gratefully acknowledge Compagnia di San Paolo for financial support to NanoFacility Piemonte at INRIM. This work has been partly carried out under project NANO-X jointly approved and funded by University of Torino and Compagnia di San Paolo. V.B. and M.T. acknowledge partial support from the "Departments of Excellence" (L. 232/2016) grant, funded by the Italian Ministry of Education, University and Research (MIUR).

Conflicts of interest

There are no conflicts to declare.


Luciano Porretta

Daniele Catanzaro

Bjarni V Halldórsson

Bernard Fortz

Abstract A point-interval (I v , p v ) is a pair constituted by an interval I v of R and a point p v ∈ I v . A graph G = (V, E) is a Max-Point-Tolerance (MPT) graph if each vertex v ∈ V can be mapped to a point-interval in such a way that (u, v) is an edge of G iff I u ∩ I v ⊇ {p u , p v }. MPT graphs constitute a superclass of interval graphs and naturally arise in genetic analysis as a way to represent specific relationships among DNA fragments extracted from a population of individuals. One of the most important applications of MPT graphs concerns the search for an association between major human diseases and chromosome regions from patients that exhibit loss of heterozygosity events. This task can be formulated as a minimum cost clique cover problem in a MPT graph and gives rise to a N P-hard combinatorial optimization problem known in the literature as the Parsimonious Loss of Heterozygosity Problem (PLOHP). In this article, we investigate ways to speed up the best known exact solution algorithm for the PLOHP as well as techniques to enlarge the size of the instances that can be optimally solved. In particular, we present a Branch&Price algorithm for the PLOHP and we develop a number of preprocessing techniques and decomposition strategies to dramatically reduce the size of its instances. Computational experiments show that the proposed approach is 10-30x faster than previous approaches described in the literature, and suggest new directions for the development of future exact solution approaches that may prove of fundamental assistance in practice.

Introduction

Interval graphs (namely, intersection graphs of intervals on a line [START_REF] Booth | Testing for the consecutive ones property, interval graphs, and graph planarity using PQ-tree algorithms[END_REF][START_REF] Fishburn | Interval orders and interval graphs: A study of partially ordered sets[END_REF][START_REF] Fulkerson | Incidence matrices and interval graphs[END_REF]) constitute an important class of graphs in computer science and discrete mathematics both for their significance in practical applications [START_REF] Diestel | Graph Theory[END_REF] and because classical N Phard combinatorial optimization problems, such as the maximum clique problem [START_REF] Fulkerson | Incidence matrices and interval graphs[END_REF][START_REF] Garey | Computers and Intractability: A guide to the theory of NP-Completeness[END_REF], the weighted independent set problem [START_REF] Frank | Some polynomial algorithms for certain graphs and hypergraphs[END_REF][START_REF] Garey | Computers and Intractability: A guide to the theory of NP-Completeness[END_REF] and the coloring problem [START_REF] Garey | Computers and Intractability: A guide to the theory of NP-Completeness[END_REF][START_REF] Golumbic | Algorithmic graph theory and perfect graphs[END_REF], can be solved in polynomial time in this particular class of graphs [START_REF] Fishburn | Interval orders and interval graphs: A study of partially ordered sets[END_REF][START_REF] Fulkerson | Incidence matrices and interval graphs[END_REF].

In the last three decades, numerous research efforts have been carried out to extend the properties of interval graphs to more general classes of graphs [START_REF] Asinowski | Vertex intersection graphs of paths on a grid[END_REF][START_REF] Corneil | Extensions of permutation and interval graphs[END_REF][START_REF] Diestel | Graph Theory[END_REF][START_REF] Gavril | A recognition algorithm for the intersection graphs of paths in trees[END_REF][START_REF] Golumbic | Algorithmic graph theory and perfect graphs[END_REF][START_REF] Golumbic | A generalization of interval graphs with tolerances[END_REF]. In this context, an important generalization is constituted by tolerance graphs, first introduced in [START_REF] Golumbic | A generalization of interval graphs with tolerances[END_REF]20]. Specifically, a graph G = (V, E) is a tolerance graph (namely, a max-tolerance graph [20]) if every vertex v ∈ V can be associated with both an interval I v ⊆ R and a tolerance value t v ∈ R in such a way that (u, v) is an edge of G iff |I u ∩ I v | ≥ max{t u , t v }. In other words, G is a tolerance graph if each pair of intervals can "tolerate" a non-empty intersection (without forming an edge) as long as the length of this intersection is strictly smaller than the maximum of the corresponding tolerance values. As for interval graphs, solving the maximum clique problem on a tolerance graph can be performed in polynomial time both in the weighted and the unweighted case [START_REF] Kaufmann | Maxtolerance graphs as intersection graphs: Cliques, cycles and recognition[END_REF]. In contrast, the recognition problem has been shown to be N P-complete [START_REF] Kaufmann | Maxtolerance graphs as intersection graphs: Cliques, cycles and recognition[END_REF] and little is known concerning the complexities of other classical combinatorial optimization problems in this class.

Another generalization relevant to this work is represented by a variant of tolerance graphs, called Max-Point-Tolerance (MPT) graphs, recently introduced in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and characterized in [START_REF] Catanzaro | Max point-tolerance graphs[END_REF]. Specifically, a graph G = (V, E) is a MPT graph if each vertex v ∈ V can be mapped to a point-interval, i.e., a pair (I v , p v ) constituted by an interval 1), in such a way that (u, v) is an edge of G if and only if

I v = (l v , r v ) of R and a point p v ∈ I v (see Figure
I u ∩ I v = (l u , r u ) ∩ (l v , r v ) ⊇ {p u , p v }.
In other words, G is a MPT graph if each pair of intervals can "tolerate" a nonempty intersection (without forming an edge) as long as their distinct points are not contained in this intersection. As an example, Figure 2 shows the MPT graph associated to the point-intervals shown in Figure 1.

As for interval graphs, classical combinatorial optimization problems over MPT graphs can be solved in polynomial time, including the maximum clique problem [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and the maximum weight independent set problem [START_REF] Catanzaro | Max point-tolerance graphs[END_REF]. In contrast, the coloring problem and the minimum cost clique cover problem [START_REF] Garey | Computers and Intractability: A guide to the theory of NP-Completeness[END_REF] have been shown to be N P-hard [START_REF] Catanzaro | Max point-tolerance graphs[END_REF]. The latter problem is particularly relevant to this article due to its remarkable importance in practical applications. Specifically, MPT graphs naturally arise in genome-wide association studies as a way to represent relationships of loss of heterozygosity among DNA fragments extracted from a patient cohort affected by a common genetic disease [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF][START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF][START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF][START_REF] Stefansson | Large recurrent microdeletions associated with schizophrenia[END_REF]. In such applications, an interval I represents the maximal boundary on a chromosome region from a patient that may carry a loss of heterozygosity; the point p ∈ I instead represents a site in the considered region that shows evidence for a loss of heterozygosity [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF][START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF]. Catanzaro et al [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] showed that searching across the genomes extracted from the patient cohort for massive shared loss of heterozygosity events that can be associated with a given genetic disease is equivalent to solving a particular version of the minimum cost clique cover problem stated as follows:

C * = arg min C⊆C(G) C∈C f (C)
where C(G) is the set of all maximal cliques and singletons in G. To simplify the notation, we define f (v) = f ({v}), for all v ∈ V .

The N P-hardness of the PLOHP as well as the need to solve practical instances of the problem have justified in recent years the development of enumerative search methods and approximation algorithms such as those described in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF][START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF]. In particular, Halldorsson et al [START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF] proposed both an exact algorithm and a greedy constructive heuristic to solve the problem. The exact algorithm is based on a brute force enumeration of all of the possible clique covers of G and proved to be unpractical even for small instances of the PLOHP [START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF]. The constructive heuristic instead builds a solution to the PLOHP by finding first the maximal cliques of G and by greedily combining them to obtain a cover of G. As a drawback, the approximate approach may largely overestimate the loss of heterozygosity events across the considered genomes and lead to wrong associations [START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF]. Catanzaro et al [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] investigated possible ways to overcome the performances of the exact algorithm presented in [START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF]. In particular, the authors first proved the N P-hardness of the PLOHP and investigated some fundamental properties of MPT graphs. Subsequently, the authors proposed an Integer Linear Programming (ILP) formulation and a number of valid inequalities able to optimally solve instances containing up to 250,000 point-intervals within 5 hours of computing time. In this article, we investigate ways to improve upon the performances of the ILP formulation described in Catanzaro et al [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] s well as techniques to solve instances having size much larger than those solvable via that formulation. To this end, in Section 3 we revisit the ILP formulation described in Catanzaro et al [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and we develop a new one based on inverse projection [START_REF] Martin | Large Scale Linear and Integer Optimization: A Unified Approach[END_REF] and solvable by column generation approaches and Branch&Price techniques. In particular, we show that the pricing oracle reduces to finding a maximum node-weighted cliques in a MPT graph and we develop a polynomial time solution algorithm for this purpose. In Section 2, we introduce a number of preprocessing and decomposition strategies to dramatically reduce the size of a given instance of the problem and to divide a reduced instance into a family of independent subproblems. In Section 4, we measure the performance of the proposed algorithm on a set of artificial and biological instances of the PLOHP and we compare them versus the current state of the art presented in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF][START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF]. The results show that proposed algorithm is on average 10-30x faster than the ILP formulation described in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and enable the solution of very large practical instances of the PLOHP containing over half million point-intervals within 1 hour computing time.

Preprocessing

In this section we develop a number of preprocessing techniques both to reduce the size of an instance of the PLOHP and to decompose an instance into a family of independent subproblems having a smaller size than the original one. We will discuss the computational impact of these techniques in Section 4.

Removing Vertices

A natural approach to reduce the size of an instance of the PLOHP consists of removing point-intervals inducing isolated vertices in the MPT graph. In fact, any feasible solution to the PLOHP is characterized by identifying as singleton any isolated vertex of the MPT graph. The isolated vertices can be efficiently detected by inspecting the set I for intervals having no-compatibility between each other. For example, the point-interval I v induces an isolated vertex in V if there does not exists in I a distinct point-interval

I u such that (l u ≤ p v ≤ r u ) ∧ (l v ≤ p u ≤ r v ).
It is easy to see that the removal of the isolated vertices from an instance of the PLOHP can be performed in O(|I|(|I| -1)/2).

It is worth noting that the isolated vertices are not the only vertices that can be removed from an instance of the PLOHP. For example, given a point-interval

I v = [(l v , r v ), p v ] ∈ I,
consider the situation in which there exists another interval

I u = [(l u , r u ), p u ] ∈ I such that l u = l v , p u = p v , r u = r v .
In such a case, we say that I u is a cloned interval of I v . Let Γ I v be the set of cloned intervals of I v . Then, it is easy to see that the set of vertices induced by Γ I v can be removed from V as I v is representative of the induced equivalence class Γ I v . The identification of the cloned intervals can be performed in time O(|I|(|I| -1)/2) and can be carried out together with the removal of the isolated vertices.

Decomposing a MPT graph into connected components

The particular nature of point-intervals may give rise to MPT graphs with several connected components. For example, the point-intervals shown in Figure 1 leads to Fig. 3: An example of block decomposition of the point-intervals shown in Figure 1.

to a MPT graph with by two connected components (see Figure 2). The potential existence of connected components in a MPT graph can be exploited both to reduce the time necessary to identify a maximal clique in G as well as to decompose the problem into independent subproblems of smaller size. A possible way to quickly identify connected components in a MPT graph consists of ordering the pointintervals according to a specific precedence rule. In particular, given two distinct point-intervals I u , I v ∈ I, we say that I v precedes I u , in symbols I v ≺ I u , if one of the following cases is verified:

1. l v < l u ; 2. if l v = l u then p v < p u ; 3. if l v = l u and p v = p u then r v ≤ r u .
For example, Figure 3 shows the result of the application of the precedence rule on the point-intervals shown in Figure 1. It is easy to realize that, if the point-intervals have been sorted according to the above precedence rule, a connected component arises in G whenever there exist two I u = {(l u , r u ), p u } and

I v = {(l v , r v ), p v } such that p v ≤ l u .
The set of connected components can be identified either by an exact approach such as that described in [START_REF] Cormen | Introduction to Algorithms[END_REF] or by heuristics. As the exact approach is time consuming for very large graphs and not particularly beneficial in terms of identification of connected components, we opted for an heuristic approach able to exploit the particular structure of a MPT graph. In particular, the heuristic provided in Algorithm 1 runs in O(|I|) and exploits the sufficient condition p i < l i+1 to identify a set of nodes not connected with a given node. The algorithm takes as input the Algorithm 1 -Connected Components Finder

Input: sorted I Output: K: set of connected components of G 1: i = 0, K = ∅, K = ∅ 2: while i < |I| do 3: if p i > l i+1 then 4: K = K ∪ {I i } 5: else 6: K = K ∪ {K}, K = ∅ 7: return K
set I sorted according to the above precedence rule. Subsequently, the algorithm visits each point-interval in I and constructs a connected component B of G by collecting all of the point-intervals I i such that I i does not satisfy the condition p i ≤ l i+1 . The decomposition of G into connected components enables running in parallel several Branch&Price algorithms; this approach may vastly shorten the solution time of the overall solution algorithm for the PLOHP.

A Branch&Price algorithm for the PLOHP

In this section we describe an exact solution algorithm for the PLOHP based on Branch&Price. Throughout the section we assume, without loss of generality, that the cost function f in the PLOHP is a size-defined submodular set function, i.e., it is such that f (S) = ψ(|S|), for any subset S ⊆ V and for some ψ : [0 . . . |V |] → R + 0 . This characterization of f is commonly assumed in loss of heterozygosity studies [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and allows to solve the PLOHP in polynomial time when the MPT graph is reducible to an interval-graph [START_REF] Dion | Clique partitioning of interval graphs with submodular costs on the cliques[END_REF].

3.1 An integer linear programming formulation for the PLOHP Consider a set I of point-intervals and let G = (V, E) be the MPT graph induced by I. Then, the following proposition holds [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF]:

Proposition. A MPT graph contains at most |I|(|I| -1)/2 maximal cliques.
As the number of maximal cliques in a MPT graph is polynomially bounded and a polynomial time algorithm to determine all of them exists (see [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF]), a polynomial sized formulation for the PLOHP can be obtained as follows. Consider a vertex v ∈ V and let

C v = {C ∈ C(G) : v ∈ C}. Let x v be a decision variable equal to 1 if vertex v ∈ V forms
a singleton in a solution to the problem and 0 otherwise. Similarly, let y C be a decision variable equal to 1 if the maximal clique C ∈ C(G), with |C| ≥ 2, is selected in a solution to the problem and 0 otherwise. Then, a valid ILP formulation for the PLHOP is:

Catanzaro-Labbé-Halldorsson' ILP Formulation (CLHF). z BIP = min C∈C(G) |C|≥2 f (C)y C + v∈V f (v)x v (1a) s.t. C∈C v (G) |C|≥2 y C + x v ≥ 1 ∀ v ∈ V (1b) x v ∈ {0, 1} ∀ v ∈ V (1c) y C ∈ {0, 1} ∀ C ∈ C(G), |C| ≥ 2. ( 1d
)
The objective function (1a) minimizes the cost of finding a clique cover of G, and constraints (1b) impose that, for each vertex v ∈ V , the solution contains either a clique or a singleton in C(G) that covers v. Note that, as the overall number of maximal cliques in a MPT graph is polynomially bounded, the CLHF is polynomial-sized. However, when the instances of the PLOHP are very large (e.g., they include more than 250,000 point-intervals), the use of the CLHF may become unpractical due to the large number of y variables. A possible approach to overcome such a limitation consists of reformulating the CLHF by inverse projection [START_REF] Martin | Large Scale Linear and Integer Optimization: A Unified Approach[END_REF]. To this end, denote

C k (G) as the subset of C(G) that includes k ≤ |C(G)| maximal cliques in G and set C k v (G) = {C ∈ C k (G) : v ∈ C}. Now,
z RMP = min C∈C k (G) |C|≥2 f (C)y C + v∈V f (v)x v (2a) s.t. C∈C k v (G) |C|≥2 y C + x v ≥ 1 ∀ v ∈ V (2b) x v ≥ 0 ∀ v ∈ V (2c) y C ≥ 0 ∀ C ∈ C k (G), |C| ≥ 2. ( 2d
)
Note that in the RMP constraints x v ≤ 1, for all v ∈ V , and y C ≤ 1, for all C ∈ C k (G), |C| ≥ 2, have been omitted as they are redundant. The RMP can be solved by standard linear programming techniques. To this end, denote µ v as the dual variables associated with constraints (2b). The dual of the RMP is characterized by the following dual constraints:

µ v ≤ f (v) ∀ v ∈ V (3a)
v∈V :v∈C

µ v ≤ f (C) ∀C ∈ C k (G) : |C| ≥ 2. (3b)
A variable with negative reduced cost to be added to the RMP corresponds to a dual constraint violated by the current optimal solution. Since all of the x v variables are present in the RMP, constraints (3a) will never be violated. In contrast, constraints (3b) are violated if

∃ Ĉ ∈ C(G) : | Ĉ| ≥ 2 and v∈V :v∈ Ĉ µ v > f ( Ĉ). (4)
Checking whether (4) holds in the current optimal solution to the RMP involves solving a maximum node-weighted clique problem in G with weights {µ}. This task can be performed e.g., by means of Algorithm 3 that is described in Section 3.2. Algorithm 2 implements the column generation technique by using Algorithm 3 as pricing oracle. Specifically, Algorithm 2 takes G as an input and iterates the following steps: line 2 solves the RMP; line 3 gets the dual values {µ} and associates a weight µ v to each vertex v ∈ V ; line 5 calls the oracle and searches in C(G) for a maximum node-weighted clique Ĉ satisfying (4); if such a clique exists, line 7 adds variable y Ĉ (and its corresponding column) to the RMP and iterates lines 1-8; if, there are no violated constraints (2b), the solution to the RMP is provably optimal and Algorithm 2 returns the optimal value to the RMP. We observed that adding many cliques C having an overall weight larger than f (C) together with the most violated node-weighted clique at each iteration of the pricing oracle may

Algorithm 2 -Column Generation

Input: G Output: Optimal solution to the RMP 1: repeat 2: Solve the RMP 3:

Get the dual vector µ associated to (2b) 4:

Set weight µv on each vertex v ∈ V 5:

Ĉ =PricingOracle(G,µ)

6: if Ĉ = ∅ then 7:
Add y Ĉ and its column to the RMP 8: until Ĉ = ∅ 9: return z * RMP decrease the running time of Algorithm 2. Section 4 explores, via computational experiments, the problem of deciding how many cliques should be added in a generic instance of the PLOPH.

Finding a maximum node-weighted clique in a MPT graph

A critical step in Algorithm 2 consists of finding a maximum node-weighted clique in G. Although this problem is N P-hard for general graphs, it can be solved in polynomial time when the graph is MPT and the objective function is sizedefined submodular [START_REF] Catanzaro | The pure parsimony haplotyping problem: overview and computational advances[END_REF]. In fact, it is worth noting that, as G is a MPT, any subset of vertices that forms a maximal clique C in G is characterized by having two vertices whose corresponding points, say p l and p r , are the furthest to the left and to the right, respectively, in the set of point-intervals induced by C. Each vertex v in the clique is connected to these two vertices and its corresponding point-interval is such that

p v ∈ [p l , p r ] ⊆ [l v , r v ].
Hence, a clique can be defined by the leftmost and rightmost vertices, respectively. Then, a possible approach to enumerate the polynomial bounded number of cliques in G consists of picking any pair of adjacent vertices in G (i.e., distinct point-intervals in I, say (I u , p u ) and (I v , p v ), such that

I u ∩ I v = (l u , r u ) ∩ (l v , r v ) ⊇ {p u , p v })
and searching for the largest set of points that fall inside [p i = min{p u , p v }, p j = max{p u , p v }], i.e., the largest set of pointintervals I k that satisfy condition (p i ≤ p k ≤ p j ) ∧ (l k < p i ) ∧ (r k > p j ). Algorithm 3 outlines the pseudo code necessary to perform the overall task. We note that the complexity of Algorithm 3 is O(|I| 3 ).

Initializing the columns of the RMP

In order to initialize the RMP we first search for a subset C k (G) ⊂ C(G) covering all of the vertices in G and such that any clique in C k (G) has cardinality greater than or equal to 2. A possible approach to carry out this task is outlined in Algorithm 4.

In particular, the algorithm takes it as input the set I in the non-decreasing order described in Section 2.2 to rapidly identify the maximal cliques. Subsequently it iterates the following procedure: for each point-interval I k ∈ I (i) find a maximal clique C that cover (the vertex induced by) I k and (ii) remove from I the pointintervals (whose induced vertices are) in C. The procedure stops when I = ∅. The complexity of Algorithm 4 is O(|I| 3 ).

Algorithm 3 -Pricing Oracle

Input: G, µ Output: C : a maximum node-weighted clique in G 1: set C = ∅, max = 0, get I from G 2: for all I i , I j ∈ I : i < j and adjacent(I i , I j ) do 3: set cost = 0 4: for all I k ∈ I do 5:

if

(p i ≤ p k ≤ p j ) ∧ (l k < p i ) ∧ (r k > p j ) then 6: C = C ∪ {p k } 7: for all v ∈ V : v ∈ C do 8: cost = (cost + µv) 9:
if cost > max then 10: max = cost 11:

C = C 12: return C

Branching rules

In order to find the optimal solution to the CLHF, we embody the previously described column generation approach in a Branch-&-Bound algorithm, in which the branching is performed on both x and y variables. Specifically, concerning the x variables, the algorithm branches according to the most fractional variable rule (see [START_REF] Nemhauser | Optimization[END_REF]). The same rule also applies to the y variables with the exception that a y variable that has been set to 0 in a previous node of the search tree cannot be set to 1 in one of its descendants. To ensure the application of this rule we perform the following steps: (1) we store all the maximal clique generated in each node; (2) during the column generation process, we filter, from among all of the maximal cliques, the ones associated to the y variables that have been set to 0 in the parent node; (3) we propagate the filtered set to all the descendant nodes.

Computational experiments

In this section we analyze the performance of the Branch&Price algorithm to solve instances of the parsimonious loss of heterozygosity problem. Our experiments

Algorithm 4 -C k (G)-Generator Input: Set of intervals I Output: C k (G) 1: C = ∅, Q = I 2:
for all I i , I j ∈ I : i < j and adjacent(I i , I j ) and Q = ∅ do 3:

C = ∅ 4:
for all I k ∈ I do 5:

if compatible(I i , I j , I k ) then 6: were motivated by a number of goals, namely: to compare the performance of the Branch&Price algorithm with the ones obtained by the CLHF [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF], which currently is the best exact algorithm for the PLOHP; to evaluate the benefits obtained by using the presolving and block decomposition strategies previously described; and finally, to allow the analysis of larger data sets with respect to the ones currently analyzed. As in Catanzaro et al [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and Halldorsson et al [START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF], we emphasize that our experiments simply aim to evaluate the computational performance of the exact algorithms. We neither attempt to study the efficiency of the Branch&Price algorithm to predict LOH events across the genomes of a population of individuals nor to compare its accuracy versus LOH predictors that do not use the parsimony criterion. This analysis has been performed by [START_REF] Conrad | A high-resolution survey of deletion polymorphism in the human genome[END_REF][START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF][START_REF] Speed | Identification of Deletion Polymorphisms from Haplotypes[END_REF] and [START_REF] Mccarroll | Integrated detection and population-genetic analysis of snps and copy number variation[END_REF] and we refer the interested reader to these articles for further information.

Q = Q \ I k 7: C = C ∪ I k 8: C = C ∪ C 9:

Implementation

We implemented the CLHF in ANSI C++ by using FICO Xpress 7.6, Optimizer libraries v26.01.04 (64-bit Hyper capacity). This implementation complies with the methodology described in Catanzaro et al [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF]. We implemented the Branch&Price algorithm in ANSI C++ by using SCIP Optimization Suite 3.1.0 [START_REF] Achterberg | Scip: solving constraint integer programs[END_REF] to handle the column generation and the Branch&Price routines. Moreover, we used FICO Xpress Optimizer as linear programming solver. The experiments have been performed on an Intel Core i7-4930K CPU, 3.40GHz, equipped with 64 GByte RAM and operating system Ubuntu release 12.10 (kernel Linux 3.5.0-41-generic). During the runtime of the CLHF we activated Xpress automatic cuts, Xpress presolving strategy, and Xpress primal heuristic to generate the first upper bound for the PLOPH. Finally, similar to Halldórsson et al. [START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF] and Catanzaro et al.

[5], we set the maximum runtime for the Formulations to 5 hours. The executable codes used in the experiments can be downloaded together with the instances at http://homepages.ulb.ac.be/ ~lporrett/PLOHP/Code.tar.bz2.

Calibrating the pricing oracle

Determining the number of variables to be added in the RMP at each iteration is crucial for the performances of the proposed Branch&Price. To perform this task, we generated a learning set of instances of the PLOHP constituted by 100 random instances. We generated such a set by using the same approach and parameters discussed in [START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF]. Each random instance has between 35725 and 66670 point-intervals and between 16915 and 52522 cliques on the real segment [START_REF] Achterberg | Scip: solving constraint integer programs[END_REF]1000]. Table 1 reports on the main characteristics of this set. We introduced a parametrized number of violated cliques to be added by the pricing oracle at runtime. Such a number is a fraction of the overall number of maximal cliques contained in a random instance and we set it to be equal to one of the following seven different percentages: 1%, 3%, 5%, 7%, 10%, 15% and 20%. Figure 4 shows the box-and-whisker plot histogram of the solution time taken by the Branch&Price algorithm to exactly solve the learning set of instances of the PLOHP when considering different percentages of violated cliques. A box shows the range between the 25% and the 75% quantile of the data. The median of the data is indicate by a bar. The whiskers extend to the most extreme data point which is no more than 1.5 times the interquartile range from the box. The figure shows that the median of the solution time taken by the Branch&Price algorithm to exactly solve the learning set of instances of the PLOHP is minimum when adding in the RMP no more than 1% of violated cliques. To exclude the presence of statistical equivalence between the contained in a given instance. The last three columns report on the number of Connected components, cloned intervals and singletons contained in a given instance. Table 3 shows that the instances "gens-x α ", x ∈ {1, 5} are characterized by a number of connected components ranging from 29 to 37, a number of singletons ranging from 0 to 5 and over 70% of cloned point-intervals. The instances "gensx β ", x ∈ {1, 5} are characterized by a number of connected components ranging from 33 to 40, a number of singletons ranging from 3 to 7 and over 70% of cloned point-intervals. Finally, the instances "gens-x γ ", x ∈ {1, 5} are characterized by a number of connected components ranging from 40 to 61, a number of singletons ranging from 2 to 13 and over 80% of cloned point-intervals. We also considered a set of 5 artificial datasets of the PLOHP, hereafter denoted as "Sim-1", "Sim-2", "Sim-3", "Sim-4" and "Sim-5", each containing 10 random instances of the problem and mainly differing from one another both by the length and the number of the point-intervals contained in each of the corresponding instances. We generated the instances in each dataset by first fixing the real segment [START_REF] Achterberg | Scip: solving constraint integer programs[END_REF]6000] as a reference and by generating, by means of the Mersenne Twister library [START_REF] Matsumoto | Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator[END_REF], a random number of triplets 4 summarizes the characteristics of the considered datasets. Finally, in order to evaluate the performance of the Branch&Price algorithm also on biological instances of the PLOHP, we considered a biological instance provided by deCode Genetics and hereafter called "Bio22". This instance is constituted by 20022 point-intervals on the real segment [START_REF] Achterberg | Scip: solving constraint integer programs[END_REF]6000], each having minimum length 2 and maximum length 6000. The pointintervals in "Bio22" refers to genetic fragments from chromosome 22 extracted from a population of 18360 individuals and consists of 6000 single nucleotide polymorphisms. All of the instances used in our experiments can be downloaded at http://homepages.ulb.ac.be/ ~lporrett/PLOHP/Datasets.tar.bz2.

(l v , p v , r v ) such that 1 ≤ l v ≤ p v ≤ r v ≤ 6000. Table

Computational performances

Table 5 shows the performances of the CLHF and the Branch&Price in solving the instances derived from the 5 real instances described in [START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF] with and without the presolving strategies. In particular, the table is constituted by three main sections, each one corresponding to the datasets "gens-x α ", "gens-x β " and "gens-x γ ", x ∈ {1, 5}, previously described. Moreover, the columns are divided in four sections, each of them showing the solution time (expressed in seconds) and the number of cliques needed to optimally solve a given instance of the PLOHP. Specifically, Table 5: Computational performances comparison of the CLHF and the Branch&Price algorithm with and without presolving strategies the first pair of columns is related to the execution of the CLHF. The second pair of columns is related to the execution of the Branch&Price. The third pair of columns is related to the execution of the CLHF using the presolving strategies. Finally, the fourth pair of columns are related to the execution of the Branch&Price using the presolving strategies. Table 5 shows the benefits of using the presolving strategies. Specifically, we have observed the removal of at least 70% of cloned intervals (hence x v variables in the CLHF) in each of the considered instances. Moreover, the detection of the connected components in the considered instances enabled the decomposition of the corresponding induced MPT graphs into at least 29 connected components (i.e., independent subproblems); the combination of the two strategies allowed a speed up of about 18x of the solution times taken by the CLHF. Table 5 also shows that, independently of the presence (or absence) of the presolving strategies, the Branch&Price is always faster than the CLHF. This fact alone justifies and shows the importance of using column generation when dealing with very large instances of the PLOHP. In this context, it is worth noting that, although the Branch&Price is faster than the CLHF, it also proves unable to solve the instances "gens-x γ ", x ∈ {1, 5}, within the considered time limit. This fact is mainly due to the pricing oracle, which proves particularly time-consuming in gens-x γ .

4.5 Computational performances on larger datasets of the PLOHP Table 6 shows the performance of the CLHF and the Branch&Price in solving the "Sim" instances when using the presolving strategies. Both the CLHF and the Branch&Price algorithm without presolving strategies proved unable to solve any instance in this set. The table reports on the mean and standard deviation of (i) the solution time (expressed in seconds) necessary to solve a random instance of the PLOHP; (ii) the gap (expressed in percentage), i.e., the difference between the optimal value to a given instance of the PLOHP in a specific dataset and the objective function value of the linear programming relaxation at the root node of the respective search tree, divided by the optimal value; and (iii) the nodes explored in the search tree. In general, Table 6 shows that the CLHF and the Branch&Price, combined with presolving strategies, is able to solve each simulated instance in about half an hour, except for the parameter set "Sim-5" that required on average 26 minutes for the Branch&Price while the CLHF exceeds the 64GB of memory available for the experiments. The table also shows that the solution times necessary to solve the artificial instances of the PLOHP are longer than the corresponding ones necessary to solve any other instance considered in this article. This fact is possibly due to the larger number of point-intervals present in the artificial instances which in turn increases the number of cliques in the corresponding MPT graphs.

Table 7 reports on the minimum, maximum and average number of connected components, intervals, singletons, edges and cloned nodes in the considered simulated instances. Table 8 reports on the minimum, maximum and average time (expressed in seconds) to filter both the singletons and the cloned intervals as well as to identify the connected components. Table 9 shows the performance of the Branch&Price in solving "Bio22" instance using the presolving strategies. The table reports on (i) the solution time (expressed in seconds) necessary to solve the biological instance of the PLOHP; (ii) the gap (expressed in percentage); (iii) the nodes explored in the search tree; and (iv) the overall number of connected components, intervals, cloned intervals and singletons, respectively.

Table 9 shows that the proposed Branch&Price, combined with presolving strategies, is able to solve "Bio22" in 664.40 seconds. This instance is characterized by having a number of point-intervals smaller than those relative to the instances considered in Table 3 as well as by a smaller number of cloned ones. Filtering the cloned point-intervals in "Bio22" reduces the number of x v variables by the 20% as opposed to 70% in Table 5. However, the relative higher number of connected components in "Bio22" allows to decompose the problem into a bigger family of small independent subproblems; this fact justifies the remarkable speed up in terms of solution time.

Conclusion

The Parsimonious Loss of heterozygosity Problem (PLOHP) is a N P-hard combinatorial optimization problem consisting of solving the minimum cost clique cover problem on a MPT graph. The optimal solutions to the PLOHP have a remarkable importance in practice, as they enable the association of major human diseases with chromosome regions from patients that exhibit loss of heterozygosity events. As genome-wide association studies usually involve the analysis of massive amount Table 8: The minimum, the maximum and the average time (expressed in seconds) both to filter the singletons and the cloned intervals, and to identify the connected components in Sim1, Sim2, Sim3, Sim4, Sim5.

of data, practical instances of the PLOHP may have very large sizes. This fact precludes the use of exact solution approaches to the PLOHP such as those described in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF][START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF]. In this article, we have investigated ways (i) to overcome this limit, (ii) to speed up the solution time of the best known exact solution algorithm for the PLOHP, and (iii) to enlarge the size of the instances that can be optimally solved within the same time limit. We developed a number of preprocessing strategies able to both dramatically reduce the size of the instances analyzed and decompose the reduced instances into a collection of small independent subproblems. These strategies allowed to speed up of about 8 times the solution times of the best known exact solution algorithm for the PLOHP as well as enabled the resolution, within the same time limit, of instances previously too large to be tackled. Moreover, we combined such preprocessing strategies with a Branch&Price algorithm based on column generation techniques able to exploit the combinatorial properties of the MPT graphs. overall performances of the proposed algorithm proved to be up to 30 times faster than previous approaches described in the literature; such performances could be further enhanced, by parallelizing the way in which the Branch&Price algorithm handle the collection of small independent subproblems provided by the preprocessing strategies. This result helps to illustrate the power of ILP approaches to tackle problems arising from genome-wide association studies and characterized by complex sets of constraints. Hopefully, future research efforts will provide new insights on the combinatorics of MPT graphs and will enable the analysis of even larger genomic fragments.

A

  Branch&Price Algorithm for the Minimum Cost Clique Cover Problem in Max-Point Tolerance Graphs L. Porretta, • D. Catanzaro, • B. V. Halldórsson, • B. Fortz

Fig. 1 :

 1 Fig. 1: An example of point-intervals Fig. 2: The MPT graph induced by the point-intervals shown in Figure 1.

  consider the following Restricted Master Problem (RMP) associated with the linear programming relaxation of the CLHF: the CLHF -The Restricted Master Problem (RMP).

Fig. 4 :

 4 Fig. 4: BoxPlot of the solution time taken by the Branch&Price algorithm to exactly solve the learning set of instances of the PLOHP when considering different percentages of violated cliques.

Table 1 :

 1 Summary of number of nodes, edges and cliques contained in the learning set of instances of the PLOHP used to calibrate the pricing oracle.

		return C
	10: function compatible(Iu, Iv, Iz)
	11:	return (pu ≤ pz ≤ pv) ∧ (lz < pu) ∧ (rz > pv)

Table 2 :

 2 Wilcoxon signed-rank test on solution time for each pair of percentages

	Set of	Point-interval	Max. number of Min. number of Average number of
	Instances	Length	Point-intervals	Point-intervals	Point-intervals
	Sim-1	2	23963	23777	23851 ±52.79
	Sim-2	5	24012	23823	23950.5±58.03
	Sim-4	7	23615	23454	23553.5±50.67
	Sim-3	9	23777	23803	23828 ±42.55
	Sim-5	random in [1,9]	29020	28699	28904 ±94.11

Table 4 :

 4 Parameter sets used to generate random instances of the PLOHP.

Table 6 :

 6 Computational performances comparison of the CLHF+Preprocessing and the Branch&Price algorithm + Preprocessing on artificial instances of the PLOHP.

			CLHF			Branch&Price
	Instance	Time	Gap	Nodes	Time	Gap	Nodes
		(sec.)	(%)		(sec.)	(%)
	Sim-1	1 217.96±41.14	0.017±0.18	1.012±0.20	577.77±18.24	0.003±0.13	1.003±0.12
	Sim-2	1 200.47±39.02	0.009±0.17	1.010±0.18	552.44±40.46	0.004±0.12	1.004±0.11
	Sim-3	1 200.24±48.61	0.007±0.13	1.010±0.19	575.27±32.15	0.001±0.05	1.001±0.06
	Sim-4	1 162.67±94.31	0.007±0.14	1.012±0.26	509.01±84.75	0.003±0.09	1.003±0.08
	Sim-5	OOM±OOM	OOM±OOM OOM±OOM	969.56±424.29 0.002±0.09	1.001±0.06

Table 7 :

 7 Statistics of simulated instances Sim1, Sim2, Sim3, Sim4, Sim5.

			Sim1			Sim2
		Min	Max	Average	Min	Max	Average
	Singleton	0.43	0.52	0.44	0.43	0.52	0.43
	Cloned Inter.	0.42	0.43	0.42	0.42	0.50	0.42
	Conn. Comp.	6.2•10 -5	1.3•10 -5	6.2•10 -5	6.1•10 -5	1.7•10 -5	6.2•10 -5
			Sim3			Sim4
		Min	Max	Average	Min	Max	Average
	Singleton	0.43	0.52	0.43	0.41	0.42	0.41
	Cloned Inter.	0.42	0.50	0.42	0.40	0.41	0.42
	Conn. Comp.	6.1•10 -5	6.2•10 -5	6.2•10 -5	6.0•10 -5	1.3•10 -5	6.1•10 -5
			Sim5
		Min	Max	Average
	Singleton	0.73	0.78	0.74
	Cloned Inter.	0.65	0.66	0.65
	Conn. Comp.	7.7•10 -5	7.9•10 -5	7.8•10 -5

Acknowledgements Part of this work has been developed when L. Porretta and D. Catanzaro were visiting both Reykjavik University and DeCODE Genetics. L. Porretta acknowledges support from the Belgian National Fund for Scientific Research (FRS-FNRS) who funded his short-term research visit. D. Catanzaro acknowledges support from the D. Catanzaro acknowledges support from the FRS-FNRS via the grant "Crédit de Recherche" ref. S/25-MCF/OL J.0026.17, the Université Catholique de Louvain via the "Fonds Spéciaux de Recherche" (FSR) 2017-2019, the Fondation Louvain via the grant "Le numérique au service de l'humain", and the computational resources provided by the "Consortium des Équipements de Calcul Intensif" (C ÉCI), funded by the FRS-FNRS via the grant ref. 2.5020.11. The authors thank the anonymous referees, whose insightful suggestions helped to improve the manuscript.

 [START_REF] Fay | Wilcoxon-mann-whitney or t-test? On assumptions for hypothesis tests and multiple interpretations of decision rules[END_REF]. Table 2 reports on the p-values so obtained. Specifically, the presence of a value smaller then 0.05 in a cell means that the sample in the row is statistically significant smaller than the sample in the column. In contrast, the presence of a value bigger than 0.05 in a cell means that the sample in the row is statistically significant bigger than the sample in the column. In particular, Table 2 shows that adding the 1% of the violated cliques always provide a smaller solution time compared to the others approaches. For this reason, we decided to use this setup in applying the Branch&Price algorithm to the remaining experiments.

Benchmark instances

In order to compare the performance of the Branch&Price algorithm with the ones obtained by the CLHF [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF], we considered the instances of the PLOHP described in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] and [START_REF] Halldorsson | The Clark phaseable sample size problem: Long-range phasing and loss of heterozygosity in GWASlark phaseable sample size problem: long-range phasing and loss of heterozygosity in GWAS[END_REF]. These instances, hereafter denoted as "gens-1", "gens-2", "gens-3", "gens-4" and "gens-5", are characterized by having 645572, 645552, 645510, 645543 and 659885 point-intervals on the real segment [START_REF] Achterberg | Scip: solving constraint integer programs[END_REF]3575], respectively, each having minimum length 1 and maximum length 9. Catanzaro et al observed that the ILP formulation presented in [START_REF] Catanzaro | An integer programming formulation of the parsimonious loss of heterozygosity problem[END_REF] proved unable to solve any of these instances within 1 hour of computing time. Hence, the authors also considered smaller instances of the PLOHP obtained from the given ones by extracting from each "gens-x", x ∈ {1, 5}, the point-intervals contained in the real segments α = [1, 1500], β = [1, 2000] and γ = [1, 3575], respectively. The authors obtained this way 15 instances of the PLOHP, that are hereafter grouped into three datasets denoted as "gens-x α ", "gens-x β " and "gens-x γ ", x ∈ {1, 5}, respectively. Table 3 shows the main characteristics of the 15 instances. In particular, the table is composed by three main sections, each one corresponding to the datasets "gens-x α ", "gens-x β " and "gensx γ ", x ∈ {1, 5}, previously described. The first three columns of