Teresa Mele

Daniele Generali

Stephen Fox

Maria Pia Brizzi

Alessandra Bersiga

Manuela Milani

Giovanni Allevi

Simone Bonardi

Sergio Aguggini

Marco Volante

Teresa Mele¹ ⁻²

Daniele Generali¹

Manuela Milani¹

Giovanni Allevi¹

Simone Bonardi¹

Sergio Aguggini¹

Luigi Dogliotti²

Alberto Bottini¹

Adrian Harris

Alfredo Berruti²

Anti-angiogenic effect of tamoxifen combined with epirubicin in breast cancer patients

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Chemotherapy and endocrine therapy are the standard treatments of breast cancer with positive hormonal receptors and their association is known to be more efficacious than either treatment alone in terms of disease free survival (DFS) and overall survival (OS) [START_REF] Carlson | Adjuvant endocrine therapy in hormone receptor-positive postmenopausal breast cancer: evolution of NCCN, ASCO, and St Gallen recommendations[END_REF] [2].

Chemotherapy followed by endocrine therapy is now recognized as the optimal method to combine the two treatment strategies. A multicenter prospective American trial comparing FEC with concomitant Tamoxifen versus sequential Tamoxifen administration, demonstrated that concomitant administration is detrimental in terms of DFS (INT 100 trial [START_REF] Albain | Adjuvant chemohormonal therapy for primary breast cancer should be sequential instead of concurrent: initial results from Intergroup Trial 0100 (SWOG-8814)[END_REF]).

Whilst acknowledging that chemotherapy and endocrine therapy have different mechanisms of action, pharmacological targets and toxicity patterns, it is still possible that concomitant administration could be helpful since tamoxifen has effects besides inhibition of proliferation. Numerous animal models of breast cancer have previously demonstrated a pro-effect of estrogens and an anti-angiogenic effect of tamoxifen in vivo [START_REF] Garvin | Effects of oestradiol and tamoxifen on VEGF, soluble VEGFR-1, and VEGFR-2 in breast cancer and endothelial cells[END_REF]; [START_REF] Haran | Tamoxifen enhances cell death in implanted MCF7 breast cancer by inhibiting endothelium growth[END_REF]; [START_REF] Lindner | Effects of tamoxifen and interferon-beta or the combination on tumor-induced angiogenesis[END_REF]; [START_REF] Guo | An antiangiogenic urokinase-derived peptide combined with tamoxifen decreases tumor growth and metastasis in a syngeneic model of breast cancer[END_REF]; [START_REF] Garvin | Tamoxifen Inhibits Secretion of Vascular Endothelial Growth Factor in Breast Cancer in Vivo[END_REF]; [START_REF] Elkin | An angiogenic switch in breast cancer involves estrogen and soluble vascular endothelial growth factor receptor 1[END_REF]. Estradiol increases extracellular levels of vascular endothelial growth factor (VEGF) while tamoxifen inhibits the secretion of VEGF in breast cancer in vivo [START_REF] Garvin | Tamoxifen Inhibits Secretion of Vascular Endothelial Growth Factor in Breast Cancer in Vivo[END_REF]. These data suggest that VEGF expression in breast is at least in part modulated by estrogen receptor, and VEGF gene is a target of for estrogen receptor [START_REF] Applanat | Vascular endothelial growth factor is a target gene for estrogen receptor and contributes to breast cancer progression[END_REF].

VEGF-A, a key mediator of tumour angiogenesis, is known to exert its angiogenic effects via two tyrosine kinase receptors, VEGFR-1 (flt-1) and VEGFR-2 (flk-1/KDR) [START_REF] Kerbel | Tumor angiogenesis[END_REF]; [START_REF] Ferrara | The biology of VEGF and its receptors[END_REF]; [START_REF] Ferrara | Angiogenesis as a therapeutic target[END_REF]; [START_REF] Shibuya | Differential roles of vascular endothelial growth factor receptor-1 and receptor-2 in angiogenesis[END_REF]; [START_REF] Rydén | Tumor specific VEGF-A and VEGFR2/KDR protein are co-expressed in breast cancer[END_REF]; [START_REF] Neufeld | Vascular endothelial growth factor (VEGF) and its receptors[END_REF]. These receptors regulate physiological as well as pathological angiogenesis. VEGFR2 signalling after its phosphorylation plays a crucial role in VEGF-A signalling mediating cell proliferation, survival, migration, and actin reorganization through a variety of signalling cascades including MEK/ERK and PI3-kinase/AKT [START_REF] Kerbel | Tumor angiogenesis[END_REF] [START_REF] Banerjee | Mechanisms of Disease: angiogenesis and the management of breast cancer[END_REF]. To our knowledge it is unknown whether the concomitant administration of tamoxifen with cytotoxic therapy might result in a greater anti-angiogenic activity.

Primary systemic antineoplastic therapy is an excellent method to explore the effects of treatment on tumor biology and represents the ideal model to assess drug interactions within tumoral microenvironment. In the present study, we evaluated VEGF and VGFR2 tissue expression before and after treatment in a series of breast cancer patients enrolled in a prospective randomized trial comparing single agent epirubicin versus epirubicin plus tamoxifen as primary systemic therapies [START_REF] Bottini | Cytotoxic and antiproliferative activity of the single agent epirubicin versus epirubicin plus tamoxifen as primary chemotherapy in human breast cancer: a single-institution phase III trial[END_REF]. The primary aim of the study was to evaluate the effect of study treatments on these two markers, secondary aims were: a) to correlate angiogenic markers with clinical and biological features, b) to correlate the changes of angiogenic markers with tumor response, c) to evaluate the prognostic role of VEGF and VEGFR2 assessed either at baseline or at post chemotherapy residual histology.

PATIENTS AND METHODS

Patients

Patients with T2-4 N0-1 breast cancer were recruited to a randomized trial comparing single agent epirubicin (EPI arm) versus epirubicin plus tamoxifen (EPI-TAM arm) as the primary systemic treatment [START_REF] Bottini | Cytotoxic and antiproliferative activity of the single agent epirubicin versus epirubicin plus tamoxifen as primary chemotherapy in human breast cancer: a single-institution phase III trial[END_REF]. They were required to have an Eastern Cooperative Oncology Group performance status < 2, adequate bone marrow reserve (WBC count > 3.5 X 109/L, platelets > 100 X 109/L, and haemoglobin > 10 g/dl), hepatic function (AST, ALT, bilirubin, and alkaline phosphatase levels <1.25 the upper limit of normal value) and renal function (serum creatinine <1.25 the upper limit of normal value). Patients were accrued from January 1997 to December 2001. The study was approved by the Institutional Investigations Committee. All patients gave written informed consent to the diagnostic procedures, the proposed treatment, and the biological evaluations. Twohundred and eleven patients were enrolled, 105 were randomized to receive epirubicin alone, and 106 were randomized to receive epirubicin plus tamoxifen. On first presentation, an incision biopsy was done on each patient and a small tissue sample (0.5-0.8 cm) was removed. Chemotherapy was started within 2 days of diagnosis. Patients in the EPI arm received 60 mg/m2 of epirubicin (Farmorubicina, Pharmacia, Milan, Italy) by slow i.v. push on days 1 and 2; whereas patients on the EPI-TAM arm received 60 mg/m2 of epirubicin by slow i.v. push on days 1 and 2 and 30 mg of tamoxifen (Kessar, Pharmacia) daily. Epirubicin injections were repeated every 21 days for three or four cycles before definitive surgery, whereas tamoxifen was given continuously until definitive surgery. All patients postoperatively received four cycles of the CMF regimen [i.v. cyclophosphamide (600 mg/m2), i.v. methotrexate (40 mg/m2), and i.v. 5-fluorouracil (600 mg/m2) on days 1 and 8, every 28 days]. Patients with estrogen receptor (ER) positive primary tumor in both treatment arms received tamoxifen (20 mg, i.e., lower than the primary dose) starting after surgery, up to progression or for a maximum of 5 years. The median follow up of patients was 53 months (range 13-95).

Treatment evaluation

Each month the size of the primary tumor and the size of the axillary lymph nodes, when appreciable, were measured by the same clinician using a caliper. Response was assessed before definitive surgery by the clinical measurement of the changes in the product of the two largest diameters recorded in two sequential evaluations. According to WHO Criteria, tumor progression was defined as an increase of at least 25% in tumor size; stable disease as an increase of less than 25% or a reduction less than 50%; partial response as tumor shrinkage greater than 50%; and complete response as the complete disappearance of any clinical sign of disease. Pathological complete response was defined as the complete absence of neoplastic cell either in the breast or in the axillary lymph nodes.

Surgery was performed after primary chemotherapy and clinical reassessment (quadrantectomy or modified radical mastectomy, together with axilllary node dissection, when indicated). If subjected to quadrantectomy, the patients underwent irradiation of residual breast.

Histopathologic and immunohistochemistry.

Tumour grade was evaluated using the Nottingham prognostic index [START_REF] Andrew | The Nottingham Prognostic Index for Invasive Carcinoma of the Breast[END_REF].

Immunohistochemichal evaluation was done on paraffin embedded tumor samples obtained at diagnosis and at definitive surgery. Bcl-2, p53, ER, PgR, HER2 and Ki67 staining were done at the Pathology Unit of Azienda Ospedialiera Istituti Ospitalieri of Cremona, Italy as described elsewhere [START_REF] Bottini | p53 but not bcl-2 immunostaining is predictive of poor complete response to primary chemotherapy in breast cancer patients[END_REF]. All staining was scored by counting the number of positively stained cells and expressed as a percentage of the total tumor cells (at least 1000) counted across several representative fields of the section using a standard light microscope equipped with a 10 x 10 square graticule. Reproducibility of counting was assessed by a second investigator rescoring 10 slides. The relative intensity of ER and PgR staining was assessed in a semiquantitative fashion, incorporating both the intensity and distribution of specific staining. A value (HSCORE) was derived from the sum of the percentages of positive-stained epithelial cells multiplied by the weighted intensity of staining. Specimens were deemed receptor positive if the HSCORE was greater than 100 [START_REF] Bottini | p53 but not bcl-2 immunostaining is predictive of poor complete response to primary chemotherapy in breast cancer patients[END_REF].

Immunohistochemistry for angiogenesis markers (HIF1α (hypoxia inducible factor 1-α) and Ca IX (carbonic anydrase IX), VEGF, VEGFR2) was performed on 5-µ sections of tissue microarray containing two 1-mm diameter cores taken from selected morphologically representative tumor regions from the incisional biopsy and from tumor remaining at definitive surgery. Quality control was assessed on each block by hematoxylin and eosin staining. The Envision HRP kit (Dako; Cambridgeshire, United Kingdom) system was used for subsequent visualization [START_REF] Generali | Phosphorylated ER , HIF-1 , and MAPK Signaling As Predictors of Primary Endocrine Treatment Response and Resistance in Patients With Breast Cancer[END_REF], [START_REF] Fox | CITED4 inhibits hypoxia-activated transcription in cancer cells, and its cytoplasmic location in breast cancer is associated with elevated expression of tumor cell hypoxia-inducible factor 1alpha[END_REF].

The first step consisted on an antigen-retrieval procedure, by heating a tissue section in a citrate buffer. The primary antibodies were: VEGF (VGI) (Oxford University) dilution 1 : 4, overnight incubation at room temperature, KDR (34a) (Oxford University) dilution 1 : 2, overnight incubation at room temperature. All sections had a negative control slide (no primary antibody) of an adjacent section to preclude nonspecific staining. Positive controls included breast carcinomas known to exhibit high levels of each marker. A single pathologist, blinded to patient outcome and to the origin of the samples, used a semiquantitative method. Intensity was semi-quantitively assessed: 0 (no staining), 1 (weak staining), 2 (moderate staining), or 3 (strong staining) for VEGF, VEGFR2 (Figure 2).

The immunohistochemical analysis of the angiogenesis markers was performed at the Weatherall Institute of Molecular Medicine of John Radcliffe Hospital of Oxford, UK.

Statistical methods

For statistical computations VEGF and VEGFR2 expression were considered either as discrete variables or dichotomized as follows: low levels, score 0,1, high levels, score 2,3.

Comparison of categorized variables was done by χ² or χ² for trend test when indicated.

Mann-Withney U test or Wilkoxon rank sum test for non parametric data were used, when indicated, to compare continuous variables. Spearman test for non parametric data was employed to analyze the correlations between variables. Disease free survival (DFS) and overall survival (OS) curves were estimated using the Kaplan-Meier method and compared using the log rank test.

Cox proportional hazard models were used to estimate hazard ratios (HR) for the role of VEGF and VEGF2R expression in predicting disease progression or death. Cox models and logistic regression models were used to assess the presence of heterogeneity in the effect of VEGF and VEGF2R expression on disease progression and disease response, respectively, in patient subgroups defined by treatment randomization, by including in the model the appropriate treatment/covariate interaction terms. This procedure is a test of the homogeneity of the hazard ratios associated with VEGF and VEGF2R between strata defined by treatment arm.

Missing data were dealt with by excluding patients from particular analyses if their files did not contain data on the required variables. All P values reported are the result of twosided tests. P values of less than 0.05 were considered to indicate statistical significance.

Statistical computation was done by SPSS for Windows software (version 16.0).

RESULTS

Patient characteristics.

191 out of 211 (90.52%) patients prospectively enrolled in the trial were evaluable for VEGF, and VEGFR2 expression. Patient's characteristics and expression levels of each parameter are shown in Table 1. 91 patients (47.6%) were randomized to the EPI arm and 100 patients (52.4%) were randomized to the EPI-TAM arm respectively. As outlined in the Consort diagram (Figure 1), VEGF was evaluated at baseline excision in 160 patients (83,7%) and VEGFR2 in 171 (89,5%), respectively. The corresponding evaluation at residual histology was 143 (74,8%) for VEGF and 137 (71,7 %) for VEGFR2.

Relationship between VEGF, and VEGFR2 expression and clinico-pathological variables.

In a univariate analysis, VEGF at baseline significantly correlated with HIF1-α (p<0.0001), whereas it failed to show any relationship with T-N status, grading, hormonal receptor status, p53, her-2, bcl2, Ki67 and CA IX (all p>0.05).(Table 2). Baseline VEGFR2 failed to show any significant relationship with the clinical and biological parameters considered (all p>0.05). A significant but weak relationship was found between VEGF and VEGFR2 expression (Spearman R 0.55, p<0.01).

VEGF and VEGFR2 expression and prediction of response to treatment.

Among the 191/192 patients with evaluable histological samples, 35 (18.3%) had a complete clinical response, 112 (58.6%) had a partial clinical response and 43 (22.5%) had no response. VEGF and VEGFR2 expression showed no association with overall clinical response and clinical complete response in all cases and or in each treatment arms (data not shown).

Effect of treatment on VEGF and VEGFR2 expression.

The effect of epirubicin plus/minus tamoxifen on neoangiogenesis marker expression was assessed in patients with matched tumor samples before and after treatment.

VEGF expression significantly decreased after treatment in overall patients (p<0.02) (data not shown). Stratifying patients according to treatment arm, VEGF expression significantly decreased in the EPI-TAM arm (p <0.01) whereas it did not significantly change in the EPI arm (p=0.54) (Table 3). As a consequence of these changes, while baseline VEGF expression did not differ between the two arms (chi square for trend p= 0.83), VEGF at residual histology was significantly higher in EPI patients as opposed to EPI-TAM patients (chi square for trend p= 0.035).

VEGFR2 expression, conversely, significantly increased after treatment in all patients (p<0.002) (data not shown), and in patients the EPI arm and EPI-TAM arm (p<0.02 and p=0.03 in, respectively) (Table3). Baseline VEGFR2 expression did not differ between the two arms (chi square for trend p=0.115), whereas VEGFR2 expression at residual histology was significantly higher in patients randomized in the EPI arm as opposed to those randomized in the EPI-TAM arm (p = 0.045).

A similar pattern of VEGF and VEGFR2 expression before and after treatments was also observed in patients with estrogen receptor positive tumors only. Unfortunately the number of patients with ER negative tumors is too few to assess whether in this subset there is a different behaviour of angiogenic markers or not.

Relationship between changes in VEGF and VEGFR2 expression and clinical response.

VEGF expression in overall patients decreased in 44 patients (37.3%) while no change or increase was observed in 74 patients out of 118 (62.7%).

As shown in Table 4, no significant difference was observed in the distribution of VEGF decrease according to disease response, in all patients. Similar results were observed in both treatment arms.

VEGFR2 expression in all patients decreased in 27 out of 121 patients (22.3%) whilst there was no change or an increase in 94 (77.7%). The proportion of patients showing decreased VEGF2R expression was greater in responding patients as opposed to nonresponders (p=0.02). This was statistically significant in the EPI-TAM arm (p=0.05) and not in the EPI arm (p=0.26) (Table 4), but no difference between the 2 arms was observed at the interaction test (p=0.54).

VEGF and VEGFR2 expression and disease outcome.

VEGF expression at baseline was significantly associated with a shorter disease free survival (p=0.047), whereas VEGFR2 showed no association (p=0.65) (Figure 3a and3b).

Conversely, high VEGFR2 expression (p=0.032) but not high VEGF (P=0.66) at post chemotherapy residual histology, was significantly predictive of poor disease free survival (Figures 3c and3d).

Dividing patients according to the treatment arms, the negative prognostic role of VEGF expression at baseline conditions was evident in patients randomized in the EPI arm but not in those randomized in the EPI-TAM arm (interaction test, p<0.05) (Figure 4). No interaction between treatment arms was observed for VEGFR2 at baseline and VEGF and VEGFR assessed at post chemotherapy residual histology.

Both VEGF and VEGFR2, either at baseline condition or at residual tumor histology, did not show any relationship with overall survival (data not shown).

Discussion

The present study aimed to investigate prospectively the effects of tamoxifen administered in association with epirubicin on angiogenic markers in a randomized trial of primary systemic therapy. Although the relationship between VEGF and tamoxifen has previously been explored in various other clinical settings, the effects of chemotherapy plus tamoxifen on angiogenesis and a comparison of these effects with chemotherapy alone have not been studied. Matched tumor samples from the neoadjuvant prospective trial of epirubicin plus minus tamoxifen has allowed an unparalleled opportunity to analyse changes in tissue expression of VEGF and VEGFR2.

These results clearly demonstrate that epirubicin administered at the conventional maximum dose of 120 mg/m2 every 21 days has no anti-angiogenic effect since there was no change in VEGF expression after treatment. However, the addition of tamoxifen to epirubicin resulted in a significant anti-angiogenic activity with a reduction in VEGF expression. Furthermore, VEGFR2 expression at post chemotherapy residual histology, were significantly lower in patients randomized to receive epirubicin plus tamoxifen compared with patients who received epirubicin alone. The contrasting changes in VEGF and VEGFR2 expression before and after treatment could constitute, a feed back loop with loss of VEGF ligand, relating in part to the increased VEGFR2 expression at residual tumor histology.

VEGF and VEGFR2 assessed at baseline were not predictive of subsequent clinical response and failed to demonstrate a predictive role of angiogenesis markers for tumor response to chemotherapy as already described [START_REF] Bottini | Changes in Microvessel Density As Assessed by CD34 Antibodies after Primary Chemotherapy in Human Breast Cancer[END_REF][24] [START_REF] Baar | A vasculature-targeting regimen of preoperative docetaxel with or without bevacizumab for locally advanced breast cancer: impact on angiogenic biomarkers[END_REF]. Treatment-induced-decrease in tissue VEGF expression also did not correlate with tumor response either in overall cases or stratifying patients according to treatment arm. However, a decrease in VEGFR2 was more frequently associated with disease response than no change or increase. In this trial the addition of tamoxifen to epirubicin resulted in a greater antitumor activity than epirubicin alone [START_REF] Bottini | Cytotoxic and antiproliferative activity of the single agent epirubicin versus epirubicin plus tamoxifen as primary chemotherapy in human breast cancer: a single-institution phase III trial[END_REF], and it is possible that these results at least in part are due to a greater anti-angiogenic activity of the combination therapy.

VEGF expression has been repeatedly correlated with poor outcome of breast cancer patients. Our data confirm previously published papers [26] [27] [28] [29] [30] [31] [32]

since high expression of VEGF was predictive of shorter DFS. VEGF expression at residual histology, however, was not associated with prognosis. Conversely, VEGFR2 expression at post treatment residual histology was significantly associated with short DFS but did not have any prognostic effect when assessed at baseline. The different behavior of these two markers before and after treatment (i.e a decreasing trend for VEGF and increasing trend for VEGFR) might account for the different prognostic significance of VEGF and VEGFR2 observed at baseline and post treatment.

Interestingly, dividing patients according to treatment arm the negative prognostic value of elevated VEGF expression at baseline was evident in the EPI patients but not in the EPI-TAM patients. These data suggest the reduction in VEGF induced by tamoxifen may have modified tumour growth in those randomized to the EPI-TAM arm. Conversely no interaction with the treatment arm was observed as far VEGFR2 expression either at baseline or after treatment. It is reported that tamoxifen has both an estrogen and an estrogen-independent effect on angiogenesis. Tamoxifen inhibits angiogenesis in a chick egg chorioallantoic membrane model. This inhibition is not altered by the presence of excess estrogens, suggesting that the mechanism is independent of tamoxifen's effect on the estrogen receptor [START_REF] Haran | Tamoxifen enhances cell death in implanted MCF7 breast cancer by inhibiting endothelium growth[END_REF]. This interesting issue could not be explored in our series, since, due to the low number of patients with ER negative tumors, we were unable to evaluate the anti-angiogenic effect of tamoxifen in this subset.

Although a trend with overall survival was observed, due to the low number of events observed in this trial, both markers failed to be significantly associated with overall survival.

In conclusion, the addition of tamoxifen to epirubicin resulted in an anti-angiogeneic activity and supports, at least in part, the greater activity of the combination observed in the clinical trial [START_REF] Bottini | Cytotoxic and antiproliferative activity of the single agent epirubicin versus epirubicin plus tamoxifen as primary chemotherapy in human breast cancer: a single-institution phase III trial[END_REF]. This result seems to be opposite to the INT100 data showing a negative interaction of tamoxifen with chemotherapy. The different chemotherapy employed in the two studies, FEC vs single agent epirubicin, the greater epirubicin dose employed in our study and the possible interpherence of TAM in the metabolism of cyclophosphamide in the INT100 study could have accounted for these discrepancies.

Nevertheless, this study provides further support for a link between VEGF and the estrogen receptor pathways suggesting that further exploration of the interaction between chemotherapy and endocrine therapy administered concomitantly is warranted. It will be of interest to evaluate whether aromatase inhibitors have similar effects.

Legends to Figures

Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure 1: Consort diagram

Figure 1 .

 1 Figure 1. Consort diagram

Table 1

 1 Patient characteristics and biological markers at baseline

	Patients	91	100
		Epi	Epi + Tam
	Clinical Characteristics

Table 2 .

 2 Distribution of prognostic variables and immuno-histochemical features according to VEGF and VEGFR2 expression

	VEGF	0	1	2	3	p
	T₁₋₂	26/32(81.25%)	32/39(82.05%)	27/37(72.97%)	39/52(75.00%)	0.358
	T₃₋₄	6/32(18.75%)	7/39(17.95%)	10/37(27.03%)	13/52(25.00%)
	N₀	20/32 (62.50%)	25/39(64.10%)	15/37(40.54%)	29/52(55.77%)	0.282
	N₁	12/32(37.50%)	14/39(35.90%)	22/37(59.46%)	23/52(44.23%)
	G2	7/32(21.88%)	16/39(41.03%)	5/36(13.89%)	14/50(28.00%)	0.771
	G3	25/32 (78.13%)	23/39(58.97%)	31/36(86.11%)	36/50(72.00%)
	ER -	4/32 (12.50%)	8/39 (20.51%)	6/37 (16.22%)	12/51(23.53%)	0.296
	ER +	28/32(87.50%)	31/39(79.49%)	31/37(83.78%)	39/51(76.47%)
	PgR -	15/32(46.88%)	23/39(58.97%)	15/37(40.54%)	27/51(52.94%)	0.981
	PgR +	17/32(53.13%)	16/39(41.03%)	22/37(59.46%)	24/51(47.06%)
	p53	20/32(62.50%)	18/39(46.15%)	18/36(50.00%)	23/52(44.23%)	0.176
	her-2	6/32(18.75%)	10/39(25.64%)	9/37(24.32%)	16/52(30.77%)	0.256
	bcl2	24/31(77.42%)	29/39(74.36%)	32/37(86.49%)	34/52(65.38%)	0.304
	HIF1-α	20/30(66.67%)	27/39(69.23%)	33/37(89.19%)	50/52(96.15%)	0.00006
	CA IX	6/28 (21.43%)	7/38 (18.42%)	10/36(27.78%)	15/50(30.00%)	0.239
	Topo2A	1/21(4.76%)	3/23(13.04%)	16/27(59.26%)	27/46(58.70%)	0.000001
	Ki67	21.53	22.05	24.62	19.44	0.64
	95 % IC	(16.12 -26.94)	(15.40 -28.70)	(17.45 -31.79)	(14.89 -23.99)
	VEGFR2	0	1	2	3	p
	T₁₋₂	11/13 (84.62%)	45/57 (78.95%)	30/40 (75%)	44/61 (72.13%)	0.257
	T₃₋₄	2/13 (15.38%)	12/57 (21.05%)	10/40 (25%)	17/61 (27.87%)
	N₀	7/13 (53.85%)	34/57 (59.65%)	23/40 (57.5%)	34/61 (55.74%)	0.828
	N₁	6/13 (46.15%)	23/57 (40.35%)	17/40 (42.5%)	27/61 (44.26%)
	G2	5/13 (38.46%)	11/53 (20.75%)	10/40 (25%)	17/61 (27.87%)	0.91
	G3	8/13 (61.54%)	42/53 (79.25%)	30/40 (75%)	44/61 (72.13%)
	ER -	2/13 (15.38%)	13/57 (22.81%)	8/40 (20%)	10/60 (16.67%)	0.619
	ER +	11/13 (84.62%)	44/57 (77.19%)	32/40 (80%)	50/60 (83.33%)
	PgR -	5/13 (38.46%)	31/57 (54.39%)	24/40 (60%)	27/60 (45%)	0.729
	PgR +	8/13 (61.54%)	26/57 (45.61%)	16/40 (40%)	33/60 (55%)
	p53	7/13 (53.85%)	36/57 (63.16%)	14/39 (35.9%)	28/61 (45.9%)	0.089
	her-2	3/13 (23.08%)	15/57 (26.32%)	10/40 (25%)	17/61 (27.87%)	0.754
	bcl2	10/13 (76.92%)	40/56 (71.43%)	30/40 (75%)	46/61 (75.41%)	0.782
	HIF1-α	7/11 (63.64%)	44/54 (81.48%)	35/39 (89.74%)	51/60 (85%)	0.171
	CA IX	4/12 (33.33%)	12/53 (22.64%)	14/39 (35.9%)	10/58 (17.24%)	0.327
	Topo2A	1/8 (12.5%)	6/32 (18.75%)	14/30(46.67%)	27/54 (50%)	0.002
	Ki67	25	22.82	21.55	21.72	0.29
	95 % IC	(8.19 -41.80)	(18.34 -27.30)	(14.04 -29.05)	(17.34 -26.10)

Table 3 .

 3 Effect of treatment on VEGF and VEGFR2 expression

	EPI				EPI-TAM
	Baseline	Post-treatment	p *	Baseline	Post-treatment	p *

Table 4 .

 4 Relationship between changes in VEGF and VEGFR2 expression and clinical response.

		NR	PR	CR	p
	VEGF decrease
	Overall 11/29 (37.93%) 26/76 (34.21%) 7/13 (53.85%) 0.521
	EPI arm	7/20 (35.0%) 12/37 (32.43%)	3/6 (50.0%) 0.699
	EPI-TAM arm	4/9 (44.44%) 14/39 (35.90%) 4/7 (57.14%) 0.685
	VEGFR2 decrease
	Overall	2/30 (6.67%) 20/76 (26.32%) 5/15 (33.33%) 0.020
	EPI arm*	2/21 (9.52%)	9/36 (25.0%)	1/5 (20.0%) 0.260
	EPI-TAM arm*	0/9 (0%) 11/40 (27.50%) 4/10 (40.0%) 0.05
	*Interaction test p=0.54


R Wawryk
email: r.wawryk@int.pan.wroc.pl

J Mucha

H Misiorek

Z Henkie

Thermal conductivity of USb 2 and UBi 2 single crystals

Keywords: anisotropy, thermal transport user supplied): uranium dipnictides, thermal conductivity, anisotropy uranium dipnictides, thermal conductivity, anisotropy

The thermal conductivity (κ) of single crystals of tetragonal uniaxial antiferromagnets USb 2 (T N = 202 K) and UBi 2 (T N = 180.8 K) has been measured along the a-axis (κ a ) over the temperature range from 0.5 K to 300 K and along the c-axis (κ c ) from 0.5 K to 70 K. The as-grown samples have RRR values of about 500 -600 and 100 -150 for UBi 2 and USb 2 , respectively. The anisotropy of the thermal conductivity (κ a (T)/κ c (T) ~ 5) and the low-T Lorenz ratios are discuss in relation to Fermi surface topology for both compounds.

§ 1. INTRODUCTION

Uranium compounds display a wide variety of properties that are related to the behaviour of 5f electrons. These phenomena vary from one uranium compound to another depending on a number of different factors including crystal structure, crystal electric field (CEF) effects, electronic configuration, etc. Thus, in order to gain a better understanding of the intrinsic nature of 5f electron systems, it is useful to examine the evolution of physical properties in series of compounds. For this reason, the series of uranium dipnictides UX 2 (X=P, As, Sb, or Bi) has received considerable attention. These materials become antiferromagnets below room temperature and have ordered magnetic moments and effective paramagnetic moments as large as about 2 µ B and 3.1 µ B [1,2,3], respectively. Earlier measurements of the magnetic properties of UX 2 compounds were interpreted in terms of localised 5f moments of U 4+ cations in the crystal electric field of the surrounding pnictogen anions [4]. More recent data addressing hyperfine interactions in the UX 2 series, studied by Mössbauer spectroscopy [5], were interpreted in terms of hybridisation between the 5f and conduction electrons. The later approach was consistent with the results from angleresolved photoemission spectroscopy for USb 2 [6], which revealed a narrow, 5fderived band with small dispersion at the Fermi energy. The topology of the Fermi surface of the UX 2 series, when determined either experimentally by de Haas-van Alphen (dHvA) and Shubnikov-de Haas (SdH) measurements [7,8] or by means of ab initio calculations [9], both predict a significant two-dimensional character for the electronic properties. The two dimensional behaviour is thought to vary through the entire series and to reach the clearest difference between USb 2 and UBi 2 .

The uranium diantimonide and dibismuthide compounds crystallize in the tetragonal structure of anti-Cu 2 Sb type (space group P4/nmm or reveal antiferromagnetic ordering of the uranium magnetic moments, with values of 1.88 µ B for USb 2 and 2.1 µ B for UBi 2 below the Néel temperatures, T N of 202 K and 180.8 K, respectively. The moments are parallel to the c-axis and create ferromagnetic planes which alternate perpendicular to the c-axis with the sequence (↑↓↓↑) in USb 2 and (↑↓) in UBi 2 . This indicates a higher flattening of magnetic Brillouin zone (MBZ) in USb 2 [8] as compared to that in UBi 2 . Thus the ratios k c /k a, where k c is the wave vector of the MBZ along the c-axis and k a is the wave vector along the a-axis, for USb 2 and UBi 2 are equal to 0.24 and 0.5, respectively.

The Fermi surfaces (FS) and paramagnetic Brillouin zone (PBZ) are basically the same in USb 2 and UBi 2 . Roughly speaking the FS consists of two types of branches. The first type of FS branch is a spherical surface and the second is a slightly corrugated cylinder parallel to the c-axis [8]. Magnetic ordering in UBi 2 does not change the BZ and hence the FS. The magnetic ordering in USb 2 reduces twice the MBZ and the spherical FS crosses the MBZ and transforms to the two cylinders. The same happens with the corrugated cylinder, which also transforms to the two cylinders. Thus the flattening of the MBZ in USb 2 enhances the two-dimensionality. Consequently, its resistivity anisotropy, which is several times lower than that of UBi 2 at room temperature, becomes several times higher at low-T [11]. Here we extend the transport examination to studies of the heat conductivity in USb 2 and UBi 2 single crystals. § 2. EXPERIMENTAL DETAILS The "molten metal solution" method was used to synthesize plate-like single crystals of USb 2 and UBi 2 [12]. From the as-grown crystals, it was possible to prepare rectangular specimens with typical dimensions of about 6.5×1.1×1.3 mm 3 (UBi 2 ) and 4.0×1.6×0.7 mm 3 (USb 2 ) along the a×b×c axes for measurements of the thermal conductivity along the a-axis, κ a (T), and specimens of about 0.7×0.7×0.6 mm 3 (UBi 2 ) and 0.6×0.4×1.1 mm 3 (USb 2 ) along the a×b×c axes for measurements of the thermal conductivity along the c-axis, κ c (T). The thermal conductivity was determined with total uncertainty less then 2.5 % by the stationary heat flux method in the temperature range 5 -300 K and setting the temperature gradient along the specimen's length from the range 0.1 -0.5 K. The experimental setup and the measurement procedure have been described in detail in ref. [13,14]. In the range 0.5 K -80 K the κ(T) was measured by the modified set-up for the thermoelectric power measurements [15]. § 3. RESULTS OF THE THERMAL CONDUCTIVITY OF USb 2 AND UBi 2

The a-axis thermal conductivity data for USb 2 (specimen 2S a ) and UBi 2 (specimen 2B a ) measured in the temperature range 5 K < T < 300 K are presented in Fig. 1. The values of κ a (16.7 and 17.8 W m -1 K -1 at 300 K for USb 2 and UBi 2, respectively) decrease with temperature down to T N where κ a (T) traces a shallow valley. Visible change of slope of the dκ a /dT vs. T dependence at T N , as shown in the inset for both compounds, relates this behaviour to the antiferromagnetic ordering temperature. The straight lines following the experimental data and constrained to cross at T N are used to highlight the anomaly. Though the κ a (T) anomalies at T N are detectable, they are drastically weaker than the anomalies seen in the electrical resistivities along the a-axes shown by two lower curves in the inset, which are due to conduction electron scattering by phonons and fluctuations of magnetic moments [11]. UBi 2 (specimen 3B c ) and USb 2 (specimen 3S c ) measured along the caxis. The inset presents the derivative dκ a /dT for the specimen's 2B a of UBi 2 and 2S a of USb 2 around the magnetic transition temperature T N and compares them with the ρ(T) data [11] shown by two lower curves.

For temperature below T N , κ a (T) passes through a maximum κ a max at temperature T max , which is dependant on the purity of the specimen, and finally proceeds to zero with decreasing T. In order to better demonstrate the impurity effect, we have determined the low-T behaviour of κ a (T) for several additional specimens prepared from crystals from different batches of USb 2 , -specimen (spec.) 1S a -and UBi 2 -spec. 1B a . Finally the specimens prepared for measurements of the thermal conductivity along the c-axis are denoted 3S c and 3B c , respectively. The electrical resistivity, ρ, of each specimen was also examined and residual resistivity ρ 0 , for T → 0, as well as the residual resistivity ratio, RRR = ρ 300 /ρ 4.2 were determined; here ρ 300 and ρ 4.2 are the resistivities at 300 K and at 4.2 K, respectively. It is well know that higher RRR values often indicate high purity. We can characterize specimens with sets of numbers (RRR):(κ max in W m -1 K -1 ):(T max in K) which are for USb 2 specimens: 1S a -( 97 We use a metallic approach [START_REF] Berman | Thermal Conduction in Solids[END_REF] to analyse the thermal transport in USb 2 and UBi 2 and assume validity of the Wiedemann-Franz-Lorenz law and electronic thermal resistivity W e = ρ/LT where L is the Lorenz ratio which, under certain conditions, is equal to the standard Lorenz number L 0 = 2.45 × 10 -8 W Ω K -2 . We assume that for the lowest temperatures W e = 1/κ e exp is the sum of two components. One is the thermal resistivity due to the scattering of electrons by defects, (1) The power n for scattering by phonons is expected to be 2 at low temperatures and a plot of W e T against T 3 should give something like a straight line. However, deviation from that expected behaviour is often observed. As a matter of fact, this is the case for USb 2 and UBi 2 . The W e T data for four specimens of these compounds, plotted against T 2.73 in the inset of Fig. 2, give something like a straight line. More detailed analysis shows that the widest temperature range of linearity is obtained for n selected for particular samples. It is observed that n decreases from 1.8 to 1.5 when W e increases from the lowest values, as observed for the 2S a specimen, to the highest values for 3B c . This is illustrated in main part of Fig. 2 where log (WT -A) data are plotted vs. log T for A values assuring the best linear behaviours in the temperature range of about 10 ─ 50 K. These values of A point out κ el = TL 0 /ρ approaching or even exceeding κ exp which, however, is expected to be higher then κ el . Such higher values than that expected for A obtained in this approach results presumably from lattice contribution, not negligible in the considered temperature range. This is discussed in the § 4.

W e d = A/T = ρ 0 /L 0 T. The log T 1.0 1.2 1.4 1.6 1.8 2.0 log (WT-A) -3.0 -2.0 -1.0 0.0 1.0 1B a 2B a 1S a 2S a 3B c 3S c WT = A + BT n n = 2.5 n = 2.8
Finally we have assumed that the heat transport at the lowest-T (0.5 -2 K) is entirely controlled by defects. Then κ e exp = T/A comes from eq. ( 1), allowing us to determine A values as equal to reciprocal of slope of the κ exp vs. T dependence. Such values of A determined at the lowest temperature point to L 0 (in 10 -8 W Ω K -2 ) equal to 4.1 and 72 for USb 2 along the a-and c-axis and 0.95 and 10.6 for UBi 2 also along the a-and c-axis, respectively. The values of κ el = TL 0 /ρ calculated for all these Lorenz numbers are compared in Fig. 3 to corresponding experimental data of the total thermal conductivity in USb 2 and UBi 2 and lattice conductivity, κ ph + κ mag . It is seen that beyond the lowest temperatures the lattice conductivity clearly exceeds the electronic conductivity in the all cases presented in Fig. 3 . Electronic and phonon/magnon contributions to the total thermal conductivity of USb 2 and UBi 2 single crystals in the a and c directions. The following values of L 0 were taken to evaluation of electronic contribution: 4.1×10 -8 , 9.5 ×10 -9 , 72 ×10 -8 , and 10.6 ×10 -8 W Ω K -2 for the samples 2S a , 2B a , 3S c , and 3B c , respectively. § 4. DISCUSSION AND SUMMARY Aoki et al. [8] revealed by dHvA and SdH experiments that the spherical α and cylindrical β branches of the FS in UBi 2 occupies 9.9 % and 4.8 % BZ volume, respectively. The magnetic unit cell contains two molecules of UBi 2 and it is expected In order to explore a possible correlation between the determined Lorenz number and the topology of the FS, the latter will be characterized, for simplicity, by the difference between its cross section as appropriate for the direction of the thermal conductivity examination and the cross section of spherical surface of a simple metal. For the a-axis thermal conductivity of UBi 2 , the essential factor will be the shape of the FS cross sections by planes perpendicular k c . All cylindrical FS branches (2β) are parallel to k c and their cross sections are circles, like the cross section of the spherical α FS branch. For this case L 0 = 0.95 ×10 -8 W Ω K -2 was found. This value is lower than the standard L 0 value and can be compared with L 0 = 1×10 -8 W Ω K -2 , as found for ferromagnetic nickel [START_REF] White | [END_REF].

By comparison to UBi 2 , there are 20 % more carriers in USb 2 , which form two times as many FS branches. In spite of this difference, the FS are all corrugated cylinders parallel to k c . Thus the cross sections of the USb 2 FS with planes perpendicular to k c are composed of six small circles (high curvature) so it is farther from spherical surface of simple metal than in the case of UBi 2 . For thermal conductivity along the a-axis of USb 2 the value of L 0 = 4.1×10 -8 W Ω K -2 has been determined. This is in turn higher than the standard L 0 value.

The cylindrical FS branches in UBi 2 can contribute to the transport coefficient along the c-axis only as a result of the corrugations in these branches. Any of the cross section with planes parallel to the k c axis does not resemble cross sections of a spherical surface. The only circular FS cross section may come from the spherical FS α branch, but in spite of this the overall similarity to a spherical cross section is reduced and an increased value of L 0 = 10.6×10 -8 W Ω K -2 was determined for the caxis thermal conductivity in UBi 2 . Finally, L 0 increases up to 72×10 -8 W Ω K -2 in the case of USb 2 for the c-axis where all FS branches are corrugated cylinders parallel to k c .

Based on a simple theoretical treatment, (see ref. [START_REF] Berman | Thermal Conduction in Solids[END_REF]) the lattice thermal conductivity should be about one-third of the electronic component at high temperature if only phonon-electron scattering is important. Calculations done for copper show that the lattice contribution gives a very small fraction to the total thermal conductivity at high and low temperatures. It will suffice to remark that at high temperature the ratio κ/σ T is close to the standard value L 0 . This is obviously not the case of our compounds, which show layered chemical structure and are low carrier density systems. Instead, we make a comparison to the thermal conductivity of YBa 2 Cu 3 O 7-y as examined by Hagen et al. [18]. This material is a low carrier density system with a layered chemical structure and similar thermal conductivity anisotropy in the normal state where κ ab /κ c ≈ 5 -6. For this compound, it was found that the standard L 0 indicates that the electronic contribution is 55 % of the total κ ab and 1 -2 % of κ c. Thus, it was argued that detailed analysis of the latter casts serious doubt about the validity of applying the Bloch-Boltzmann theory to the transport behaviour along the c-axis. In the case of USb 2 and UBi 2 the low temperature L 0 behaviour points towards lower contribution of the electronic thermal conductivity to the total conductivity as compared to that determined by standard L 0 in the case of the a-axis thermal conductivity of UBi 2 -10 % -and the higher one in the case κ el for the a-axis of USb 2 -36 % -and κ el for the c-axis of USb 2 -38 % -and UBi 2 -23 %. In the latter case the κ el contribution is nearly 30 times higher than that calculated for standard L 0 .

Electronic thermal resistivity due to the scattering of electrons by phonons should show quadratic temperature dependence at low-T. In our case, we observe a thermal resistivity term W ∼ T n with n decreasing from 1.8 for specimens with the lowest thermal resistivities to 1.5 for specimens with the highest resistivities, on this temperature range (see Fig. 2). The total thermal resistivity seems to be the main parameter controlling the decrease of n. Presumably, n is an effective parameter, which decreases with an additional scattering mechanism and gives an additional contribution to the total thermal resistivity. At these temperatures, we first consider a lattice term for which the temperature dependant contribution varies like T ─2 and increases with increasing ρ 0 . The lattice thermal conductivity would be about onethird of the electronic component at high temperature if only phonon-electron scattering was important in determining both contributions [START_REF] Berman | Thermal Conduction in Solids[END_REF]. As the examined compounds are antiferromagnets, we should also consider the possible scattering of electrons by magnetic moments, for which the contribution to the electronic thermal conductivity seems to be well described by κ el (T) = L 0 T/ρ(T). Therefore it is surprising that the κ a (T) data for both compounds show weak anomalies at T N as compared with ρ a (T) (inset to Fig. 1). There also is a weak and weakly temperature dependent anisotropy of the thermal conductivity. The ratio κ a (T)/κ c (T) at 80 K (maximum of ρ c (T) for USb 2 ) for the best specimens of USb 2 and UBi 2 is equal to 5.1 and 5.4, respectively, and does not change by more than 30 % with temperature down to 20 K. This result can be compared with the anisotropy of the electrical resistivity, as shown for the both compounds by the lower inset in Fig. 4. To show high specimen sensitivity of another electronic transport coefficient, we present thermoelectric power for specimen 2S a showing RRR 150 and compare it with the previously examined specimen 3S a [11]. The increase of RRR enhances the phonon drag peak but lowers the thermoelectric power in the paramagnetic region. One might think that the weak manifestation of the aforementioned peculiarities of electrical resistivity behaviour in the thermal conductivity are possibly due to their signatures being overwhelmed by another scattering mechanism which dominates the total thermal conductivity. temperature for USb 2 single crystals along the a-axis. The squares denote the data for 2S a sample, while the dotted line represents the data for 3S a sample measured earlier [11]. Curve 1 in the upper inset displays the equation: S(T) = aT + b/T with the fitting parameters a and b equal to 0.205 µV K -2 and 285 µV, respectively. The straight line 2 and curve 3 present the diffusion and phonon/magnon drag contributions to the total thermoelectric power, respectively. The lower inset presents the ratio (ρ c -ρ 0 )/(ρ a -ρ 0 ) for USb 2 and UBi 2 single crystals taken from ref. [11].

We conclude that the electronic thermal transport analysis at low temperatures leads to an estimation of various values of the Lorenz number L 0 and gives a more realistic electronic contribution κ el to the total thermal conductivity for both compounds UBi 2 and USb 2 . The lack of the spherical FS branches in USb 2 and flat MBZ influence on the anisotropy of the electronic thermal conductivity κ el a /κ el c that is not so large as in electrical resistivity and thermoelectric power. The anisotropy of the total thermal conductivity is comparable with the other layered structures e.g. mentioned above YBCO.

Figure 1 .

 1 Figure 1. Thermal conductivity of UBi 2 (specimens 1B a and 2B a ), and USb 2 (specimens 1S a and 2S a ) single crystals vs. temperature measured along the a-axis and κ(T) data for

  ):(58):(22); 2S a -(150):(92):(24); 3S c -(2.4):(12.8):(26) and for UBi 2 specimens: 1B a -(498):(26):(24); 2B a -(598):(41):(20); 3B c -(130):(8.7):(18.5).

Figure 2 .

 2 Figure 2. The plot log (WT -A) vs. log T for spec. 1B a , 2B a , 3B c , 1S a , 2S a , and 3S c . The solid lines represent eq. WT = A + BT n for n = 2.5 and 2.8. The inset shows WT as a function T 2.73 .

  component is the resistivity due to the scattering of electrons by phonons/magnons W e p = BT n as the examined compounds undergo to antiferromagnetic transitions. Thus, W e = 1/κ e exp = W e d + W e p = A/T + BT n .

Figure 3 .

 3 Figure 3.Electronic and phonon/magnon contributions to the total thermal conductivity of USb 2 and UBi 2 single crystals in the a and c directions. The following values of L 0 were taken to evaluation of electronic contribution: 4.1×10 -8 , 9.5 ×10 -9 , 72 ×10 -8 , and 10.6 ×10 -8 W Ω K -2 for the samples 2S a , 2B a , 3S c , and 3B c , respectively.

Figure 4 .

 4 Figure 4. Thermoelectric power vs. temperature for USb 2 single crystals along the a-axis. The squares denote the data for 2S a sample, while the dotted line represents the data for 3S a sample measured earlier[11]. Curve 1 in the upper inset displays the equation: S(T) = aT + b/T with the fitting parameters a and b equal to 0.205 µV K -2 and 285 µV, respectively. The straight line 2 and curve 3 present the diffusion and phonon/magnon drag contributions to the total thermoelectric power, respectively. The lower inset

  to be a compensated metal. Aoki et al. assumed that the spherical α branch is an electron FS and it is compensated with two cylindrical β branches that are the hole Fermi surfaces. This results in either electron and hole carrier concentrations in UBi 2 to be equal 1.1 × 10 21 carrier/cm 3 (2×0.099 carrier/4.44 2 ×8.91×10 -24 cm 3 ) which is approximately 75 times lower than the concentration of conduction electrons in copper. The magnetic unit cell of USb 2 contains four molecules of USb 2 and is also expected to be a compensated metal. Its two cylindrical FS, α and δ, originating from spherical FS in the PBZ occupies 16.8 % and 5.4 % of volumes of the magnetic BZ, respectively, and there are two electron FS branches. The two other cylindrical FS type branches, ε and γ, respectively, originating from cylindrical FS branch in PBZ, which occupy 7.9 % and 3.3 % of the MBZ. Aoki et al. assumed that 2ε and 2γ are hole FS branches occupying 22 % of MBZ compensate electron α and δ FS branches thus occupying also 22 % of the MBZ. This results in either electron and hole carrier concentrations in USb 2 to be equal 1.4 × 10 21 carrier/cm 3 (2×0.22 carrier/2×4.27 2 ×8.75×10 -24 cm 3 ) which is approximately 61 times lower than the concentration of conduction electrons in copper.

	F o r
	P
	e
	e r
	R
	e
	κ tot κ el κ ph + κ mag	v i e w
		O n l
		y

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

ACKNOWLEDGEMENTS

The authors are indebted to Dr. R.E. Baumbach for his assistance at the manuscript editing.


Julien Bruchon
email: bruchon@emse.fr

Daniel Pino Muñoz

François Valdivieso

Sylvain Drapier

Guillaume Pacquaut

-François Valdivieso

3D simulation of the matter transport by surface diffusion within a level-set context

Keywords: Surface diffusion, level-set method, curvature, Laplace-Beltrami operator, sintering process MOTS-CLÉS : Diffusion surfacique, méthode level-set, courbure, opérateur de Laplace-Beltrami, procédé de frittage

Introduction

The general context of this work is the direct simulation of the formation as well as of the subsequent evolution of a powder compact microstructure during a sintering process. Such a process consists in consolidating a metallic or ceramic powder compact, by thermally activating diffusion phenomena to a temperature lower than the melting temperature of the material. When considering a pressure-free sintering process (see [START_REF] Ashby | A first report on sintering diagrams[END_REF] for further details), the surface diffusion is the path of matter transport which is first activated (see Figure 1). Indeed, in order to minimize the surface energy, the matter flows over the grain free boundary towards the maximum (in absolute value) of curvature. The location of the grain centres remains unchanged under this transport path, that is, the porosity does not decrease (no shrinkage). Due to this first diffusion, stresses appear in the material, giving rise to the boundary diffusion and to the volume diffusion, with a subsequent decrease of the global porosity. This paper deals with the matter transport by the surface diffusion from the grain surface towards the neck formed between two grains. This situation does not involve the mechanical response of the material, since the diffusion is described only by a geometrical criterion: the surface Laplacian of the curvature. An outline of the paper is as follows. The physical description of the surface diffusion is first detailed in Section 2. Section 3 is devoted to the level-set strategy developed to solve the surface diffusion equations. Finally, numerical simulations are shown in Section 4.

Surface diffusion equations

Let Ω g be a set of grains, and let S g be the free boundary of this set (see Figure 1). Regarding the literature about the ceramic sintering process modelling [START_REF] Ashby | A first report on sintering diagrams[END_REF][START_REF] Bouvard | Deformation of interparticle necks by Diffusion-Controlled Creep[END_REF], the matter flow by surface diffusion is characterized by a flux j s along the free surface S g , driven by chemical potential gradients. In turn, these gradients depend on the gradient of the mean curvature K:

j s = δ s D s γ s kT ∇ s K [1]
where T is the absolute temperature, k is the Boltzmann's constant, D s is the surface diffusion coefficient, δ s is the thickness in which the diffusion occurs and γ s is the sur-face free energy. In expression [1], operator ∇ s denotes the surface gradient, defined as the tangential component of the gradient

∇ s K = ∇ K -(∇ K • n)n
where n is the outward-pointing unit vector normal to the free surface S g . The surface gradient can easily be rewritten by introducing P , the projection matrix onto the plane tangent to S g :

∇ s K = P ∇ K with P = I -n ⊗ n [2]
where I is the identity matrix.

The surface flux j s results in the deposition or removal of material, which gives rise to a displacement rate, assumed to be normal to the surface. This surface diffusion velocity is then written v s = v n n. The mass balance between, on one hand the matter which flows over the free surface S g , and on the other hand the normal displacement of this surface, leads to the following expression of the surface velocity [START_REF] Bouvard | Deformation of interparticle necks by Diffusion-Controlled Creep[END_REF]:

v s = Ω(∇ s •j s )n = C 0 (∆ s K)n [3]
where Ω is the molar volume of the material and with C 0 = δ s D s Ωγ s kT which is assumed to be constant in this approach. It means that the temperature is assumed to be at least uniform in the computational volume. In fact, all the simulations presented in the following correspond to isothermal cases. The operator ∆ s , the "surface Laplacian" operator is the so-called Laplace-Beltrami operator.

Finally, let us remark that the surface velocity given by Equation [3] preserves the volume. Indeed, if |Ω g | denotes the grain volume measure, its variation during the process is expressed by

d|Ω g | dt = S g v s • n dS = S g ∇ s •j s dS
Since the flux is continuous and S g is a closed surface, this variation vanishes, and d|Ω g |/dt = 0. Consequently, the grain volume remains constant. This point, as well as the absence of shrinkage, will be a criterion to assert the relevancy of our simulations.

Level-set formulation

Equation [3] is discretized by using a finite element method. The simulation of the matter transport by surface diffusion requires only a surface description of the grains. However, the final objective of this work is to take into account the different paths of the matter transport, including in volume. That is why the numerical strategy proposed in this paper is based on a volume description of the grains, i.e., the whole computational domain is meshed. Furthermore, in order to deal implicitly with the topological changes arising in the grain surface (formation of the neck), an Eulerian approach is chosen to describe the grain surface changes. Hence, the computational domain Ω ⊂ IR 3 (the unit cube) is made up of two parts: the grain set S g , and the surrounding media Ω a (the air, for example), Ω = Ω g ∪ Ω a . An additional function, the level-set function, φ : Ω → IR, is then required to describe the grain free boundary. This point is detailed in the next section.

Description of the grains

For any time t, and any point x ∈ Ω, the value of the level-set function φ(x, t) is positive if x ∈ Ω a (t), and is negative if x ∈ Ω g (t). Consequently the grain free boundary is defined by the zero isosurface of φ:

S g (t) = {x ∈ Ω ; φ(x, t) = 0}
This level-set function φ, which has to be smooth enough in the vicinity of S g , is defined as following, according to [START_REF] Coupez | Réinitialisation convective et locale des fonctions Level Set pour le mouvement de surfaces et d'interfaces[END_REF][START_REF] Bernacki | Level set framework for the numerical modelling of primary recrystallization in polycrystalline materials[END_REF]. At the initial time t = 0, φ is initialized by the grain initial state φ 0 , φ(x, 0) = φ 0 (x), with

φ 0 (x) =            2E π sin( πd(x, S g (0)) 2E ) if |d(x, S g (0))| ≤ E + 2E π if d(x, S g (0)) ≥ E - 2E π if d(x, S g (0)) ≤ -E [4]
where d(x, S g (0)) is the signed distance from the point x to the initial free surface S g (0). The numerical parameter E, called the interface width, is chosen equal to 1.5 × mesh size in the presented simulations. Note that, since sin x ≈ x when |x| 1, φ 0 (x) ≈ d(x, S g (0)) in the vicinity of the interface S g (0). Hence, the level-set function φ corresponds to the signed distance function to the interface in the neighbourhood of this interface. Furthermore, regarding [4], φ is extended by a constant value (positive or negative) equal to ±2E/π outside the narrow band [-E, +E] around the interface S g , and is differentiable through the boundaries of this band, i.e., at each x such that d(x, S g ) = ±E.

At any time t > 0, φ is solution, in a neighbourhood of the interface, of the "classical" advection equation,

∂φ ∂t + v s • ∇φ = 0, [5]
where v s is given by Equation [3], while outside this neighbourhood, a renormalization procedure is applied by enforcing

∇ φ = 1 -( πφ 2E ) 2
This approach is similar to the classical renormalization technique consisting in enforcing ∇ φ = 1, but applied to a function φ of the form [4] (indeed, d dx sin(x) = cos(x) = 1 -sin 2 (x)). Hence, the advection and the renormalization steps are performed by solving one single equation as described in [START_REF] Coupez | Réinitialisation convective et locale des fonctions Level Set pour le mouvement de surfaces et d'interfaces[END_REF].

Curvature and surface velocity level-set formulation

By definition, the value of φ(x, t) gives the position of a point x at a time t with respect to the grain free boundary. Furthermore, and it is a key point for the surface diffusion simulation, both the unit vector n φ normal to the interface and the mean curvature K φ of this interface can be expressed in a level-set way by the two following relations [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF][START_REF] Osher | Level Set Methods: An Overview and Some Recent Results[END_REF]:

n φ = ∇ φ ∇ φ [6]
and,

K φ = ∇ • ∇ φ ∇ φ [7]
where • denotes the Euclidian norm in IR 3 . Hence, if the level-set function is known, then the curvature can be calculated, at least theoretically. Note that expressions [6] and [7] define the normal vector and the mean curvature in all the computational domain Ω, and not only over the free surface {φ = 0}. Of course, n φ and K φ vanish outside the narrow band [-E, E] around the interface, while they correspond to the usual normal vector n and curvature K in the vicinity of the interface (Bruchon et al., 2009a). Following [START_REF] Burger | A level-set approach to anisotropic flows with curvature regularization[END_REF], the surface diffusion velocity [3] can be rewritten in a level-set form as

v s = v n ∇ φ ∇ φ [8]
and,

v n = C 0 1 ∇ φ ∇ •( ∇ φ P φ ∇ K φ ) [9]
where the projection matrix [2] is now defined by

P φ = I - ∇ φ ∇ φ ⊗ ∇ φ ∇ φ
The velocity [8] is defined in the whole computational domain, and corresponds, in the vicinity of the zero level set of φ, to the surface diffusion velocity. The differential operators (gradient and divergence) involved in Equation [9] are now expressed explicitly with the cartesian coordinates (x, y, z), and can therefore be computed within the context of an Eulerian description of the grains.

Curvature -normal velocity stabilized mixed system

Formulation [9] is discretized by using a finite element approach. The computational domain Ω is discretized by a simplex mesh T h (Ω), that is, a set of tetrahedrons. The unknowns φ, K and v n are approximated by φ h , K h and v nh , respectively, chosen as being continuous and piecewise linear over Ω. Furthermore, the partition

0 = t 0 < t 1 < • • • < t Θ = Θ of the time interval [0, Θ] is introduced: φ t h , K t h and v n t
h denote the finite element approximations evaluated at the time t. Regarding Section 3.2, all seems easy: if φ t h is assumed to be known, then K t h and v n t h could be directly computed by introducing Equations [6], [7] and [9] into the finite element formulation. However, three difficulties appear and must be understood. First, the velocity v nh depends on the fourth-order spatial derivative of φ h . Since φ h is piecewise linear, the gradient ∇ φ h is piecewise constant, and the nth-order derivative ∇ (n) φ h is identically equal to zero when n ≥ 2. Equations [7] and [9] must therefore be considered in a weak sense (see Equation [12] and [START_REF] Burger | A level-set approach to anisotropic flows with curvature regularization[END_REF][START_REF] Bänsch | A finite element method for surface diffusion: the parametric case[END_REF]). The second difficulty is the nonlinear coupling between the level-set function φ, the curvature and the surface diffusion velocity. We have chosen to treat the velocity in an explicit way in the transport equation [5]. In other words, φ t+∆t h is computed by solving Equation [5] with v s t h , the velocity evaluated on the configuration described by φ t h . Finally, the last difficulty is that the fully explicit scheme which consists in computing successively K t h and v n t h , assuming that φ t h is known, and then to transport φ h from time t to time t + ∆t, leads to numerical oscillations of the level-set function [START_REF] Bruchon | Modélisation et simulation du transport de matière par diffusion surfacique à l'aide d'une approche Level-Set[END_REF].

To overcome these difficulties, the numerical method proposed in this section consists in building a system, the unknowns of which are the curvature K t h and the normal velocity v n t h (which is equal to C 0 ∆ s K h regarding Equation [3]). In order to introduce implicitly a regularization term in this formulation, the following first-order Taylor's expansion is considered:

φ t+ 1 2 h = def. φ t h + ∂φ t h ∂t ∆t = φ t+∆t h + o(∆t) Hence, φ t+ 1 2 h
is a first-order approximation of φ t+∆t h . Since the level-set function φ is solution of the transport Equation [5], ∂φ/∂t = -v s • ∇ φ, and the previous relation can be turned into

φ t+ 1 2 h = φ t h -v s t h • ∇ φ t h ∆t
Finally, taking into account expression [8] leads to

φ t+ 1 2 h = φ t h -v n t h ∇ φ t h ∆t [10]
The system with the unknowns K/ ∆ s K can now be constructed by considering

φ t+ 1 2 h instead of φ t h in relation [7]
giving the curvature of the discretized level-set function. Hence, at a time t, assuming that φ t h is known, the system in curvature K t h / normal velocity v n t h is given by

K t h + ∇ • ∆t A ∇ v n t h = ∇ • 1 A ∇ φ t h v n t h ∇ φ t h -C 0 ∇ •( ∇ φ t h P φ t h ∇ K t h ) = 0 [11]
The key point of this approach is that using φ

t+ 1 2 h
instead of φ t h to compute the curvature, induces "naturally" an additional term in the left-hand side of the first equation of [11]. This term, coupling K and v n , is a regularization term of the form ∆ v nh , with the regularization parameter ∆t/A. The parameter denoted by A should be equal to ∇ φ t+ 1 2 h . However, to avoid to deal with this nonlinear term, A is simply taken equal to ∇ φ t h -∆t ∇ v n t-∆t h

. Furthermore, the term ∇φ t which appears in Equation [10] has been taken identically equal to 1 in this curvature expression.

Since the effective computation of K h and v nh is carried out by using a finite element discretization of [11], the mixed weak formulation of [11] has to be written: At time t, assuming φ t h known, find

(K t h , v n t h ) ∈ V h × V h solution of Ω K h ψ h dΩ -∆t Ω 1 A ∇ v n t h • ∇ ψ h dΩ = - Ω 1 A ∇ φ t h • ∇ ψ h dΩ Ω ∇ φ t h v n t h ψ h dΩ + Ω C 0 ∇ φ t h (P φ t h ∇ K t h ) • ∇ ψ h dΩ = 0 for all ψ h ∈ V h [12]
Note that one single type of weighting functions ψ h has been used for the curvature equation and for the velocity one, since both curvature and velocity belong to the same functional space V h , the space of the continuous functions which are piecewise linear. Furthermore, previous Equation [12] does not require the enforcement of Dirichlet conditions: since φ h is constant far from the free boundary, K h and v nh vanish over the mesh boundary.

To summarize the previous developments, after an initialization step consisting in calculating the initial level-set function by Equation [4], our time-stepping strategy consists, at time t and assuming φ t h is known, in the computation of the surface velocity by [12] and [8], and in the subsequent transport of the level-set function by solving Equation [5]. A mesh adaptation strategy is also used to obtain an accurate description of the grain free boundary (Bruchon et al., 2009a).

Numerical simulations

This section presents several direct simulations of change in free surface by surface diffusion, using the numerical strategy developed previously. The presented developments have been implemented in the CIMLIB finite element library. This C++ library, which is highly parallel, is developed at Centre de Mise en Forme des Matériaux (Mines ParisTech, UMR CNRS 7635) by the team of T. Coupez [START_REF] Digonnet | Object-oriented programming for fast-and-easy development of parallel applications in forming processes simulation[END_REF].

Change in two grains of the same size by surface diffusion

Figure 2 shows the temporal change in the zero isosurface of the level-set function corresponding to two spherical grains of radius R = 0.2. These grains are initially in contact (see Figure 2(a)), though the discretization with a mesh does not allow the grains to be perfectly tangent). Despite the initial "roughness" of the contact area between the grains, and due to the matter diffusion, this surface becomes quickly smooth. The grain volume is well-preserved, and as expected, no shrinkage phenomenon occurs. The radius x of the circular contact area is called the neck radius. Theoretical models, based on geometrical assumptions and established within the context of the sintering process modelling [START_REF] Rahaman | Ceramic processing and sintering[END_REF], predict that the adimensional neck radius x/R behaves as the power 1/7 of the time:

x(t) R = 56C 0 R 4 t 1/7 = 1.78t 1/7 [13]
where t is an adimensional time, defined by t = C 0 R 4 t. Figure 3 shows, in logarithmic scale, the temporal change in the adimensional radius x/R, obtained by direct simulation for different values of the grain radius, ranging from 0.1 to 2.5. Accordingly to the analytical model, the simulation shows that the change in x/R versus t , does not depend on the grain radius R, and behaves as t 1/7 (curve mentioned as "Simulation, 1/7" in Figure 3). The only difference between the analytical model and the simulation results, is the coefficient 1.78 provided by the analytical model (curve referred to as "Analytical model, 1/7"). Our simulations provide a coefficient approximatively equal to 1.3, whatever the time step (ranging from 10 -6 to 10 -2 ) and the mesh size (ranging from 10 -4 to 10 -2 , with isotropic or anisotropic remeshing). In Equation ( 13), the two key parameters which characterize the surface diffusion from the grain surface towards the neck, are, on one hand the power 1/7, and on the other hand, the power 4 applied to R in the denominator [START_REF] Rahaman | Ceramic processing and sintering[END_REF]. The differences between the simulations and the analytical model may be explained by the geometrical approximations made on the curvature in the analytical model.

Change in a small cluster of grains by surface diffusion

The sintering by surface diffusion of 20 grains (see Figure 4) is investigated in this section. The grain radii vary randomly uniformly between 0.08 and 0.2, while the grain position is initially randomly uniformly distributed, with the only restrictions that there is no isolated grain, and that two grains in contact must be tangent, as shown in Figure 4(a). The computational domain is the unit-cube, discretized with an unstructured mesh made up of 1,320,000 tetrahedrons (227,000 nodes). This mesh is adapted dynamically, in order to be refined in the vicinity of the free boundary, according to the strategy developed in (Bruchon et al., 2009a) (the mesh size ranges from 0.01 to 0.08). As previously, the contact area between two grains are not initially well-defined. However, it can be seen when comparing Figures 4(a) and 4(b) that the diffusion smoothes quickly these surfaces. Neck formation is then observed, and the cluster evolves towards an equilibrium state. As expected, no shrinkage phenomenon occurs (the dimensions of the bounding box containing the grains remain constant along the simulation), and the grain volume is well preserved during the simulation as shown in Figure 5. The CPU time of this computation involving 160 time steps, is of 165 minutes by using a parallel computing strategy on four cores (Intel Xeon 2.2 GHz processors).

Conclusion

Within the general context of the sintering process simulation, a level-set formulation of the surface diffusion problem has been proposed. The main difficulty induced by this approach, the dependence of the surface velocity on the fourth-order spatial derivative of the level-set function, has been overcome by considering a mixed system in curvature / surface Laplacian of the curvature. This formulation induces implicitly a stabilisation term ∆v n with an associated stabilisation parameter depending on the time step. The resulting discretized formulation is shown to be stable in all our simulations. The accuracy of this approach has been proved by investigating the surface diffusion between two grains and by comparing, in this case, the simulation with an analytical model. Furthermore, the ability of our numerical strategy, to describe the changes occuring in complex geometries, has been shown by detailing a simulation involving 20 grains. The next step in the development of this work is to describe the join separating two particles, in order to model the matter transport by grain boundary diffusion.

Figure 1 .

 1 Figure 1. Matter transport between two identical grains

Figure 2 .Figure 3 .

 23 Figure 2. Change in the free boundary {φ h = 0} due to the surface diffusion between two grains of equal size (radius R = 0.1)

Figure 4 .Figure 5 .

 45 Figure 4. Surface diffusion simulation: change in the free boundary {φ h = 0} of 20 grains of radii randomly uniformly distributed between 0.08 and 0.2 (C 0 = 10 -7 , ∆t = 10 -2 )

EJCM -19/2010. Giens 2009

EJCM -19/2010. Giens 2009

ANNEXE POUR LE SERVICE FABRICATION