Etienne Deborde
email: etienne.deborde@chru-strasbourg.fr

Benjamin Dubourg

Sofiane Bejar

Anne-Claire Brehin

Sébastien Normant

Paul Michelin

Jean-Nicolas Dacher

Differentiation between Fabry disease and hypertrophic cardiomyopathy with cardiac T1 mapping

Keywords: Cardiac magnetic resonance imaging (MRI), Fabry disease, Left ventricular hypertrophy, Hypertrophic cardiomyopathy Fabry disease LVH: left ventricular hypertrophy ERT: enzyme replacement therapy CMR: cardiac magnetic resonance LVM: left ventricular mass LGE: late gadolinium enhancement LV: left ventricle MRI: magnetic resonance imaging HCM: hypertrophic cardiomyopathy CI: confidence interval SD: standard deviation

Purpose.

To evaluate the potential of non-contrast myocardial T1 mapping on cardiovascular magnetic resonance examination (CMR) in differentiating patients with Fabry disease (FD) from those with hypertrophic cardiomyopathy (HCM) and healthy control subjects.

Materials and methods.

Seventeen patients with FD (8 men, 9 women; mean age, 48 years ± 18; [range: 19-73 years]; 53% with left ventricular hypertrophy [LVH]) were matched with 36 patients with hypertrophic cardiomyopathy (HCM) (22 men, 14 women; mean age, 57 years ± 16 [SD] [range: 22-85 years]) and 70 healthy control subjects (34 men, 36 women; mean age, 38 years ± 15 [SD]; [range: 18-65 years]). Cardiac T1 mapping was performed using the modified Look-Locker inversion (MOLLI ® ) sequence on a 1.5-T magnet. T1 values were calculated, on midventricular section, for septal left ventricular segments (S8-S9) and all mid-ventricular ones (global T1 values; S7-S12). Statistical analysis included unpaired Mann-Whitney test, receiver operating characteristic curve and likelihood ratios. Results. Septal native T1 values were significantly decreased in patients with AFD (889 ± 61 [SD] ms; range: 784-980 ms) compared to those with HCM (995 ± 48 [SD] ms; range: 935-1125 ms) (P<0.001) and versus healthy controls (965 ± 29 [SD] ms; range: 910-1028 ms) (P<0.001). Global native T1 values were also significantly decreased in patients with AFD (891 ± 49 [SD] ms; range 794-970 ms) compared to those with HCM (995 ± 34 [SD] ms; range: 952-1086 ms) (P<0.001) and versus healthy controls (966 ± 27 [SD] ms; range: 920-1042 ms) (P<0.001). A septal left ventricular native T1 cutoff value of 940 ms could distinguish AFD from HCM with 88% sensitivity (95% CI: 73-100%) and 92% specificity (95% CI: 83-100%

). Positive likelihood ratio was 11, negative likelihood ratio was 0.12.

Compared to controls, the same threshold could distinguish AFD with 88% sensitivity (95% CI: 73-100%) and 86% specificity (95% CI: 78-94%). Positive likelihood ratio was 6.3, negative likelihood ratio was 0.14. T1 value was abnormal in 4 of 8 (50%) of FD patients who did not have LVH.

Conclusion.

Native T1 values are significantly lower in patients with AFD compared with HCM and healthy volunteers.

Introduction

Fabry disease (FD) is a rare multisystem X-linked lysosomal storage disorder of glycosphingolipid metabolism due to complete or partial deficiency of the lysosomal enzyme α-galactosidase A (α-Gal A) caused by mutations in the GLA gene (located on Xq 22) [START_REF] Linhart | The heart in Anderson-Fabry disease and other lysosomal storage disorders[END_REF].

More than 600 pathogenic mutations have been identified, most of them being private, occurring in individual families [START_REF] Clarke | Narrative review: Fabry disease[END_REF]. Disease severity in patients with FD is inversely related to α-Gal A activity. As a result, progressive intracellular accumulation of glycosphingolipids occurs in multiple organs. Premature death is the consequence of cardiac failure, arrhythmia, stroke and renal failure [START_REF] Linhart | The heart in Anderson-Fabry disease and other lysosomal storage disorders[END_REF]. Progressive left ventricular hypertrophy (LVH) and myocardial fibrosis may induce systolic/diastolic dysfunction, valvular heart disease, arrhythmias and heart failure [START_REF] Clarke | Narrative review: Fabry disease[END_REF].

The incidence of FD has been estimated 1 in 40,000 to 1 in 117,000 live births for males [START_REF] Linhart | The heart in Anderson-Fabry disease and other lysosomal storage disorders[END_REF]. Screening of patient populations at risk of FD (i.e., those with cardiomyopathy [START_REF] Sachdev | Prevalence of Anderson-Fabry disease in male patients with late onset hypertrophic cardiomyopathy[END_REF][START_REF] Monserrat | Prevalence of Fabry disease in a cohort of 508 unrelated patients with hypertrophic cardiomyopathy[END_REF], stroke [START_REF] Rolfs | Acute cerebrovascular disease in the young: the stroke in young Fabry patients study[END_REF], or renal failure [START_REF] Wallin | Fabry disease: results of the first UK hemodialysis screening study[END_REF]), has identified a low but not negligible prevalence of FD in the general symptomatic population. Recently, studies in high-risk patient cohorts suggested that FD was much more common and accounted for at least 3% of unexplained LVH in middle aged men [START_REF] Linhart | The heart in Anderson-Fabry disease and other lysosomal storage disorders[END_REF][START_REF] Mehta | Fabry disease defined: baseline clinical manifestations of 366 patients in the Fabry Outcome Survey[END_REF]. The diagnosis of FD is often delayed by as much as 14 years in males and 16 years in females although the first signs of FD, such as acroparaesthesia, usually appear in childhood in affected patients.

The assessment of cardiac involvement in patients with FD is crucial because early enzyme replacement therapy (ERT) is likely to postpone the progression of LVH and development of myocardial fibrosis and maintain functional capacities [START_REF] Weidemann | Improvement of cardiac function during enzyme replacement therapy in patients with Fabry disease: a prospective strain rate imaging study[END_REF][START_REF] Weidemann | Long-term effects of enzyme replacement therapy on Fabry cardiomyopathy: evidence for a better outcome with early treatment[END_REF]. Cardiac magnetic resonance (CMR) examination is now the gold standard imaging technique to assess myocardial anatomy, regional and global function and approach tissue characterization. CMR examination has emerged as a key imaging modality to provide both quantitative and qualitative assessment of cardiomyopathies [START_REF] Karamitsos | The role of cardiovascular magnetic resonance imaging in heart failure[END_REF][START_REF] Wong | Association between extracellular matrix expansion quantified by cardiovascular magnetic resonance and short-term mortality[END_REF][START_REF] Klem | Prognostic value of routine cardiac magnetic resonance assessment of left ventricular ejection fraction and myocardial damage: an international, multicenter study[END_REF]. In the context of FD, CMR examination helps identify and quantify increased left ventricular mass (LVM) and the presence of scarring on late images obtained after intravenous administration of a gadolinium chelate (i.e., "late gadolinium enhancement" [LGE]). However, those anomalies are neither sensitive nor specific of FD [START_REF] Hughes | Effects of enzyme replacement therapy on the cardiomyopathy of Anderson-Fabry disease: a randomised, double-blind, placebo-controlled clinical trial of agalsidase alfa[END_REF]. As glycosphingolipids are supposed to decrease the T1 longitudinal relaxation time, we hypothesized that native myocardial T1 values would be significantly lower in patients with FD compared to those with hypertrophic cardiomyopathy (HCM) or healthy subjects. In this regard, Sado et al. have shown that native myocardial T1 values were substantially lower in patients with FD compared to those with other causes of LVH [START_REF] Sado | Identification and assessment of Anderson-Fabry disease by cardiovascular magnetic resonance noncontrast myocardial T1 mapping[END_REF].

The goals of this study were to evaluate the potential of non-contrast myocardial T1 mapping for differentiating patients with FD from those with HCM and healthy controls.

Materials and Methods

Patient selection

Seventeen patients with clinically and genetically confirmed FD, who underwent a CMR examination from May 2013 to March 2018 were retrospectively included. Four of them were receiving ERT. Patients with usual magnetic resonance imaging (MRI) contraindication such as pregnancy, or legal protection were excluded. Thirty-six patients with HCM were also included in the same period of time. The diagnosis of HCM was defined as a LVH (wall thickness > 15 mm or 13 mm in the context of documented familial disease) associated with a non-dilated left ventricle (LV) and required the absence of another cardiac or systemic disease [START_REF] Maron | American College of Cardiology/European Society of Cardiology clinical expert consensus document on hypertrophic cardiomyopathy: a report of the American College of Cardiology Foundation Task Force on Clinical Expert Consensus Documents and the European Society of Cardiology Committee for Practice Guidelines[END_REF]. The diagnosis of HCM was established based on international guidelines [START_REF] Gersh | ACCF/AHA guideline for the diagnosis and treatment of hypertrophic cardiomyopathy: a report of the American College of Cardiology Foundation/American Heart Association Task Force on Practice Guidelines. Developed in collaboration with the American Association for Thoracic Surgery, American Society of Echocardiography, American Society of Nuclear Cardiology, Heart Failure Society of America, Heart Rhythm Society, Society for Cardiovascular Angiography and Interventions, and Society of Thoracic Surgeons[END_REF]. Exclusion criteria included documented coronary artery disease, hypertension, aortic valvular disease, and significant renal dysfunction with estimated glomerular filtration rate < 30 mL/min/1.73m 2 and usual contraindications to MRI examination.

Concomitantly, 70 asymptomatic healthy control participants with no documented history of cardiovascular disease and normal CMR findings formed a control group. Healthy control participants were enrolled from the "Heart Map Study", a study performed in our center to determine the normal T1, T2 and T2* values of the myocardium in a cohort of 175 healthy volunteers (ClinicalTrials.gov Identifier: NCT01712581) [START_REF] Bentatou | Distribution of left ventricular trabeculation across age and gender in 140 healthy Caucasian subjects on MR imaging[END_REF]. All healthy volunteers gave written informed consent to participate in the study. The institutional review board approved this retrospective study, and informed consent was waived. Clinical and demographic data were obtained from the electronic patient records.

No patients were excluded. The final cohort consisted of 17 patients with genepositive FD; there were 8 men and 9 women (mean age, 48 ± 18 [SD] years; range: 19-73 years) (Table 1). Four patients with FD (24%) were receiving ERT. Thirty-six patients with HCM were included; there were 22 men and 14 women (mean age, 57 ± 16 [SD] years; range: 22-85 years). Finally, 70 healthy volunteers were included; there were 34 men and 36 women (mean age, 38 ± 15 [SD] years; range: 18-65 years). Figure 1 shows patients inclusion.

CMR protocol

CMR examinations were performed on a 1.5 T Avanto ® unit (Siemens Healthineers) with a cardiac 32 channel-coil for signal reception and electrocardiographic (ECG) gating.

CMR examinations included cardiac function assessment, native T1 mapping and LGE images.

Cardiac function assessment was performed using cine-MRI sequences obtained with retrospectively ECG-gated balanced steady-state-free precession (SSFP) sequences in long axis and short axis providing full coverage of the LV. The following parameters were used: TR, 36 ms; TE, 1.5 ms; flip angle, 60 °; field of view (FOV), 360x420 mm; matrix size, 216x256; slice thickness, 6 mm; interslice gap, 3mm; cardiac phases per cardiac cycle, 20-25.

Mid ventricular native T1 mapping was performed in the short axis plane (S7-12) in all patients. Modified Look-Locker inversion recovery (MOLLI ® ) technique was used to evaluate T1 values. Inversion groupings acquisition scheme 3(3)3(3)5 was applied [START_REF] Messroghli | Modified Look-Locker inversion recovery (MOLLI) for high-resolution T1 mapping of the heart[END_REF]. For patients with AFD, three short-axis sections were acquired at the base, mid-ventricle, and apex before intravenous administration of gadobutrol (Gadovist ® , Bayer Healthcare) at a dose of 0.1 mmol/kg of body weight [START_REF] Kellman | T1-and extracellular volume mapping in the heart: estimation of error maps and the influence of noise on precision[END_REF]. Imaging parameters were as follows: TR, 3.7 ms; TE, 1.1 ms; flip angle, 35°; matrix acquisition size, 256 x 128; FOV, 360 x 270 mm; slice thickness, 8 mm; TI, 75-750 ms; trigger delay, 300 ms; and views per segment, 24. Images were processed with a curve fitting technique to generate T1 maps.

LGE sequences were performed 10 minutes after the intravenous administration of gadobutrol to identify regional fibrosis using a T1-weighted inversion recovery gradient echo technique. Imaging parameters were as follows: TR, 7.1 ms; TE, 3.1 ms; slice thickness, 8 mm; matrix acquisition size, 256 x 192; with an inversion time individually determined to null the myocardial signal. LGE sequences were performed using standard long-axis and short-axis views of the LV. Healthy volunteers were not injected.

Image analysis

All CMR images were analysed using SyngoVia ® (Siemens Healthineers) by one observer (E.D.) with an experience of two years in cardiac MRI.

Evaluation of LV function and wall thickness

LV short-axis epicardial and endocardial borders were semi-automatically contoured at enddiastole and end-systole for determining end-diastolic and end-systolic volumes, LV ejection fraction and LVM. The papillary muscles were included as part of the LV cavity volume.

Maximum end-diastolic LV wall thickness was manually determined on short axis cine steady-state free precession images. LVH was considered when LVM indexed to body surface area was above 85 g/m 2 in men and 81 g/m 2 in women [START_REF] Kawel-Boehm | Normal values for cardiovascular magnetic resonance in adults and children[END_REF].

Evaluation of myocardial T1 values

For each short-axis image, and for each myocardial segment, a region of interest (ROI) was manually drawn within the LV myocardium to calculate T1 values, excluding papillary muscles. Care was taken to avoid LV cavity blood pool and potential partial volume effects at the endocardial and epicardial borders. A total of six segments per patient were analysed. T1 values comparisons were established on the short axis midventricular section as it was the only available one in healthy volunteers. Segments 8 and 9 of the American Heart Association were considered for estimating septal T1 values while global T1 values were inferred from segments 7-12. As midwall myocardial MOLLI ® T1 was shown to be reproducible, with intra-individual, intra-and inter-MRI variabilit ≤ 2%) [START_REF] Piechnik | Normal variation of magnetic resonance T1 relaxation times in the human population at 1.5 T using ShMOLLI[END_REF], we did not perform reproducibility testing. The relatively long duration of T1 mapping sequences occasionally led to some undesired breathing motion, which generated artifacts and non-evaluable segments. Less than 1% of all segments were excluded and these segments were not used for the analysis. An analysis with short axis slice used in the study is described in Figure 2.

Evaluation of LGE images

The LGE images were visually examined for the presence of regional fibrosis.

LGE was considered present only if identified in two orthogonal views and excluding RV insertion points [START_REF] Maron | Clinical utility of cardiovascular magnetic resonance in hypertrophic cardiomyopathy[END_REF]. Presence or absence of LGE was mentioned for each myocardial segment. The analysis of LGE sequences was performed independently from the T1 analyses. Fourteen patients with FD received an intravenous administration of gadolinium chelate; of these, one patient had not interpretable CMR because of artifacts so that 13 CMR examinations could be analyzed for LGE findings.

Statistical analysis

Quantitative variables were expressed as mean ± standard deviation (SD) and ranges.

Qualitative variables were expressed as raw numbers, proportions and percentages.

Comparison of continuous variables was performed using unpaired Mann-Whitney test.

Comparison of proportions was done with chi-squared (χ2) test. Receiver operating characteristic curves (ROC) were created to evaluate the ability of native T1 values to discriminate between the three groups. Areas under the receiver operating characteristic curve (AUC) for native T1 values were tested for equality [START_REF] Delong | Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach[END_REF]. Bivariate correlations were searched for using Pearson correlation test, with r < 0.3 indicating none or very weak correlation, 0.3 < r < 0.5 a weak correlation, 0.5 < r < 0.7 a moderate correlation and r > 0.

Results

Study population

A total of 123 patients were evaluated during the study period, including 17 patients with FD, 70 asymptomatic control subjects and 36 patients with HCM. The 4 patients with FD (24%) under ERT were not excluded from the analysis. Patients' characteristics are presented in Table 1.

CMR imaging findings

CMR imaging findings are shown in

Relationship between T1 values and functional parameters in patients with AFD

When considering the whole FD group, we observed a strong positive correlation between septal native T1 values and non-septal native T1 values (r = 0.85; P < 0.001). Septal native T1 values did not differ significantly between men (886 ± 59 [SD] ms) and women (893 ± 67 [SD] ms) (P = 0.61) (Fig. 5). Fifty percent of patients with FD and no LVH had T1 values under the limit of 907 ms (normal values -2SD). T1 values were weakly and inversely correlated with indexed LVM (r=-0.483; P = 0.049) and was also moderately and inversely correlated with the maximal LV wall thickness (r = -0.572; P = 0.016). Septal native T1 values did not differ significantly between patients with FD whether they were on ERT or not (respectively 872 ms ± 54 [SD] vs. 895 ms ± 64 [SD]) (P = 0.395). Normal or even elevated T1 in the LV inferolateral wall was shown in some patients with FD (Fig. 2). This variation was observed in patients with (n=5; 38 %) or without (n=8; 62 %) LGE (respectively 978 ms ± 22 [SD] vs. 917 ms ± 52 [SD]) (P = 0.009) (Table 2).

Discussion

Differentiating FD from HCM on conventional MRI remains challenging. Our study confirmed that myocardial native T1 mapping helps differentiate between patients with FD and healthy subjects and those with HCM. At 1.5-T and using MOLLI ® , a septal native T1 women and 1120 ms in men [START_REF] Thompson | T1 mapping with cardiovascular MRI is highly sensitive for Fabry disease independent of hypertrophy and sex[END_REF].

By comparison to MOLLI ® and ShMOLLI ® , Roujol et al. reported that SASHA ® and SAPPHIRE ® had better accuracy, lower precision, and similar reproducibility. They found that MOLLI ® and ShMOLLI ® led to underestimating myocardial T1 values in both phantom and healthy subject studies. These underestimations were found to be more prominent with higher T1 values [START_REF] Roujol | Accuracy, precision, and reproducibility of four T1 mapping sequences: a head-to-head comparison of MOLLI, ShMOLLI, SASHA, and SAPPHIRE[END_REF]. In our study we chose to use MOLLI ® for several reasons; first, we had experience with the technique. Secondly, we had carried out a study in a population of 175 healthy volunteers (NCT01712581) [START_REF] Bentatou | Distribution of left ventricular trabeculation across age and gender in 140 healthy Caucasian subjects on MR imaging[END_REF]. Finally, the midwall myocardial native T1 drawn from MOLLI ® sequence was shown to be reproducible with an intra-individual, intra-and inter-MRI variability ≤ 2% [START_REF] Piechnik | Normal variation of magnetic resonance T1 relaxation times in the human population at 1.5 T using ShMOLLI[END_REF].

A large-scale study including more than 1,200 subjects and using ShMOLLI ® reported that patients with FD had significantly lower native T1 (863 ± 23 [SD] ms) compared to patients with normal CMR examination or the normal reference myocardium in patients with other cardiovascular disease [START_REF] Nacif | Myocardial T1 mapping with MRI: comparison of look-locker and MOLLI sequences[END_REF]. Another recent study reported that patients with FD and LVH had reduced native T1 values and more positive circumferential strain compared to those without LVH. They showed that the basal septum could be the most appropriate region for measuring native T1 values because of a smaller standard deviation. They also showed that basal septal T1 value was superior to circumferential strain in discriminating among patients with LVH [START_REF] Liu | Measurement of myocardial native T1 in cardiovascular diseases and norm in 1291 subjects[END_REF].

In a previous publication that compared native myocardial T1 values at 3.0 T between patients with FD and those with HCM, a septal native T1 cutoff point value of 1220 ms was defined [START_REF] Karur | Use of myocardial T1 mapping at 3.0 T to differentiate Anderson-Fabry disease from hypertrophic cardiomyopathy[END_REF][START_REF] Weingärtner | Myocardial T1-mapping at 3T using saturation-recovery: reference values, precision and comparison with MOLLI[END_REF][START_REF] Dabir | Reference values for healthy human myocardium using a T1 mapping methodology: results from the International T1 Multicenter cardiovascular magnetic resonance study[END_REF]. This cutoff value allowed distinguishing patients with FD from those with HCM with high sensitivity (29 of 30 [97 %]), specificity (28 of 30 [93 %]), and accuracy (57 of 60 [95 %]). Although global T1 values had slightly better reproducibility compared with septal LV native T1 values, septal T1 values were recommended in the setting of patients with suspected FD, because of better discriminating capabilities and ease of measurements [START_REF] Karur | Use of myocardial T1 mapping at 3.0 T to differentiate Anderson-Fabry disease from hypertrophic cardiomyopathy[END_REF]. Regarding the potential of evaluating global vs. septal T1 values, our results did not bring definite conclusion although the global AUCs appeared more discriminating than septal ones. The values of T1 did not differ significantly between men and women in our study as well as in data from the literature [START_REF] Niemann | Differences in Fabry cardiomyopathy between female and male patients: consequences for diagnostic assessment[END_REF].

The decrease in myocardial T1 values in patients with FD could be explained by the accumulation of glycosphingolipids in myocytes [START_REF] Eng | Safety and efficacy of recombinant human alpha-galactosidase a replacement therapy in Fabry's disease[END_REF]. On the other hand, we noted a pseudonormalization of T1 values in the lateral-inferior basal segment in patients with a focal fibrous scar visualized on the LGE sequences. A significant percentage of patients with FD having renal dysfunction, this paradoxical increase in T1 values could make it possible to identify focal fibrosis without contrast injection. Although this could be of interest in this population at risk of gadolinium-induced complications, further studies are needed to confirm this trend.

Septal T1 values do not seem to be affected by the development of fibrosis; previous autopsy work found very little amount of fibrosis in the septum in 3 cases of severe cardiac involvement in FD [START_REF] Sado | Identification and assessment of Anderson-Fabry disease by cardiovascular magnetic resonance noncontrast myocardial T1 mapping[END_REF][START_REF] Thompson | T1 mapping with cardiovascular MRI is highly sensitive for Fabry disease independent of hypertrophy and sex[END_REF][START_REF] Sheppard | A detailed pathologic examination of heart tissue from three older patients with Anderson-Fabry disease on enzyme replacement therapy[END_REF]. We observed that in patients with concentric hypertrophy, the fibrous scar on LGE images has most often an inferior-lateral-basal location and this was recently demonstrated by Deva et al. in 39 patients [START_REF] Deva | Cardiovascular magnetic resonance demonstration of the spectrum of morphological phenotypes and patterns of myocardial scarring in Anderson-Fabry disease[END_REF].

As reported by Wilson et al., the study of myocardial strain in patients with FD using cine sequences could be a useful addition to T1 mapping [START_REF] Wilson | Comparison of native T1, strain, and traditional measures of cardiovascular structure and function by cardiac magnetic resonance imaging in patients With Anderson-Fabry disease[END_REF]. Nordin et al. suggested that inflammation could contribute to LGE in patients with FD [START_REF] Nordin | Cardiac Fabry disease with late gadolinium enhancement is a chronic inflammatory cardiomyopathy[END_REF]. T2 mapping and troponins could be assessed to evaluate inflammation, stratify the risk and monitor therapies. It has been shown that in the absence of myocardial fibrosis, ERT could improve myocardial morphology, function and exercise capacity [START_REF] Weidemann | Long-term effects of enzyme replacement therapy on Fabry cardiomyopathy: evidence for a better outcome with early treatment[END_REF]. Today, as there is no available morphological or functional biomarker of treatment effectiveness, we are performing a prospective multicentric research (ClinicalTrials.gov Identifier: NCT02956954) in order to try to address this question.

Our study has several limitations. First, it is retrospective single center study that included a limited number of patients with FD. However, FD is a rare disease and considering the significantly lower native T1 values compared to control subjects, only a small sample size was sufficient to detect significant changes. Previous sample size calculations have been performed showing that 7 patients were required to detect a 40 ms difference [START_REF] Bhuva | Reproducibility of native T1 mapping using ShMOLLI and MOLLI -implications for sample size calculation[END_REF].

Considering the small number of patients with FD, we could not divide them into sub groups (e.g., patients without hypertrophy or patients with ERT). Also, we did not measure the extracellular volume in our population. Deva et al. reported that in the case of non-concentric hypertrophy, the topography of LGE was different and most often multiple [START_REF] Deva | Cardiovascular magnetic resonance demonstration of the spectrum of morphological phenotypes and patterns of myocardial scarring in Anderson-Fabry disease[END_REF]. We were unable to reproduce this result in our study as the vast majority of patients with FD had concentric hypertrophy. Reproducibility of T1 values measurements was not assessed because previous studies had shown it to be excellent [START_REF] Piechnik | Normal variation of magnetic resonance T1 relaxation times in the human population at 1.5 T using ShMOLLI[END_REF][START_REF] Knobelsdorff-Brenkenhoff | Myocardial T1 and T2 mapping at 3T: reference values, influencing factors and implications[END_REF]. A single operator performed all ROI tracings accordingly with rules of good practice. Multivariate analysis (age, sex and heart rate) was not possible considering the limited number of patients. A previous study showed that the heart rate was the only physiological factor that had a small effect on myocardial T1 values (-6 ms/10 bpm) [START_REF] Thompson | T1 mapping with cardiovascular MRI is highly sensitive for Fabry disease independent of hypertrophy and sex[END_REF]. This study involved only the MOLLI ® technique while at least 4 other T1 mapping schemes have been assessed [START_REF] Fontana | Comparison of T1 mapping techniques for ECV quantification: histological validation and reproducibility of ShMOLLI versus multibreath-hold T1 quantification[END_REF][START_REF] Or | Role of transcytolemmal water-exchange in magnetic resonance measurements of diffuse myocardial fibrosis in hypertensive heart disease[END_REF]. However, in different Look-Locker sequences, T1 values were shown to correlate in clinical studies [START_REF] Nacif | Myocardial T1 mapping with MRI: comparison of look-locker and MOLLI sequences[END_REF]. In accordance with a recent study, our results may potentially only apply to our equipment given the interequipment variability [START_REF] Brunelle | Variability induced by the MR imager in dynamic contrast-enhanced imaging of the prostate[END_REF].

In conclusion, native CMR imaging T1 values at 1.5 T are reduced in FD patients compared with HCM patients and healthy volunteers; therefore, this imaging biomarker could be instrumental in the early identification of cardiac involvement in patients with FD. Finding a low T1 value should be considered as a possible hint of FD and should refer the patient to a specialist in order to confirm the diagnosis.

Figure 1 .

 1 Figure 1. Flowchart of the study.

Figure 2 .

 2 Figure 2. 50 year-old woman with Fabry disease. A, Short axis T1 map reveals decreased T1 value in the septum (915 ms) and pseudo-normalized T1 value in the infero lateral wall (1006 ms). B, Late gadolinium enhancement image shows typical scar in the infero-lateral segment (arrow) indicating fibrosis.

Figure 3 .

 3 Figure 3. Graph shows mean septal left ventricular T1 values (in ms) in the three study groups. Healthy indicates healthy control subjects; HCM indicates hypertrophic cardiomyopathy; FD indicates Fabry disease.

Figure 4 .

 4 Figure 4. Graph shows mean global left ventricular T1 values (in ms) in the three study groups. Healthy indicates healthy control subjects; HCM indicates hypertrophic cardiomyopathy; FD indicates Fabry disease.

Figure 5 .

 5 Figure 5. Graph shows mean global left ventricular T1 values (in ms) in men and women with Fabry disease.

Table 2

 2

	. Nine of 17 patients (53 %) with FD and

Table 1 .

 1 Baseline characteristics in patients with Anderson-Fabry disease, patients with hypertrophic cardiomyopathy and healthy volunteers.

Table 2 .

 2 Cardiac MR imaging findings.

Acknowledgements

The authors express their thanks to CMR technologists of Rouen University Hospital whose expertise and dedication to excellent image quality made this research possible. We also thank the Statistics Department for their helpful input.

Disclosure of interest

No potential conflict of interest relevant to this article was reported.

Conflicts of interest

The authors have no conflicts of interest to disclose in relation with this study. * Only 14 patients of the AFD Group received an injection of gadolinium chelate, one examination was not interpretable.

Variables


Sergey N Solodukhin
email: sergey.solodukhin@lmpt.univ-tours.fr

Logarithmic terms in entropy of Schwarzschild black holes in higher loops

The Bekenstein-Hawking (BH) entropy is expected to be modified by certain correction terms in the quantum loop expansion. As is well known the logarithmic terms in the entropy of black holes appear as a one-loop addition to the classical BH entropy. In this note we study the further modifications of the logarithmic terms in the entropy of the Schwarzschild black holes due to higher quantum loops: up to three loops in a general renormalizable theory of gauge fields, scalars and fermions and two loops in quantum gravity. For a large class of field multiplets (including that of the Standard Model) that include graviton and for a certain range in the values of the couplings these modifications manifest themselves in cooling down the black holes at later stages of evaporation and, respectively, in increasing the life time of the black hole. If this picture persists to even higher loops, then the small black holes formed in the early stages of the cosmic evolution do not evaporate completely by now as is predicted in the standard picture. Instead, their long-lived (Planckian mass) remnants are present in abundance in today's Universe.

Since the works of Bekenstein and Hawking [1] it is generally accepted that the black holes, considered to be exact solutions to General Relativity, are characterized by a certain entropy proportional to the area of the horizon. The Schwarzschild black hole is a solution that is described by a single dimensionful parameter, mass M or horizon radius r + = 2GM . Classically, it has the entropy

S BH = A + 4G = 4π M 2 M 2 P L , (1)
where A + = 4πr 2 + is horizon area and we introduced the Planck mass, M P L = 1/ √ G, and the temperature T -1 H = ∂ M S BH = 8πM/M 2 P L . This temperature becomes unbounded as soon as mass M of the black hole decreases. So that the evaporation rate

dM dt ∼ -T 4 H A + (2)
accelerates for small black holes as dM/dt ∼ -M 4 P L /M 2 and the black hole evaporates in finite time t BH ∼ M 3 0 /M 4 P L , where M 0 is the initial mass.

At the quantum level, when matter and/or gravity, are quantized the classical formula (1) gets modified. One way to think about this is to consider the entanglement entropy of the quantum fields. This entropy requires a UV regulator to be properly defined. The entropy contains both the area law similar to (1) [2] and a logarithmic term, which for a Schwarzschild black hole takes a simple form [3] - [6],

S ent = A + 48π 2 + s 0 ln r + , (3)
where the pre-factor in front of the logarithmic term is

s 0 = 1 45 (N 0 + 7 2 N 1/2 -13N 1 - 233 4 N 3/2 + 212N 2 + 91N A ) (4)
for a multiplet of N 0 scalars, N 1/2 Dirac fermions, N 1 vector fields, N 2 particles of spin 2 (gravitons), N 3/2 fields of spin 3/2 (gravitino) and N A rank 2 antisymmetric tensor fields. It should be noted that the existing technique allows us to compute only the UV divergent parts of the entropy. The dependence on r + in the logarithmic term then comes from the two facts: that the entropy is a dimensionless quantity and that there is only one dimensionful parameter, r + that characterises the geometry. So that the entanglement entropy may only be a function of the combination (r + / ). The other remark is that the scalar field contribution in (4) to logarithmic term is the same for any value of non-minimal coupling ξ provided the latter is introduced as ξφ2 R. The quantum entropy (3) can be considered as a one-loop correction to the classical entropy (1). In the total entropy, S BH + S ent , the UV divergences then can be absorbed in the renormalization1 of the Newton constant, 1/G ren = 1/G + 1/12π 2 , and the R 2 couplings in the gravitational action, see [9]. So that one finds for the total black hole entropy expressed in terms of the renormalized couplings

S 1-loop = 4π M 2 M 2 P L + s 0 ln M/µ , (5)
where M P L is defined with respect to the renormalized G ren and we omitted the irrelevant constants. Respectively, the temperature is modified,

T -1 = 8πM/M 2 P L + s 0 /M . (6)
The modification becomes important when the mass of the black hole approaches the Planck mass. Extrapolating this approximation to those masses one finds that the evaporation scenario depends on the sign of s 0 . If s 0 = -q < 0 then there exists a certain mass M min = √ 8πqM P L for which the entropy as function of mass develops a minimum and the temperature becomes infinite (just as for M = 0 in the classical case). The black hole evaporates down to this minimal mass in a finite time t ∼ (M 0 -M min ) 5 /(M 2 0 M 4 P L ). Formally extended to region where M < M min the temperature becomes negative there.

On the other hand, if s 0 > 0 then the entropy is monotonic function of mass while the temperature develops a maximum at M 2 * = s 0 M 2 P L /8π , T max = M P L / √ 32πs 0 . Extrapolating the formulas below this mass the temperature decreases and the black hole cools down so that T ∼ M for small M . The evaporation rate (2) then slows down and the black hole evaporates in infinite time. Notice that in the case when s 0 1 the critical mass M * can be well above the Planck scale. For masses M > M * the specific heat is negative in a similar way as it follows from the Bekenstein-Hawking entropy [1] and, irrespectively of black holes, is in fact valid for a wide class of gravitationally coupled systems [START_REF] Thirring | Systems with negative specific heat[END_REF]. However, for smaller masses M < M * the specific heat becomes positive 2 .

The entropy (5) is not automatically positive function of mass that may signal of the breakdown of the one-loop approximation in the region where the entropy becomes negative. The entropy is positive (and the approximation is reliable) for M > µ. For sufficiently small µ (or large values of s 0 ) the threshold point lies well below the critical point M * so that the most of the interesting modifications discussed above happen well above the point where the approximation breaks down.

The other way to obtain entropy (5) is to consider the effective gravitational action,

W gr = - 1 16πG R + W Q (7)
where the quantum part W Q is represented as an expansion in powers of the Riemann curvature and its covariant derivatives. It contains both local and non-local terms. In the one-loop approximation W Q = 1 2 ln det D , where D is the operator that governs the quadratic perturbations. The entropy then arises as a response of the gravitational action to a small angle deficit δ = 2π(1 -α) at the horizon Σ, S = (α∂ α -1)W gr | α=1 , see [4] for more details on this formalism. The conical singularity manifests itself in the singular terms in the Riemann curvature concentrated at the horizon Σ, see [14],

R µναβ = 2π(1 -α)((n µ n α )(n ν n β ) -(n µ n β )(n ν n α )δ Σ + . . . , (8)
where . . . stand for the regular terms in the curvature and n a µ , a = 1, 2 is a pair of vectors normal to Σ, (n µ n ν ) = 2 a=1 n a µ n a ν . Discussing the entropy of the Schwarzschild black holes, for which the Ricci tensor identically vanishes, we need to look only at the terms that contain the Riemann tensor. Focusing on the quadratic terms one finds

W Q = s 0 64π 2 R µναβ R µναβ ln + . . . . (9)
Applying formula (8) to ( 9) and using that Σ R µναβ n a µ n b ν n a α n b β = 8π one arrives at (3). Notice that ln r + term is then restored by the dimensionality arguments. This term, in fact, comes from the non-local (UV finite) part of the quantum action W Q .

Notice that formula ( 5) is valid for any, not necessarily conformal, massless field. However, if computed in a 4d CFT the logarithmic term takes the form [8]

s 0 = 64π 2 (C -A) (10)
and is, thus, related to the conformal charges A and C that appear in the conformal anomaly,

< T >= AE 4 -CW 2
, where E 4 is the Euler density and W 2 is the square of the Weyl tensor. Notice that in certain theories s 0 vanishes. This is so, for instance, for N = 4 super-Yang-Mills theory.

In the entanglement entropy of massive fields there appear new terms, both UV divergent and finite, that are due to mass m. However, most of them go away in the total entropy expressed in terms of the renormalized Newton constant. Such terms appear already in Minkowski spacetime, see [10],

S ent = c s A + 48π 2 ds s 2 e -ms 2 , (11)
where for a scalar c 0 = 1 and for a Dirac fermion c 1/2 = 2. The respective renormalization of the Newton constant goes as

1 G ren = 1 G + c s A + 12π 2 ds s 2 e -ms 2 . ( 12
)
In a curved space-time there appear more terms in the entropy that depend on m. However, the logarithmic term in the entropy remains unchanged. This is due to the fact that in the logarithmically UV divergent terms the mass m may appear only in a combination m 2 R that contributes to the renormalization of the Newton constant that is already taken into account in [START_REF] Thirring | Systems with negative specific heat[END_REF]. There, however, may appear new UV finite terms in the entropy of the form f (r + m) which grow slower than logarithm 3 . Exact form of those terms is not yet known and is to be determined. In our analysis we shall ignore these terms concentrating our attention on the logarithmic terms only.

The story described up to this point is what we have at the level of one loop, when only the interaction of the matter with the background gravitational field is taken into account. The sign of the logarithmic term in ( 5) is essential when we discuss the evaporation of the black hole since the logarithmic term becomes important at the later stages of the evaporation when mass of the black hole becomes comparable to the Planck mass. For negative values of s 0 the evaporating black hole reaches, in a finite time, the stage of infinite temperature, just as in the classical case. The further evolution is difficult to predict within the present analysis. Since this configuration corresponds to the minimum of the entropy the system would probably tend to absorb the energy rather than emit it. Conceptually, this situation is not an improvement over the classical case.

On the other hand, in the case of positive s 0 the black hole evaporation is much less violent. As soon as the black hole reaches some maximal temperature its mass then further decreases and the hole starts to quickly cool down. As a result, in this case, the black holes never evaporate completely, at the later times they are present as the long-lived sub-Planckian remnants.

As we see it from (4) the sign of s 0 is not a priori definite, it depends on the multiplet of fields existing in Nature. Although many fields contribute positively to s 0 some of them, as vector fields, contribute negatively. s 0 can be computed in the Standard Model which contains 24 fermions, 12 gauge bosons and 4 scalars (a complex doublet). Provided one graviton is added to this multiplet eq. ( 4) gives us a positive value, s 0,SM = 16/5. Note that without the graviton this value would be negative, -68/45.

The goal of this note is to investigate how the one-loop approximation ( 5) is modified when the interactions are turned on and the respective higher quantum loops are taken into account.

By interactions here we mean both the interaction of the matter fields and the purely gravitational interaction when particles of spin 2 are considered. In the case of gravity the coupling constant is dimensionful and we expect there to appear corrections to (5) of the form M 2 P L /M 2 . Identifying those correction terms one expands the domain of validity of the approximation.

Our strategy is first to look at the logarithmic UV divergent terms in the effective action, then compute the respective UV terms in the entropy using (8) and finally restore the dependence on r + using the dimensionality arguments, in the same line as in one loop, (3) and (5). We base our further analysis on the results available in the literature: the instrumental works of Jack and Osborn [20] and Jack [21] - [22] (see also a recent paper [24] where some typos in the previous publications have been corrected) for the matter fields and the classical papers by Goroff and Sagnotti and by van de Ven [29] for gravity. (In the case of gravity the one-loop UV terms were computed in [28] that results in the graviton contribution in (4).) We use these earlier results and compute the corresponding modifications in the entropy. Earlier, it was suggested by Sen [7] that the logarithmic term is a one-loop effect and it does not get renormalized in higher loops. This is valid for a theory with a dimensionful coupling constant such as quantum gravity. In the interacting quantum field theories with dimensionless couplings one expects the possible corrections to the one-loop result and, indeed, as we show below these corrections do appear in the higher loops.

Before proceeding one technical remark is that we use here a regularisation, such as heat kernel or Pauli-Villars, in which the UV regulator is dimensionful. In the literature one quite often uses the dimensional regularisation, d = (4 -d). The conversion rule between the two regularisations is the following: 1 d = -ln . Earlier works on entanglement entropy in the interacting field systems include [15] - [19].

Below we consider some examples of interacting quantum field theories.

General renormalizable scalar field theories. Consider a multi-component scalar field φ i , i = 1, . . . , N 0 with a general Lagrangian

L = 1 2 (∇φ) 2 + V (φ) , (13)
where the potential takes the form

V (φ) = 1 4! λ ijkl φ i φ j φ k φ l + 1 3! g ijk φ i φ j φ k + 1 2 m ij φ i φ j + 1 2 ξ ij φ i φ j R + . . . . (14)
The corresponding UV divergences up to 3 loops in this theory have been calculated by Jack and Osborn [20] in 1984. Leaving aside the terms linear in Ricci scalar R, that as we discussed above will contribute to the renormalization of the Newton constant we are interested in terms which are quadratic in the Riemann tensor. Those terms in the theory (13) appear in one loop, the corresponding contribution to the entropy is given by the scalar field part in (5) and in three loops 4 , as was shown by Jack and Osborn 5 , W

Q = 1 (4π) 6 µ -d d 1 2592 λ ijkl λ ijkl (R αβµν R αβµν -2R µν R µν + 1 3 R 2 ) + . . . , (3)
where we keep only terms quadratic in curvature. Notice that ( 15) represents the renormalization of the Weyl-square term. Applying (8) we compute the corresponding contribution to the entropy of the Schwarzschild black hole and find that s 0 in the three loop approximation reads

s 0 (λ) = N 0 45 - 1 (4π) 4 1 648 λ ijkl λ ijkl . ( 16
)
This value of s 0 is smaller than the one for the free scalar fields.

Gauge theory coupled to fermions and scalars. A system of a non-abelian gauge field A = A a µ t a coupled to Dirac fermions ψ and real scalars φ carrying the representations of a simple gauge group G is described by the Lagrangian

L = 1 4g 2 Tr(F µν F µν ) + 1 2 (Dφ) T Dφ + ψD /ψ + 1 2 ψY i ψφ i , (17)
where

F µν = ∂ µ A ν -∂ ν A µ + [A µ ,
A ν ] and one included the Yukawa couplings Y i . The two loop computation in this theory was performed in [21], [22], [23] (some typos have been corrected in [24], see also [25] for further developments and [26] for a non-perturbative CFT consideration).

As in (15) the appropriate UV divergence comes from the renormalization of Weyl-square term, W

Q = 1 (4π) 4 µ -d d g 2 ( 2 9 C A 2 - 7 72 C ψ 2 - 1 18 C φ 2 )N 1 + 1 96 TrY 2 R αβµν R αβµν + . . . , (2)
where C A 2 (C ψ 2 and C φ 2 ) is Casimir operator for gauge fields (respectively for fermions and scalars). In these notations the beta-function for the gauge coupling, [27]. So that one finds in the two loop approximation that

β g = g 3 (4π) 2 (-11 3 C A 2 + 2 3 C ψ 2 + 1 6 C φ 2 )
s 0 (g, Y ) = s f ree 0 - 1 16π 2 g 2 ( 2 9 C A 2 - 7 72 C ψ 2 - 1 18 C φ 2 )N 1 + 1 96 TrY 2 . ( 19
)
We see that the contributions of the gauge fields and of the Yukawa couplings tend to decrease the value of s 0 while the contributions of the matter fields (fermions and scalars) increase it.

In the QCD sector C A 2 = 3 and C ψ 2 = n F (number of flavors) and no scalars one has that

2 9 C A 2 -7 72 C ψ 2 = 48-7n F 72 . It is negative if n F ≥ 7.
In the Standard Model n F = 6 and, hence, 4 In two loops the UV term W

(2) Q contains terms at most linear in R , g ijk g ijk R , that contributes to the renormalization of the Newton constant and does not effect the logarithmic term. 5 Notice that our definition for W Q differs by sign from the one used in [20].

s 0 (g) < s f ree 0 . This observation and that the Yukawa couplings and λ-couplings in ( 16) make the negative contributions to s 0 appear to indicate that in the Standard Model the value of s 0 is smaller than in the free theory. To see whether this value remains positive (provided we included the graviton's contribution in one loop) requires a more careful analysis including the actual values of all couplings involved6 .

Quantum Gravity. The one-loop calculation of UV divergences in quantum gravity has been first performed by 't Hooft and Veltman and later confirmed by various authors [28]. For pure gravity the divergent term in the effective action is

W (1) Q = - 1 (4π) 2 µ -d d ( 53 45 E 4 + 7 20 R 2 µν + 1 120 R 2 ) , (20)
where E 4 is the Euler density. This term results in graviton's contribution to logarithmic term s 0 presented in (5).

The two loop computation in quantum gravity was performed by Goroff and Sagnotti and later by van de Ven [29],

W (2) Q = - 209 2880 µ -d d κ (4π) 4 R µναβ R αβσρ R σρ µν , (21)
where κ = 32π 2 G. We keep here the term that contains only the Riemann tensor. This term is gauge independent and not vanishes on-shell. In fact, in two loops there may appear more terms than just (21), in particular the products of two Riemann tensors and either Ricci tensor or Ricci scalar (see paper of Goroff and Sagnotti for some of such terms). The complete information on the exact numerical factors with which those terms appear in the effective action is not available in the literature. These terms, however, appear to depend on the gauge and are unlikely to make any contributions to physical quantities such as the entropy.

Concentrating our attention only on the Riemann tensor we see that the one-loop and twoloop results (20) and (21) have the same sign. So that the two-loop correction due to gravitons to s 0 comes with the same sign as the one-loop result. Using (8) one finds that on a conical manifold

κ Mα R µναβ R αβσρ R σρ µν = 12κπ(1 -α) Σ R abσρ R σρ ab + • • • = 192π 2 (1 -α) κ r 2 + + . . . , (22)
where . . . stands for terms quadratic in (1 -α) and r + = 2GM is the horizon radius. With this relation we find

S(M ) = 4π M 2 M 2 P L + s(M ) ln M µ , s(M ) = s 0 + σ M 2 P L M 2 , σ = 209 480 ( 23
)
for the entropy of the Schwarzschild black hole in the two-loop approximation. In agreement with the arguments of [7] there is no higher loop gravity corrections to value of s 0 . However, the two loop computation results in a new σ -term which is negligible for masses much larger than the Planck mass but it gives an important positive contribution for M ∼ M P L . This σ -term, combined with the standard s 0 -term, can still be regarded as a "logarithmic term" keeping in mind that the coefficient in front of the logarithm now becomes a function of the mass. In the higher gravitational loops there will appear higher powers of M 2 P L /M 2 in the function s(M ).

Modified temperature and evaporation. With the entropy as function of the mass given by eq.( 23) one finds the modified temperature

M P L T -1 = 1 x 3/2 (8πx 2 + s 0 x -σ ln x λ ) , x = M 2 /M 2 P L , λ = µ 2 e M 2 P L (24)
In the complete analysis of the temperature one should study the dependence on three parameters: s 0 , σ and λ. Here we just mention that the presence of the σ -term with σ > 0 makes the temperature positive for any M even if s 0 is negative. Indeed, for -s 0 = q > 0 one finds that the function in the right hand side of ( 24) is everywhere positive if two conditions are satisfied: 1) q < 16πλ and 2)

q 2
32π 2 ln (16πλ/q) < σ . Thus, for sufficiently small values of q the temperature remains positive function of M and it develops a maximum below which the evaporating black hole cools down as soon as its mass decreases. Extrapolating (2) to the process of evaporation with the temperature (24) one finds that at later times the mass of the evaporating black hole falls off as t -1/13 that is a much slower rate compared to t -1/5 in the case of vanishing σ .

We conclude with some remarks.

1. Strictly speaking, the extrapolation of the 1-and 2-loop results for the entropy to values of mass M approaching the Planck mass goes beyond the validity of the loop approximation and, thus, is outside the domain of applicability. However, the same and even for more reasons could (and should) be said about the applicability of the classical BH formula for the entropy. BH formula is not valid for small black holes of mass M ∼ M P L and should be replaced with a better one. We believe that the higher loop results, presented here and yet approximative, show us a certain tendency in the behaviour of S(M ) for small M which, likely, will become even stronger in the higher loops when more terms in the function s(M ) will be available. The higher loop contributions will appear as the higher order terms in the series with respect to M 2 P L /M 2 . It is possible that, provided all loops are included, these terms will sum up to something non-analytic at M = 0.

2. The other limitation of our results comes from the fact that here the quantum corrections to the entropy are computed in a (still classical) Ricci flat geometry. In a fully consistent consideration one would have to also analyse the modification of the geometry itself due to the quantum corrections. On one hand, this is a more difficult problem. On the other, in a nonperturbative analysis of the quantum corrected Einstein equations the black hole may appear to be replaced with a rather different (and yet longer-lived) object as the consideration of [30] shows. These are different indications that the existence of the long-lived compact objects in the complete theory is a likely possibility.

3. The analysis based on the use of the BH entropy plays an important role in the story of primordial black holes (for a review see, for instance, [31]). Namely, it imposes a lower limit on the mass of the primordial black holes since the small ones have enough time to evaporate completely until the present days. Provided one uses the modified version of the entropy with the logarithmic terms present this conclusion is no more in place: the black hole evaporation slows down as soon as mass decreases below the critical mass M * ∼ √ s 0 M P L that leads to much longer life time of the black holes. In this picture the small primordial black holes do not evaporate completely. Instead, all of them are present in today's Universe as some long-lived Planckian compact objects. From a phenomenological point of view this is a particular case of a more general conclusion of [32] about the much longer life time of black holes for a wide class of possible modifications of the Hawking's formula for the temperature as function of mass M .

There are certain subtleties[10] due to the presence of the non-minimal coupling[11]. These are not, however, essential in the present discussion.

In the critical point M = M * the specific heat diverges that makes it quite attractive to interpret this point as a second order phase transition (for a discussion of singular behavior of specific heat in the second order transition see, for instance,[START_REF] Abe | Singularity of Specific Heat in the Second Order Phase Transition[END_REF]). This is an interesting direction to explore that we leave for a separate study.

For instance, one anticipates an exponential term of the form e -r 2 + m 2 which for a large range of parameters (M, m) leads to smaller effects than the logarithmic term.

In N = 4 SU (N ) SYM the C -charge is not expected to renormalise and, hence, the shift in s 0(19) is expected to vanish. I thank K. Skenderis for this remark.

Acknowledgements

I would like to thank B. Pioline for a discussion on the logarithmic terms, K. Krasnov for discussions on the loop calculations in the background field formalism, and, especially, M. Shaposhnikov for warm hospitality in his group at EPFL and many inspiring discussions during the initial stages of this project. The present work is supported in part by the ERC-AdG-2015 grant 694896. I also acknowledge the useful communications with I. Jack regarding his papers and with K. Skenderis on renormalization of Weyl-square term in 4d CFTs and many insightful conversations with my colleagues in Tours: M. Chernodub, N. Mohammedi and S. Nicolis.


Dong Ding
email: dong.ding@utc.fr

Abdellatif Ouahsine
email: ouahsine@utc.fr

Weixuan Xiao

Peng Du

Numerical study of ballast-flight caused by dropping snow / ice blocks in high-speed railways using Discontinuous Deformation Analysis (DDA)

Keywords: High speed railway, ballast flight, snow-dropping, numerical modeling, discontinuous deformation analysis

Frozen snow / ice blocks drop at high speeds from train causing the ballast to fly up and damage the car body. Thus, in this paper, we propose a numerical model based on the discontinuous deformation analysis (DDA) method to study the ballast flight caused by dropping snow / ice blocks in high-speed railways and to analyze the dynamic behavior of ballast particles during their collision with a snow / ice block. The validation of the proposed model is done by comparing the numerical results with the theoretical and the experimental ones. The numerical results show that the velocity, shape and incident angle of snow / ice block play an important role in the ballast flight. Specifically, the number and the maximum displacement of ballast particles increase as the train speed increases and the incident angle greatly affects the movement direction of ballast particles. The shape of the ice block affects the amount and extent of ballast flight.

Introduction

The phenomenon of ballast flight is one of the major problems in high speed ballasted track, which has resulted in major maintenance costs and safety concerns. Flying ballast particles may hit the rail, the train body or the passengers through stations. Furthermore, small particles of ballast may come to rest between the railhead and the wheels of rail vehicles, which cause substantial local bending damage to the rail [START_REF] Quinn | A full-scale experimental and modelling study of ballast flight under high-speed trains[END_REF]. The 5 aerodynamic effect is commonly regarded as the main cause of ballast flight [START_REF] Sanz-Andres | The initiation of rotational motion of a lying object caused by wind gusts[END_REF][START_REF] Premoli | Ballast flight under high-speed trains: Wind tunnel full-scale experimental tests[END_REF]. However, as counted by M. R. Saat et al. [START_REF] Saat | Identification of high-speed rail ballast flight risk factors and risk mitigation strategies[END_REF], about 50% ballast flight incidents occur in snow condition. This is because high-speed trains run through snow zone and blow up the snow which sticks to the underfloor equipment and freezes rapidly into ice, then, the frozen snow and ice drops at high speeds from trains due to temperature changes, train vibration and heat from the brakes, causing the ballast to fly up [START_REF] Loponen | Studies of snow-dropping from a train on a turnout due to dynamic excitations[END_REF][START_REF] Kawashima | Experimental studies on ballast-flying phenomenon caused by dropping of accreted snow/ice from high-speed trains[END_REF], as shown in Fig. 1. The effects 10 of ballast flight caused by dropping snow/ice may lead to serious catastrophic consequences due to the high initial velocity, large mass, and chain reaction of the falling snow/ice. Snow / ice dropping and influencing factors of ballast flight have been studied in many recent studies.

Loponen et al [START_REF] Loponen | Studies of snow-dropping from a train on a turnout due to dynamic excitations[END_REF] studied the amount of the excitation required to drop snow from the train underframe. [START_REF] Kawashima | Experimental studies on ballast-flying phenomenon caused by dropping of accreted snow/ice from high-speed trains[END_REF] carried out an experimental study using air cannon tests to investigate the ballast-flight phenomenon caused by the dropping of accreted ice, and gave a relationship between the number of flying ballast stones and the mass, shape, speed, and angle of the ice. Furthermore, computational fluid dynamics (CFD) methods have been applied to the study of the flow between the train underbody and the track bed around the bogie area and its impact on the ballast flight [START_REF] Khayrullina | CFD simulation of train aerodynamics: traininduced wind conditions at an underground railroad passenger platform[END_REF][START_REF] Zhu | Flow between the train underbody and trackbed around the bogie area and its impact on ballast flight[END_REF][START_REF] Paz | Numerical methodology for evaluating the effect of sleepers in the underbody flow of a high-speed train[END_REF]. Xie et al [START_REF] Xie | Study of snow accumulation on a high-speed train's bogies based on the discrete phase model[END_REF] The numerical result showed that the snow particles accumulated and moved on the train bogies by the high-speed air. However, to the best of our knowledge, there are no computational techniques focused on the dynamic impact behavior between the ice block and the ballasted track bed which is the key reason resulting in the ballast flight during winter. In this study, the discontinuous deformation analysis (DDA) method is used to determine the displacement of irregular ballast particles and the collision between ice block and ballast particles. DDA is a kind of discrete element method (DEM), developed by G.H. Shi in 1989 for application of rock fracture mechanics and geotechnical or structural problems [START_REF] Shi | Generalization of two-dimensional discontinuous deformation analysis for forward modelling[END_REF][START_REF] Shi | Discontinuous deformation analysis: a new numerical model for the statics and dynamics of deformable block structures[END_REF]. Ishikawa T. et al. first used DDA for simulating deformation behavior of ballast [START_REF] Ishikawa | DDA applied to deformation analysis of coarse granular materials (ballast)[END_REF], which showed good performance in the study of coarse granular materials.

The main purpose of this paper is to determine the influence factors of ballast flight and to analyze the dynamic behavior of ballast particles during their collision with a snow / ice block. In the DDA method, the kinematic conditions of the contact surfaces are enforced through the penalty method in order to avoid the inter-penetration between blocks. The performance of the ballast flight is evaluated by considering the velocity, the shape and the incident angle of snow / ice blocks. The numerical results show that these factors influence significantly the ballast particles dynamics and also their flight as well as the collision between ballast particles.

Theoretical background and governing equations

As the computation domain is composed of discrete elements, we use the Discontinuous Deformation Analysis (DDA) method which is based on the principle of minimum potential energy. The DDA method does not have the meshing procedure of blocks and therefore, no refinement is needed to improve the quality of the calculated solution, which has the advantage of reducing the computation time [START_REF] Ji | Impacts of ship movement on the sediment transport in shipping channel[END_REF][START_REF] Linde | Three-dimensional numerical simulation of ship resistance in restricted waterways: Effect of ship sinkage and channel restriction[END_REF].

In the DDA method, the displacement (u, v) at any point (x, y) of a block i can be represented by six variables: two translations( u c , v c ) of the block gravity center (x c ,y c ) in x and y directions, a rotation γ c around (x c ,y c ), and two normal and a shear strains (ε xx , ε yy , ε xy ); therefore, the variables vector associated with the block i is given:

D i = u c , v c , γ c , ε xx , ε yy , ε xy T (1)
According to the first-order expression of any point (x, y), the displacement (u, v) for an individual block i can be written as:

U i = u v = T i D i (2)
where

T i =   1 0 -(y -y c ) (x -x c ) 0 (y -y c )/2 0 1 (x -x c ) 0 (y -y c ) (x -x c )/2   (3)
The strain of the block i can be expressed by the relationship between the strain and displacement:

ε i = LU i (4)
Where

L =      ∂ ∂x 0 0 ∂ ∂y 1 2 ∂ ∂y 1 2 ∂ ∂x     
is the differential operator matrix for 2D problem.

Substituting Eq. (2) into Eq. ( 4), we get:

ε i = LT i D i = BD i (5)
Assuming that the deformation is elastic and linear, the stress tensor is written as follows:

σ i = Eε i = EBD i ( 6
)
Where E is the elastic matrix of deformation planes and

B =      0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1      .
The total potential energy Π p of the block i, defined as the the sum of the elastic strain energy Π e , initial stress potential energy Π σ0 , body force potential energy Π b , and inertial energy Π i , is given by:

Π p = Π = Π e + Π σ0 + Π b + Π i = Ai 1 2 ε T i σ i dA i + Ai ε T i σ 0 dA i - Ai U T i f b dA i + Ai U T i m Di dA i (7)
Π p = Ai 1 2 ε T i Eε i dA i + Ai ε T i σ 0 dA i -D T i T T i Ai f b dA i - Ai mT i Di dA i (8)
where f b denotes the body forces applied on a block i, and m is the block mass per unit area. σ 0 is the initial stress of the block.

Substituting Eq. ( 5) and (6) into Eq. ( 8), the total potential of a system of N blocks is expressed as:

Π np = N 1 D T i M Di + 1 2 D T i KD i -D T i f e (9)
where M = Ai mT T i T i dA i is the mass matrix, K= Ai B T EBdA i is the stiffness matrix, f e = Ai (T T i f b -B T σ 0 )dA i is the external forces matrix. According to the minimized potential energy, the block system equations of motion can be represented in the compact form:

∂Π p ∂D i = 0 ⇒ M D + KD = F (10)
Using the Newmark-β method to solve the kinematic equation of motion Eq. ( 10), the displacement and velocity can be approximated as:

D n+1 = D n + ∆t Ḋn + ∆t 2 2 (1 -2β) Dn + 2β Dn+1 Ḋn+1 = Ḋn + ∆t (1 -γ) Dn + γ Dn+1 ( 11
)
where D and Ḋ are the acceleration and the initial matrix of each time step, while β and γ are velocity and acceleration weighting parameters respectively. Specifications for β and γ parameters lead to various time integration schemes. In the present investigation, we assume that β = 1/2 and γ = 1, which is an implicit scheme that gives a stable and accurate computed solution [START_REF] Doolin | Time integration in discontinuous deformation analysis[END_REF][START_REF] Khanh | Assessment of ship manoeuvrability by using a coupling between a nonlinear transient manoeuvring model and mathematical programming techniques[END_REF][START_REF] Cai | Moving immersed boundary method[END_REF]. Substituting Eq. (11) into

Eq. ( 10), we get the equations of motion:

(K + 2M ∆t 2 )D n+1 = F + 2M ∆t Ḋn ( 12
)
The final motion equations can be written in compact form as:

KD = F (13)
where K= 2M ∆t 2 + K, and F = F + 2M ∆t Ḋn . Note that in the present study we used an implicit scheme, thus the explicit damping is assumed to be zero. Furthermore, the DDA method is used in both static and dynamic analyses. The only difference is that in static calculations, the velocity at the beginning of each time step is assumed to be zero, while the dynamic calculation inherits the velocity of the previous time step which is contributed by the inertia force.

At every time steps, the simultaneous equilibrium of N blocks Eq. ( 13) can be written as:

           K 11 K 12 K 13 . . . K 1n K 21 K 22 K 23 . . . K 2n K 31 K 32 K 33 . . . K 3n . . . . . . . . . . . . . . . K n1 K n2 K n3 . . . K nn                       D 1 D 2 D 3 . . . D n            =            F 1 F 2 F 3 . . . F n            (14)
where F i and D i are the 6×1 matrices of forces and displacements of block i. The stiffness sub-matrices Thereafter, it should be noted that in the present investigations, in order to get the contact springs between blocks and avoid inter-penetration between blocks, the surface contact constraints are enforced by the penalty method [START_REF] Kaidi | Stability of breakwaters under hydrodynamic loading using a coupled DDA/FEM approach[END_REF], which is equivalent to placing a spring between the two blocks, by considering the contact of two blocks i and j, the penetration distance δ, of point P 1 into edge P 2 P 3 of block j, as shown in Fig. 3. This method does not increase the number of governing equations and the solution can easily be obtained by adding a contact sub-matrix to the stiffness matrix [START_REF] Shi | Discontinuous deformation analysis: a new numerical model for the statics and dynamics of deformable block structures[END_REF][START_REF] Jing | Formulation of discontinuous deformation analysis (DDA):an implicit discrete element model for block systems[END_REF]. The strain energy of the contact spring is:

K ij can
Π c = 1 2 k n δ 2 n + 1 2 k t δ 2 t ( 15
)
where k n and k t are the normal and tangential stiffness, while δ n and δ t are the penetration overlap in normal and tangential direction. The contact strain energy should also be minimized; therefore, the penalty number k n and k t are assumed to be large, which can achieve an accurate result.

Validation of the DDA method

Case1: Frictionless impact -Free fall

In order to examine the accuracy of impact in DDA, the progress of free fall and elastic rebound was simulated. The simulation consisted of 1 m× 1 m block falling 5.0 m onto a 1 m× 10 m base block with four fixed points, as shown in Fig. 4(a). We have an initial velocity of 0 m/s, and we assume an acceleration due to the gravity of 10 m/s 2 . The velocity of fall and rebound can be written as:

     V f all = V 0 + gt 0 ≤ t ≤ 1 V rebound = V t -gt 1 < t ≤ V t g (16)
where V f all and V rebound are the velocity during the fall and rebound. V t is the velocity when the block reaches the impact plane.

The velocity of free fall and rebound for four various contact stiffness values was simulated using the DDA method. The theoretical value of free fall and elastic rebound velocity after the frictionless impact was calculated by Eq. ( 16). and the numerical results are shown in Fig. 4(b), where we assume that the contact stiffness is related to the Young modulus of block material [START_REF] Doolin | Displacement accuracy of discontinuous deformation analysis method applied to sliding block[END_REF]. The results (Fig. 4(b)) show that the greater the contact stiffness is, the less the impact damping will be. It shows that the best agreement between theoretical and numerical results is met when the contact stiffness value is around 10 9 N/m. It is worth noting that for smaller contact stiffness, the penetration overlap is large and then the impact between blocks becomes inelastic. Comparison between the theoretical value (Eq. ( 16)) and DDA results for different contact stiffness

Case2: Dynamic behavior of ballast bed after collision. Comparison with experimental results

In this second test, the dynamic behavior of ballasted track bed at the time of collision with an ice block was verified according to the air cannon experimental results by KawashimaI et al [START_REF] Kawashima | Experimental studies on ballast-flying phenomenon caused by dropping of accreted snow/ice from high-speed trains[END_REF]. The behavior of a 6 600 mm × 105 mm ice block with a velocity of 80 km/h was simulated, and the computed results were compared with the experimental observations. The DDA simulations exhibited phenomenon of ballast flight similar to the ones observed experimentally, as shown in Fig. 5. The number of ballast particles thrown 330 mm (distance from the ballast to the train underframe) or higher N 33 for two different mass was calculated by the current DDA model. Due to the shape of the ice block in the air cannon test is cylindrical, its mass can be calculated by:

M ice = ρπ( d 2 ) 2 H ( 17
)
where ρ = 0.9340 × 10 -6 kg/mm 3 is the density of ice [START_REF] Petrovic | Review mechanical properties of ice and snow[END_REF], and d = 105 mm is the diameter of the cylinder.

The mass of the ice block depends on the height of the ice cylinder H. The results from both experimental and numerical simulations showed that the number of the flying ballast particles increases as the velocity and the mass of ice block increase, as shown in Fig. 6. The numerical results agree well with the experimental results.

Simulation results

Velocity of ice block

Cross section

The ballast shoulder has a width of 500 mm and a height of 100 mm, whereas the ballast bed has a thickness of 350 mm (see Fig. 7). Based on the European standard: aggregates for railway ballast (13450, 2002) [START_REF] Okonta | Effect of grading category on the roundness of degraded and abraded railway quartzites[END_REF], the gradation of ballast particles is from 31.5 mm to 63 mm, where the ballast particles have an irregular shape and are completely compacted. The ice block is assumed to be a circular shape with a diameter of 100 mm, which falls from a height of 600 mm on the ballast bed, where the falling distance is approximately equal to the height of the rail plus half of the height of the bogie. Adding the unit mass, elastic modulus, Poisson ratio, friction angle, etc. to the stiffness matrix in accordance to the physical characteristics of the ballast and the ice [START_REF] Petrovic | Review mechanical properties of ice and snow[END_REF][START_REF] Suhr | Parametrisation of a DEM model for railway ballast under different load cases[END_REF], as shown in Table 1. shows that the ice block after the collision first rebounded and then rolled. The displacement of ice block increases with the initial vertical velocity. When the vertical velocity exceeds 5 m/s, the ice block hits the rail, or even bounce and a secondary strike on the track bed happens when the velocity is more than 10 m/s. For the ballast after the collision, the higher the initial vertical velocity, the higher rebound of marked ballast is ejected, the severity of ballast bed responses is generated, and the greater the number of ballast particles is ejected. From Fig. 9, we can get the same results, the number of flying ballasts and flying height of the ballast increase with the initial vertical velocity.

Fig. 10 shows the displacement of the ice block and a marked ballast. In general, larger initial velocity induces greater displacement. However, unlike in the case of small and moderate initial vertical velocities (0 m/s and 5 m/s), during the collision process at initial vertical velocity of 10 m/s, there is a collision between ice block and rail, which may reduce significantly the displacement of ice block after the impact (see Fig. 10a). It follows that the ballast trajectories and displacement show that its movement may be divided into three phases. Taking the initial vertical velocity of 5 m/s as an example (see Fig. 10(b) ), the first phase sited between t= 0 s and t= 0.11 s corresponds to a stable state of the ballast which remains motionless before its impact by the ice blocks. The second phase sited between t= 0.11 s and t= 0.34 s corresponds to the impact and flying processes in which the ballast undergoes a rebound following a parabolic curve. The third phase sited after t= 0.34 s corresponds to the marked ballast fly-back. It should be mentioned that the severity and the duration of every phase strongly depend on the initial vertical velocity. The increased initial vertical velocity induces greater displacement of the ballast, and longer flight duration, which results in a more serious consequence.

Longitudinal section

The longitudinal section is more realistic since it considers the real velocity of the ice block dropping from high-speed trains. The high-speed railway is commonly defined by the maximum running speed exceeds 200 km/h. However, the running speed of high-speed railways in China, Japan and France exceeds 300 km/h and even the maximum test speed of the eastern line in France is 574.8 km/h [START_REF] Jing | High-speed railway ballast flight mechanism analysis and risk management-a literature review[END_REF]. In addition, many ultrahigh-speed railways are under development and under construction. In order to analyze the initial dropping velocity of the ice block, the initial collision position was set on the surface of the ballast bed. An initial vertical velocity of 3.5 m/s which corresponded to the falling height of 600 mm, the initial longitudinal velocity equal to the train running speed, and four initial longitudinal velocities from 100 km/h to 400 km/h with an interval of 100 km/h were used.

The trajectory of the ice block and marked ballast, as well as the phenomenon of ballast bed after collision, is represented in Fig. 11. It is clear that the track bed was damaged after the collision, the ballast particles were hit and displaced, and the surface of the track bed showed different degrees of depression. The severity of the damage of the track bed had a positive relationship with the initial longitudinal velocity. For the trajectory of the ice block, the ice block quickly flew upward and forward after the collision. For the trajectory of marked ballast, we know that the ballast flew to the front direction. The height and distance of the ballast flight increased with the initial longitudinal velocity. When the dropping longitudinal velocity is more than 300 km/h, the marked ballast may collide with the train underframe since the maximum height of the marked ballast after the collision is more than 330 mm (distance from the ballast bed to train underframe) [START_REF] Kawashima | Experimental studies on ballast-flying phenomenon caused by dropping of accreted snow/ice from high-speed trains[END_REF]. If the ice and ballast particles reach the train underframe, they may accelerate significantly due to the collision with the train or may strike off more ice /snow blocks, causing more serious consequences.

Indeed, the effects of ballast projection may vary from depend on where the particle lands. By considering the ballast particles flying height H, the impact risk is classified into three categories:

• If H ≤ 184 mm(height of UIC-60 rail), flying ballast particles have the possibility to hit the rail;

• If 184 mm < H ≤ 330 mm, flying ballast particles have the possibility to hit the wheels of the train ;

• If H >330 mm, flying ballast particles will certainly hit the bottom of the train, and probably the rail and the wheels of the train.

The displacement of the marked ballast and the ice block was calculated, as shown in Fig. 12. Results

show that the displacement of marked ballast is linear with increasing initial longitudinal velocity at t=0.5

s. The displacement of marked ballast of 400 km/h is about four times of that 100 km/h. Significantly, the initial longitudinal velocity of the ice block directly affects the displacement of ice block and ballast particles. From Fig. 13, the numerical results show that the number of flying ballast particles and the vertical displacement of marked ballast increases from 100 km/h to 400 km/h. When the longitudinal velocity of ice block is in the range of 200 km/h to 300 km/h, the results exhibit a peak of 33 % and 47 % of the total number of flying ballast and maximum flying height of marked ballast. Above 300 km/h, this process evolution increases gradually until a plateau is reached. This non-uniform evolution may be related to the gradation of ballast composed of different particle sizes. As the longitudinal velocity increases, most of the small size and a part of the middle size of ballast particles in the impacted area fly above the bed surface.

According to the above simulation results, the initial longitudinal velocity of the snow / ice block which is directly dependent on the running speed of the train, has a great impact on the ballast flight. Reducing the speed of trains is the simplest and most effective mitigation strategy. However, this method may defeat the main purpose of high-speed lines. The national center of operation of SNCF in France proposed a time schedule which divided train speed into three kinds: acceptable, tolerable and unacceptable. The weather forecast, mechanical simulation, and safety approach were taken into consideration.

Incident angle of ice/snow

Regarding the angle of incidence of the ice block, in the longitudinal section, this angle can correspond to different attachment position of the ice, as shown in Fig. 14. The incident angle is not an independent variable, while it is related to the train speed and the distance of drop as well as the ice dropping position; therefore, four angles (0 o , 30 o , 60 o , and 90 o ) were investigated in the simulation. The incident angles were controlled by the velocity of the longitudinal and vertical directions. The total velocity was set to 20 m/s.

Fig. 15 shows the trajectory of the ice block and the marked ballast as well as the response of the track bed after the collision at different ejection angles. Based on these simulation results, it can be seen that the incident angle is also one of the main factors affecting ballast flight. When the incident angle is 0 o , the

Shape of ice block

In the third simulation, circular, triangular and square ice blocks with a diameter/side length of 100 mm were placed 600 mm above the ballast bed surface at the center of the track. With the initial vertical velocity of 10 m/s, the impact results are shown in Fig. 17.

Three shapes of ice block trajectories are significantly different. The triangular ice intrudes the ballast bed directly after colliding, while the square ice bounces off. The trajectory of the circular bounced and hit the trackbed twice. Furthermore, according to the response of the ballast bed after collision, the triangular ice ejected the greatest number of particles and caused the greatest height of ejection, on account of its contact area with the track bed being the smallest, resulting in maximum pressure. The square ice has the widest and the deepest impact among the three, due to its larger contact area and greater mass. The results of the displacement of the ice block and the marked ballast are shown in Fig. 18. After the collision, the triangular ice has the smallest displacement and the square ice has the largest rebound. The triangular ice causes the maximum displacement of ballast particles, followed by the circular ice and then the square ice.

An investigation by the Japan Railway Technical Research Institute showed that the density of the packing snow and ice at the bottom of a train is from 150 to 900 kg/m 3 and the maximum weight of snow / ice that may drop is about 15 kg. Furthermore, the shape of the snow / ice block which drops from the train is various and irregular [START_REF] Navikas | Modelling of snow cover thickness influence on the railway construction temperature regime under variable weather conditions[END_REF]. Therefore, the ballast flight caused by different shapes of snow / ice block is very complex and serious.

Conclusions

A numerical model based on the discontinuous deformation analysis (DDA) method, was proposed to study the dynamic behavior of ballast stones and their collision with a snow / ice block. This study took into account the shapes of the ice blocks and the contacts between ballast particles, where we assumed that contact constraints were imposed through the penalty method. The ballast flight induced by the dropping snow / ice with some variations in intensity depends on the velocity, the incident angle, and shapes of the ice blocks. The main findings derived from the numerical simulation may be summarized as follows:

• The velocity of the snow / ice block, which directly depends on the running speed of the train and the position of the ice attached, has a great impact on the ballast flight. In the longitudinal-section, the number of flying ballast particles and their displacement increases from 100 km/h to 400 km/h. When the longitudinal velocity of the ice block is in the range of 200 km/h to 300 km/h, the results exhibit a peak of 33 % and 47 % of the total number of flying ballast and maximum of flying height of marked ballast. When the dropping longitudinal velocity is more than 300 km/h, the marked ballast may collide with the train underframe since the maximum height of the marked ballast after the collision is more than 330 mm. Hence, setting the maximum operating speed according to the weather conditions is an effective measure to reduce serious consequences.

• The angle of impact greatly influences the direction of the movement of the ballast particles. The angle of incoming snow / ice blocks, which indeed are the impacting projectiles, is closely related to the location of a detaching projectile. Hence, it is worthwhile that, in order to significantly reduce the ballast flight caused by the melting of snow / ice blocks, some measures should be taken to prevent snow accumulating in the train and snow settlement.

• The shape of the snow / ice block affects the extent of ballast flight. The triangular ice block intrudes the ballast bed directly after collision, while the square ice bounces off. The triangular ice ejects the greatest number of particles and caused the greatest height of ejection. The square ice has the widest and the deepest impact among the three.

Also, it is worth noting that very little research has been done on ice block and ballast particle breakage from the point of view of aerodynamic interactions and numerical modeling . The breakage process is very important as it may indicate the violence of the impact after the flying ballast collision step. Thus, the analysis of the breakage process by the DDA method is still an open problem which should be undertaken in future studies.

Figure 1 :

 1 Figure 1: Schematic illustration of the mechanism and process of ballast flight caused by dropping ice block.

used a 3 -

 3 D numerical model based on the coupling between Navier-Stokes equations based model and a discrete phase model to investigate the flow field in the presence of snow accumulation in the train underframe.

Figure 2 :

 2 Figure 2: The structure of stiffness matrix in the case of three blocks:(a) Location of three blocks; (b) Structure of stiffness matrix

Figure 3 :

 3 Figure 3: Interaction between two contacting blocks

Figure 4 :

 4 Figure 4: Falling block for testing frictionless impact and elastic rebound:(a) Schematic diagram of free fall motion; (b)

Figure 5 :

 5 Figure 5: Comparison between the numerical and experimental responses [6] of the blocks for the DDA method.

Figure 6 :

 6 Figure 6: Comparison between the numerical and experimental values of N 33 [6] for the DDA.

Figure 7 :

 7 Figure 7: Track bed model

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Computed ballast flight depending on the initial vertical velocity: (a) Vy= 0 m/s (free fall); (b) Vy= 5 m/s; (c) Vy= 10 m/s

Figure 11 :Figure 12 :Figure 13 :Figure 14 :

 11121314 Figure 11: Relationship between ballast flight and longitudinal velocity: (a) Vx= 100 km/h; (b) Vx= 200 km/h; (c) Vx= 300 km/h; (d) Vx= 400 km/h

Figure 15 :Figure 16 :

 1516 Figure 15: Ballast flight at different incident ejection angles: (a) α = 0 o ; (b) α = 30 o ; (c) α = 60 o ; (d) α = 90 o

Figure 17 :Figure 18 :

 1718 Figure 17: Relationship between ballast flight and ice block shape: (a) Triangle; (b) Circle; (c) Square

Acknowledgements

The study is supported by the China Scholarship Council.