F Gatuingt

L Snozzi

J F Molinari

Numerical determination of the tensile response and the dissipated fracture energy of concrete: role of the meso-structure and influence of the loading rate

At the mesoscopic scale concrete can be considered as a mix of coarse aggregates with a mortar paste matrix. In this paper we investigate numerically the influence of aggregates arrangements and loading rate on the tensile response of concrete. Each coarse aggregate is assumed to be circular with six different radiuses following the aggregates size distribution of real gravel. Rateindependent cohesive elements are used to model failure within the mesostructure. Our results show that the spatial distribution of heterogeneities does not influence the peak strength, while it changes the post-peak macroscopic response. This implies that our specimen size is large enough for strength computation but that larger mesostructures should be considered to obtain fully reliable toughness predictions. While, the cohesive approach is able to capture the transition from one macro-crack in quasi-static to multiple micro-cracks in fast dynamics, which increases the dissipated fracture energy, our results suggest that the full extent of the high-rate strengthening of concrete observed experimentally for loading rates greater than ε = 1/s cannot be captured with rate independent constitutive laws.

Introduction

Reinforced concrete is one of the main materials used in constructions. Industrial or civil structures could be subjected to dynamic loads and these loads have to be taken into account in the design process. For example, structures that involve public safety have to be designed to resist not only to the static loading but also to the dynamic loading produced by extreme conditions such as earthquakes, explosions or accidental collisions. To achieve the design and analyses of these structures, it is important to investigate the dynamic mechanical response of concrete. Several macroscopic models for concrete have been developed [START_REF] Ottosen | Constitutive model for short-time loading of concrete[END_REF][START_REF] Sluys | Wave propagation, localisation and dispersion in softening solids[END_REF][START_REF] Dubé | Rate dependent damage model for concrete in dynamics[END_REF][START_REF] Sercombe | Viscous hardening plasticity for concrete in high-rate dynamics[END_REF][START_REF] Bažant | Fracturing rate effect and creep in microplane model for dynamics[END_REF][START_REF] Gatuingt | Coupled damage and plasticity modelling in transient dynamic analysis of concrete[END_REF][START_REF] Hervé | On numerical implementation of a coupled rate dependent damage-plasticity constitutive model for concrete in application to high-rate dynamics[END_REF][START_REF] Pedersen | An analysis of dynamic fracture in concrete with a continuum visco-elastic visco-plastic damage model[END_REF][START_REF] Gatuingt | Anisotropic 3d delay-damage model to stimulate concrete structures[END_REF][START_REF] Pontiroli | Predicting concrete behaviour from quasi-static loading to hypervelocity impact An overview of the PRM model[END_REF][START_REF] Desmorat | Nonstandard thermodynamics framework for robust computations with induced anisotropic damage[END_REF] to perform non-linear numerical analysis of such problems. In these models a rate effect in tension has often been introduced to represent the experimental data for loading strain rates exceeding 1 /s [START_REF] Weerheijm | The rate dependency of concrete in tension -New data for wet, normal and dry conditions. DYMAT 2009: 9th Inter[END_REF][START_REF] Erzar | An Experimental Method to Determine the Tensile Strength of Concrete at High Rates of Strain[END_REF]. Indeed, it is well known that for loading rates ranging from static (10 -4 /s) to intermediate (10 -1 -0 /s), a moderate rate effect in tensile strength is observed. Beyond a loading strain rate of 1 /s a very steep strength increase occurs.

Figure 1 shows a set of static, Split Hopkinson Bar (SHB) and modified SHB tests. The derived data on strength and fracture energy of theses tests are summarized in Table 1. These experimental Figure 1: Experimental stress-displacement curves for static, SHB and Modified SHB tests (from [START_REF] Weerheijm | Research developments and experimental data on dynamic concrete behaviour[END_REF]).

data show a rate effect on strength, and on the dissipated fracture energy, in the high loading rate regime. Table 1 also shows that for the higher loading velocity the dynamic increase factor is more important for the dissipated fracture energy than for the tensile strength. In addition, some authors have shown experimentally that the deformation capability of concrete under uniaxial tensile loading increases due to a rate effect [START_REF] Kormeling | Strain rate effects on steel fibre concrete in uniaxial tension[END_REF][START_REF] Oh | Fracture behavior of concrete under high rates of loading[END_REF][START_REF] Rossi | Effect of loading rate on the strength of concrete subjected to uniaxial tension[END_REF]. As previously mentioned, it is usual to introduce in numerical models this dynamic increase factor directly in the macroscopic behaviour as shown in Figure 2. This has to be done with care as it is not clear how important the dissipation at the structural level is. Indeed, the increase of the fracture energy with the strain rate may be caused, in the experiments, by a structural effect (at the concrete sample level) instead of a material one's [START_REF] Van Doormaal | Experimental and numerical determination of the dynamic fracture energy of concrete[END_REF][START_REF] Vegt | Failure mechanisms of concrete under impact loading[END_REF].

Stress (MPa)

Figure 2: Uniaxial tension response for different loading rates of a classical visco-damage model (from [START_REF] Chambart | Endommagement anisotrope et comportement dynamique des structures en béton armé jusqu'à la ruine[END_REF]).

Multiple cracking and micro-inertia effects in the fracture process zone are thought as the main factors for the rate dependency of concrete in the high dynamic regime. However, the influence of the internal microstructure of concrete and the mechanisms that lead to different crack patterns when varying the loading rate, remain open questions. In order to give some answers to these questions, concrete has to be considered as a heterogeneous material where the nature of the heterogeneity depends essentially on the scale of observation. If one takes the scale of sand's grain, concrete can be considered as a biphasic material made of aggregates of different sizes randomly distributed in a mortar matrix. Thus, the failure of a concrete sample in tension is related to processes that take place at the so-called mesoscale. In order to investigate the role of the meso-structure on the dynamic tensile response of concrete, one can use a numerical approach. This is achieved in this paper using a 2D finite element mesoscopic description of concrete (aggregates + matrix) with cohesive capability. This method has proven its efficiency on numerical simulations of fracture of brittle materials [START_REF] Levy | Dynamic fragmentation of a ring: Predictable fragment mass distributions[END_REF] and has been extended to concrete like materials [START_REF] Snozzi | Influence of the meso-structure in dynamic fracture simulation of concrete under tensile loading[END_REF]. This paper extends the work conducted by Snozzi et al. [START_REF] Snozzi | Influence of the meso-structure in dynamic fracture simulation of concrete under tensile loading[END_REF]. In the first section, we propose a more realistic description of the aggregates heterogeneities in a concrete. We generate different meshes based on five meso-structure descriptions of a concrete mix. In section 3, we describe the numerical tensile tests and investigate the mesh objectivity of the observed response. The effect of the stiffness of the inclusions (an additional complement to [START_REF] Snozzi | Influence of the meso-structure in dynamic fracture simulation of concrete under tensile loading[END_REF]) on the macroscopic tensile response is also discussed. In sections 4 and 5, the main aspects of this work are described. These are the study of the influence of the meso-structure on the tensile strain-stress curve and on the dynamic increase factors for the tensile strength and dissipated fracture energy.

Meso-scale model

The concrete constitutive behavior can be formulated either at the macro-scale or at the meso-scale. At the macro scale, the ingredients that characterize concrete's heterogeneity are not represented and one considers it as a homogeneous material. Nevertheless, it is a material whose heterogeneities (aggregates) can have sizes of the order of tens of millimeters, which influences fracture properties (the crack path is for instance clearly dependent on these heterogeneities). The approach used in this paper therefore considers concrete as a bi-phasic material: aggregates in a mortar paste, which is commonly called a meso-scale approach (as proposed in [START_REF] Wittmann | Structure of concrete with respect to crack formation[END_REF]). It allows an explicit representation of some concrete components, which enables reducing the number of model parameters, and a description of the interactions between matrix and inclusions. In the literature one can find several meso-scale models for concrete. They can be divided in two main computational classes. A first family is represented by lattice models (for instance [START_REF] Schlangen | Simple lattice model for numerical simulation of fracture of concrete materials and structures[END_REF][START_REF] Bolander | Fracture analysis using spring networks with random geometry[END_REF]), where the continuum is replaced by a system of discrete particles and the mechanical properties of the lattice beams aim to represent the concrete meso-structure [START_REF] Kozicki | Effect of aggregate structure on fracture process in concrete using 2d lattice model[END_REF][START_REF] Grassl | Meso-scale approach to modelling the fracture process zone of concrete subjected to uniaxial tension[END_REF][START_REF] Cusatis | Strain-rate effects on concrete behavior[END_REF]. The second class resorts to the finite-element approach, in which concrete is usually represented as a biphasic material, made of a mixture of aggregates embedded in a matrix phase with an interfacial transition zone (ITZ) between them [START_REF] Roelfstra | Le béton numerique[END_REF][START_REF] Stankowski | Numerical simulation of progressive failure in particle composites[END_REF][START_REF] Lopez | Microstructural analysis of concrete fracture using interface elements. application to various concretes[END_REF][START_REF] Caballero | 3d meso-mechanical numerical analysis of concrete using interface elements[END_REF][START_REF] Carol | Micromechanical analysis of quasi-brittle materials using fracturebased interface elements[END_REF][START_REF] Snozzi | Influence of the meso-structure in dynamic fracture simulation of concrete under tensile loading[END_REF].

From those works it results that the meso-scale level of material observation is suitable to capture the main characteristics of the overall mechanical behavior of concrete: micro-cracking, fracture initiation and propagation, coalescence and localization. Moreover the constitutive equations might be relaxed when comparing with a macro-scale formulation and thus the number of parameters may be reduced.

In our model, the aggregates, the mortar and the interfacial transition zone (ITZ) will be represented as separate phases with different material properties. In order to study the specific effect of the meso-structure on the dynamic tensile response, no random fields have been introduced on the properties. A finite element approach with cohesive capability is then used to simulate the cracks opening and propagation.

Meso-structure modelisation

The macroscopic response of a heterogeneous material at the mesoscopic scale depends not only on the properties of the different phases but also on their geometric distribution. The aim of this paper is to see the influence of the heterogeneities on the mechanical response in tension. We have chosen constant mechanical properties and focus only on the spatial distribution of the heterogeneities.

Real concrete is based on randomly distributed aggregates in a cement paste matrix. Figure 3a shows particle size distributions for real sand and gravel used in a real concrete. These were obtained by a sieving method [START_REF]Tests for geometrical properties of aggregates -Part 1 : determination of particle size distribution -Sieving method[END_REF]. The aggregates called sand can be included either in the heterogeneities or in the cement paste depending on the fineness of the aggregates description in the meso-scale model. In our description, we chose to describe explicitly only the aggregates with a diameter larger than 4 mm and smaller than 25 mm with six classes to capture roughly the aggregates size distribution (see Figure 3b). The smallest aggregates are then taken into account in the homogeneous matrix of mortar. The size of the generated specimens is 100×100 mm 2 . As in [START_REF] N'guyen | A mesoscopic approach for a better understanding of the transition from diffuse damage to localized damage[END_REF], knowing the aggregate size distribution, the number of perfect circular (in 2D) aggregates of each class can be calculated. To avoid an artificial increase of the volume matrix due to the boundary effect, the heterogeneities are placed in a larger sample with a final cut of all parts outside of the concrete specimen. Figure 4 shows five images of the meso-structures obtained using this method with the same aggregate size distribution of Figure 3b.

A 2D finite element mesh is then generated from the meso-structure images using the PPM2OOF [START_REF] Langer | The ppm2oof manual[END_REF] public domain software created at the National Institute of Standards and Technology (NIST). This software superimposes a mesh on any image. The simplest meshing algorithm tries to divide each image pixel into two triangular elements. In our approach, we resort to an adaptive meshing algorithm available in PPM2OOF. In this case, the mesh that is created represents the user's compromise between the number of elements and the homogeneity of each element with respect to its underlying pixel properties and the production of equilateral triangular elements.

Figure 5 shows a zoom on four different meshes obtained with the same image of the first mesostructure (Meso 1) of Figure 4. In mesh A, the average size of the elements is 1 mm (20000 elements), while in the mesh B the size is 0.5 mm (80000 elements) and in mesh C we have 180000 elements of approximately 0.3 mm. A fourth, very fine mesh (thereafter called D), has also been generated with 320000 elements (0.25 mm) in order to study the convergence of the numerical results. One can observe that the definition of the contour of the aggregates and the real surface of paste is better when the size of the mesh decreases. To illustrate this point, in the proposed zoom, one can see that the mortar in between two grains can disappear for the coarse meshes while it is well represented with the finer one.

Notice that this mesh generation is possible since the cohesive element methodology here does not require an a priori definition of the possible fracture planes [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF][START_REF] Roelfstra | Le béton numerique[END_REF][START_REF] Stankowski | Numerical simulation of progressive failure in particle composites[END_REF][START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF][START_REF] Carol | Micromechanical analysis of quasi-brittle materials using fracturebased interface elements[END_REF][START_REF] Caballero | 3d meso-mechanical numerical analysis of concrete using interface elements[END_REF].

The cohesive element method

The cohesive element method allows us to model dynamic crack propagation and damage in a brittle material like concrete. The fracture process is described by the cohesive approach (introduced by Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblat [START_REF] Barenblatt | The mathematical theory of equilibrium of cracks in brittle fracture[END_REF] in the 1960s) as a separation process occurring at the crack tip in a small region of material called cohesive zone. This can be introduced into a standard finite element environment using interface elements with null thickness and with a fracture-based constitutive law. We assume that the bulk material outside the cohesive zone remains elastic. In our case the crack path is not known a priori and all lines in the mesh are considered as a potential crack path. During the simulation, the stress at the interface between two adjacent continuum elements is computed and compared to the fracture criterion at the end of every time step. The interfacial stress, σ, is calculated averaging stresses of the adjacent Gauss points of the two continuum elements. If the inter-element stress exceeds the critical stress value, the nodes located at the inter-element boundary are doubled, the two elements are topologically disconnected and a cohesive element is inserted (see Figure 6). After the nodal disconnection, the interfacial stress starts being controlled by the traction separation law implemented in the cohesive element (Figure 7).

The coupling of the above mentioned fracture criterion together with an initially rigid cohesive law allows capturing the initiation of new cracks without an a priori definition of the possible cracks. Nevertheless, the cracks are constrained to propagate following the inter-element boundaries and the fineness of the mesh can affect it. The law we use is the linear irreversible softening law proposed by Camacho and Ortiz [START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF]. The authors hypothesize the existence of a free potential energy Φ, which depends on one effective scalar displacement which has the following form:

δ = ∆ 2 n + β 2 ∆ 2 s (1)
where -→ ∆ = (∆ n , ∆ s ) is the relative displacement vector and the parameter β accounts for the coupling between normal and tangential displacements. The value of β has to be estimated (e.g. by correlating experimental results with numerical simulations [START_REF] Ruiz | Three-dimensional cohesive modeling of dynamic mixed-mode fracture[END_REF]) but it has been shown by [START_REF] Snozzi | Influence of the meso-structure in dynamic fracture simulation of concrete under tensile loading[END_REF] that the peak strength is only slightly affected by β for a tensile loading. One will keep its value equal to one in this study.

The differentiation of the free potential energy with respect to the opening displacement leads to the cohesive tractions law:

-→ t = ∂Φ ∂ -→ ∆ = t δ (∆ n -→ n + β 2 ∆ s -→ s ) (2)
where t represents an effective cohesive traction. This traction in case of crack opening is given by:

t = f ct 1 - δ δ c for δ = δ max and δ > 0 (3)
Where f ct represents the local material strength and δ c represents the effective relative displacement beyond which complete decohesion occurs. Whereas for crack closure or reopening (δ smaller than δ max ) the functional form is assumed to have the form:

t = δ δ max t max for δ < δ max ( 4
)
where t max is the value of the effective traction when δ is equal to δ max , in which is stored the maximal effective opening displacement attained up to the moment. Moreover δ max also accounts for the irreversibility of the law allowing successive loading, unloading and reloading. The evolution of the linear decreasing law is graphically shown in Figure 7. Note that the definition of f ct and δ c implicitly establishes the existence of an effective fracture energy G c , which corresponds to the area under the curve of Figure 7:

G c = 1 2 f ct δ c (5)
Partially damaged cohesive elements have dissipated an energy W < G c . As the present work is focused on tensile uniaxial loading, the fracture energy of the different phases of our concrete is identified to the experimental one for Mode I, i.e. G c = G I c .

3 Numerical simulation of direct tensile test in dynamics

Initial and boundary conditions

The aim of the present study is to analyze the dynamic tensile response of a concrete specimen (2D plane strain). The specimen is loaded under displacement control with an imposed strain rate ε, all the nodes of the finite element mesh which are located on the upper (respectively lower) boundary are forced to move at a constant velocity V 0y = V 0 (respectively V 0y = -V 0 ) with V 0x = 0:

δ δ c δ max G c f ct t max
V 0 = εh 2 ( 6
)
where h is the height of the studied specimen. To avoid stress wave propagation and an early fracture near the boundaries, all the nodes of the finite element mesh are prescribed an initial velocity in accordance to their vertical position y [START_REF] Miller | Modeling and simulation of dynamic fragmentation in brittle materials[END_REF] as illustrated in Figure 8:

V y (y) = 2V 0 h y (7) h v y v 0 -v 0 y v 0 v 0 x y Figure 8: Velocity initial condition
No initial transversal velocity is applied on the specimen (V x = 0). This means that due to Poisson effect, the stress state is not fully uniaxial (note that this is also the case in real experiments). Indeed, it is well known that the rate effect in compression is mostly due to a structural confinement and the main aim of this tension test is to try to understand if we have the same origin for the experimental tensile rate effect.

The finite element code used for this study is based on a Newmark explicit time integration scheme. This scheme is conditionally stable so that the time step has to be smaller than a critical value usually computed as follows (Courant-Friedrichs-Lewy condition):

∆t crit = Ne min e=1 l e c ( 8
)
where c is the pressure wave speed in the phase of the material of the element e and l e is the size of this element (in 2D the radius of the larger inner circle). The stable time step has to be chosen equal to the smallest value of the critical time step -over all the elements -multiplied by a security coefficient α (equal to 0.1 in our simulation). For example, for the meshes shown on Figure 5, one obtains for the computational time step ∆t = 1.838 10 -9 s and 7.882 10 -10 s for the meshes A and D respectively.

To obtain the stress-strain curves, we define the macroscopic stress as the boundary reaction force F y divided by the initial width of the specimen and the macroscopic strain as the displacement U y of the boundary divided by the initial height h of the specimen.

Material parameters and mesh objectivity

As shown in section 2 our numerical model is based on three phases: the mortar paste, the aggregates and an ITZ. For these phases we have to identify the elastic and rupture parameters of equation 5.

In the present work, we used the parameters given in Table 2. The relationship between the cohesive element law and G c introduces a length scale -called the cohesive zone length l z -into the material description. Rice and Palmer [START_REF] Palmer | The growth of slip surfaces in the progressive failure of over-consolidated clay[END_REF] calculated this length for linear elastic fracture mechanics in Mode I plane strain crack loading:

l z = 9 32 πEG c (1 -ν 2 )f 2 ct (9)
This length has an important influence around the crack tip in numerical simulations. Typically one assumes that for mesh independency, the cohesive zone length has to be smaller than the specimen size and should contain at least three to four finite elements. For the material parameters of Table 2, the smallest and more constraining (or penalizing) cohesive zone length is of the order of 40 mm for the aggregates (see Table 3). It is important to notice that the definition of l z is not clear for the ITZ. Indeed, in our representation the ITZ is not meshed (null thickness) with the consequence that the Young modulus is not known in this region. Nevertheless, based on Equation 9 and the Mode I assumption, all the meshes of Figure 5 should give the same results due to an average element size lower than 1 mm. In a previous work (Gatuingt et al. [START_REF] Gatuingt | Influence of the meso-structure in numerical simulation of dynamic tensile fracture of concrete[END_REF]) with other parameters for the aggregates, the cohesive zone length was equal to 4mm and we showed that the mesh A was too coarse in this case. 9 where the macroscopic tensile responses obtained for the different mesh refinements and two loading rates are shown. One can see on this figure that the differences between the four meshes are very slight. Only the mesh A gives something slightly different than the others for ε = 10/s. For this mesh, the element size (≈ 1 mm) is a little bit too coarse to have a good representation of the cement paste between two close aggregates (see Figure 5). From the second mesh size (mesh B) the results are converged for all the loading rates applied on the numerical specimen. This is confirmed with the Figure 10 where the crack path is plotted for all the meshes. One can see that in this case the crack path is not strongly mesh dependent. Indeed, it is well known that for a classical concrete (our study) the crack path in tension is mainly governed by the aggregate debonding caused by the low value of f ct for the ITZ. A numerical and an experimental study of the impact of the ITZ strength on the global strength are given in [START_REF] Pedersen | Mesoscopic modelling of concrete under impact[END_REF][START_REF] Rocco | Effect of aggregate shape on the mechanical properties of a simple concrete[END_REF] respectively. In the rest of the paper we choose for all the meso-structures studied a mesh refinement similar to the mesh B. This mesh size is small enough for the definition of l z and to describe the smallest aggregates and large enough to induce an acceptable CPU time for parametric studies.

Influence of the specimen size

In order to see the influence of the specimen size on our results, we generate a large sample in which we cut a smaller one in order to have the same aggregate/paste ratio and the same aggregate 11). The large sample has a surface S L = 0.15 × 0.15 cm 2 while for the small one S S = 0.1 × 0.1 cm 2 . The consequence is that we increase the width of the sample of 50% and the surface of 225%. Figure 12 and Figure 13 show the results obtained for two different loading rates ( ε=1/s and ε=100/s). We can notice that the crack path obtained is very similar for the two sample sizes. This result is confirmed with the macroscopic stress-strain responses of the specimens. Nevertheless we can observe that the peak strength is a little bit smaller for the large sample compared to the small one. This size effect in quasi-brittle materials such as concrete is a well known phenomenon and there are a number of experimental and theoretical studies (see for example [START_REF] Reinhardt | Similitude of brittle fracture of structural concrete[END_REF][START_REF] Bažant | Size effect in blunt fracture: Concrete, rock, metal[END_REF][START_REF] Ožbolt | The size effect in concrete structures[END_REF]) which confirm its existence. Moreover in Figure 13, the post-peak response is clearly different for the two sample sizes. This result is not surprising because the post-peak response is more governed by the fracture energy. This is confirmed in Figure 14a where the evolution of the dissipated fracture energy (per unity in the sample depth) is plotted for different strain rates. As expected, for a larger sample we obtain a higher dissipated energy. If we assume an horizontal single crack, we can divide this energy by the width of the specimen and plot a normalized dissipated fracture energy (Figure 14b). In this figure we can see that the assumption of a single crack is not inaccurate for the lower strain rate ( ε = 1/s) where a "single" crack appears (Figure 12). For higher strain rates ( ε = 100/s), due to the multiple cracks this assumption is no longer true and the normalized dissipated fracture energy strongly depends on the specimen size. The larger the sample, the bigger is the Dynamic Increase Factor (DIF). In our case, we obtain a DIF equal to 8 for the smaller sample and equal to 10 for the larger one. We can conclude that in high loading rate regimes the numerical (and experimental) results are affected by this structural effect and it is not easy to properly quantify the true fracture energy values (the size of the sample is important).

To choose the sample we have to used, we can notice that the size of the numerical "small" specimen is bigger than the experimental one (cylinder with a diameter of 74 mm [START_REF] Weerheijm | Research developments and experimental data on dynamic concrete behaviour[END_REF]) in which case a smaller DIF value should be obtained. We recall that the numerical DIF is 8 for the small sample whereas the experimental DIF is 12. As the aim of this work is not to study how we can reproduce the size effect with our model (even if the global trend seems adequate) we decided to keep the smaller size of the specimen for the subsequent sections.

Role of the meso-structure on the dynamic tensile response

The meso-structures used in this part are those presented in Figure 4 with a mesh refinement similar to the mesh B in order to have a quasi-independence of the chosen mesh on the numerical results.

As previously mentioned, the numerical model does not contain any random field to represent the scattering often observed on the macroscopic response of concrete. This scattering can be important in experimental tests for the same initial concrete mix but tested on several specimens. The main differences are then due to the different aggregate arrangements and to the macro-porosity distribution. In this section, we study only the influence of spatial distribution of aggregates and not the role of the macro-porosity. This porosity has probably a great influence on the crack initiation and ultimately on the crack path but in order to clearly separate the influence of aggregates, it will be omitted in our simulations.

The results of the strain-stress curves obtained for our five numerical concrete specimens in tension for several strain rates are presented in Figure 15. One can see on this figure that the microstructure has a slight influence (less than 4 %) on the peak strength value of our concrete. The macroscopic tensile strength is mainly governed in our case by the ITZ strength and not by the meso-structure. This can be best explained by examining the failure process in quasi-static for which we can observe a gradual transition from diffuse micro-cracking to strain localization and finally to a macroscopic crack. Figure 16 shows this evolution to a unique macroscopic crack for ε = 0.1 /s. In the upper row of Figure 16, which was obtained at the peak stress just before strain localization, one can remark that for the five meso-structures, the micro-cracks are diffuse and that their numbers are almost independent of the aggregates spatial distribution. This explains the slight influence on the peak strength value observed at the macroscopic scale. Conversely, on the lower row of Figure 16, which is extracted from the softening part, we notice that the crack path is strongly dependent on the mesostructure. This is why one obtains a stronger difference (50 %) on the stress value for ε = 3.10 -4 at the end of the loading for ε = 0.1 /s.

Even if the failure process is not really the same for high strain rates (see next section), one can observe on Figure 15 that the meso-structure affects only the post-peak macroscopic response of our numerical concrete. It is interesting to notice that the influence of the meso-structure does not play the same role for all loading rates. For example the meso-structure 2 has the smallest strain at complete failure of the specimen for ε = 10 /s but this is not the case for the other loading rates. Therefore, while heterogeneities play in general an important role, no conclusion can be made on how their influence can be extrapolated from one loading rate to another. Experimental results on dynamic tension tests show a high rate sensitivity of tensile strength [START_REF] Vegt | Modelling of impact behaviour of concrete -an experimental approach[END_REF][START_REF] Erzar | An Experimental Method to Determine the Tensile Strength of Concrete at High Rates of Strain[END_REF]. Figure 17 presents results obtained in literature. If one considers these experimental results as a material behavior at the macro scale, this leads to models which are difficult to relate to physical mechanisms that occur during failure.

Nevertheless one can see from Figure 17 that two distinct modes appear with regards to the increase of the tensile strength with respect to the strain rate. In the first modeε < 1 /s -the dynamic resistance increase is probably due to the presence of water in the material. The free water in the micro-pores is assumed to exhibit the so-called Stefan-effect causing a strengthening effect in concrete with increasing loading rate [START_REF] Rossi | The dynamic behaviour of concrete: influence of free water[END_REF]. This Stefan-effect is the phenomenon that occurs when a viscous liquid is trapped between two plates that are separated quickly, causing a reaction force on the plates that is proportional to the velocity of separation. Cadoni et al. [START_REF] Cadoni | Strain rate effect on the tensile behaviour of concrete at different relative humidity levels[END_REF] give a different explanation for the influence of the moisture content. Their interpretation is based on the principle of wave propagation in concrete. When a pore is not filled with water, it will locally reflect the incoming stress wave. The multiple reflections of all pores together can cause a considerable increase Meso in stress. When a stress wave meets a pore that is filled with liquid, the reflected stress is not big enough to locally provoke the damage of the material. Therefore the wet concrete will exhibit a more pronounced rate effect than the dry concrete. This interpretation only gives an explanation of the difference between wet and dry concrete and does not explain the increase in strength between static and dynamic loading. At a larger scale (the constitutive level), Ozbolt et al. [START_REF] Ožbolt | Microplane model for concrete with relaxed kinematic constraint[END_REF] believe that the rate dependency consists of two parts: a rate dependency related to the formation of the micro-cracks, and a rate dependency due to the creep of concrete between the micro-cracks. Then, the influence of inertia forces on the rate effect is not part of the constitutive law, however, this effect is automatically accounted for in dynamic analysis in which the constitutive law interacts with inertia forces. Figure 17: Experimental rate effect in tension (from [START_REF] Erzar | An Experimental Method to Determine the Tensile Strength of Concrete at High Rates of Strain[END_REF])

In this work we will attempt to show if micro-inertia effects in the fracture process zone are sufficient to explain the rate dependency of concrete and the high fracture energy dissipation in a high rate regime. That is why the model used for the cohesive traction law of Figure 7 is not rate dependent. Moreover the concrete moisture content is usually assumed not to be dominant in the high strain rate regime. Thus, we will not consider it in this study. Figure 18 shows the numerical rate effect in tension obtained for three meso-structures with our assumptions. One can see on this figure that the computed increase with strain rate is slight concerning the tensile strength. We obtain here a dynamic increase factor -equal to the ratio of the static versus the dynamic strengths -of 1.3 for ε = 100 /s while it is equal to approximately 3 in the experiments. This result shows that the assumption of no rate effect in the cohesive traction law is not realistic. The micro-inertia effects in the fracture process zone are not sufficient to explain the rate dependency of concrete in tension even at the highest loading rate velocities. One has to take into account a rate dependency at the material level, for example by linking the cohesive strength f ct to the strain rate.

Figure 19 depicts the evolution of the dissipated fracture energy as function of the macroscopic strain of the specimen for different loading rates and heterogeneities. One can see on this figure that, as already observed in Figure 18, the dissipated fracture energy strongly depends on the loading rate even with a rate independent local fracture energy (see Equation 5 and Figure 7). Here we obtain a ratio G f (dyn) /G f (stat) equal to approximately 8 -depending on the meso-structure -for ε=100 /s. In the experimental data (see Table 1), with a Young's modulus of the concrete approximately equal to 40 GPa we obtain a ratio of 12 for ε=40 /s (compared to quasi-statics). We note that even if the dynamic increase factor for the tensile strength is not large enough, it is less the case for the dissipated fracture energy. Moreover, and as expected, the dissipated fracture energy depends on the meso-structure (Figure 19a vs Figure 19b) due to the differences observed in Figure 15 on the post-peak curves. According to the structural effect hypothesis, the explanation of higher dissipated fracture energy for high loading rates resides in a more diffuse micro-cracking. In the higher loading rate regime, the initiation of the micro-cracks are similar to that observed for a slow regime (see Figure 16 upper row) around the peak stress of the macroscopic strain-stress curve, but multiple cracks appear in the softening regime. This effect can be seen on Figure 20 where the final aspect of the concrete specimen based on the meso-structure Meso 1 are compared for low and high loading rates. This behavior only explains that the differences between low and high loading rates are not important for the peak stress but only for the total dissipated fracture energy, which is in contradiction with experimental observations. Moreover in Figure 20b/ we can notice that, contrary to experimental observations, we do not have aggregates failure at high rates. The origin of this phenomenon is not clear (several assumptions could be removed to better understand this, including, spherical shape of the aggregates, rate independency of the ITZ, choice of the material parameters).

On the contrary, the material view point sustains that the increase of toughness could be an effect of internal material "hardening". One can argue that this mechanism should be related to the chosen modeling scale. i.e. when considering a large modeling scale (macroscale), the toughness parameter of the cohesive approach could be a function of the loading rate, as the mesh size cannot be small enough to capture all micro-cracking mechanisms. Whereas, decreasing the scale should entail a rate independent toughness parameter. This part is not really clear on Figure 21, which shows the evolution of the total fracture energy dissipated for mesh B and mesh D of the meso-structure called Meso 1. Indeed, for the mesh D the average size of the elements is 0.25 mm and one observes a small increase in the dissipated fracture energy in comparison with the mesh B (0.5 mm) with no real effect on the peak strength. Moreover with the mesh D, the finite element size is close to the average size of the entrained air bubbles (micro-pores) [START_REF] Mehta | Concrete structure, properties, and materials[END_REF], the larger heterogeneities at this scale. So it could be pertinent to keep this scale and to incorporate a material rate effect in the cohesive traction law to represent the Stefan effect [START_REF] Rossi | The dynamic behaviour of concrete: influence of free water[END_REF] at the paste level. Snozzi et a. [START_REF] Snozzi | Influence of the meso-structure in dynamic fracture simulation of concrete under tensile loading[END_REF] have shown that choosing a rate dependent cohesive traction law -where the cohesive strength is a material constant value, while the critical opening displacement is a function of the crack opening rate -does not involve a difference in the peak strength values and dissipated fracture energy in comparison with a rate-independent law. They attribute this to a size effect due to the physical dimension of the RVE. Indeed, for the concrete phases parameters and for the loading rates considered in this paper, we observe that the numerical results are more controlled by the strength parameters (strength controlled regime) than by the toughness.

To improve our simulations and the representativeness of our numerical concrete, the challenge will be to find a rate dependent cohesive traction law (at least for the f ct value) that will not increase too much the total dissipated fracture energy and the material ductility . One can for example use a function closed to the one proposed by Bažant and Li [START_REF] Bažant | Cohesive crack with rate-dependent opening and viscoelasticity: I. mathematical model and scaling[END_REF]. But an open question remains. Should one use a rate dependent function for only the ITZ, mortar, or aggregates or for all the parts of the meso-structure? This point will be difficult to answer without some specific experiments on the different constituents of a real concrete.

Conclusion

In this paper we proposed a meso-scopic model for the analysis of dynamic tensile failure of concrete. This model is based on a 2D finite element description with cohesive capability of a mix of aggregates larger than 4 mm in a mortar paste matrix. The distribution of the coarse aggregates is assumed to follow a realistic aggregates size distribution. Five meso-structures are generated with this distribution. The influence of the heterogeneous meso-structure of concrete and the loading rate on the tensile response and the dissipated fracture energy are studied. With our specimen size we observe a small impact of the aggregates arrangement on the tensile strength. This is because for our model parameters, the size of our specimen is large enough for a strength analysis and thus we can consider that a "material" response (independent of aggregates arrangement) is properly captured. This is of course based on the strong assumption of a perfect mix with no macro-porosity. Nevertheless when the post-peak response and the dissipated fracture energy are considered, we show that the mesostructure has an influence (as well as the sample size) and then the "homogeneous" regime has not been obtained. It is important to notice that the definition of homogeneity depends not only on the material description but also on the phenomenon observed.

The second main interest of the paper is the study of the dynamic loading rate effects on the tensile strength and on the dissipated fracture energy with our approach. The cohesive law used is independent of the local or global strain rate to investigate if the experimentally observed dynamic increase factor is due to a material effect or a structural (inertial) effect. With regards to tensile strength, the numerical results exhibit a small rate effect that is not coherent with the experiments. Interestingly, for the dissipated energy the numerical results are close to the experimental ones. This implies that the tensile strength increase is probably mostly due to a material effect (moisture, viscoelasticity, . . . ) while the dissipated fracture energy observed in the post-peak part is more due to the increase of the number of micro-cracks with the loading rate. With adding a small rate dependency in the cohesive law could permit to have a realistic virtual material to test the influence of other experimental configurations.

Figure 4 :

 4 Figure 4: Different images of the meso-structures generated with the same aggregates size distribution.

Figure 5 :

 5 Figure 5: Different mesh sizes generated with PPM2OOF on the same meso-structure.

Figure 6 :

 6 Figure 6: Schematic of the nodal disconnection between two continuum elements in which the interfacial stress has exceeded the critical stress

Figure 7 :

 7 Figure 7: Representation of the linear cohesive law.

Figure 9 :

 9 Figure 9: Influence of the mesh refinement on the stress-strain curves: a/ ε = 0.1/s b/ ε = 10/s

Figure 10 :

 10 Figure 10: Influence of the mesh refinement on the crack path for two loading rates.

Figure 11 :

 11 Figure 11: Large vs small size of the specimen

Figure 12 :Figure 13 :

 1213 Figure 12: Influence of the specimen size on the fracture path: a/ ε = 1/s b/ ε = 100/s

Figure 14 :

 14 Figure 14: Influence of the specimen size on the dissipated fracture energy: a/ total energy, b/ normalized energy

Figure 15 :

 15 Figure 15: Influence of the meso-structure on the stress-strain curves for different strain rates: a/ ε = 0.1 /s, b/ ε = 1 /s, c/ ε = 10 /s, d/ ε = 100 /s

Figure 16 :

 16 Figure 16: Influence of the meso-structure on the crack path for different meso-structures at ε = 0.1 /s: upper row at ε = 7.10 -5 for ≈ the peak stress (disp. × 2000), lower row at ε = 3.10 -4 (disp. × 100)

Figure 18 :

 18 Figure 18: Numerical rate effect in tension for three meso-structures: a/ Meso 1, b/ Meso 3, c/ Meso 5

Figure 19 :

 19 Figure 19: Influence of the meso-structure on the dissipated fracture energy for different strain rates: a/ Meso 1, b/ Meso 5

Figure 20 :

 20 Figure 20: Influence of the loading rate on the final cracking for the meso-structure Meso 1 (disp. × 100): a/ ε = 0.1 /s, b/ ε = 100 /s

Figure 21 :

 21 Figure 21: Influence of the mesh refinement on the dissipated fracture energy for ε = 100 /s

Table 1 :

 1 Strength, fracture energy and dynamic/static ratio for concrete at different loading rates (from[START_REF] Weerheijm | Research developments and experimental data on dynamic concrete behaviour[END_REF]) Load Rate Tens. Strength Fract. Energy f t(dyn) /f t(stat) G f (dyn) /G f (stat)

		σ (GP a/s)	f t (MP a)	G f (N/m)
	Static	10 -4	3.3	120	1	1
	SHB	39	5.58	120	1.7	1
	Mod. SHB	1685	17.	1505	5.2	12.5

Table 2 :

 2 Material properties used for the mesh sensitivity study

	Aggregate Mortar paste ITZ

Table 3 :

 3 Cohesive length zone obtained with material properties of Table2

	Aggregate Mortar paste ITZ

Contrary to[START_REF] Snozzi | Influence of the meso-structure in dynamic fracture simulation of concrete under tensile loading[END_REF] this value is estimated from G c of the ITZ with δ c equal to the value obtained for the mortar. With this assumption one obtains a f ct smaller in the ITZ than in the mortar, which is commonly accepted for a classical concrete

Acknowledgments

This material is based on the work supported by the Swiss National Foundation under Grant no. 200021 122046/1. We want to thank the ENAC school at EPFL for supporting the sabbatical visit of the first author.


Mathieu Beau
email: mbeau@stp.dias.ie

Baptiste Savoie
email: baptiste.savoie@gmail.com.

Gaussian decay for a difference of traces of the Schrödinger semigroup associated to the isotropic harmonic oscillator

Keywords: 35J10, 47D08, 81Q10, 81Q15, 82B10 Quantum harmonic oscillator, Gibbs semigroups, Mehler's formula, Duhamel-like formula, Geometric perturbation theory

This paper deals with the derivation of a sharp estimate on the difference of traces of the one-parameter Schrödinger semigroup associated to the quantum isotropic harmonic oscillator. Denoting by H∞,κ the self-adjoint realization in

L ) e -tH L,κ , t > 0 has for L sufficiently large a Gaussian decay in L. Furthermore, the estimate that we derive is sharp in the two following senses: its behavior when t ↓ 0 is similar to the one given by Tr L 2 (R d ) e -tH∞,κ = (2 sinh( κ 2 t)) -d and the exponential decay in t arising from Tr L 2 (R d ) e -tH∞,κ when t ↑ ∞ is preserved. For illustrative purposes, we give a simple application within the framework of quantum statistical mechanics.

1 Introduction and the main result. Within the framework of statistical mechanics, the thermodynamic description of large macroscopic equilibrium systems is obtained by considering the thermodynamic (i.e. infinite volume) limit of finite systems described by suitable Gibbs ensembles (e.g. micro-canonical, canonical, grand-canonical). Whenever the thermodynamic limit exists, the surface effects disappear and one is left with the bulk properties. If the limit only depends on the intensive parameters then the system has the extensive property, i.e. the thermodynamic quantities are asymptotically proportional to the system size. Investigating the existence of the thermodynamic limit is part of fundamental topics addressed in rigorous statistical mechanics. It is usually a non-trivial problem which is even more difficult that the interactions involved are singular. The existence of the thermodynamic limit depends on the nature of the interactions, and therefore, may not exist. There is a wide literature on this topic, we refer to [START_REF] Ruelle | Statistical Mechanics, Rigorous results[END_REF] for an overview. Whenever the limit exists, one has to prove in addition that it is independent of the boundary conditions and of the sequence of confining boxes considered. A related topic consists in estimating the boundary contributions. This allows one to answer the question: how negligible is the error made by replacing the finitevolume statistical quantities with the corresponding thermodynamic limit?

One of the most widespread models in physics is the isotropic quantum harmonic oscillator. Turning to the statistical description of d-dimensionnal (d ∈ {1, 2, 3}) harmonically trapped gases, the finite-volume quantities in the grand-canonical ensemble can be written under certain conditions in terms of the trace of the one-parameter 'finite-volume' semigroup (the 'inverse temperature' plays the role of parameter) generated by the Hamiltonian (1.1) with Dirichlet boundary conditions. Without loss of generality, we here consider cubic confining boxes. For such systems, the existence of the conventional thermodynamic limit is trivial since the one-parameter 'infinitevolume' semigroup generated by the Hamiltonian (1.2) on the whole space is trace class. This stems from the discrete nature of the spectrum of (1.2) along with the explicit knowledge of its eigenvalues, see (1.3)-(1.5) below. The error made by replacing the trace of the finite-volume semigroup with the trace of the corresponding infinite-volume semigroup is expected to be negligible on the basis of an heuristical argument: the confinement arising from the harmonic potential overwhelms the confining effects of the 'walls' for sufficiently large boxes. From the foregoing, a natural question arises: how small is this error?

The purpose of this paper is to estimate the remainder term in the asymptotic expansion of the trace of the finite-volume semigroup generated by the Hamiltonian in (1.1) in the large-volume limit. We prove that for any dimension d ∈ {1, 2, 3} the remainder term has a Gaussian decay in the characteristic length L of the cubic confining box for L sufficiently large, see Theorem 1.1 below. To our best knowledge, there is no proof of such a result in literature. A series of remarks placed below Theorem 1.1 discuss how the estimate is sharp. Its proof relies on a geometric perturbation theory applied to bounded operators on separable Hilbert spaces. Such a method has been originally developed to treat approximations of resolvent operators in [START_REF] Cornean | The Faraday effect revisited: Thermodynamic limit[END_REF], see also [START_REF] Cornean | Sharp trace asymptotics for a class of 2D-magnetic operators[END_REF][START_REF] Savoie | On the atomic orbital magnetism: A rigorous derivation of the Larmor and Van Vleck contributions[END_REF] for further applications. Our paper extends the method to treat approximations of semigroups.

For the sake of completeness, we mention that an asymptotic expansion in the high-temperature regime of the trace of the finite-volume semigroup (generated by operators of type (1.1)) was derived in [START_REF] Van Den Berg | Bounds on Green's functions of second-order differential equations[END_REF]Thm. 3]. The method is essentially based on the Feynman-Kac formula and Wiener estimates. However, no information on the large-volume behavior can be inferred from. Note that [START_REF] Van Den Berg | Bounds on Green's functions of second-order differential equations[END_REF]Thm. 3] was applied to investigate the thermodynamics of quantum gases confined by weak external potentials, see [START_REF] Van Den Berg | On the free Boson gas in a weak external potential[END_REF] for further details.

The setting and the main result.

For any d ∈ {1, 2, 3} and

L ∈ (0, ∞), denote Λ d L := {x ∈ R d : -L 2 < x l < L 2 , l = 1, . . . , d} and |Λ d L | its Lebesgue-measure. On C ∞ 0 (Λ d L )
, define ∀κ > 0 the family of operators:

H L,κ := 1 2 (-i∇ x ) 2 + 1 2 κ 2 |x| 2 . (1.1)
By standard arguments, ∀κ > 0 (1.1) extends to a family of self-adjoint and bounded from below operators ∀L ∈ (0, ∞), denoted again by H L,κ , with domain

D(H L,κ ) = W 1,2 0 (Λ d L ) ∩ W 2,2 (Λ d L ).

Obviously this definition corresponds to choose Dirichlet boundary conditions on the boundary

∂Λ d L . Since the inclusion W 1,2 0 (Λ d L ) ֒→ L 2 (Λ d L )
is compact, then ∀κ > 0 H L,κ has a purely discrete spectrum with an accumulation point at infinity. When Λ d L fills the whole space (i.e. when L ↑ ∞), define ∀κ > 0 on C ∞ 0 (R d ) the operator:

H ∞,κ := 1 2 (-i∇ x ) 2 + 1 2 κ 2 |x| 2 . (1.2)
From [14, Thm. X.28], ∀κ > 0 (1.2) is essentially self-adjoint and its self-adjoint extension, denoted again by H ∞,κ , is semi-bounded. By [15, Thm. XIII.16], the spectrum of H ∞,κ is purely discrete with eigenvalues increasing to infinity. From the one-dimensional problem, the eigenvalues and eigenfunctions of the multidimensional case can be written down explicitly. The eigenvalues of the one-dimensional problem are all simple (i.e. non-degenerate) and given by, see e.g. [4, Sec. 1.8]:

ǫ (s) ∞,κ := κ s + 1 2 , s ∈ N. (1.3)
The corresponding eigenfunctions which form an orthonormal basis in L 2 (R), read as:

∀x ∈ R, φ (s) ∞,κ (x) := 1 √ 2 s s! κ π 1 4 e -κ 2 x 2 H s √ κx , s ∈ N, (1.4)
where H s , s ∈ N are the Hermite polynomials defined by: H s (x) := (-1) s e x 2 d s dx s (e -x 2 ), ∀x ∈ R. The eigenvalues and eigenfunctions of the multidimensional case (i.e. d = 2, 3) are respectively related to those of the one-dimensional case by:

E (s) ∞,κ := d j=1 ǫ (sj ) ∞,κ = κ d j=1 s j + 1 2 , s = {s j } d j=1 ∈ N d , (1.5)
ψ (s) ∞,κ (x) := d j=1 φ (sj ) ∞,κ (x j ), x = {x j } d j=1 ∈ R d . (1.6)
From (1.3)-(1.5) and by the use of the min-max principle, one has for any L ∈ (0, ∞):

∀κ > 0, inf σ (H L,κ ) ≥ inf σ (H ∞,κ ) = E (0) ∞,κ = dǫ (0) ∞,κ > 0, ǫ (0) ∞,κ := κ 2 .
Let us turn to the one-parameter strongly-continuous semigroup (the so-called C 0 -semigroup in the Hille-Phillips terminology [START_REF] Hille | Functional Analysis and Semi-Groups, Revised and Expanded Edition[END_REF]) generated by the operators introduced above. At finite-volume, it is defined ∀L ∈ (0, ∞) and ∀κ > 0 by {G L,κ (t) := e -tHL,κ :

L 2 (Λ d L ) → L 2 (Λ d L )} t≥0 .
It is a self-adjoint and positive operator on L 2 (Λ d L ) by the spectral theorem and the functional calculus, see e.g. [START_REF] Reed | Methods of Modern Mathematical Physics, II : Fourier Analysis, Self-Adjointness[END_REF]Sec. X.8]. The same hold true for the one-parameter semigroup on the whole space {G ∞,κ (t) := e -tH∞,κ :

L 2 (R d ) → L 2 (R d )} t≥0 . Moreover ∀0 < L ≤ ∞, ∀κ > 0 and ∀t > 0, G L,κ (t) is a Gibbs semigroup, i.e. G L,κ (t) (resp. G ∞,κ (t)) belongs to the Banach space of trace-class operators on L 2 (Λ d L ) (resp. L 2 (R d ))
, see e.g. [START_REF] Uhlenbrock | Perturbation of statistical semigroups in quantum statistical mechanics[END_REF][START_REF] Angelescu | On the perturbation of Gibbs semigroups[END_REF]. A basic feature is the monotonicity property for the finite-volume trace, see Lemma 4.4 in Sec. 4:

∀L ∈ (0, ∞), Tr L 2 (Λ d L ) {G L,κ (t)} ≤ Tr L 2 (R d ) {G ∞,κ (t)} = 2 sinh κ 2 t -d , κ > 0, t > 0.
Our main result is the following sharp estimate on the difference of traces of the semigroups:

Theorem 1.1. For any d ∈ {1, 2, 3} there exists a constant C d > 0 and ∀0 < κ 0 < 1 there exists a real L κ0 > 0 s.t. ∀L ∈ [L κ0 , ∞), ∀κ ∈ [κ 0 , ∞) and ∀t > 0:

Tr L 2 (Λ d L ) {G L,κ (t)} -Tr L 2 (R d ) {G ∞,κ (t)} ≤ C d 1 + √ κ (1 + κ) d (1 + t) 3(d+ 1 2 ) Tr L 2 (R d ) {G ∞,κ (t)} e -κ 32 L 2 4 tanh( κ 2 t) . (1.7)
Remark 1.2. The estimate is sharp in the following senses. Its behavior when t ↓ 0 is given by the term Tr L 2 (R d ) {G ∞,κ (t)} = (2 sinh( κ 2 t)) -d . We recall that sinh(x) ∼ x when x ↓ 0. Moreover, the r.h.s. of (1.7) has an exponential decay in t when t ↑ ∞ arising from Tr L 2 (R d ) {G ∞,κ (t)}. Note that one can get rid of the polynomial growth in t appearing in the r.h.s. of (1.7) via (4.13).

Remark 1.3. The L in Theorem 1.1 is independent of κ whenever κ ∈ [1, ∞). If 0 < κ < 1, then the estimate holds for L large enough chosen accordingly (i.e. L ≥ cste × κ -1 2 ).
Remark 1.4. In (1.7), the powers on the factors (1 + √ κ), (1 + κ), (1 + t) and the constant appearing in the argument of the exponential can be optimized.

An application in quantum statistical mechanics.

Consider a d-dimensional ideal quantum gas composed of a large number of non-relativistic spin-0 identical particles confined in the box Λ d L and trapped in an isotropic harmonic potential. Such a system is considered to figure out the Bose-Einstein condensation phenomenon created by cold alkali atom gases in magnetic-optical trap, see e.g. [START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF]Chap. 10] and references therein. Within the one-body approximation, the dynamics of a single Boson is determined by (1.1). Suppose that the system is at equilibrium with a thermal and particles bath. In the grandcanonical ensemble, let (β, z, |Λ d L |) be the external parameters. Here, β := (k B T ) -1 > 0 is the 'inverse' temperature (k B stands for the Boltzmann constant) and z = e βµ the fugacity (µ is the chemical potential). The finite-volume single-particle partition function is defined as, see e.g. [START_REF] Ruelle | Statistical Mechanics, Rigorous results[END_REF]:

Φ L,κ (β) := Tr L 2 (Λ d L ) {G L,κ (β)} , β > 0. (1.8)
The grand-canonical average number of particles at finite-volume is related to (1.8) by, see [START_REF] Beau | Rigorous investigation of the reduced density matrix for the ideal Bose gas in harmonic traps by a loop-gas-like approach[END_REF]: (HL,κ) .

N L,κ (β, z) := ∞ l=1 z l Φ L,κ (lβ), β > 0, z ∈ 0, e β inf σ
(1.9) Theorem 1.1 allows one to get the large-volume behavior of the quantities in (1.8)-(1.9). Indeed, one gets ∀0 < κ 1 < κ 2 < ∞, ∀0 < β 1 < β 2 < ∞ and for any compact subset K ⊂ (0, e β1E (0) ∞,κ 1 ):

Φ ∞,κ (β) := lim L↑∞ Φ L,κ (β) = Tr L 2 (R d ) {G ∞,κ (β)} = e -βE (0) ∞,κ (1 -e -βκ ) d , N ∞,κ (β, z) := lim L↑∞ N L,κ (β, z) = ∞ l=1 z l Φ ∞,κ (lβ), uniformly in (κ, β, z) ∈ [κ 1 , κ 2 ] × [β 1 , β 2 ] × K.
Moreover, one has the following asymptotics:

Φ L,κ (β) = Φ ∞,κ (β) + O e -cL 2 , N L,κ (β, z) = N ∞,κ (β, z) + O e -cL 2 ,
for some L-independent constant c = c(κ, β) > 0. We emphasize that the upper bound in (1.7) plays a crucial to prove the thermodynamic limit of (1.9) for any z ∈ (0, e βE (0) ∞,κ ), see [3, Sec. A].

2 Proof of Theorem 1.1.

The starting-point consists in rewriting the difference between the traces involving the difference between the semigroup integral kernels. We refer the reader to Sec. 4 in which we have collected some basic properties on the semigroup kernel. Since ∀L ∈ (0, ∞] and ∀κ > 0, {G L,κ (t)} t>0 is a Gibbs semigroup with a jointly continuous integral kernel G

(d) L,κ (• , • ; t) : R d × R d → C, then: Tr L 2 (R d ) {G ∞,κ (t)} -Tr L 2 (Λ d L ) {G L,κ (t)} = Y (d) L,κ (t) + Z (d) L,κ (t),
with ∀d ∈ {1, 2, 3}, ∀L ∈ (0, ∞), ∀κ > 0 and ∀t > 0:

Y (d) L,κ (t) := Λ d L G (d) ∞,κ (x, x; t) -G (d) L,κ (x, x; t) dx, (2.1)
Z (d) L,κ (t) := R d \Λ d L G (d) ∞,κ (x, x; t) dx. (2.2)
Here, we used [START_REF] Cornean | On the magnetization of a charged Bose gas in the canonical ensemble[END_REF]Prop. 9]. Note that ∀κ > 0, ∀t > 0 the kernel G

(d) ∞,κ (• , • t)
is explicitly known and it is given by Mehler's formula, see (4.3)-(4.4). It is derived from (1.4)-(1.6) and (1.5).

Next, it remains to estimate each one of the above quantity. For the quantity in (2.2):

Lemma 2.1. ∀d ∈ {1, 2, 3}, ∀L ∈ (0, ∞), ∀κ > 0 and ∀t > 0:

Z (d) L,κ (t) ≤ 2 sinh κ 2 t -d e -dκ L 2 4 tanh( κ 2 t) .
Proof. Let β > 0 and κ > 0 be fixed. Because of (4.4), it is enough to treat only the case of d = 1. From (4.3) and by setting x = y, one has by direct computations:

∀L ∈ (0, ∞), ∀t > 0, Z (d=1) L,κ (t) = erfc κ tanh κ 2 t L 2 2 sinh(κt) tanh κ 2 t
, where erfc denotes the complementary error function, see e.g. [1, Eq. (7.1.2)]. From the Chernoff inequality which reads as erfc(α) ≤ e -α 2 ∀α ≥ 0, along with the identity:

sinh(αt) = 2 sinh α 2 t cosh α 2 t , ∀α ≥ 0, ∀t > 0, (2.3)
one arrives at:

∀L ∈ (0, ∞), ∀t > 0, Z (d=1) L,κ (t) ≤ 2 sinh κ 2 t -1 e -κ L 2 4 tanh( κ 2 t) .
As for the the quantity defined in (2.1), we establish the following estimate:

Proposition 2.2. For any d ∈ {1, 2, 3}, there exists a constant C d > 0 and ∀0 < κ 0 < 1 there exists a real

L κ0 > 0 s.t. ∀L ∈ [L κ0 , ∞), ∀κ ∈ [κ 0 , ∞) and ∀t > 0: Y (d) L,κ (t) ≤ C d 1 + √ κ (1 + κ) d (1 + t) 3(d+ 1 2 ) 2 sinh κ 2 t -d e -κ 32 L 2 4 tanh( κ 2 t) . (2.4)
Theorem 1.1 follows from Lemma 2.1 and Proposition 2.2 together. The rest of this section is devoted to the proof of Proposition 2.2.

Proof of Proposition 2.2.

In view of (2.1), the first step consists in writing an expression for the difference between the two semigroup kernels. It is contained in the following lemma: Lemma 2.3. ∀L ∈ (0, ∞), ∀κ > 0 and ∀t > 0:

∀(x, y) ∈ Λ 2 L , G (1) ∞,κ (x, y; t) -G (1)
L,κ (x, y; t) = - 1 2 t 0 G (1) ∞,κ (x, - L 2 ; s) ∂ z G (1)
L,κ (-

L 2 , y; t -s) -G (1) ∞,κ (x, L 2 ; s) ∂ z G (1) L,κ ( L 2 , y; t -s) ds, (2.5)
and in the case of d = 2, 3, for any (x, y) ∈ Λ 2d L :

G (d) ∞,κ (x, y; t) -G (d) L,κ (x, y; t) = - 1 2 t 0 ∂Λ d L G (d) ∞,κ (x, z; s) n z • ∇ z G (d) L,κ (z, y; t -s) dσ(z) ds, (2.6
) where dσ(z) denotes the measure on ∂Λ d L and n z the outer normal to ∂Λ d L at z. The proof of Lemma 2.3 in the case of d = 3 can be found in [START_REF] Briet | Diamagnetic expansions for perfect quantum gases[END_REF]Lem. 4.2], see also [START_REF] Briet | Diamagnetic expansions for perfect quantum gases II: Uniform bounds[END_REF]. Since the generalization to d = 1, 2 can be easily obtained by similar arguments, we do not give any proof.

Recall that the kernel G (d)

∞,κ is explicitly known and given by Mehler's formula. In view of (2.1) along with the expressions (2.5)-(2.6), the actual problem comes down to derive a sufficiently sharp estimate on the gradient of the finite-volume semigroup kernel allowing us to bring out a Gaussian decay in L for the quantity in (2.1). It is contained in the following proposition: Proposition 2.4. For any d ∈ {1, 2, 3}, there exists a constant C d > 0 and ∀0 < κ 0 < 1 there exists a

L κ0 > 0 s.t. ∀L ∈ [L κ0 , ∞), ∀κ ∈ [κ 0 , ∞), ∀(x, y) ∈ Λ 2d
L and ∀t > 0:

∇ x G (d) L,κ (x, y; t) ≤ C d P (d) ∞,κ (x, y; t) + R (d) L,κ (x, y; t) , (2.7)
P (d) ∞,κ (x, y; t) := (1 + √ κ)(1 + t) 5 2 coth κ 2 t G (d) ∞,κ (x, y; t, 8); (2.8) R (d) L,κ (x, y; t) := κ d 2 (1 + κ) d 2 t d-1 2 (sinh(κt)) d 2 (1 + t) 5d 2 +1 e -κ 16 L 2 4 tanh( κ 2 t) G (d)
∞,0 (x, y; 4t).

(2.9)

Here, G (d) ∞,κ (• , • ; t, γ), κ > 0 and γ > 0 is defined in (4.7) and G (d) ∞,0 (• , • ; t) in (4.2)-(4.4
). Remark 2.5. The L in Proposition 2.4 can be chosen uniformly in κ whenever κ ≥ 1. If 0 < κ < 1 then the estimate holds for L large enough chosen accordingly (i.e. L ≥ cste/ √ κ).

Note that Proposition 2.4 actually contains the key-estimate of this paper; its proof is placed in Sec. 2.2. We mention that the derivation of such an estimate leans on a Duhamel-like formula for the finite-volume semigroup G L,κ (t), L ∈ (0, ∞) obtained via a geometric perturbation theory.

Proof of Proposition 2.2. Let 0 < κ 0 < 1 be fixed. Denote ς L = ±L/2. Start with the case of d = 1. In view of (2.1), (2.5) and (2.7), we need to estimate ∀L ∈

[L κ0 , ∞) and ∀κ ∈ [κ 0 , ∞): ∀t > 0, Y (d=1),1 L,κ (t) := 1 2 t 0 Λ 1 L G (d=1) ∞,κ (x, ς L ; s, 1)P (d=1) ∞,κ (ς L , x; t -s) dx ds, (2.10)
Y (d=1),2 L,κ (t) := 1 2 t 0 Λ 1 L G (d=1) ∞,κ (x, ς L ; s, 1)R (d=1) L,κ (ς L , x; t -s) dx ds. (2.11)
Here, we have commuted the two integrals; this will be justified by what follows. We first estimate (2.10). In view of (4.3) and (2.8), then from (4.16) for any L ∈ [L κ0 , ∞), κ ∈ [κ 0 , ∞) and t > 0:

Y (d=1),1 L,κ (t) ≤ C √ κ 1 + √ κ (1 + t) 5 2 1 2 sinh(κt) e -κ 8 L 2 4 tanh( κ 2 t) t 0 coth κ 2 (t -s) ds, (2.12)
for some constant C > 0. By using the upper bound of the following inequality:

1 α ≤ coth(α) := 1 tanh(α) ≤ 1 + α α , α > 0, (2.13)
along with the inequality:

1 2 sinh(κt) = tanh κ 2 t 2 sinh(κt) tanh κ 2 t = tanh κ 2 t 2 sinh κ 2 t ≤ 1 2 sinh κ 2 t , (2.14)
justified by (2.3), then there exists another constant

C > 0 s.t. ∀L ∈ [L κ0 , ∞) and ∀κ ∈ [κ 0 , ∞): ∀t > 0, Y (d=1),1 L,κ (t) ≤ C 1 + √ κ √ 1 + κ (1 + t) 7 2 2 sinh κ 2 t e -κ 8 L 2 4 tanh( κ 2 t) . (2.15)
Next, we estimate (2.11). From (2.9) and ( 4.

3), one has ∀L ∈ [L κ0 , ∞), ∀κ ∈ [κ 0 , ∞) and ∀t > 0: Y (d=1),2 L,κ (t) ≤ √ κ √ 1 + κ(1 + t) 7 2 t 0 1 sinh(κ(t -s)) e -κ 16 L 2 4 tanh( κ 2 (t-s)) ds × Λ 1 L G (d=1) ∞,κ (x, ς L ; s, 1)G (d=1) ∞,0 (ς L , x; 4(t -s)) dx.
For the following, we need to make appear from the integration over Λ 1 L a Gaussian decay in L while having the argument s. To do so, let us remark that on R 2d , d ∈ {1, 2, 3} one has:

∀s > 0, G (d) ∞,κ (x, y; s, 1) ≤ e -κ 4 (|x| 2 +|y| 2 ) tanh( κ 2 s) G (d) ∞,κ (x, y; s, 2). (2.16)
To get (2.16), we expanded in (4.3) the squares and used that 2ab ≤ (a 2 + b 2 ) with the fact that coth(α) -tanh(α) ≥ 0 ∀α > 0. Now, from (2.16) followed by the lower bound:

tanh(αs) + tanh(α(t -s)) = tanh(αt) {1 + tanh(αs) tanh(α(t -s))} ≥ tanh(αt), ∀α ≥ 0, ∀t > s > 0, (2.17)
then by using the upper bound in the first inequality of (4.8) along with (4.15), one has ∀L ∈

[L κ0 , ∞) and ∀κ ∈ [κ 0 , ∞): ∀t > 0, Y (d=1),2 L,κ (t) ≤ Cκ √ 1 + κ (1 + t) 7 2 √ t e -κ 16 L 2 4 tanh( κ 2 t) t 0 √ s sinh(κs) sinh(κ(t -s))
ds, for some constant C > 0. It remains to use successively the identity:

coth(αs) + coth(α(t -s)) = sinh(αt) sinh(αs) sinh(α(t -s))
, ∀α > 0, ∀t > s > 0, (2.18) followed by (2.13) and (2.14). This yields:

∀t > 0, Y (d=1),2 L,κ (t) ≤ C √ κ(1 + κ) √ t (1 + t) 4 2 sinh κ 2 t e -κ 16 L 2 4 tanh( κ 2 t) , (2.19)
for another constant C > 0. Gathering (2.15)-(2.19) together, we get (2.4) in the case of d = 1. Subsequently, we turn to the case of d = 2. Since (2.6) is made up of four terms, then the same holds for the quantity in (2.1). Since these terms have the same structure, it is enough to treat only one of them. In view of (2.7), we need to estimate ∀L ∈ [L κ0 , ∞), ∀κ ∈ [κ 0 , ∞) and ∀t > 0:

Y (d=2),1 L,κ (t) := 1 2 t 0 Λ 2 L Λ 1 L G (d=2) ∞,κ (x, (z 1 , ς L ); s, 1)P (d=2) ∞,κ ((z 1 , ς L ), x; t -s) dz 1 dx ds, (2.20)
Y (d=2),2 L,κ (t) := 1 2 t 0 Λ 2 L Λ 1 L G (d=2) ∞,κ (x, (z 1 , ς L ); s, 1)R (d=2) L,κ ((z 1 , ς L ), x; t -s) dz 1 dx ds. (2.21)
The strategy consists in using the property (4.4) in order to use the results stated in the case of d = 1. Let us first estimate the quantity in (2.20). In view of (4.4) and (2.8), then from (4.16):

Λ 2 L Λ 1 L G (d=2) ∞,κ (x, (z 1 , ς L ); s, 1)G (d=2) ∞,κ ((z 1 , ς L ), x; t -s, 8) dz 1 dx ≤ C R 1 G (d=1) ∞,κ (x 1 , x 1 ; t, 8) dx 1 R 1 G (d=1) ∞,κ (x 2 , ς L ; s, 1)G (d=1) ∞,κ (ς L , x 2 ; t -s, 8) dx 2 ,
for some constant C > 0. From (4.19), the first integral in the above r.h.s. is nothing but the trace (multiplied by a constant). Then for any L ∈ [L κ0 , ∞), κ ∈ [κ 0 , ∞) and t > 0, we arrive at:

Y (d=2),1 L,κ (t) ≤ C 1 + √ κ (1 + t) 5 2 2 sinh κ 2 t × t 0 coth κ 2 (t -s) ds R 1 G (d=1) ∞,κ (x 2 , ς L ; s, 1)G (d=1) ∞,κ (ς L , x 2 ; t -s, 8) dx 2 ,
for another C > 0. The integral w.r.t. s has been estimated in the case of d = 1, see (2.12). Therefore, it remains to mimic the arguments leading to (2.15) to conclude. Next, we estimate (2.21). In view of (4.3)-(4.4) and (2.9), from the first upper bound in (4.8) along with (4.16):

Λ 2 L Λ 1 L G (d=2) ∞,κ (x, (z 1 , ς L ); s, 1)G (d=2) ∞,0 ((z 1 , ς L ), x; 4(t -s)) dz 1 dx ≤ C κ sinh(κs) √ s × Λ 1 L G (d=1) ∞,0 (x 1 , x 1 ; 4t) dx 1 R 1 G (d=1) ∞,κ (x 2 , ς L ; s, 1)G (d=1) ∞,0 (ς L , x 2 ; 4(t -s)) dx 2 ,
for some constant C > 0. Note that the integrand in the first integral of the above r.h.s. is nothing but a constant. This will make appear a factor L, but we will get rid of it at the end. Ergo, in view of (2.9) and (4.3), there exists another

C > 0 s.t. ∀L ∈ [L κ0 , ∞), ∀κ ∈ [κ 0 , ∞) and ∀t > 0: Y (d=2),2 L,κ (t) ≤ Cκ 3 2 (1 + κ)L (1 + t) 6 √ t t 0 √ s √ t -s sinh(κs) sinh(κ(t -s)) e -κ 16 L 2 4 tanh( κ 2 (t-s)) ds × R 1 G (d=1) ∞,κ (x 2 , ς L ; s, 1)G (d=1) ∞,0 (ς L , x 2 ; 4(t -s)) dx 2 .
The rest of the proof mimics the strategy used for the case of d = 1. By using the upper bound in the first inequality of (4.8) along with (4.15), one has ∀L ∈ [L κ0 , ∞), ∀κ ∈ [κ 0 , ∞) and ∀t > 0:

Y (d=2),2 L,κ (t) ≤ Cκ 2 (1+κ)L (1 + t) 6 t t 0 s √ t -s sinh(κs) sinh(κ(t -s)) e -κ 16 L 2 4 tanh( κ 2 (t-s)) e -κ 4 L 2 4 tanh( κ 2 s) ds.
By using successively (2.17), (2.18) and (2.50) one straightforwardly gets:

Y (d=2),2 L,κ (t) ≤ Cκ(1 + κ) 2 L √ t (1 + t) 7 2 sinh(κt) e -κ 16
L 2 4 tanh( κ 2 t) ,
for another L-independent C > 0. It remains to use (4.13) to get rid of the L-factor:

L 2 sinh(κt) e -κ 32 L 2 4 tanh( κ 2 t) ≤ C √ κ 1 2 sinh(κt) tanh κ 2 t ≤ C √ κ 1 2 sinh κ 2 t 2 .
Gathering the above estimates together, one arrives ∀L ∈ [L κ0 , ∞), ∀κ ∈ [κ 0 , ∞) and ∀t > 0 at:

Y (d=2),2 L,κ (t) ≤ C √ κ(1 + κ) 2 √ t (1 + t) 7 2 sinh κ 2 t 2 e -κ 32 L 2 4 tanh( κ 2 t) ,
for another constant C > 0. The case of d = 3 can be deduced by similar arguments.

Proof of Proposition 2.4.

As it was previously mentioned, Proposition 2.4 contains the key-estimate to prove Theorem 1.1. The proof leans on an approximation of the finite-volume semigroup operator via a geometric perturbation theory. Although this method had been originally developed for the resolvent operators, see [START_REF] Cornean | The Faraday effect revisited: Thermodynamic limit[END_REF] and also [START_REF] Cornean | Sharp trace asymptotics for a class of 2D-magnetic operators[END_REF][START_REF] Savoie | On the atomic orbital magnetism: A rigorous derivation of the Larmor and Van Vleck contributions[END_REF], below we extend the method to the semigroup operators.

An approximation via a geometric perturbation theory.

The key-idea consists in isolating in Λ d L the region close to the boundary from the bulk where the semigroup G ∞,κ (t) will act. The underlying difficulty is to keep a good control of the remainder terms arising from this approximation. This will be achieved by using well-chosen cutoff functions.

For any 0 < η < 1, 0 < ϑ ≤ 1000, d ∈ {1, 2, 3} and L ∈ (0, ∞) define:

Θ L,η (ϑ) := x ∈ Λ d L : dist x, ∂Λ d L ≤ ϑL η . (2.22)
For L sufficiently large, Θ L,η (ϑ) models a 'thin' compact subset of Λ d L near the boundary with Lebesgue-measure |Θ L,η (ϑ)| of order O(L (d-1)+η ). For any 0 < η < 1, let

L 0 = L 0 (η) ≥ 1 s.t. Θ L0,η (1000) Λ d L0 , L 0 -L η 0 ≥ L 0 / √ 2, (2.23)
and L 0 large enough. Let us now introduce some well-chosen families of smooth cutoff functions. Let f L,η and f c L,η , L ∈ [L 0 (η), ∞) be a partition of the unity of Λ d L satisfying:

f L,η + f c L,η = 1 on Λ d L ; Supp (f L,η ) ⊂ Λ d L \ Θ L,η 1 16 , f L,η = 1 if x ∈ Λ d L \ Θ L,η 1 8 , 0 ≤ f L,η ≤ 1; Supp f c L,η ⊂ Θ L,η 1 8 , f c L,η = 1 if x ∈ Θ L,η 1 16 .
Moreover, there exists a constant C > 0 s.t.

∀L ≥ L 0 (η), D σ f L,η ∞ ≤ CL -|σ|η , ∀|σ| ≤ 2, |σ| = σ 1 + • • • + σ d .
Also, let fL,η and fL,η , L ∈ [L 0 (η), ∞) satisfying:

Supp fL,η ⊂ Λ d L \ Θ L,η 1 64 , fL,η = 1 if x ∈ Λ d L \ Θ L,η 1 32 , 0 ≤ fL,η ≤ 1; Supp fL,η ⊂ Θ L,η 1 2 , fL,η = 1 if x ∈ Θ L,η 1 4 , 0 ≤ fL,η ≤ 1.
Moreover, there exists another constant C > 0 s.t.

∀L ≥ L 0 (η), max D σ fL,η ∞ , D σ fL,η ∞ ≤ CL -|σ|η , ∀|σ| ≤ 2.
With these properties, one straightforwardly gets:

fL,η f L,η = f L,η ; (2.24) dist Supp D σ fL,η , Supp (D τ f L,η ) ≥ CL η , ∀1 ≤ |σ| ≤ 2, ∀0 ≤ |τ | ≤ 2; (2.25) fL,η f c L,η = f c L,η ; (2.26) dist Supp D σ fL,η , Supp D τ f c L,η ≥ CL η , ∀1 ≤ |σ| ≤ 2, ∀0 ≤ |τ | ≤ 2, (2.27)
for some L-independent constant C > 0.

Afterwards, let us define

∀0 < η < 1, ∀L ∈ [L 0 (η), ∞) (see (2.23)) and ∀κ > 0 on C ∞ 0 (Λ d L ): h L,κ,η := 1 2 (-i∇ x ) 2 + 1 2 κ 2 V L,η (x), V L,η (x) := |x| 2 , if x ∈ Supp fL,η , 1 4 (L -L η ) 2 , otherwise.
(2.28)

By standard arguments, (2.28) extends to a family of self-adjoint and semi-bounded operators

∀L ∈ [L 0 (η), ∞), denoted again by h L,κ,η , with domain D(h L,κ,η ) = W 1,2 0 (Λ d L ) ∩ W 2,2 (Λ d L ). For any 0 < η < 1, L ∈ [L 0 (η), ∞) and κ > 0, let {g L,κ,η (t) := e -thL,κ,η : L 2 (Λ d L ) → L 2 (Λ d L )
} t≥0 be the strongly-continuous one-parameter semigroup generated by h L,κ,η . It is an integral operator with an integral kernel jointly continuous in (x, y, t)

∈ Λ d L × Λ d L × (0, ∞).
We denote it by g

(d) L,κ,η . Next, introduce ∀0 < η < 1, ∀L ∈ [L 0 (η), ∞) and ∀κ > 0 the following operators on L 2 (Λ d L ): ∀t > 0, G L,κ,η (t) := fL,η G ∞,κ (t)f L,η + fL,η g L,κ,η (t)f c L,η , (2.29)
W L,κ,η (t) := - 1 2 ∆ fL,η + i ∇ fL,η • (-i∇) G ∞,κ (t)f L,η + - 1 2 ∆ fL,η + i ∇ fL,η • (-i∇) g L,κ,η (t)f c L,η .
(2.30)

Sometimes, we will use the shorthand notations:

∀t > 0, G (p) L,κ,η (t) := fL,η G ∞,κ (t)f L,η , G (r)
L,κ,η (t) := fL,η g L,κ,η (t)f c L,η .

(2.31)

The main result of this paragraph is the following Duhamel-like formula:

Proposition 2.6. ∀d ∈ {1, 2, 3}, ∀0 < η < 1, ∀L ∈ [L 0 (η), ∞) and ∀κ > 0, it takes place in the bounded operators sense on L 2 (Λ d L ):

∀t > 0, G L,κ (t) = G L,κ,η (t) - t 0 G L,κ (t -s)W L,κ,η (s) ds. (2.32)
The proof of Proposition 2.6 can be found in Sec. 2.2.3; it is essentially based on the application of [START_REF] Cornean | On the magnetization of a charged Bose gas in the canonical ensemble[END_REF]Prop. 3] taking into account the features of the cutoff functions introduced previously.

Remark 2.7. One can derive the following upper bounds on the operator norms. ∀d ∈ {1, 2, 3} there exist two constants C d , c > 0 s.t. ∀0 < η < 1, ∀L ∈ [L 0 (η), ∞), ∀κ > 0 and ∀t > 0:

G L,κ,η (t) ≤ G (p) L,κ,η (t) + G (r) L,κ,η (t) ≤ (cosh(κt)) -d 2 + C d e -κ 2 16 L 2 t , (2.33)
W L,κ,η (t) ≤ C d √ 1 + κ √ 1 + t √ t e -c L 2η t 1 + (1 + t) d-

Conclusion of the proof.

The starting-point in the proof of Proposition 2.4 is the Duhamel-like formula in (2.32). Taking its adjoint, one has ∀d ∈ {1, 2, 3}, ∀0 < η < 1, ∀L ∈ [L 0 (η), ∞) (see (2.23)) and ∀κ > 0 on L 2 (Λ d L ):

∀t > 0, G L,κ (t) = G * L,κ,η (t) - t 0 W * L,κ,η (s)G L,κ (t -s) ds, (2.35)
where the adjoint operator of G L,κ,η (t) and W L,κ,η (t) reads respectively as, see (2.29)-(2.30):

G * L,κ,η (t) = f L,η G ∞,κ (t) fL,η + f c L,η g L,κ,η (t) fL,η , (2.36)
W * L,κ,η (t) = -f L,η G ∞,κ (t) 1 2 ∆ fL,η + if L,η {(-i∇) G ∞,κ (t) -[(-i∇) , G ∞,κ (t)]} ∇ fL,η + -f c L,η g L,κ,η (t) 1 2 ∆ fL,η + if c L,η {(-i∇) g L,κ,η (t) -[(-i∇) , g L,κ,η (t)]} ∇ fL,η . (2.37)
Here, [• , • ] denotes the usual commutator, and in the bounded operators sense:

[(-i∇) , G ∞,κ (t)] = - t 0 G ∞,κ (t -s) [(-i∇) , H ∞,κ ] G ∞,κ (s) ds, (2.38) [(-i∇) , g L,κ,η (t)] = - t 0 g L,κ,η (t -s) [(-i∇) , h L,κ,η ] g L,κ,η (s) ds.
(2.39) Writing (2.35) in the kernels sense, it follows the identity:

∀(x, y) ∈ Λ 2d L , ∀t > 0, ∇ x G (d) L,κ (x, y; t) = ∇ x G * L,κ,η (d) (x, y; t) - t 0 Λ d L ∇ x W * L,κ,η (
i) ∀0 < η < 1, ∀L ∈ [L 0 (η), ∞), ∀κ > 0, ∀(x, y) ∈ Λ 2d
L and ∀t > 0:

∇ x G * L,κ,η (d) (x, y; t) ≤ C d P (d) ∞,κ,η (x, y; t) + R (d)
L,κ,η (x, y; t) ,

P (d) ∞,κ,η (x, y; t) := 1 + √ κ coth κ 2 t G (d) ∞,κ (x, y; t, 2), (2.41) R (d) L,κ,η (x, y; t) := (1 + t) d √ t e -κ 2 8 L 2 4 t G (d)
∞,0 (x, y; 2t).

(2.42)

(ii). ∀ 1 4 < η < 1, ∀L ∈ [L 0 (η), ∞), ∀κ > 0, ∀(x, y) ∈ Λ 2d
L and ∀t > 0:

∇ x W * L,κ,η (d) (x, y; t) ≤ C d r (d) ∞,κ,η (x, y; t) + r (d)
L,κ,η (x, y; t) , (2.43)

r (d) ∞,κ,η (x, y; t) := 1 + √ κ coth κ 2 t (1 + t)e -cκL 2η coth( κ 2 t) G (d) ∞,κ (x, y; t, 8), (2.44) r (d)
L,κ,η (x, y; t) :=

(1 + t) 2d √ t e -κ 2 16 L 2 4 t e -c L 2η t χ ΘL,η( 1 8 ) (x)G (d)
∞,0 (x, y; 4t)χ ΘL,η ( 12 ) (y).

(2.45)

Here, χ ΘL,η(ϑ) , ϑ > 0 denotes the indicator function associated with Θ L,η (ϑ) defined in (2.22).

Remark 2.9. In (ii), η has been restricted to ( 1 4 , 1) only to make the estimates more elegant. Proof of Proposition 2.4. From now on, we set η = 1 2 in the r.h.s. of (2.40). In view of the second term, (2.43) with (2.44)-(2.45) and (4.5), we need to estimate the two quantities:

Q (d) ∞,κ (x, y; t) := t 0 R d r (d) ∞,κ,η= 1 2 (x, z; s)G (d) ∞,κ (z, y; t -s, 1) dz ds, (2.46)
Q (d) L,κ (x, y; t) := t 0 R d r (d) L,κ,η= 1 2 (x, z; s)G (d)
∞,κ (z, y; t -s, 1) dz ds.

(2.47)

Let L ≥ L 0 (η = 1 2 ) defined in (2.23). We start with (2.46). From (2.44) followed by (4.16), there exist two constants c, C > 0 s.t. ∀κ > 0, ∀(x, y) ∈ Λ 2d

L and ∀t > 0:

Q (d) ∞,κ (x, y; t) ≤ C 1 + √ κ (1 + t) coth κ 2 t coth κ 2 t G (d) ∞,κ (x, y; t, 8) t 0 coth κ 2 s e -cκL coth( κ 2 s) ds.
By using (4.13) to get rid of the coth in the integrand, followed by the lower bound in (2.13) for the (artificial) denominator in the above r.h.s., then the upper bound in (2.8) follows. Let us turn to the quantity in (2.47). From (2.45), one has ∀κ > 0, ∀(x, y) ∈ Λ 2d L and ∀t > 0:

Q (d) L,κ (x, y; t) ≤ (1 + t) 2d t 0 e -κ 2 16 L 2 4 s √ s ds R d G (d) ∞,0 (x, z; 4s)χ Θ L, 1 2 ( 1 2 ) (z)G (d) ∞,κ (z, y; t -s, 1) dz.
To make appear a Gaussian decay in L from the integration over R d , we use (2.16). Here, the presence of the characteristic function in the integrand plays a crucial role. Since

∀z ∈ Θ L, 1 2 ( 1 2 ), |z| ≥ 1 2 (L - √ L) (recall that L ≥ L 0 (η = 1 2 ), L 0 (η) as in (2.23)) then ∀(x, y) ∈ Λ 2d
L and ∀t > 0: )s ≤ e -κ 2 0 32 L 2 4 s holds. By using an inequality of type:

Q (d) L,κ (x, y; t) ≤ (1 + t) 2d t 0 1 √ s e -κ 2 16 L 2 4 s e -κ 8 L 2 4 tanh( κ 2 (t-s)) ds × ΘL,η( 1 2 ) G (d) ∞,0 (x, z; 4s)G (d) ∞,κ (z
Q (d) L,κ (x, y; t) ≤ Cκ d (1 + t) 2d G (d) ∞,0 (x, y; 4t) t 0 s d-1 2 (t -s)
e -κ 2 16 L 2 4 s e -κ 8 L 2 4 tanh( κ 2 (t-s)) ≤ e -κ 8 L 2 4 [tanh( κ 2 s)+tanh( κ 2 (t-s))] ≤ e -κ 8 L 2 4 tanh( κ 2 t) , 0 < s < t, (2.49)
which is justified by the lower bound in (2.13) together with (2.17), it follows from the identity in (2.18)

that ∀L ∈ [L κ0 , ∞), ∀κ ∈ [κ 0 , ∞), ∀(x, y) ∈ Λ 2d
L and ∀t > 0:

Q (d) L,κ (x, y; t) ≤ Cκ d (1 + t) 2d (sinh(κt)) d 2 e -κ 16 L 2 4 tanh( κ 2 t) G (d) ∞,0 (x, y; 4t)× × t 0 s d-1 2 (t -s) d 2
(coth(κs))

d 2 + (coth(κ(t -s))) d 2 ds,
for some constant C > 0. Here, we used that (a + b) δ ≤ 2 δ (a δ + b δ ) ∀a, b, δ > 0. To conclude this estimate, it remains to use that there exists another constant C > 0 s.t. ∀t > 0 and ∀κ > 0:

max t 0 s d-1 2 (t -s) d 2 (coth(κs)) d 2 ds, t 0 s d-1 2 (t -s) d 2 (coth(κ(t -s))) d 2 ds ≤ C (1 + κ) d 2 κ d 2 (1 + t) d 2 t d+1 2 . (2.50)
To derive (2.9), we have to modify the upper bound in (2.42) by mimicking the method used above to make appear the singularity (sinh(κt))

∂S ∂t (t)φ + HS(t)φ = R(t)φ.
Then the following two statements are true: (i). The sequence of bounded operators (n > [1/t]):

T n (t) := t-1 n 1 n exp[-(t -s)H]R(s) ds,
converges in norm; let T (t) be its limit;

(ii). The following equality takes place on B(H ): exp(-tH) = S(t) -T (t).

Before giving the proof, we need a series of estimates related to the kernel of the semigroup generated by the operator in (2.28). The proof of the below lemma can be found in Sec. 2.2.4. Lemma 2.11. ∀d ∈ {1, 2, 3} there exists a constant

C d > 0 s.t. ∀0 < η < 1, ∀L ∈ [L 0 (η), ∞), ∀κ > 0, ∀(x, y) ∈ Λ 2d
L and ∀t > 0:

g (d) L,κ,η (x, y; t) ≤ C d e -κ 2 4 L 2 4 t G (d) ∞,0 (x, y; t), (2.51)
∇ x g (d) L,κ,η (x, y; t) ≤ C d (1 + t) d √ t e -κ 2 8 L 2 4 t G (d) ∞,0 (x, y; 2t), (2.52)
∆ x g (d) L,κ,η (x, y; t) ≤ C d (1 + t) 2d t e -κ 2 16 L 2 4 t G (d) ∞,0 (x, y; 2t). (2.53)
Proof of Proposition 2.6. The only thing we have to do is verify the assumptions of Proposition 2.10 in which G L,κ,η (t) plays the role of S(t). Let 0 < η < 1, L ∈ [L 0 (η), ∞) and κ > 0 kept fixed. (A) From (2.33), G L,κ,η (t) is uniformly bounded in t by some constant C d > 0. (B) By using that s -lim t↓0 G ∞,κ (t) = ½ and s -lim t↓0 g L,κ,η (t) = ½ in the kernels sense, then:

∀φ ∈ L 2 (Λ d L ), lim t↓0 G L,κ,η φ = fL,η f L,η + fL,η f c L,η φ = f L,η + f c L,η φ = φ,
where we used (2.24) and (2.26). Next, let us investigate the strong differentiability. From (2.31):

∀φ ∈ L 2 (Λ d L ), 1 δt G (p) L,κ,η (t + δt)φ (• ) -G (p) L,κ,η (t)φ (• ) = fL,η (• ) 1 δt R d R d G (d) ∞,κ (• , z; t)G (d) ∞,κ (z, y; δt) dz -G (d) ∞,κ (• , y; t) f L,η (y)φ(y) dy. Since G ∞,κ (t)L 2 (R d ) → D(H ∞,κ
), then the Stone theorem (in the kernels sense) provides:

lim δt↓0 1 δt G (p) L,κ,η (t + δt)φ (• ) -G (p) L,κ,η (t)φ (• ) = -fL,η (• )H ∞,κ R d G (d) ∞,κ (• , y; t)f L,η (y)φ(y) dy.
By using similar arguments to treat the contribution coming from G

L,κ,η (• ), we therefore obtain:

lim δt↓0 1 δt {G L,κ,η (t + δt)φ -G L,κ,η (t)φ} = -fL,η H ∞,κ G ∞,κ (t)f L,η φ -fL,η h L,κ,η g L,κ,η (t)f c L,η φ. (2.54) (C) Let D 0 := {φ ∈ C 1 (Λ d L ) ∩ C 2 (Λ d L ), φ| ∂Λ d L = 0, ∆φ ∈ L 2 (Λ d L
)} be the domain on which H L,κ is essentially self-adjoint. In the weak sense for any ϕ ∈ D 0 , ψ ∈ C ∞ 0 (Λ d L ) and t > 0:

l L (ϕ, ψ) := H L,κ ϕ, G L,κ,η (t)ψ L 2 (Λ d L ) = -ϕ, ∂G L,κ,η ∂t (t)ψ L 2 (Λ d L ) + ϕ, W L,κ,η (t)ψ L 2 (Λ d L ) ,
where ∂GL,κ,η ∂t (t) denotes the operator in the r.h.s. of (2.54). Note that the second equality is obtained by performing some integration by parts, and afterwards by using the following identities:

H L,κ fL,η = H ∞,κ fL,η = H ∞,κ , fL,η + fL,η H ∞,κ ,
as well as (recall that the potential V L in (2.28) satisfies V L (x) = |x| 2 on Supp( fL,η )):

H L,κ fL,η = h L,κ,η fL,η = h L,κ,η , fL,η + fL,η h L,κ,η . Since l L (ϕ, • ) is a bounded linear functional ∀ϕ ∈ D 0 , then C ∞ 0 (Λ d L ) ∋ ψ → l L (ϕ, ψ
) can be extended in a linear and bounded functional on L 2 (Λ d L ) by the B.L.T. theorem. As well, since l L (• , ψ) is a bounded linear functional ∀ψ ∈ L 2 (Λ d L ) then ϕ → l L (ϕ, ψ) can be extended on the self-adjointness domain D(H L,κ ). This means that ∀t > 0, Ran(G L,κ,η (t)) ⊂ D(H L,κ ). Hence:

ϕ, H L,κ G L,κ,η (t)ψ L 2 (Λ d L ) = -ϕ, ∂G L,κ,η ∂t (t)ψ L 2 (Λ d L ) + ϕ, W L,κ,η (t)ψ L 2 (Λ d L ) .
Finally, from (2.34)

W L,κ,η (t) ≤ Ct -1 2 ∀0 < t ≤ 1. Hence W L,κ,η (t) is integrable in t ∼ 0.

Proof of intermediary results.

Proof of Lemma 2.11. (2.51) follows from the Feynman-Kac formula in [14, Thm. X.68] together with (4.5) and the definition of the

L 0 in (2.23) leading to (L -L η ) 2 ≥ L 2 /2 ∀L ∈ [L 0 (η), ∞).
Next, let us turn to the proof of (2.52)-(2.53). To do that, let us introduce an operator of reference.

∀d ∈ {1, 2, 3}, ∀0 < η < 1, ∀L ∈ (0, ∞) and ∀κ > 0, define on C ∞ 0 (Λ d L ): hL,κ,η := 1 2 (-i∇ x ) 2 + κ 2 2 ṼL,η (x), ṼL,η (x) := 1 4 (L -L η ) 2 .
(2.55) By standard arguments, (2.55) extends to a family of self-adjoint and semi-bounded operators ∀L ∈ (0, ∞), denoted again by hL,κ,η . For any 0 < η < 1, L ∈ (0, ∞) and κ > 0, let {g L,κ,η (t) := e -t hL,κ,η :

L 2 (Λ d L ) → L 2 (Λ d L )
} t≥0 be the strongly-continuous one-parameter semigroup generated by hL,κ,η . Its integral kernel denoted by g(d) L,κ,η is explicitly known and reads as:

∀(x, y) ∈ Λ 2d L , ∀t > 0, g(d) L,κ,η (x, y; t) = e -κ 2 8 (L-L η ) 2 t G (d) L,0 (x, y; t), (2.56)
where

G (d)
L,0 is the kernel of the semigroup generated by the Dirichlet Laplacian in L 2 (Λ d L ), see (4.6). Note that (2.56) directly follows from the Feynman-Kac formula. The starting-point of the proof of (2.52)-(2.53) is a Duhamel-like formula to express the semigroup {g L,κ,η (t)} t>0 in terms of {g L,κ,η (t)} t>0 whose integral kernel is given in (2.56). Let 0 < η < 1, L ∈ [L 0 (η), ∞) (see (2.23)) and κ > 0 be fixed. In the bounded operators sense on L 2 (Λ d L ), it takes place:

∀t > 0, g L,κ,η (t) = gL,κ,η (t) - t 0 gL,κ,η (s) h L,κ,η -hL,κ,η g L,κ,η (t -s) ds, (2.57)
where we used the self-adjointness of the semigroups {g L,κ,η (t)} t≥0 , {g L,κ,η (t)} t≥0 . Proof of (2.52). From (2.57), it follows in the kernels sense:

∀(x, y) ∈ Λ 2d L , ∀t > 0, ∇ x g (d) L,κ,η (x, y; t) = ∇ x g(d) L,κ,η (x, y; t) - 1 2 q (d)
L,κ,η (x, y; t), (2.58)

q (d) L,κ,η (x, y; t) := κ 2 t 0 Λ d L ∇ x g(d) L,κ,η (x, z; s) V L,η (z) -ṼL,η (z) g (d)
L,κ,η (z, y; t -s) dz ds.

Recall that V L,η (z) -ṼL,η (z) = |z| 2 -1 4 (L -L η ) 2 on Supp( fL,η ), 0 otherwise. Let us estimate the first kernel in the r.h.s. of (2.58). From (4.11)

∀t > 0, (-i∇) g L,κ,η (t) = (-i∇) gL,κ,η (t) - t 0 gL,κ,η (s) (-i∇) h L,κ,η -hL,κ,η g L,κ,η (t -s) ds.
It follows in the kernels sense:

∀(x, y) ∈ Λ 2d L , ∀t > 0, ∆ x g (d)
L,κ,η (x, y; t) = ∆ x g(d) L,κ,η (x, y; t) -

1 2 2 l=1 u (d),l L,κ,η (x, y; t), u (d),1 L,κ,η (x, y; t) := κ 2 t 0 Λ d L ∇ x g(d) L,κ,η (x, z; s) (∇ z V L,η ) (z)g (d)
L,κ,η (z, y; t -s) dz ds,

u (d),2 L,κ,η (x, y; t) := κ 2 t 0 Λ d L ∇ x g(d) L,κ,η (x, z; s) ṼL,η (z) -V L,η (z) ∇ z g (d)
L,κ,η (z, y; t -s) dz ds.

From (4.12) and (2.56), there exists a constant C d > 0 s.t.

∀(x, y) ∈ Λ 2d L , ∀t > 0, ∆ x g(d) L,κ,η (x, y; t) ≤ C d (1 + t) d t e -κ 2 4 L 2 4 t G (d)
∞,0 (x, y; 2t).

Subsequently, by mimicking the method leading to (2.60), there exists another C d > 0 s.t.

∀(x, y) ∈ Λ 2d L , ∀t > 0, u (d),1 L,κ,η (x, y; t) ≤ C d (1 + t) d √ t e -κ 2 8 L 2 4 t G (d)
∞,0 (x, y; 2t).

By the same method again but replacing the estimate (2.51) with (2.52), we have:

u (d),2 L,κ,η (x, y; t) ≤ C d (1 + t) 2d t e -κ 2 16 L 2 4 t G (d) ∞,0 (x, y; 2t) t 0 ds √ s √ t -s .
Gathering the three above estimates together, then the proof of (2.53) is over.

Proof of Lemma 2.8. Let d ∈ {1, 2, 3}, 0 < η < 1, L ∈ [L 0 (η), ∞) and κ > 0 kept fixed. (i). From (2.36) written in the kernels sense, then ∀(x, y) ∈ Λ 2d L and ∀t > 0:

∇ x G * L,κ,η (d) (x, y; t) = (∇f L,η ) (x)G (d) ∞,κ (x, y; t) fL,η (y) + f L,η (x)∇ x G (d) ∞,κ (x, y; t) fL,η (y)+ + ∇f c L,η (x)g (d) L,κ,η (x, y; t) fL,η (y) + f c L,η (x)∇ x g (d)
L,κ,η (x, y; t) fL,η (y).

(2.41) is an upper bound for the two first kernels in the above r.h.s. obtained from (4.3)-(4.4) and (4.9). (2.42) is an upper bound for the two last kernels obtained from (2.51) and (2.52).

(ii). From (2.37) written in the kernels sense, then ∀(x, y) ∈ Λ 2d L and ∀t > 0:

∇ x W * L,κ,η (d) (x, y; t) = 4 m=1 Q (d),m
L,κ,η (x, y; t), with:

Q (d),1 L,κ,η (x, y; t) := -i (∇f L,η ) (x) [(-i∇) , G ∞,κ (t)] (x, y) ∇ fL,η (y)+ + (∇f L,η ) (x)∇ x G (d) ∞,κ (x, y; t) ∇ fL,η (y) -(∇f L,η ) (x)G (d) ∞,κ (x, y; t) 1 2 ∆ fL,η (y), (2.61)
Q (d),2 L,κ,η (x, y; t) := -if L,η (x)∇ x [(-i∇) , G ∞,κ (t)] (x, y) ∇ fL,η (y)+ + f L,η (x)∆ x G (d) ∞,κ (x, y; t) ∇ fL,η (y) -f L,η (x)∇ x G (d) ∞,κ (x, y; t) 1 2 ∆ fL,η (y), (2.62)
Q (d),3 L,κ,η (x, y; t) := -i ∇f c L,η (x) [(-i∇) , g L,κ,η (t)] (x, y) ∇ fL,η (y)+ + ∇f c L,η (x)∇ x g (d) L,κ,η (x, y; t) ∇ fL,η (y) -∇f c L,η (x)g (d)
L,κ,η (x, y; t)

1 2 ∆ fL,η (y), (2.63)
Q (d),4 L,κ,η (x, y; t) -if c L,η (x)∇ x [(-i∇) , g L,κ,η (t)] (x, y) ∇ fL,η (y)+ + f c L,η (x)∆ x g (d) L,κ,η (x, y; t) ∇ fL,η (y) -f c L,η (x)∇ x g (d)
L,κ,η (x, y; t)

1 2 ∆ fL,η (y).
(2.64)

Let us first estimate (2.61). In view of (4.9), (2.41) is clearly an upper bound for the two last terms in the r.h.s. of (2.61). For the first term in (2.61), we use (2.38) in the kernels sense. Then, there exists a constant

C d > 0 s.t. ∀(x, y) ∈ Λ 2d L and ∀t > 0: κ 2 2 t 0 R d G (d) ∞,κ (x, z; t -s) ∇ z |z| 2 G (d) ∞,κ (z, y; s) dz ds ≤ C d (κ 2 L + κ 3 2 )tG (d) ∞,κ (x
: i (∇f L,η ) (x) [(-i∇) , G ∞,κ (t)] (x, y) ∇ fL,η (y) ≤ C d L -3η (1 + L 1-η )te -cκL 2η coth( κ 2 t) G (d)
∞,κ (x, y; t, 4). Restricting to 1 > η > 1 4 , and gathering the above estimate with (2.41) together, then there exist two other constants c,

C d > 0 s.t. ∀L ∈ [L 0 (η), ∞), ∀(x, y) ∈ Λ 2d
L and ∀t > 0:

Q (d),1 L,κ,η (x, y; t) ≤ C d 1 + √ κ coth κ 2 t (1 + t)e -cκL 2η coth( κ 2 t) G (d) ∞,κ (
f L,η (x)∇ x G (d) ∞,κ (x, y; t) 1 2 ∆ fL,η (y) + f L,η (x)∆ x G (d) ∞,κ (x, y; t) ∇ fL,η (y) ≤ C d √ κ coth κ 2 t e -cκL 2η coth( κ 2 t) G (d) ∞,κ (x, y; t, 4).
Here, the property (2.25) is essential to remove a coth(κt) in the numerator of (4.10). For the first term of (2.62), we use the same reasoning than the one leading to (2.65) combined with the property (2.25). Thus, there exist two other constants c, C d > 0 s.t. ∀(x, y) ∈ Λ 2d L and ∀t > 0:

κ 2 2 f L,η (x) t 0 R d ∇ x G (d) ∞,κ (x, z; t -s) ∇ z |z| 2 G (d) ∞,κ (z, y; s) ∇ fL,η (y) dz ds ≤ C d 1 + √ κ L -4η (1 + L)(1 + t)e -cκL 2η coth( κ 2 t) G (d) ∞,κ (x, y; t, 8).
Restricting to 1 > η > 1 4 , and gathering the above estimates together, then there exist two other constants

c, C d > 0 s.t. ∀L ∈ [L 0 (η), ∞), ∀(x, y) ∈ Λ 2d
L and ∀t > 0:

Q (d),2 L,κ,η (x, y; t) ≤ C d 1 + √ κ coth κ 2 t (1 + t)e -cκL 2η coth( κ 2 t) G (d) ∞,κ (
κ 2 2 t 0 Λ d L g (d) L,κ,η (x, z; t -s) (∇ z V L,η ) (z)g (d) L,κ,η (z, y; s) dz ds ≤ C d e -κ 2 8 L 2 4 t e -c L 2η t G (d)
∞,0 (x, y; 2t).

We conclude that there exist two other constants c, 4 Appendix-The semigroup: A review of some properties.

C d > 0 s.t. ∀(x, y) ∈ Λ 2d L : ∀t > 0, Q (d),3 L,κ,η (x, y; t) ≤ C d (1 + t) d √ t e -κ 2 8 L 2 4 t e -c L 2η t G (d) ∞,0 (x, y; 4t). ( 2
Here, we collect the technical results we use throughout the paper involving the semigroup generated by H L,κ , see Sec. 1.2.

For simplicity's sake, we hereafter use the notation Λ ∞ := R. From (1.1)-(1.2), recall that:

∀L ∈ (0, ∞], H L,κ = 1 2 (-i∇ x ) 2 + 1 2 κ 2 |x| 2 in L 2 (Λ d L ), d ∈ {1, 2, 3}. (4.1)
Below, we allow the value κ = 0; in that case, H L,0 with L < ∞ is nothing but the Dirichlet Laplacian and H ∞,0 the free Laplacian on the whole space whose self-adjointness domain is W 2,2 (Λ d ∞ ). Recall some properties on the strongly continuous one-parameter semigroup {G L,κ (t) := e -tHL,κ : L If κ > 0, the one-dimensional kernel is given by the so-called Mehler formula, see [12, pp. 176]: We mention that, if κ = 0, the kernel G

∀(x, y) ∈ Λ 2 ∞ , ∀t > 0, G (
∇ x G (d

1. 1

 1 Introduction.

2 e -κ 8 L 2 4 16 L 2 4

 282162 ) sinh(κ(t -s))} d tanh( κ 2 (t-s)) ds, for another constant C > 0. Here, we artificially made appear a (sinh(κs)) d 2 under the integration w.r.t. s. This leads to the appearance of a e d 2 κs in the numerator. If κ ≥ 1, we can get rid of it via the term e -κ 2 s (for L large enough) since κ ≤ κ 2 . If 0 < κ < 1, one has to choose L large enough accordingly to κ, i.e. L ≥ cste/ √ κ. Given a κ 0 > 0, let L = L κ0 ≥ L 0 ( 1 2 ) s.t. ∀L ≥ L κ0 , the inequality e -

  . was partially supported by the Lundbeck Foundation, and the European Research Council under the European Community's Seventh Framework Program (FP7/2007-2013)/ERC grant agreement 202859. A part of this work was done while the second author was visiting DIAS-STP (Dublin), B.S. is grateful for invitation and financial support. Both authors warmly thank Horia Cornean, Tony Dorlas and Cyril Levy for helpful and stimulating discussions.

  2 (Λ d L ) → L 2 (Λ d L )} t≥0 generated by H L,κ in (4.1). We refer to [14, Sec. X.8] and[START_REF] Simon | Schrödinger semigroups[END_REF] Sec. B]. As already mentioned, ∀κ ≥ 0 and ∀L ∈ (0, ∞] it is a self-adjoint and positive operator on L 2 (Λ d L ) by the spectral theorem and the functional calculus. Moreover, since {G L,κ (t)} t>0 is bounded fromL 2 (Λ d L ) to L ∞ (Λ d L ), then it is an integral operator by the Dunford-Gelfand-Pettis theorem. Let us turn to the integral kernel of {G L,κ (t)} t>0 we denote by G(d) L,κ . ∀κ ≥ 0 and ∀L ∈ (0, ∞], G (d) L,κ is jointly continuous in (x, y, t) ∈ Λ d L × Λ d L × (0, ∞) and vanishes if x ∈ ∂Λ d L or y ∈ ∂Λ d L . When L = ∞, it is explicitly known. If κ = 0,it is the so-called heat kernel reading for d = 1 as:∀(x, y) ∈ Λ 2 ∞ , ∀t > 0, G

G

  d=1) ∞,κ (x, y; t) = κ 2π sinh(κt) e -κ 4 [(x+y) 2 tanh( κ 2 t)+(x-y) 2 coth( κ 2 t)] . (4.3) Note that the multidimensional kernel (i.e. d = 2, 3) is directly obtained from (4.2) or (4.3) by:∀κ ≥ 0, G (d) ∞,κ (x, y; t) := d j=1 (d=1) ∞,κ (x j , y j ; t), x := {x j } d j=1 , y := {y j } d j=1 .(4.4)When restricting to L ∈ (0, ∞), the mapping L → G(d)L,κ (x, y; t) is positive and monotone increasing. This leads to the following pointwise inequality which holds ∀κ ≥ 0 and ∀L ∈ (0, ∞):∀(x, y, t) ∈ Λ d L × Λ d L × (0, ∞), G(d)L,κ (x, y; t) ≤ sup L>0 G (d) L,κ (x, y; t) = G (d) ∞,κ (x, y; t). (4.5)

  is explicitly known and reads as, see[8, Eq. (4.13)]:∀(x, y) ∈ Λ 2d L , ∀t > 0, G y + 2mL) 2 2t -exp -(x + y -2mL -L) 2 2t .In view of (4.3)-(4.4), let us introduce ∀κ > 0 the new notation:∀γ > 0, G (d) ∞,κ (x, y; t, γ) := κ 2π sinh(κt) d 2 d j=1 e -κ 4γ [(xj+yj) 2 tanh( κ 2 t)+(xj-yj ) 2 coth( κ 2 t)] ,(4.7)with the convention:G (d) ∞,κ (• , • ; t) = G (d) ∞,κ(• , • ; t, 1). Here are collected all the needed estimates: Lemma 4.1. ∀d ∈ {1, 2, 3}, there exists a constant C d > 0 s.t. (i). ∀κ > 0, ∀γ > 0, ∀(x, y) ∈ Λ 2d ∞ and ∀t > 0: G (d) ∞,κ (x, y; t, γ) (x, y; γt) ≤ (2πt) -d 2 , (4.8)

  Lemma 2.8. ∀d ∈ {1, 2, 3} there exist two constants c, C d > 0 s.t.: (

d) (x, z; s)G (d) L,κ (z, y; t -s) dz ds. (2.40)

Next, we need the following lemma whose proof can be found in Sec. 2.2.4:

  Next, we extend the integration w.r.t. z to R d and we use the first upper bound in (4.8) followed by(4.15). Then, we introduce a factor s under the integral w.r.t s. Next, we successively use the lower bound in (2.13) and the upper bound cosh(α) ≤ e α ∀α ≥ 0 leading to (κs) -d

			d-1 2 s -d-1 2
	(coth(κs))	d 2 ≤ e	2 ≤ 2 . On this way, we get under the same conditions than (2.48): 2 κs (sinh(κs)) -d d

, y; t -s, 2) dz.

(2.48)

  and (2.56), there exists a constantC d > 0 s.t.where we used in the last inequality(4.15). Finally use (4.13) to get rid of the L 2 what leads to:

	q	(d) L,κ,η (x, y; t) ≤ C d	(1 + t) d √ t	e -κ 2 8	L 2 4 t G (d) ∞,0 (x, y; 2t),	(2.60)
	for another constant C d > 0. It remains to gather (2.59) and (2.60) together.
	Proof of (2.53). Starting from the below identity which holds in the bounded operators sense:
	∀t > 0, [(-i∇) , gL,κ,η (t)] = -	0	t	gL,κ,η (t -s) (-i∇) , hL,κ,η gL,κ,η (s) ds,
	then by using that [(-i∇), hL,κ,η ] = 0, one gets from (2.57) on L 2 (Λ d L ):
	∀(x, y) ∈ Λ 2d L , ∀t > 0, ∇ x	g(d) L,κ,η (x, y; t) ≤ C d	(1 + t) d √ t	e -κ 2 4	L 2 4 t G (d) ∞,0 (x, y; 2t).	(2.59)
	Subsequently, from (2.59) along with (2.51), there exists another constant C d > 0 s.t.
	∀(x, y) ∈ Λ 2d L , ∀t > 0, q	(d) L,κ,η (x, y; t) ≤ C d κ 2 L 2 (1 + t) d e -κ 2 4	L 2 4 t G (d) ∞,0 (x, y; 2t)	0	t	ds √ s	,

  The estimate in (2.45) follows by adding (2.68) and (2.69) together, then by taking into account the support of the cutoff functions introduced in Sec. 2.2.1.

						.68)
	Turning to (2.64), one can prove that there exist two other constants c, C d > 0 s.t. on Λ 2d L :
	∀t > 0, Q	(d),4 L,κ,η (x, y; t) ≤ C d	(1 + t) 2d √ t	e -κ 2 16	L 2 4 t e -c L 2η t G (d) ∞,0 (x, y; 4t).	(2.69)
	Here, we used (2.27) combined with (4.13) to get rid of a	√ t in the denominator of (2.53).

Proof. From the lower bound sinh(α) ≥ α ∀α ≥ 0 and the one in (2.13), (4.2) is an upper bound for (4.3). Hence (4.8). (4.9)-(4.10) are obtained by direct calculations. The main ingredients are: ∀µ, ν > 0, ∀x ≥ 0, x µ e -νx ≤ 2µ eν µ e -ν 2 x , (4. [START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF] and the following identity: We continue with the following lemma expressing the semigroup property in the kernels sense:

Proof. (i) follows from the semigroup property which reads as:

The proof of (ii) is based on the following explicit calculation: Then, set a 0 := tanh( κ 2 δu), b 0 := coth( κ 2 δu), c 0 := tanh( κ 2 δ(t -u)) and d 0 := coth( κ 2 δ(t -u)). From the identity in (2.18), the following one: tanh(αs) + tanh(α(t -s)) = sinh(αt) cosh(αs) cosh(α(t -s))

, ∀α ≥ 0, ∀t > s > 0, followed by (2.3), one gets:

The rest of the proof consists in using some identities involving the hyperbolic functions to simplify each one of the factor inside the exponentials in the r.h.s. of (4.17). It is (quite) easy to get:

Now, we give some estimates on the operator and trace norms of the semigroup {G L,κ (t)} t>0 . For any κ ≥ 0 and L ∈ (0, ∞], {G L,κ (t)} t>0 is a contraction semigroup, see e.g. [START_REF] Hille | Functional Analysis and Semi-Groups, Revised and Expanded Edition[END_REF]: Lemma 4.3. ∀d ∈ {1, 2, 3}, ∀κ ≥ 0 and ∀t > 0:

Proof. The first inequality follows from the fact that the semigroup {G L,κ (t)} t≥0 is increasing in L in the sense of [START_REF] Broderix | Continuity properties of Schrödinger semigroups with magnetic fields[END_REF]Eq. (2.39)]. The Shur-Holmgren criterion provides the estimate on the operator norms. When κ > 0, we used (4.17) (with c = 0 = d) along with (4.14).

Restricting to κ > 0, ∀L ∈ (0, ∞] {G L,κ (t)} t>0 is a Gibbs semigroup (i.e. trace class, see [START_REF] Angelescu | On the perturbation of Gibbs semigroups[END_REF]):

Lemma 4.4. ∀d ∈ {1, 2, 3}, ∀κ > 0 and ∀L ∈ (0, ∞], {G L,κ (t)} t>0 is a trace class operator on

∞,κ = d κ 2 , one has for any L ∈ (0, ∞):

be the Banach space of Hilbert-Schmidt and trace class operators on L 2 (Λ d L ) respectively. We start with d = 1. Let κ > 0 and t > 0 be fixed. In view of (4.3), from (4.17) (we set c = 0 = d):

∞,κ (x, y; t)

∞ , from [8, Prop. 9] it follows that:

where we used the identity (2.3). By positivity of G

}. The rest of the proof leans on the estimate (4.5) which leads to G L,κ (t) 2 I2 ≤ G ∞,κ (t) 2 I2 . Hence, ∀L ∈ (0, ∞) G L,κ (t) is also a trace class operator on L 2 (Λ 1 L ), and by mimicking the above arguments, its trace norm obeys G L,κ (t) I1 = Tr L 2 (Λ 1 L ) {G L,κ (t)} ≤ G ∞,κ (t) I1 . The case of d = 1 is done. The generalization to d = 2, 3 is straightforward due to (4.4).