William Pophillat

Peter Bayer

Esther Teyssier

Philipp Blum

Cerema Guillaume Attard
email: guillaume.attard@cerema.fr

Impact of groundwater heat pump systems on subsurface temperature under variable advection, conduction and dispersion

Keywords: pump, open-loop, numerical modelling, aquifer. Highlights

The number of worldwide installed groundwater heat pump (GWHP) systems is approaching one million units. This increase in GWHP systems creates an increasing competition for groundwater use. Thus, we have simulated the long-term thermal impacts of such open-loop geothermal systems with numerical heat transport models, by studying the role of flow and heat transport parameters. The focus is set on transient conditions evolving over decades of seasonally imbalanced heat injections in aquifers with differing ambient flow regimes. The results demonstrate the existence of plume length peaking in moderate groundwater flow velocities (0.5 -1.0 m d -1 ). The longitudinal and transverse dispersivity coefficients have a strong influence on the plume extension and on transient development.

Reducing the injected temperature difference allows the control of the maximum temperature reached in the plume but implies a considerable hydraulic influence. Comparisons with numerical simulations, using an average of the thermal load, shows that the flow velocity and dispersity coefficients strongly condition the relevance of such simplified approaches. In addition, if satisfactory results are obtained for seasonal or long-term estimations, those for intermediate lengths of time are insufficient. This study contributes to improving the overall understanding of key elements involved in the thermal impact and the development of GWHP systems.

Introduction

Groundwater heat pumps (GWHP, Figure 1a) represent a common technology in the field of lowenthalpy geothermal energy utilization. They rely on the presence of groundwater at shallow depths, which is pumped from one or more wells to extract energy for space heating, air conditioning, industrial heat supply, as well as other applications [START_REF] Rafferty | Design aspects of commercial open-loop heat pump systems[END_REF]. The used cooled or warmed water is reinjected into a separate aquifer and ideally, does not interact with the production well. Often, groundwater is also utilized for cooling purposes and therefore water warmer than the originally sourced water is reinjected into the ground. According to Hähnlein et al. [START_REF] Hähnlein | Sustainability and policy for the thermal use of shallow geothermal energy[END_REF] these systems may be further distinguished based on the mode of use. A special category are those applications with the purpose of temporal heat storage such as aquifer thermal storage (ATES) systems [START_REF] Fleuchaus | Worldwide application of aquifer thermal energy storage-A review[END_REF][START_REF] Bloemendal | ATES systems in aquifers with high ambient groundwater flow velocity[END_REF][START_REF] Bloemendal | Analysis of the impact of storage conditions on the thermal recovery efficiency of low-temperature ATES systems[END_REF]. In the present study however, we focus on GWHP systems that do not consider heat storage. This is mostly the case when groundwater flow is significant and local heat storage is not feasible.

Utilization of groundwater as a heat carrier characterises GWHP systems as open-loop. However, most heat pump-based geothermal systems, are closed-loop systems with vertical borehole heat exchangers [START_REF]Transparency-Market-Research[END_REF][START_REF] Bayer | Greenhouse gas emission savings of ground source heat pump systems in Europe: a review[END_REF] (Figure 1b). In urban environments with productive aquifers, the open-loop GWHPs are growing in numbers. Around every seventh geothermal heat pump is estimated to be an open-loop system [START_REF]Transparency-Market-Research[END_REF]. In combination with reported global installed capacity of all geothermal heat pumps in 2015 [START_REF] Lund | Direct utilization of geothermal energy 2015 worldwide review[END_REF], and assuming continued marked growth rates of 10%, the equivalent number of currently installed 12 kW GWHP units has reached more than 800,000 worldwide. As these installations are concentrated primarily in the major cities of countries including USA, China, Sweden, Switzerland, Germany and France, the management of the concurrent use of shallow groundwater resources for various energy supply applications becomes increasingly essential [START_REF] Hähnlein | Sustainability and policy for the thermal use of shallow geothermal energy[END_REF][START_REF] Ferguson | Unfinished business in geothermal energy[END_REF][START_REF] Rivera | Increased Temperature In Urban Ground As Source Of Sustainable Energy[END_REF]. Regulating multiple adjacent systems, each with slowly evolving thermal effects on the subsurface, is a challenging task. The elementary tools for predicting and managing GWHPs are models that simulate the physical processes in the ground. In recent years, European cities in particular have refined their management strategies of thermal groundwater use through the application of numerical groundwater flow and transport models. Prominent examples of their use are in London [START_REF] Fry | Lessons from London: regulation of open-loop ground source heat pumps in central London[END_REF][START_REF] Arthur | Modelling large ground source cooling systems in the Chalk aquifer of central London[END_REF], Zürich, Basel in Switzerland [START_REF] Epting | The thermal impact of subsurface building structures on urban groundwater resources-A paradigmatic example[END_REF][START_REF] Mueller | Combining monitoring and modelling tools as a basis for city-scale concepts for a sustainable thermal management of urban groundwater resources[END_REF], Zaragoza in Spain [START_REF] Epting | Development of concepts for the management of thermal resources in urban areas-Assessment of transferability from the Basel (Switzerland) and Zaragoza (Spain) case studies[END_REF][START_REF] Maya | An upscaling procedure for the optimal implementation of openloop geothermal energy systems into hydrogeological models[END_REF], and Turin in Italy [START_REF] Sciacovelli | Multi-scale modeling of the environmental impact and energy performance of open-loop groundwater heat pumps in urban areas[END_REF][START_REF] Bucci | Shallow groundwater temperature in the Turin area (NW Italy): vertical distribution and anthropogenic effects[END_REF][START_REF] Russo | Groundwater Heat Pump (GWHP) system modeling and Thermal Affected Zone (TAZ) prediction reliability: Influence of temporal variations in flow discharge and injection temperature[END_REF] among others. Moreover, Yu et al. [START_REF] Yu | Numerical study based on one-year monitoring data of groundwater-source heat pumps primarily for heating: a case in Tangshan, China[END_REF] reported numerical GWHP simulations for the city of Tangshan, Liang et al. [START_REF] Liang | Modeling and performance evaluation of shallow ground water heat pumps in Beijing plain, China[END_REF] for the Beijing plain and Zhou and Zhou [START_REF] Zhou | Simulation of thermal transport in aquifer: a GWHP system in Chengdu, China[END_REF] for Chengdu, China. Rodríguez et al. [START_REF] Rodríguez | District heating by Groundwater Heat Pump in a social integration housing project located in Puente Alto[END_REF] inspected GWHP systems for district heating in Puente Alto, Chile. For decades, numerical codes have been used for examination of GWHP thermal impacts in US American cities (e.g. [START_REF] Warner | Thermal Impact of Residential Ground-Water Heat Pumps[END_REF]).

However, often data and technical capacity are inadequate for a detailed simulation. Therefore, models are often calibrated with a limited amount of data and are based on previous experience with similar systems already in use. In such cases, simplification of numerical and analytical models were proposed for a more straightforward analysis of potential thermal breakthrough between well doublets consisting of extraction and injection wells [START_REF] Lippmann | Ground-Water Use for Cooling: Associated Aquifer Temperature Changes[END_REF][START_REF] Gringarten | A theoretical study of heat extraction from aquifers with uniform regional flow[END_REF][START_REF] Casasso | Modelling thermal recycling occurring in groundwater heat pumps (GWHPs)[END_REF][START_REF] Piga | Thermal Impact Assessment of Groundwater Heat Pumps (GWHPs): Rigorous vs. Simplified Models[END_REF][START_REF] Milnes | Assessing the impact of thermal feedback and recycling in open-loop groundwater heat pump (GWHP) systems: a complementary design tool[END_REF][START_REF] Banks | Thermogeological assessment of open-loop well-doublet schemes: a review and synthesis of analytical approaches[END_REF][START_REF] Gandy | Predictive modelling of groundwater abstraction and artificial recharge of cooling water[END_REF].

From the perspective of regulators, special attention needs to be given to the potential interaction between neighbouring installations, water protection areas [START_REF] Haehnlein | International legal status of the use of shallow geothermal energy[END_REF][START_REF] Park | Minimum Well Separation for small groundwater heat pump (GWHP) systems in Korea: preliminary analysis based on regional aquifer properties[END_REF][START_REF] Pophillat | Analytical solutions for predicting thermal plumes of groundwater heat pump systems[END_REF] as well as groundwater ecosystems [START_REF] Brielmann | Shallow geothermal energy usage and its potential impacts on groundwater ecosystems[END_REF]. Some related studies have presented sensitivity and scenario analysis of crucial factors for thermal plume development through the employment of analytical and numerical models. For example, Lo Russo et al. [START_REF] Russo | Development of the thermally affected zone (TAZ) around a groundwater heat pump (GWHP) system: A sensitivity analysis[END_REF] used numerical simulations to carry out a sensitivity analysis by varying hydrogeological and thermal parameters gradually from -20 % to +20 %, with respect to their real values at a GWHP site. They compared the results from the various scenarios after one year of constant injection. Their results showed that the parameters which have the most significant influence are those related to the advective component of the ground heat flow. However, Piga et al. [START_REF] Piga | Thermal Impact Assessment of Groundwater Heat Pumps (GWHPs): Rigorous vs. Simplified Models[END_REF] employed numerical sensitivity analyses for a more extensive range of parameter values by considering variable injections and longer simulation times of fifty years. They demonstrated that among the hydraulic and thermal parameters, especially the Darcy velocity, thermal conductivity and dispersity, have a significant influence on the plume extension. For a case study, Lo Russo et al. [START_REF] Russo | Groundwater Heat Pump (GWHP) system modeling and Thermal Affected Zone (TAZ) prediction reliability: Influence of temporal variations in flow discharge and injection temperature[END_REF] investigated the role of the temporal discretization of the heat injection for the numerical estimation of the thermal impact. They highlighted the need to consider a temporal resolution that represents a temporal variation of over a year (i.e. hourly to average monthly equivalent) of the injection rate to satisfactorily estimate thermal plume evolution in the aquifer. Piga et al. [START_REF] Piga | Thermal Impact Assessment of Groundwater Heat Pumps (GWHPs): Rigorous vs. Simplified Models[END_REF] also evaluated the ability of numerical models to estimate the thermal impact between well doublets. They concluded that yearly average values can only be used to estimate the plumes in the case of heating or cooling-dominated loads. For more balanced thermal loads, using the yearly averaged load heavily overestimates the plume length and underestimates its width.

The findings from previous studies highlight the strong influence of dispersivity. However, the few studies that focus on the long-term thermal impact consider a fixed ratio between longitudinal and transverse dispersivity coefficients. Therefore, the separate influence on the plume evolution is hardly understandable. The influence of parameters that define the heat injection is also rarely studied.

Furthermore, the plume extension assessment is mostly based on the use of a fixed temperature difference threshold, for instance, 1 K, along with the potential threshold effect that must be considered.

The present work uses a sensitivity analysis, based on two-dimensional (2D) transient numerical simulations. It aims to deepen the understanding of the influence of the groundwater velocity, the longitudinal and transverse dispersivity coefficients and the parameters used to define the energy load on the long-term thermal impact due to a variable injection. We propose an approach based on the maximal plume extension reached during a cycle of injection, and the influence of the threshold chosen to define the plume boundary is discussed. Also, we specify the range of applicability of simplified averaging of the thermal load depending on the studied parameters. This simplified averaging method contributes to the overarching objective of this study to support authorities in the regulation of open GWHP systems in crowded areas, with the help of suitable simplified models.

Materials and methods

The methodology applied in this work is based on simulation of transient plumes evolving from injection wells. By systematic variation of parameter settings in different theoretical scenarios, the dominant factors for the thermal impact of seasonally operating GWHP systems are revealed. In the following passage the model scenarios are introduced, and the evaluation procedure to compare different model predictions is presented. Additionally, the specifications of the numerical models are described.

Studied scenarios

A porous, homogeneous and confined aquifer with adiabatic upper and lower boundaries is inspected. This aquifer has a porosity of n = 30 %, and hydraulic conductivity of K = 6.94 ×10 -4 m s -1 , simulating the typical characteristics of a coarse, sand aquifer [START_REF] Bayer | Three-dimensional high resolution fluvioglacial aquifer analog: Part 1: Field study[END_REF]. The thermal properties are considered homogeneous with a bulk thermal conductivity of the water-saturated aquifer material of λm = 2.24 W m -1 K -1 (water: 0.578 W m -1 K -1 , solid: 4.0 W m -1 K -1 ) and a composite heat capacity of 2.888 MJ m - 1 K -1 (water: 4.185 MJ m -1 K -1 , solid: 2.332 MJ m -1 K -1 ) [START_REF] Woodside | Thermal conductivity of porous media. I. Unconsolidated sands[END_REF]. The undisturbed temperature is Tu = 12 °C (285.15 K). In different scenarios, further settings are varied based on a reference case with an undisturbed uniform hydraulic gradient of 0.5% which is considered a moderate value for an urban aquifer [START_REF] Attard | Impact of underground structures on the flow of urban groundwater[END_REF]. The resulting groundwater flow velocity (or seepage velocity) is va = 1 m d -1 (which corresponds to a specific discharge of q = va • n = 0.33 m d -1 ). The default longitudinal and transverse dispersivity of this reference case are set αl = 5 m and αt = 0.5 m, respectively [START_REF] Stauffer | Thermal Use of Shallow Groundwater[END_REF].

In our theoretical study, the implemented GWHP system is only used for cooling during the summer season and is operated for 20 years. The thermal impact from a similar system which was used for heating would be equivalent, but with a cold plume instead of the heat plume simulated here. The assigned time-dependent function for injection is oriented at a realistic air-conditioning demand for a tertiary building located in the area of Lyon (France). In line with the results by Lo Russo et al. [START_REF] Russo | Groundwater Heat Pump (GWHP) system modeling and Thermal Affected Zone (TAZ) prediction reliability: Influence of temporal variations in flow discharge and injection temperature[END_REF] on appropriate temporal discretization, the injection is described by using a modulation function with monthly increments (Figure 2). Given that the analysis is based on a 2D approach, the various parameters describing the injected heat flux are given per unit length of the injection well. The injected heat flux by unit length of the well considered in the reference scenario is therefore constant over a month, and the value corresponds to the monthly mean of the real injected flux (Figure 1a). The corresponding water injection is characterized in the reference scenario by a constant temperature difference of ∆ = 10 K between groundwater and injected water, and a variable injection rate per unit length of the well, . During the injection period, values vary accordingly between 0.01 and 0.12 l s -1 m -1 . Per unit thickness of the aquifer, the annual total volume of water injected is 8172 m 3 and the annual total energy content is 3.42 10 11 J. Based on the reference case, three distinct influences are configured for the sensitivity analysis resulting in 20 different scenarios (Table 1). and va = 10 m d -1 . These velocities are obtained by keeping the hydraulic gradient at 0.5 % and varying the hydraulic conductivity of the aquifer from 6.94 × 10 -5 m s -1 to 6.94 × 10 -3 m s -1 .

These values correspond to materials varying from very fine sand to coarse gravel [START_REF] Spitz | A practical guide to groundwater and solute transport modeling[END_REF]. Note that groundwater flow velocity varies strongly among different field sites depending on various hydrogeological conditions. Even if two orders of magnitudes are studied for va , this range does not cover the full scope of observed velocities in practice.

(ii) Influence of the dispersivity coefficients: By revisiting field data reported in the literature, Zech et al. [START_REF] Zech | Is unique scaling of aquifer macrodispersivity supported by field data?: Is unique scaling supported by field data?[END_REF] showed that the longitudinal dispersivity, αl, and in particular its asymptotic value for an extended observation distance, depends strongly on the structural properties of the aquifer.

Based on their analysed data sample, this asymptotic value remained lower than 10 m.

Furthermore, the survey by Gelhar et al. [START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF] demonstrated that the ratio αl/αt varies for each case between 1 and 100. Because of these observations, ten scenarios are defined to evaluate the influence of longitudinal and transverse dispersivity. Six of these scenarios investigate the specific role of the transverse dispersivity. The value of αl is kept at 5 m as in the reference case, and αt is varied by the factors 0.1, 0.2, 0.5, 2, 5 and 10 (the ratio αl/αt thus varies between 1 and 100). Four scenarios inspect the role of the longitudinal dispersivity, αl. According to Zech et al.

[43], its value is kept below 10 m. The transverse dispersivity αt is set fixed at 0.5 m as in the reference case, and the αl is varied by the factors 0.1, 0.2, 0.5 and 2, i.e. the anisotropy ratio αl/αt thus varies between 1 and 20.

(iii) Influence of parameters describing the heat injection: The injected heat flux per unit length of 176 the injection well is defined by = ∆ , where Cw is the heat capacity of water. The 177 variable heat flux (i.e., heat supply by the GWHP) given in Figure 1a can thus be implemented 178 by varying or ∆ . Four scenarios are defined to scrutinize the influence of and ∆ 179 maintaining the similar heat budget. Two scenarios (∆T-5 and ∆T-2.5) focus on the formulation 180 of a fixed ∆

(5 K and 2.5 K). Accordingly, during the injection period, varies between 181 0.03 and 0.24 l s -1 m -1 and 0.06 and 0.48 l s -1 m -1 , respectively. Two further scenarios (q-0.12 and 182 q-0.24) evaluate the influence of a fixed value of (0.12 l s -1 m -1 and 0.24 l s -1 m -1 ). The 183 corresponding values of ∆ vary between 1.2 K and 9.8 K (with a mean value of 5 K), and 184 between 0.6 K and 4.9 K (with a mean value of 2.5 K). 185

Description Abbreviation qinj (l s -1 m -1 ) ΔTinj (K) K (m s -1 ) αl (m) αt (m)

Reference

Ref.

[0.01 ; 0.12] 10 6.94 × 10 -4 0.5

Influence of the groundwater flow velocity

Reference value × 0. Fixed injection rate ∆ = 5 K q-0.12 0.12 [1.2 ; 9.8] 6.94 × 10 -4 5 0.5 ∆ = 2.5 K q-0.24 0.24 [0.6 ; 4.9] 6.94 × 10 -4 5 0.5

Fixed parameters i (-) n (-) Tu (K) λw (W m -1 K -1 ) λs (W m -1 K -1 ) Cw (J m -3 K -1 ) Cs (J m -3 K -1 )

Numerical simulation

Groundwater flow and heat transport are simulated using the commercial code FEFLOW [START_REF] Diersch | FEFLOW: finite element modeling of flow, mass and heat transport in porous and fractured media[END_REF]. This software is based on the finite element method and was applied for heat transport in related studies [START_REF] Epting | The thermal impact of subsurface building structures on urban groundwater resources-A paradigmatic example[END_REF][START_REF] Russo | Groundwater Heat Pump (GWHP) system modeling and Thermal Affected Zone (TAZ) prediction reliability: Influence of temporal variations in flow discharge and injection temperature[END_REF][START_REF] Attard | Deterministic modeling of the impact of underground structures on urban groundwater temperature[END_REF][START_REF] Wagner | Numerical sensitivity study of thermal response tests[END_REF][START_REF] Hecht-Méndez | Evaluating MT3DMS for Heat Transport Simulation of Closed Geothermal Systems[END_REF][START_REF] Zhu | Groundwater temperature evolution in the subsurface urban heat island of Cologne, Germany[END_REF][START_REF] Dehkordi | Effect of thermal-hydrogeological and borehole heat exchanger properties on performance and impact of vertical closed-loop geothermal heat pump systems[END_REF]. The 2D numerical model domain used for all scenarios is illustrated in Figure 2. The domain covers a rectangle of 3000 m by 300 m. Only the injection well is simulated assuming that the production well is located far away from it and that there is no hydraulic interference among them both. The injection well is located at 150 m from the east and north boundaries. A horizontal mesh is applied using a triangular algorithm. The size of the elements ranges from several centimetres near the injection well to about three meters at the domain boundaries.

Figure 3 shows the boundary conditions (BCs) applied. Concerning the flow problem, Dirichlet BCs (fixed hydraulic head) are assigned to the west and east boundaries, resulting in a forced flow from west to east (Figure 3a). The water injection is simulated by assigning a variable nodal source BC to the well node. Concerning the heat transport problem, Dirichlet BCs (fixed temperature) are set at the east boundary to fix the temperature of the fluid entering through these boundaries. The hot water injection is simulated by implementing a variable heat nodal source BC to the well node and using the divergent form of the heat transport equation [START_REF] Diersch | About the difference between the convective form and the divergence form of the transport equation in FEFLOW White Papers[END_REF]. As we focus on a 2D approach, heat conduction into confining layers of the aquifer is not considered. Consequently, all the results are representative for GWHP systems with a sufficiently long well screen. For each scenario, 20 annual cycles of transient injection are simulated. An annual cycle starts at the beginning of the injection phase (May 1 st ) and lasts 365 days. The hydrogeological and thermal parameters used for the calculations are listed in Table 1.

In order to compare results with a simplifying seasonal average of the energy loads, the models described above are also used with continuous fixed heat injection. The yearly average of the thermal load is used with both a one-year and twenty-year simulation time. We also assess the ability to estimate seasonal impact (longer than five months) by using the seasonal average of the thermal load .

Indicators

In the reference case the injection, and thus the shape of the thermal plume, varies throughout an annual cycle. However, after a sufficiently long time at a given location, the plume evolution during a cycle follows the same pattern each year, i.e. the annual trend of the thermal impact is identical. The time it takes to reach this dynamic stability is hereafter called "stabilization time" (St). This is defined as the minimal number of years after which becomes a periodic (annual) function, meaning that for all points of the studied area and for all > St, we get , , = , , + 1 year .

In our analysis, we focus on the area where the thermal impact exceeds 1 K (1 K-plume), which is commonly used to describe the thermal impact or the so-called thermally affected zone (TAZ) caused by shallow geothermal systems [START_REF] Piga | Thermal Impact Assessment of Groundwater Heat Pumps (GWHPs): Rigorous vs. Simplified Models[END_REF][START_REF] Russo | Development of the thermally affected zone (TAZ) around a groundwater heat pump (GWHP) system: A sensitivity analysis[END_REF][START_REF] Russo | Groundwater Heat Pump (GWHP) system modeling and Thermal Affected Zone (TAZ) prediction reliability: Influence of temporal variations in flow discharge and injection temperature[END_REF][START_REF] Molina-Giraldo | Evaluating the influence of thermal dispersion on temperature plumes from geothermal systems using analytical solutions[END_REF][START_REF] Molina-Giraldo | A moving finite line source model to simulate borehole heat exchangers with groundwater advection[END_REF][START_REF] Alcaraz | Advection and dispersion heat transport mechanisms in the quantification of shallow geothermal resources and associated environmental impacts[END_REF]. Aside from this, two other thresholds (0.5 K and 2 K) are used to discuss the role of this threshold on the results. For time-dependent management and especially for regulations, maximal effects are of special interest [START_REF] Hähnlein | Sustainability and policy for the thermal use of shallow geothermal energy[END_REF]. The maximal temperature field is thus determined by computing the maximum temperature reached at each point of the domain since the beginning of the simulation:

#$% , , & = max ) Є [, ;) . ] , , (1)
Based on this field, the 1 K-max-plume (0.5 K-max-plume and 2 K-max-plume) is defined as the area where the thermal impact exceeds 1 K (0.5 K and 2 K, respectively) at least once since the beginning of the simulation. It expands each new annual cycle and stabilizes when dynamic stability is reached.

Subsequently, we focus on the maximum plume at the end of the simulation ( & ), which characterizes the maximal impact reached over the total period of 20 years. For illustration, Figure 4 shows the stabilized shape of 1 K-plumes during the 20 th annual cycle and the corresponding 1 K-max-plume. The plumes extension is analysed based on focusing on the length and maximal width of the plume downgradient from the injection well. Comparisons between scenarios are made by calculating the relative difference between the maximum plume length, width and stabilization time values according to:

01 2 = 100 • 242 567 2 567 % (2)
where X denotes the considered parameter (length, width or stabilization time) value and XREF its reference value (i.e. of the reference scenario).

Results

Simulation of the reference scenario

The final plume of the reference scenario is visualized in Figure 4b. Ambient groundwater flow moves the thermal plume in the downstream direction. Anisotropic dispersion deforms the plumes and expands them in the longitudinal direction, while together mechanical dispersion and diffusion control continuous spreading. This results in a gradual decrease of the maximum temperature in the downgradient of the well. The 1 K-max-plume length and width are 623 m and 46 m respectively, and it takes six annual cycles to reach the dynamic stability.

However, the thermal effect of the GWHP reaches much further than the 1 K-plume boundary. A crucial determinant for the plume extension is the chosen temperature threshold. To illustrate the role of the threshold, an animation is added in the supplementary material (SI_ref.avi). It illustrates the transient evolution of the 2 K-, 1 K-and 0.5 K-plumes during eight years of injection for reference conditions.

Influence of groundwater flow velocity

When the groundwater flow velocity is low (va = 0.1 m d -1 ), the injection induces a higher local hydraulic gradient around the injection well than for the reference case. It yields a wide 1 K-maxplume characterized by a much higher maximum width (>250 %) and a slightly smaller length. In fact, for flow velocities lower than 0.5 m d -1 , the dynamic stability is not reached after 20 cycles.

Furthermore, the 1 K-max-plumes in Figure 5 delineate the state of this simulation period, but not steady state conditions.

When increasing the flow velocity in the further scenarios, the local hydraulic effect of the injection well reduces and the hydraulic regime in the aquifer is increasingly dominated by the ambient flow regime. Additionally, higher flow velocities induce more mechanical dispersion, both in the longitudinal and in the transversal direction. The plume gets gradually thinner when the flow velocity increases. Also, the 1 K-max-plume length grows and reaches a peak value for a flow velocity Computing the Peclet number for each of these scenarios shows that they are always characterized by an advection dominated heat transport. Results show the existence of a limit velocity value (between 0.5 and 1 m d -1 ) beyond which a sharp decline of the 1 K-max-plume length (Fig. 5b), area and energy content (Fig. 5c) are observed. As illustrated in Figure 5a, the plume area declines at higher flow velocity as a consequence of growing longitudinal spreading. At lower flow velocity, the effect of dispersive spreading diminishes, yielding a maximum of plume length and surface at va = 0.5 m d -1 . As smaller flow velocity and thus advection attenuates downgradient propagation of the plume, plume length and area decrease again. Under such conditions, conduction has a significant effect on the heat transport but obviously is less influential than advection and dispersion. In fact, the energy content within the 1K-max-plume rises when lowering the flow velocity. This reflects that the plume may be deformed by advection, but its energy content is mostly linked to the dispersive spreading of the plume. Increase of ambient flow velocity thus yields a variable plume size depending mainly on advection and dispersion, whereas the energy content is governed by dispersion. Small plumes in case of small velocity thus have relatively high temperatures. The associated enhanced lateral conduction, however, does not substantially influence the energy content of the plume. The reversal points in the curves shown in Figures 5b,c depend on the temperature threshold. The higher the threshold, the smaller the velocity value corresponding to the reversal point (Figure 6).

Influence of dispersivity coefficients

We keep the reference settings and focus on the effect of the transverse and longitudinal dispersivity coefficients. We thus consider the variation of the dispersion in each direction separately. The influence of dispersivity values on 1 K-max plumes is shown in Figure 7. Also, two Figures (S2 and The influence of the longitudinal dispersivity coefficient on the 1K-max-plume extension is shown by Figure 7a andc. In these scenarios, the transverse dispersivity coefficient is kept at αt = 0.5 m. Results

show that the lowest αl considered (αl = 0.5 m) causes a long 1 K-max-plume, which is 99 % longer than for the reference case and has a higher stabilization time. When αl increases, the 1 K-max-plume length gradually decreases. For αl ≥ 2.5 m, the relative influence of longitudinal dispersivity coefficient variation becomes less pronounced (Figure 7c). It also has a reduced influence on the 1 Kmax-plume width and on the stabilization time (Figure 7c).

The influence of the transverse dispersivity coefficient on the 1K-max-plume extension is shown by Figure 7b andd. In these scenarios, the longitudinal dispersivity coefficient is kept at αl = 5 m. The lowest αt yields a thin and long 1 K-max-plume (Figure 7b), which is 165 % longer than for the reference case and has a high stabilization time. When αt increases, the 1 K-max-plume width also increases and its length greatly decreases (Figure 7b). The stabilization time also decreases. The relative influence of αt on plume length, width and stabilization time is significantly more pronounced than for αl. For the considered ranges, the 1 K-max-plume lengths are between 177 m and 1655 m for variable αt, and between 598 m and 1242 m when changing αl, indicating a considerable influence of transverse dispersivity. Complementary information of the influence of the dispersivity coefficients can be obtained by focusing on the fate of the plumes generated by subsequent injection phases (Figure S2). Increasing αl (or αt) (i.e. increasing the plume spreading in longitudinal (or transverse) direction) conducts to accelerate the decrease of the energy concentration into the plume and consequently, to reduce the distance from which its maximal temperature passes below 1 K (Figure S2). Again, the temperature threshold is elementary for the derived ranges. As can be seen from Figure 8a, the higher the threshold, the higher the relative influence of αl is, and by contrast, the lower the relative influence of αt is.

Characteristics of the heat injection

Figure 9 depicts the results for changing the heat injection by modifying the temperature difference ∆ and the injection rate, . When reducing ∆ we increased to keep the same power of injection. The stronger hydraulic effect of injection generates more radial spreading and thus a greater plume width. Therefore, the lower the value of ∆ , the shorter, but the wider the 1K-max-plume is, especially close to the injection well. For instance, for a fixed ∆ of 2.5 K the 1K max plume is almost 30% shorter and 50% wider than in the reference scenario.

Using a fixed of 0.12 l s -1 m -1 implies setting a ∆ that varies between 1.2 K and 9.8 K with a mean value of 5 K. As can be seen from Figure 9, it yields an 1K-max-plume that lies between those generated by injection with constant ∆ of 5 K and 10 K. Also, a fixed of 0.24 l s -1 m -1 implies setting a ∆ that varies between 0.6 K and 4.9 K with a mean value of 2.5 K. It yields a 1K-maxplume extension between those from plumes generated by injection with constant ∆ of 2.5 K and 5 K. Therefore, when ∆ varies during the injection (between its minimal value ∆ # and its maximal ∆ #$% with an average value ∆ ), the 1K-plume extension is between those generated by an identical heat injection with a fixed temperature difference ∆ = ∆ #$% and ∆ = ∆ .

Seasonal averaging of the injection

Figure 10 shows the relative difference on the 1 K-max-plume extension obtained by considering the variable thermal load or its seasonal average. The simplified model estimates are evaluated for the various scenarios (except q-0.12 and q-0.24) after five months (using the seasonal average of the load = ), one year and twenty years (using the yearly average of the load, > ).

Estimating the seasonal impact (at five months) using a constant load = yields satisfactory results (relative error < 20 %), except for cases when the groundwater flow velocity > 2 m d -1 (see Figure 9a).

Overall, the plume width is underestimated and its length overestimated. This resulting trend is reinforced by the increase of the flow velocity and reflects the increase of relative importance of the injection peak on the plume development.

As seen in Figure 10, using the yearly average of the thermal load > to estimate the thermal impact after a year yields less accurate results. The quality of the results increases with the time until the dynamic stability state is reached. However, determining the long-term impact by means of > yields satisfactory results only for a groundwater flow velocity lower than 1 m d -1 , a longitudinal dispersivity coefficient higher than 5 m, a transverse dispersivity higher lower than 0.5 m and a temperature difference higher than 5 K. Outside this range, the quality of the estimation dramatically decreases.

Discussion

The influence of the background velocity has been assessed by using an extensive range of groundwater flow velocities representative of realistic hydrogeological conditions. Our results confirm previous findings that show the strong influence of the background velocity on the seasonal plume shape and its long-term evolution [START_REF] Piga | Thermal Impact Assessment of Groundwater Heat Pumps (GWHPs): Rigorous vs. Simplified Models[END_REF][START_REF] Russo | Development of the thermally affected zone (TAZ) around a groundwater heat pump (GWHP) system: A sensitivity analysis[END_REF]. We found a similar trend for the maximal plume width.

However, in contrast to previous studies, our results show the existence of a peak for the maximal length.

The present study provides the first insight into the individual influence of the longitudinal and transverse dispersivity coefficients on the thermal plume caused by a long-term variable injection of hot water. By varying both longitudinal and transverse coefficients separately within realistic ranges, we demonstrate their strong influence. Nevertheless, the longitudinal dispersivity does not have a substantial influence on the plume width. The influence of a relative increase of the longitudinal dispersivity on the plume length decreases for values over 2.5 m, which is less pronounced for the transverse dispersivity. Overall, the latter has a significantly stronger influence on both plume extension and stabilization time. As shown in previous studies, αt is the main driver for mixing with pristine groundwater [START_REF] Zhu | Groundwater temperature evolution in the subsurface urban heat island of Cologne, Germany[END_REF]. By using a simultaneous variation of both parameters with a constant ratio of 10 between longitudinal and transverse dispersivity, Piga et al. [START_REF] Piga | Thermal Impact Assessment of Groundwater Heat Pumps (GWHPs): Rigorous vs. Simplified Models[END_REF] found that varying the longitudinal dispersivity below 2 m is negligible for the plume extension. We do not observe such a trend for a separate variation of the parameters and thus further studies should be done to deepen the understanding of the influence of the degree of anisotropy in mechanical dispersion.

For all presented scenario groups, results highlight the crucial role of the temperature threshold chosen to qualify and understand the impact of the key parameters. On the one hand, the resulting physical parameters' influence may be different. On the other hand, these parameters influence the plume spread and therefore the distance between the computed plume isotherms. When considering the operational management of a large number of installations, the suitability of referring to arbitrary thresholds should be reconsidered. An example to improve this would be at the least to inspect plume fringes within threshold ranges.

The maximal annual temperature field as proposed in this work, facilitates the understanding of the parameters' role by providing a representation that is easier to interpret than a transient plume. In addition, it avoids underestimations of the plume extension which could occur when considering instantaneous plumes in time and it is therefore suitable for operational purposes. In practice however, it may be of interest to monitor and manage the seasonal or long-term transient plume evolution of GWHP systems in detail.

This work is intended to focus on the specific role of various parameters. Given this objective, it is based on the use of two-dimensional (2D) numerical models. If the influence of vertical heat exchange cannot be neglected in the assessment of long-term thermal impacts [START_REF] Piga | Thermal Impact Assessment of Groundwater Heat Pumps (GWHPs): Rigorous vs. Simplified Models[END_REF], they are expected to attenuate the plume extension and the sensitivity of the plume evolution on the parameters studies here. Further analyses should be done with three-dimensional (3D) models to scrutinize vertical heat loss and also to examine the role of further processes such as hydraulic recharge from rainfall.

Conclusions

The number of groundwater heat pumps (GWHP) systems has sharply increased during this last decade potentially leading to cumulative thermal impacts in several urban settings. However, these systems generate thermal impacts that strongly depend on their size, as well as on the present hydrogeological conditions. As these thermal impacts are likely to affect the efficiency of neighbouring installations, a good understanding of the role of various site-specific parameters on the thermal plume evolution is required to set up adapted management strategies. Our results deepen the knowledge of the role of the flow velocity, the longitudinal and transverse dispersity and the parameters used to describe the heat injection.

We demonstrate the strong influence of the groundwater flow velocity on both plume extension as well as the time it takes to reach a dynamic stability state. We also highlight the non-monotonous relation between the velocity and the plume length and width, with a peak that occurs for moderate flow velocity (0.5 -1.0 m d -1 ). By contrast, the plume width gradually decreases when the velocity increases and the overall lowest extension and stabilization time is obtained for high flow velocities (>> 1 m d -1 ). The dispersivity influence has been inspected by separately varying the longitudinal and transverse coefficients in a broad, realistic range. The longitudinal dispersivity has no significant influence on the plume width and stabilization time. Its influence on the length is strong for low values, and it decreases remarkably for a values over 2.5 m. In comparison, the transverse dispersivity has a more pronounced influence on the plume length, width and on the stabilization time.

Several scenarios were used to inspect the influence of the parameters defining the variable injection mode of a GWHP system. Our results show that if the maximal temperature difference reached within the plume can be kept under a threshold (by adjusting ∆ ), the hydraulic impact will be higher. For a threshold below ∆ the plume will be shorter but wider close to the well compared to when obtained with an identical power injection using a higher ∆ . Therefore, such a strategy could be used to avoid thermal interference with neighbouring installations in urban areas.

Simplified models considering the seasonal average of the thermal load have been tested in order to specify their range of applicability. The groundwater flow velocity, as well as the dispersity coefficients, strongly condition the relevance of the seasonal averaging. Also, if satisfactory results are obtained for seasonal or long-term estimations, results for mid-lengths of time are weaker. If such models cannot precisely determine the actual shape of the plume, in their scope of applicability they can yield satisfactory estimates of the maximal impact within a short simulation time. This makes such models suitable for supporting GWHP management in the field. However, considering the relative errors obtained, they should only be used for first-tier assessment of the thermal impact. Additionally, the applicability of analytical solutions such as those proposed by Pophillat et al. [START_REF] Pophillat | Analytical solutions for predicting thermal plumes of groundwater heat pump systems[END_REF] should be scrutinized for long-term assessment.

This study highlights the advantages of considering the annual maximal temperature field to represent the thermal plume generated by a variable injection. Our results also demonstrate the considerable influence of the given temperature threshold to delineate the plume boundary. Accordingly, spatial management of GWHP systems should not rely on arbitrary thresholds but should compare the role of different thresholds on the computed isotherms, or ideally should consider the full temperature field that is evolving from GWHP systems.

Figure 1 .

 1 Figure 1. Scheme of closed and open shallow geothermal systems: a) ground source heat pump (GSHP) and b) ground water heat pump (GWHP) systems.

Figure 2 :

 2 Figure 2: (a) Injected heat flux qh of reference scenario; (b) corresponding injection rate qinj and temperature difference ΔTinj.

  (i) Influence of the groundwater flow velocity: Six scenarios are introduced to evaluate the influence of the ambient groundwater flow velocity. They are implemented by varying the flow velocity of the reference case (va = 1 m d -1 ) by the factors 0.1, 0.2, 0.5, 2, 5 and 10. The corresponding scenarios are therefore defined by a flow velocity varying between va = 0.1 m d -1

Figure 3 :

 3 Figure 3: Numerical model geometry, injection well position (nodal source) and boundary conditions (BCs).

Figure 4 :

 4 Figure 4: The 20 th annual cycle of the reference scenario: 1 K-max-plume and 1 K-plume for a) day 153 (end of the injection phase), and b) day 365 (end of the phase without injection).

  between va = 0.5 m d -1 and va = 1 m d -1 . For a flow velocity of 0.5 m d -1 , the maximum plume length reaches 1463 m. It drops to 623 m for va = 1 m d -1 and for higher velocities the length of the 1 K-maxplume and the stabilization time gradually decrease. In comparison, the relative shrinking of the plume width is less pronounced. For a very high groundwater flow velocity of 10 m d -1 , the plume is very small. The high dispersion associated with this velocity results in a large spreading of the injected heat and the maximum temperature quickly falls below 1 K following the completion of a heat injection phase.

Figure 5 :

 5 Figure 5: Influence of the groundwater flow velocity: a) 1 K-max-plumes after 20 years; b) relative difference with the reference scenario depending on the flow velocity (x-axis is log scale); c) 1K-max-plume energy content and area after 20 years depending on the flow velocity (x-axis is log scale).

Figure 6 :

 6 Figure 6: Absolute lengths of the 0.5 K, 1 K and 2 K-max-plumes depending on the groundwater flow velocity (x-axis is log scale).

S3) were added

  as supplementary material to understand the individual role played by αl and αt. They show the 1K plumes (instead of 1K-max-plumes) at the end of each month of the 20 th annual cycle of injection for the scenarios related to the longitudinal and transverse dispersivity coefficients.

Figure 7 :

 7 Figure 7: Influence of the dispersity coefficient. 1 K-max-plumes after 20 years for various values of αl (a) and αt (b). Relative difference compared to the reference scenario for various values of αl (c) and αt (d) (horizontal axis is log scale).

Figure 8 :

 8 Figure 8: Influence of the threshold choice on the dispersivity coefficients influence evaluated as relative difference with the reference scenario (x-axis is log scale): a) longitudinal dispersivity; b) transverse dispersivity (length values for transverse dispersivity below 0.5 m exceed the domain size).

Figure 9 :

 9 Figure 9: Influence of the characteristics of the heat injection: 1 K-max-plumes after 20 years for scenarios considering a fixed ∆8 9:; (a) and a fixed < 9:; (b); Relative difference with the reference scenario for scenarios considering a fixed ∆8 9:; (c) and a fixed < 9:; (d).

Figure 1 :

 1 Figure 1: Relative difference on the 1K-max-plume extension computed using a variable thermal load (reference) and a constant load: a) Influence of the groundwater flow velocity (x-axis is log scale); b) influence of the longitudinal dispersivity (x-axis is log scale); c) influence of the transverse dispersivity (x-axis is log scale); d) influence of the injected temperature difference.

Table 1 : Definitions and parameter specifications of studied scenarios. Bold values denote those 186 who are varied with respect to the settings of the reference case. Additionally, those physical 187 parameter values are listed, which are kept constant in all scenarios. 188 Scenarios Parameter values and ranges

 1

Acknowledgements

This work was supported by the German Research Foundation (project number BA2850/3-1) and by the French Ministry of Ecological and Solidarity Transition. We thank Kristie Ussher and Enzio Schoffer for language corrections. Finally, we also would like to thank the two anonymous reviewers for the constructive and useful comments.

Conflicts of interest

The authors declare that there are no conflicts of interest.

Nomenclature

Symbol

Variable Unit


Giancarlo D'ambrosio

David Greynat

Marc Knecht

Matching long and short distances in the form factors for K → π +

At order O(αG F ), the amplitudes for the decays K → π + -involve a form factor given by the matrix element of the time-ordered product of the electromagnetic current with the four-quark operators describing weak non-leptonic neutralcurrent transitions between a kaon and a pion. The short-distance behaviour of this time-ordered product, when considered at order O(α s ) in the perturbative expansion of QCD, involves terms linear and quadratic in the logarithm of the Euclidean momentum transfer squared. It is shown how one can exactly match these short-distance features using a dispersive representation of the form factor, with an absorptive part given by an infinite sum of zero-width resonances following a Regge-type spectrum. Some phenomenology-related issues are briefly discussed.

Introduction

The amplitude for a ∆S = +1 neutral-current transition K(k) → π(p) + (p + ) -(p -), = e, µ, where (K, π) ∈ {(K + , π + ), (K 0 , π 0 )}, takes the form [1-3]

A K→π + - (s) = -e 2 × ū(p -)γ ρ v(p + ) × 1 s × i d 4 x π(p)|T { j ρ (0)L ∆S =1 non-lept (x)}|K(k) + - G F √ 2 V us V ud × C 7V (ν) 4πα × s × π(p)|( sγ ρ d)(0)|K(k) (1.1) = -e 2 × ū(p -)( k + p)v(p + ) × W Kπ (z) 16π 2 M 2 K .
Here s = (kp) 2 = (p + + p -) 2 = zM 2 K is the square of the di-lepton invariant mass, M K denotes the kaon mass [for our present purposes, it is not necessary to distinguish between the masses of charged and neutral kaons], e is the electric charge, G F the Fermi constant, and V ud , V us are elements of the CKM matrix. Each matrix element occurring in the first expression has to be evaluated in QCD with three active flavours. In particular, the electromagnetic current j ρ is made up from the contributions of the u, d and s quarks, j ρ (x) = q=u,d,s e q ( qγ ρ q)(x), e u = + Email addresses: gdambros@na.infn.it (Giancarlo D'Ambrosio), david.greynat@gmail.com (David Greynat), marc.knecht@cpt.univ-mrs.fr (Marc Knecht) Furthermore, L ∆S =1 non-lept (x) represents the effective Lagrangian for ∆S = +1 transitions [4][5][6][7][8][9], and involves the two current-current four-quark operators, as well as the QCD penguin operators, modulated by the appropriate Wilson coefficients

L ∆S =1 non-lept (x) = - G F √ 2 V us V ud × 6 I=1 C I (ν)Q I (x; ν). (1.3)
The above representation of A K→π + -holds at order O(αG F ), α = e 2 /4π denoting the fine-structure constant, and with a corresponding form factor for each (K, π) channel, that is for the CP conserving transitions K + → π + + - and K S → π 0 + -, but also for the direct-CP violating part of the amplitude for the decay K L → π 0 + -. The purpose of the perhaps somewhat unfamiliar first expression of A K→π + -given in Eq. (1.1) above is to explicitly display the two main components of the weak transition form factor defined by the second expression in Eq. (1.1),

W Kπ (z) = W LD Kπ (z; ν) + W SD Kπ (z; ν). (1.4)
The first part is dominated by long-distance contributions, the second one is generated at short distances, at the electroweak scale or beyond in case new physics sets in at even higher energies. Notice the appearance of an ultraviolet subtraction scale ν in each of the two terms in Eq. (1.4). Their sum W Kπ (z) appears in the amplitude A K→π + -entering the physical decay rate and should of course not depend on ν anymore. Explicitly, one has

s(k + p) ρ -(M 2 K -M 2 π )(k -p) ρ × W LD Kπ (z; ν) 16π 2 M 2 K = i d 4 x π(p)|T { j ρ (0)L ∆S =1
non-lept (x)}|K(k) (1.5) for the long-distance dominated part, whereas the contribution from short distances, which arises from a local term [7,[10][11][12]]

L ∆S =1 lept (x; ν) = - G F √ 2 V * us V ud C 7V (ν)Q 7V (x), (1.6)
with

Q 7V = ( si d i ) V-A ( ¯ ) V , reads W SD Kπ (z; ν) 16π 2 M 2 K = -- G F √ 2 V us V ud × C 7V (ν) 4πα × C Kπ f + (s)
(1.7) in terms one of the form factor f + (s) describing the matrix element of the ∆S = +1 neutral current

π(p)| s γ ρ d (0)|K(k) = C Kπ (k + p) ρ f Kπ + (s) + (k -p) ρ f Kπ -(s) .
(1.8) In these last formulas, C Kπ denotes a Clebsch-Gordan coefficient, chosen such that f + (0) = 1 in the flavour-S U(3) limit. In both Eqs. (1.5) and (1.8) the terms proportional to (kp) ρ do not contribute to the amplitude A K→π + -due to the conservation of the leptonic current, (kp) ρ ū(p -)γ ρ v(p + ) = 0. For the sake of completeness, let us recall that in the standard model there is also a contribution to the short-distance part of the amplitude

A K→π + -coming from a term proportional to C 7A Q 7A (x), with Q 7A = ( si d i ) V-A ( ¯ ) A .
This term does not play any role in the present discussion, as it does not involve the short-distance scale ν and is anyway CKM suppressed. We will therefore not mention it any further.

In W SD Kπ (z; ν) the scale dependence is entirely carried by the Wilson coefficient C 7V (ν),

ν dC 7V (ν) dν = α α s (ν) 6 J=1 γ J,7V (α s ) C J (ν).
(1.9)

The anomalous dimensions γ J,7V (α s ) are known to leading [10][11][12][13] and to next-to-leading [14] orders,

γ J,7V (α s ) = γ (0) J,7V α s 4π + γ (1) J,7V α s 4π 2 + • • • . (1.10)
In W LD Kπ (z; ν), the scale dependence arises from the singular structure, at short distances, of the time-ordered product of the electromagnetic current with the ∆S = +1 effective Lagrangian [both j µ and L ∆S =+1 non-lept are finite operators]. This short-distance singularity can be studied perturbatively within the operator-product expansion (OPE) [15,16]. After renormalization [we use dimensional regularization in the MS subtraction scheme, and below q µ denotes a Euclidian momentum, with q 2 = s < 0, whose components become simultaneously large] one obtains the general structure, cf. Ref. [3],

lim q→∞ i d 4 x e iq•x T { j µ (x)L ∆S =+1 non-lept (0)} = [q µ q ρ -q 2 η µρ ]× sγ ρ (1 -γ 5 )d (1.11) × - G F √ 2 V us V ud × 1 4π 6 I=1 C I (ν)ξ I (α s ; ν 2 /q 2 ) + O(q), with ξ I (α s ; ν 2 /s) = p≥0 p+1 r=0
ξ I pr α p s (ν) ln r (-ν 2 /s).

(1.12)

The subleading terms in the OPE have finite coefficients, and do therefore not depend on the renormalization scale ν.

Consequently [recall that z ≡ s/M 2 K ], ν dW LD Kπ (z; ν) dν = -16π 2 M 2 K C Kπ f + (s) × - G F √ 2 V us V ud × 1 4π ν d dν 6 I=1 C I (ν)ξ I (α s ; ν 2 /s). (1.13)
From Eqs. (1.11) and (1.12) one infers that in the Eucliean region z → -∞ the form factor W LD Kπ (z; ν) behaves asymptotically like powers of ln(-z) times powers of the strong coupling α s .

In the present Letter, we wish to discuss how the above short-distance behaviour can be reproduced, up to the order O(α s ), by a model involving an infinite number of equally-spaced [in mass squared] zero-width resonances. Models of this kind were considered in various contexts in the past, see for instance the articles [START_REF] Shifman | At the frontier of particle physics[END_REF][START_REF] Golterman | [END_REF][19] and references therein. Such Regge-type models find their justification in the properties of the QCD spectrum in the limit of a large number of colours [20][21][22]. More recently quite efficient methods, based on the Converse Mapping Theorem [23] and the notion of harmonic sums, were developed in order to handle such models [24][25][26]. In the case at hand, the general idea put forward in Ref. [3] consists in decomposing the long-distance dominated part of the form factor as a sum

W LD Kπ (z) = W ππ Kπ (z) + W res Kπ (z; ν). (1.14)
The first term describes the contribution from the resonant P-wave two-pion intermediate state to W Kπ (z). It is constructed upon assuming that it is given by an unsubtracted dispersion integral. The absorptive part consists of the two-pion spectral density ρ ππ Kπ (s), and is obtained upon inserting a two-pion intermediate state in the representation of the form factor given in Eq. (1.1). This contribution is not relevant for the discussion of the short-distance properties, and we refer the interested reader to Ref. [3] for details. The second term describes the contributions from the intermediate states with higher thresholds as a sum of zero-width resonances. The weight with which each resonance contributes must be chosen such as to reproduce the short-distance behaviour given in Eqs. (1.11) and (1.12). In Ref. [3], this matching has been achieved at order O(α 0 s ), which involves only a constant term and a term linear in ln(-z). Here we wish to extend this matching to the order O(α s ), where we also encounter a term quadratic in ln(-z).

The remainder of this Letter is thus organized as follows. In section 2 we determine the required coefficients ξ I 01 , ξ I 11 and ξ I 12 , relying on a renormalization-group argument given in Ref. [3]. For paedagogical reasons, we then review (Section 3) the matching at order O(α 0 s ), making the discussion of Ref. [3] simpler and, hopefully, also more intuitive. The matching at order O(α s ) is then presented in Section 4. Finally, we discuss some consequences and features of our results (Section 5).

Determination of the short-distance coefficients

From the perspective of perturbative QCD, the matrix element that defines the long-distance dominated part W LD Kπ (z; ν) is described by the Feyman diagrams of the type shown on Fig. 1. This is obviously not a realistic description of W LD Kπ (z; ν), except at short distances, when the momentum transfered to the electromagnetic current becomes large in the space-like or Euclidean domain, s → -∞. The leading contribution to the OPE of the electromagnetic current and the four-quark operators describing ∆S = +1 transitions involves the neutral-current operator sγ ρ (1 -γ 5 )d, and corresponds to the diagram (a) of Fig. 1. The calculation with a bare four-quark operator Q I is straightforward, see Ref. [3], and gives

lim q→∞ i d 4 x e iq•x T { j µ (x)L ∆S =+1 non-lept (0)} = [q µ q ρ -q 2 η µρ ] × [ sγ ρ (1 -γ 5 )d](0) × - G F √ 2 V us V ud × 1 4π 6 I=1 C I ξ I 00 -ξ I 01 2 D -4 + ln -q 2 ν 2 + O(α s ) + O(q). (2.1)
with

ξ I 01 = 1 4π 8 9 × (N c , 1 , -1 , -N c , 0 , 0) , ξ I 00 = ξ I 01 3 ×    2 NDR 5 HV , (2.2)
The result for ξ I 00 depends on the scheme used to handle the Dirac matrices in D dimensions, here either naive dimensional regularization (NDR) [27] or the 't Hooft-Veltman (HV) scheme [28,29]. The divergent part is removed through the renormalization of the (bare) coupling C 7V . Notice also that at this order the Wilson coefficients are not yet running. In the presence of QCD corrections, the four-quark operators and the corresponding Wilson coefficients are renormalized. But this does not take care of all the divergences, since two external lines of the Q I operators are closed into a loop with the insertion of the electromagnetic current, see Fig. 1. These remaining divergences are again absorbed through the higher-order renormalization of C 7V . To all orders in the powers of α s , and after renormalization through minimal subtraction in the MS scheme, the leading term in the OPE then takes the form given in Eqs. (1.11) and (1.12). From the scale dependence of the Wilson coefficients C 7V (ν), given in Eq. (1.9), and C I (ν),

ν dC I (ν) dν = 6 J=1 γ J, I (α s )C J (ν), (2.3)
one infers that the total form factor in Eq. (1.4) will be scale independent provided the equation

ν dξ I (α s ; ν 2 /s) dν + 6 J=1 γ I,J (α s )ξ J (α s ; ν 2 /s) = - γ I,7V (α s ) α s (ν) (2.4) holds. At order O(α 0 s ),

this allows to recover the values of the coefficients ξ I

01 from the known one-loop anomalous dimension matrices. Including first-order QCD corrections, one obtains, after renormalization,

lim q→∞ i d 4 x e iq•x T { j µ (x)L ∆S =+1 non-lept (0)} (2.5) = [q µ q ρ -q 2 η µρ ] × sγ ρ (1 -γ 5 )d × - G F √ 2 V us V ud × 1 4π 6 I=1 C I (ν) ξ I 00 -ξ I 01 ln -q 2 ν 2 + α s (ν) ξ I 10 -ξ I 11 ln -q 2 ν 2 + ξ I 12 ln 2 -q 2 ν 2 + O(α 2 s ) + O(q).
From the above renormalization-group argument and the known two-loop anomalous dimension matrices [30][31][32][33][34], one infers the relations

ξ I 12 = 1 (4π) 2 4 27 N c - 1 N c × 0 , -8 , +11 , N f , 0 , N f , (2.6
) and

ξ I 11 = 1 (4π) 2 8 3 N c - 1 N c ×    N c 2 , -19 18 , 17 9 , 7 6 -N c 2 , 0 , 7 6 NDR N c 2 , -5 18 , 13 9 , 7 6 -N c 2 , 0 , 7 6 HV . (2.7)
The coefficients ξ I 10 cannot be obtained this way, and their determination would require a full two-loop calculation, which we will however not need to attempt here.

The matching at order O(α 0

s )

It is convenient to represent the resonance contribution to the form factor in Eq. (1.14) as a dispersive integral

W res Kπ (z; D) = C Kπ f Kπ + (zM 2 K ) 4π dx ρ res Kπ (x; D) x -zM 2 K -i0 , (3.1)
where the spectral density is constructed order by order in the expansion in powers of the strong coupling α s , to (e). The latter also form a representative subset of the diagrams contributing to the gluonic corrections at order O(α s ). The vertex represented by a circular black blob (•) on the right of diagram ( f ) corresponds to the insertion of an order O(α s ) countertem into the lowest-order diagram (a). The other external lines represent an incoming d quark and an outgoing s quark, while the loop consists of a q-q pair, where q = u, d, s.

ρ res Kπ (x; D) = ρ res;0 Kπ (x; D) + ρ res;1 Kπ (x; D) + • • • , (3.2) (a) (b) (d) (c) (f ) (e)
and is defined in D dimensions. This last point requires some explanation. Indeed, by power counting the naive superficial degree of divergence of the diagrams in Fig. 1 is two, and becomes actually zero once the Ward identity following from the conservation of the electromagnetic current is implemented. Therefore, in four dimensions the form factor satisfies a once-subtracted dispersion relation, which would thus constitute an appropriate representation to start with. However, the information from short distances at our disposal, and summarized in the preceding section, comes from calculations done within a dimensional renormalization scheme with minimal subtraction and not within a momentum subtraction one. Since we want to make direct use of this information without transforming it first into a different scheme, we choose instead to start from an unsubtracted dispersion relation in D dimensions. As will hopefully become clear in the sequel, far from being an unnecessary complication, this choice even presents some advantages in actually guiding our intuition in the process of constructing an appropriate ansatz for the spectral density ρ res;0 Kπ (x; D) or ρ res;1 Kπ (x; D). At the one-loop level, i.e. order O(α 0 s ), the resonance representation we are looking for should match, at large negative values of s, the behaviour in the same limit of the unrenormalized QCD diagram (a) of Fig. 1. In order to reproduce a logarithmic behaviour at short distances, an infinite set of resonances is required [22]. We will consider a simple Regge-type description of the resonance spectrum in terms of equally spaced [in mass squared] zero-width states. Accordingly, the general structure of ρ res;0 Kπ (x; D) is given by

ρ res;0 Kπ (x; D) = A 0 (D)(4π) 2-D 2 M 2 ν 2 MS D 2 -2 Γ 2 - D 2
× n≥1 M 2 µ (0) n (D)δ(x -nM 2 ), (3.3
dx x ρ res;0 Kπ (x; D) = A 0 (D)(4π) 2-D 2 M 2 ν 2 MS d 2 -2 ×Γ 2 - D 2 n≥1 µ (0) n (D) n . (3.5)
A 0 (4) being so far a free parameter, one may, without loss of generality, require that

n≥1 µ (0) n (D) n = 1 + (D -4)µ (0) + O((D -4) 2 ), (3.6)
with µ (0) a constant that remains unspecified for the time being. Then

dx x ρ res;0 Kπ (x; D) = A 0 (4) - 2 D -4 + ln ν 2 M 2 (3.7) -2 A 0 (4) A 0 (4) -2µ (0) + O(D -4),
with ν ≡ ν MS e -γ E /2 √ 4π the subtraction scale in the MS scheme. The constant µ (0) can actually be absorbed without loss of generality into A 0 (4), which has also not been specified so far. The remaining, subtracted, dispersive integral

wM 2 dx x ρ res;0 Kπ (x; D) x + wM 2 = A 0 (D)(4π) 2-D 2 M 2 ν 2 MS d 2 -2 (3.8) × Γ 2 - D 2 n≥1 wµ (0) n (D) n(n + w)
should then be finite as D → 4. This in turn will be the case if we require (3.11) Finally, in order to reproduce the correct matching with the short-distance behaviour, one must also require

µ (0) n (D) = (D -4) μ(0) n + O((D -4) 2 ) (3.
-2wξ(w) = ln w + Cst + • • • (3.12)
as w → +∞. At this stage, one may observe that µ (0) n ∝ (D -4) μ(0) , with μ(0) a constant, would provide a convergent series ξ(w), endowed with an asymptotic logarithmic behaviour. But it would fail to satisfy the condition (3.6). As we now show, this defect can be easily repaired. Intuitively, the task will consist in providing a convergence factor for the sum (3.6) when D < 4, but which is no longer operative for D = 4, otherwise the sum ξ(w) would converge too quickly, and would no longer exhibit a logarithmic behaviour for large values of w. Instead, it would rather behave as a constant or an inverse power of w. Consider therefore the ansatz

µ (0) n (D) = a (0) (D)n D 2 -2 . (3.13)
Then one has

n≥1 µ (0) n (D) n = a (0) (D)ζ 3 - D 2 , ( 3
a (0) (D) = 1 ζ 3 -D 2 + f (D) , (3.15)
where f (D) is an arbitrary function regular at D = 4, will lead to

n≥1 µ (0) n (D) n = 1 + f (4) 2 (D -4) + O((D -4) 2 ) (3.16)
and thus satisfy the condition (3.6) with µ (0) = f (4)/2. The condition (3.9) is then also satisfied, with μn = -1/2, so that

Γ 2 - D 2 µ n (D) = 1 + O(D -4).
(3.17)

The condition (3.10) is met as well, with

-2wξ(w) = n≥1 w n(n + w) = γ E + ψ(1 + w). (3.18)
Putting everything together, one ends up with

dx ρ res;0 Kπ (x; D) x + wM 2 = A 0 (4) - 2 D -4 + ln ν 2 M 2 -2 A 0 (4) A 0 (4) -2µ (0) -γ E + ψ(1 + w) + O(D -4). (3.19)
In this last formula, the di-gamma function ψ resums the dispersive integral

ψ(1 + w) = -γ E +M 2 w dx x 1 x + M 2 w n≥1 M 2 δ(x -nM 2 ).
(3.20) Finally, for large positive w the di-gamma function behaves as ψ(1 + w) ∼ ln w = ln(-s/M 2 ), so that the correct short-distance behaviour is also recovered, provided one takes

A 0 (4) = -16π 2 M 2 K - G F √ 2 V us V ud × 6 I=1 C I (ν)ξ I 01 (3.21)
and

A 0 (4) + µ (0) + γ E 2 A 0 (4) = +16π 2 M 2 K - G F √ 2 V us V ud × 6 I=1 C I (ν) ξ I 00 2 . (3.22)
Minimal subtraction of the divergence leads to the renormalized dispersion relation

dx ρ res;0 Kπ (x; D) x + wM 2 = -16π 2 M 2 K - G F √ 2 V us V ud × 6 I=1 C I (ν) × ξ I 01 ln ν 2 M 2 -ψ(1 + w) + ξ I 00 + O(D -4), (3.23)
with

ρ res;0 Kπ (x; D) = -16π 2 M 2 K - G F √ 2 V us V ud × 6 I=1 C I (ν) × ξ I 01 - D -4 2 ξ I 00 + ( f (4) + 2γ E )ξ I 01 × M 2 ν 2 D 2 -2 Γ 2 - D 2 (3.24) × n≥1 M 2 n D 2 -2 ζ 3 -D 2 + f (D) δ(x -nM 2 ).
We may now, in some sense, reverse-engineer the whole construction, starting directly from the representation of ρ res;0 Kπ (x; D) in D dimensions given in Eq. (3.24) above and showing that it satisfies the required properties. This will be useful in Section 4, where we will only sketch the construction of the spectral density ρ res;0 Kπ (x; D), give the result, and show that it indeed exhibits the appropriate features. Using the Mellin representation

1 1 + w n = c 1 +i∞ c 1 -i∞ du 2πi w n -u Γ(u)Γ(1 -u), (3.25)
where 0 < c 1 < 1, one has

dx ρ res;0 Kπ (x; D) x + wM 2 = -16π 2 M 2 K - G F √ 2 V us V ud × 6 I=1 C I (ν) ξ I 01 - D -4 2
× ξ I 00 + ( f (4) + 2γ E )ξ I 01 × M 2 ν 2 D 2 -2 × Γ 2 -D 2 ζ 3 -D 2 + f (D) × c 1 +i∞ c 1 -i∞ du 2πi w -u Γ(u)Γ(1 -u) ×ζ 3 - D 2 -u . (3.26)
The first singularity of the integrand lying on the left of the fundamental band c 1 = Re u ∈]0, 2 -D/2[ is a simple pole at u = 0, coming from the factor Γ(u). According to the Converse Mapping Theorem [23] one therefore has

dx ρ res;0 Kπ (x; D) x + wM 2 = w→0 -16π 2 M 2 K - G F √ 2 V us V ud × 6 I=1 C I (ν) ξ I 01 - D -4 2
× ξ I 00 + ( f (4) + 2γ E )ξ I 01 × M 2 ν 2 D 2 -2 × Γ 2 -D 2 ζ 3 -D 2 ζ 3 -D 2 + f (D) = D→4 -16π 2 M 2 K - G F √ 2 V us V ud × 6 I=1 C I (ν) -2 D -4 ξ I 01 + ξ I 00 + ξ I 01 γ E + ξ I 01 ln M 2 ν 2 + O(D -4). (3.27)
The first singularity of the integrand lying on the right of the fundamental band

c 1 = Re u ∈]0, 2 -D/2[ is a sin- gle pole occurring at u = 2 -D/2, with ζ 3 -D 2 -u ∼ -1/(u -2 + D/2), so that, according to the Converse Map- ping Theorem, dx ρ res;0 Kπ (x; D) x + wM 2 ∼ w→+∞ -16π 2 M 2 K - G F √ 2 V us V ud × 6 I=1 C I (ν) ξ I 01 - D -4 2 ξ I 00 + ( f (4) + 2γ E )ξ I 01 × M 2 ν 2 D 2 -2 × Γ 2 -D 2 ζ 3 -D 2 + f (D) × Γ 2 - D 2 Γ D 2 -1 × w D 2 -2 + O w -1 .
(3.28)

Considering next the limit D → 4, one indeed recovers the expected result:

dx ρ res;0 Kπ (x; D) x + wM 2 ∼ w→+∞ -16π 2 M 2 K - G F √ 2 V us V ud × 6 I=1 C I (ν) (3.29) × -2 D -4 ξ I 01 + ξ I 00 + ξ I 01 ln ν 2 -s + O(D -4).
Summarizing the preceding analysis, we find inded that the spectral density (3.24) reproduces the minimally subtracted dispersive integral in Eq. (3.23). Let us notice, at this stage, that the result (3.23) for D = 4 does no longer depend at all on f (D). We will return to this point in Section 5 below.

The matching at order O(α s )

In order to include the O(α s ) corrections that arise at two loops, it is necessary to start with a somewhat more involved expression of the spectral density. We again let ourselves be guided by the perturbative structure of this contribution. On the basis of the order O(α s ) diagrams shown in Fig. 1, we are led to consider as a starting point the sum of two terms,

ρ res;1 Kπ (x; D) = ρ res;1a Kπ (x; D) + ρ res;1b Kπ (x; D) (4.1)
with

ρ res;1a Kπ (x; D) = A 1a (D)(4π) 4-D M 2 ν 2 MS D-4 Γ 2 - D 2 2 × n≥1 M 2 µ (1a) n (D) δ(x -nM 2 ), ρ res;1b Kπ (x; D) = A 1b (D) D -4 (4π) 2-D 2 M 2 ν 2 MS D 2 -2 Γ 2 - D 2 × n≥1 M 2 µ (1b) n (D) δ(x -nM 2 ), (4.2)
where both functions A 1a (D) and A 1b (D) are regular at D = 4. The first term, ρ res;1a Kπ (x; D), corresponds to the genuine two-loop diagrams, like the graphs (b) to (e) of Fig. 1. The second term, ρ res;1b Kπ (x; D), stands for the diagram ( f ), i.e. the one-loop graph (a) with the insertion of a one-loop counterterm proportional to 1/(D-4). Without loss of generality, one may impose a normalization condition like (3.6) for each set of coefficients separately, and from there proceed as in Section 3. We will not go through the details of this straightforward exercise, but rather state the final result and show that it satisfies the required properties. Before that, let us make a few useful remarks.

With the normalizations set as in Eq. (3.6) for both ρ res;1a Kπ (x; D) and ρ res;1b Kπ (x; D), the corresponding dispersive integrals at w = 0 display double and single poles at D = 4. The single poles contain a contribution proportional to ln(M2 /ν 2 ), which, on general grounds [START_REF] Collins | [END_REF], is not allowed, and hence has to cancel in the sum of the two dispersive integrals. For this to happen, we need to impose the condition A 1b (4) = 4A 1a (4).

Next, ρ res;1b

Kπ (x; D) has to describe the same structure than ρ res;0 Kπ (x; D) up to an additional factor ∝ 1/(D -4) coming from the inserted counterterm. It is thus natural to consider for it the ansatz

µ (1b) n (D) = n D 2 -2 ζ 3 -D 2 + f (D) , (4.3)
where f (D) is an arbitrary function regular at D = 4. Then the resulting subtracted dispersive integral

wM 2 dx x ρ res;1b Kπ (x; D) x + wM 2 (4.4)
still contains a term proportional to ψ(1 + w)/(D -4), which has eventually to be canceled by a similar contribution from the dispersive integral involving ρ res;1a Kπ (x; D). Turning to the latter, we consider the ansatz

µ (1a) n (D) = n y(D-4) ζ (1 + 4y -yD) + g(D) , (4.5)
where g(D) is an arbitrary function regular at D = 4, and y is a so-far free parameter. Considering the dispersive integral

wM 2 dx x ρ res;1a Kπ (x; D) x + wM 2 , (4.6)
we find that the required cancellation of the unwanted pole terms takes place for the choice y = 1.

Adjusting the remaining free parameters such as to reproduce the asymptotic behaviour of Eq. (2.5), we finally arrive at the result1

ρ res;1 Kπ (x; D) = -16π 2 M 2 K - G F √ 2 V us V ud α s (ν) 6 I=1 C I (ν) × M 2 ν 2 D-4 Γ 2 - D 2 2 e γ E (4-D) ξ I 12 × n≥1 M 2 n D-4 ζ (5 -D) δ(x -nM 2 ) + M 2 ν 2 D 2 -2 Γ 2 - D 2
e γ E (4-D)/2 D -4 × 4ξ I 12 + (D -4) ξ I 11 + 2γ E ξ I 12 - (D -4) 2 2 ξ I 10 + ξ I 11 γ E -2ξ I 12 π 2 6 + γ 1 × n≥1 M 2 n D 2 -2 ζ 3 -D 2 δ(x -nM 2 ) . (4.7)
It is a straightforward exercise, making use of the techniques described in the second part of Section 3, to show that this spectral density leads to the desired properties. For instance, after minimal subtraction of the simple and double poles at D = 4, the renormalized dispersive integral reads

dx ρ res;1 Kπ (x; D) x + wM 2 = -16π 2 M 2 K - G F √ 2 V us V ud α s (ν) 6 I=1 C I (ν) ξ I 12 ψ(w) + 2ξ I 12 ln M 2 ν 2 -ξ I 11 ψ(1 + w) + ξ I 10 -2 π 2 6 + γ 1 ξ I 12 + ξ I 12 ln M 2 ν 2 -ξ I 11 ln M 2 ν 2 + O(D -4), (4.8)
with

ψ(w) = 2M 2 w dx x 1 x + M 2 w n≥1 M 2 (ln n)δ(x -nM 2 ) = 2w n≥1 ln n n(n + w) . (4.9)
We now need to establish some properties of the function ψ(w), which is defined by an absolutely convergent sum as long as w is not equal to a negative integer. Thus, ψ(0) = 0 and The first singularity on the left of the fundamental strip 0 < c 1 < 1 lies at u = 0. It consists of a simple pole coming from the factor Γ(u) of the integrand. Therefore, the Converse Mapping Theorem allows to state that ψ(w) = -2wζ (2) + • • • as w → 0, a property already established before. The first singularity one encounters on the right of the fundamental strip consists of a triple pole and a simple pole, both located at u = +1,

ψ (0) = 2 n≥1 ln n n 2 = -2ζ (2).
Γ(u)Γ(1-u)ζ (2-u) ∼ u→+1 1 (u -1) 3 + π 2 6 + γ 1 u -1 +γ 2 +O(u-1).
(4.13) The Converse Mapping Theorem then allows to conclude that

ψ(w) ∼ w→+∞ ln 2 w + 2 π 2 6 + γ 1 + • • • , (4.14)
which is precisely what is required. The function ψ(w) does not seem to be related in an obvious way to the standard sets of functions that have been studied in the mathematical literature, see e.g. [START_REF]NIST Digital Library of Mathematical Functions[END_REF] and references therein.

It is thus interesting to notice that the Mellin representation (4.12) of ψ(w), which holds for Rew > 0, can be recast into [the integral is understood as its Cauchy principal value for w real and positive]

ψ(w) = 2 ∞ 0 dx x w x -w [ψ(1 + x) + γ E ] +2[ψ(1 + w) + γ E ] ln w. (4.15)
This representation allows to extend the function ψ(w) to negative values of Rew. The integral is then regular, while the second term reproduces the poles of ψ(w) when w equals a strictly negative integer. The determination of ln w is given by w → wi , with > 0 and infinitesimal, in agreement with the prescription s → s + i . This representation also proves quite useful for the numerical evaluation of ψ(w).

Discussion

We have shown that it is possible to construct a function W res Kπ (z; ν) through an infinite sum of zero-width resonances with a Regge-type spectrum, such that its shortdistance behaviour matches, at order O(α s ), the one obtained from the leading term in the OPE. After renormalization in the MS scheme of dimensional regularisation, and up to the order O(α s ), this function reads Kπ (z; ν) to lowest order in the NDR (HV) scheme. For the ease of visualization, the location of the poles at z = nM 2 /M 2 K has been slightly shifted off the real-z axis.

[w ≡ -zM 2 K /M 2 ]

W res Kπ (z; ν) = -16π 2 M 2 K - G F √ 2 V us V ud C Kπ f Kπ + (zM 2 K ) 4π × 6 I=1 C I (ν) ξ I 00 + α s (ν)ξ I 10 + ln ν 2 M 2 ξ I 01 + α s (ν) ξ I 11 + ξ I 12 ln ν 2 M 2 -ψ(1 + w) ξ I 01 + α s (ν) ξ I 11 + 2ξ I 12 ln ν 2 M 2 +α s (ν) ξ I 12 ψ(w) - π 2 3 -2γ 1 ≡ -16π 2 M 2 K - G F √ 2 V us V ud C Kπ f Kπ + (zM 2 K ) 4π ×W res Kπ (z; ν). (5.1)
This expression involves two scales. The first one, M, is the scale of the lowest resonance in the spectrum besides the ρ(770), which is already taken into account by the contribution W ππ Kπ (z) in Eq. (1.14), see Ref. [3]. The first resonance appearing in W res Kπ (z; ν) is thus the K * (892) or the φ(1020), hence M ∼ 1 GeV. The second scale in Eq. (5.1) is the renormalization scale ν. It represents the onset of the perturbative regime of QCD. In practice, the description of the spectrum in terms of well-identified resonances extends to a few radial excitations of the resonances mentioned just above, for instance ρ(1450), ρ(1570), φ(1680)..., before they merge into the continuum. This means ν ∼ 2 GeV.

The function W res Kπ (z; ν) is shown in Fig. 2. For positive values of z it displays the expected infinite series of equally-spaced poles, whereas for negative values of z the asymptotic regime sets in rapidly. In order to draw these plots, some knowledge of the coefficients ξ I 10 , which are not fixed by the renormalization-group constraint (2.4), is needed. For the sake of illustration, we have estimated these coefficients to vary in the range -ξ I 11 ≤ ξ I 10 ≤ +ξ I 11 , taking the relation (2.2) between the lowest-order coefficients ξ I 00 and ξ I 01 as a guide. The same choice also applies to Fig. 3.

In the remainder of the Letter, we wish to address in turn two issues that we think deserve to be given some consideration, namely: i) the size of the residual dependence with respect to the short-distance scale ν and ii) some features and properties of the resonance model that we have constructed, as well as possible improvements.

Residual scale dependence

By construction, the expression of W res Kπ (z; ν) displayed in Eq. (5.1) provides a form factor that is independent of the renormalization scale ν at order O(α s ),

ν d dν W res Kπ (z; ν) + W SD Kπ (z; ν) = O(α 2 s ). (5.2)
We expect the residual scale dependence to be weaker that the one that results from the matching at lowest order only. We illustrate these changes in the case of the two constants a + and b + corresponding to the values of the form factor W Kπ (z), for (K, π) = (K + , π + ), and of its derivative at z = 0, respectively. For this purpose, we add to the sum of W res Kπ and W SD Kπ the contribution W ππ Kπ of the two-pion intermediate state evaluated in Ref. [3], see Eq. (1.14) above. This means a + | ππ = -1.58 and b + | ππ = -0.76. The improvement when going from lowest order to next-to-leading order can be appreciated from the plots shown in Fig. 3 [details on the numerical aspects can be found in the appendix]. We notice that indeed both the scale dependence and the scheme dependence become less pronounced when order O(α s ) corrections are included. Our ignorance of the coefficients ξ I 10 induces, under the conditions stated above, an uncertainty that amounts, at ν = 2 GeV, to about 5% in both a + and b + . We also observe that the corrections coming from the resonance and short-distance parts are quite small as compared to the contribution from the two-pion state.

This improvement in the control over the scheme and scale deprendences allows us to refine somewhat our evaluation of the coefficients a + and b + made in Ref. [3]. We obtain a + = -1.59 (8), b + = -0.82( 6).

(5.3)

These values are still at variance with the experimental determinations a exp + = -0.593 (11) and b exp + = -0.675(43), see Ref. [3] for a detailed discussion, especially as far as a + is concerned. On the other hand, as discussed in Ref. [3], in order to reach a definite conclusion, the contributions a + | ππ and b + | ππ need to be evaluated within a tighter framework than the one adopted there, and other exclusive contributions, in particular from two-kaon states, should eventually be accounted for explicitly, see also the discussion below. Work in this direction is on its way [37].

Properties and features of the resonance model

The resonance model that we have constructed in Sections 3 and 4 is by far not unique. Discussing the arbitrariness of the construction in full generality represents a formidable, if not impossible, task, given the fact that the corresponding dispersive integrals are only constrained to reproduce the divergences of perturbative QCD at z = 0 and the leading asymptotic behaviour from the OPE at z → -∞. The determination of sub-leading terms in the OPE might be a way to constrain the resonance model further. We will, however, not address this possibility here, but rather state a few remarks concerning the model that we have constructed above. Some of the arbitrariness of the construction is embodied in the functions A 0 (D) and f (D) in Section 3, or A 1a (D), A 1b (D) f (D) and g(D) in Section 4. These functions do not appear anymore in the final expression (5.1), partly because they are absorbed in the matching to the coefficients ξ I rp , partly because they contribute only to the terms of order O(D -4) in the renormalized dispersive integral.

These features certainly do not exhaust all the arbitrariness of the model. While the contribution of the resonances with higher masses will be constrained by the

Figure 1 :

 1 Figure 1: Diagrams contributing to the leading short-distance behaviour of the time-ordered product of the electromagnetic current (materialized by the wiggly line on the left) with one of the four-quark operators, whose insertion is shown by the crossed-circle (⊕) vertex on the right of diagram (a), of order O(α 0 s ), and of the diagrams (b)to (e). The latter also form a representative subset of the diagrams contributing to the gluonic corrections at order O(α s ). The vertex represented by a circular black blob (•) on the right of diagram ( f ) corresponds to the insertion of an order O(α s ) countertem into the lowest-order diagram (a). The other external lines represent an incoming d quark and an outgoing s quark, while the loop consists of a q-q pair, where q = u, d, s.

4 )

 4 ) with so-far unspecified functions A 0 (D) and µ(0) n (D). The additional prefactors simply account for the structure of diagram (a) of Fig.1in D dimensions. Since the one-loop divergence is already contained in the factor Γ(2 -D/2), A 0 (D) has to be regular at D = 4. It accounts for the scheme dependence, for instance in handling Dirac matrices in D dimensions, see Section 2. M denotes the mass of the lowest-lying resonance and ν MS denotes the renormalization scale in the minimal subtraction scheme. Then, with w ≡ -s/M 2 , dx ρ res;0 Kπ (x; D) x + wM 2 = A 0 (D)(4π) 2-The value of the above dispersive integral at w = 0 reproduces the divergent part of the diagram (a) of Fig.1

  Kπ (x; D) x + wM 2 = A 0 (4)(-2w)ξ(w) + O(D -4).

(4. 10 )

 10 Next, using the relation (3.25) and n≥1 n u-2 ln n = -ζ (2u), (4.11) one obtains the following Mellin representation of the function ψ(w) ψ(w) = -2w c 1 +i∞ c 1 -i∞ du 2πi w -u Γ(u)Γ(1-u)×ζ (2-u). (4.12)

Figure 2 :

 2 Figure 2: The function W res Kπ (z; ν) defined in Eq. (5.1) as a function of z for ν = 2 GeV and M = 1 Gev. The light-red (purple) band corresponds to the NDR (HV) scheme. The red (blue) line coresponds to the truncation of W resKπ (z; ν) to lowest order in the NDR (HV) scheme. For the ease of visualization, the location of the poles at z = nM 2 /M 2 K has been slightly shifted off the real-z axis.

Figure 3 :

 3 Figure 3: The variation of a + (left plot) and of b + (right plot) as a function of the short-distance renormalization scale ν. The scale dependence at lowest order corresponds to the red (blue) lines in the NDR (HV) scheme. The scale dependence after inclusion of O(α s ) corrections corresponds to the light-red (purple) bands in the NDR (HV) scheme.

Up to a factor (-1) n , the Stieltjes constants γ n give the coefficients of the Taylor expansion of the regular part of ζ(s) at s = 1, with γ 0 = γ E = 0.577216 • • • , γ 1 = -0.0728 • • • , γ

= -0.00969 • • • ; see [35, section 25.2].

Acknowledgements

One of us (M.K.) wishes to thank the INFN-Sezione di Napoli and the Universitá di Napoli Federico II for their warm hospitality. The work of M.K. has received partial support from the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the "Investissements d'Avenir" French government program managed by the ANR. The Feynman diagrams displayed in Fig. 1 have been drawn using Jaxodraw 2.1-0 [38,39].

short-distance behaviour, one might expect that the descritpion of the lower-lying resonances like K * (892) or φ(1020) [recall that the contribution from ρ(770) is already taken care of by W ππ Kπ ] provided by this model may be less realistic from a phenomenological point of view. One way to circumvent this possible drawback would be to consider additional intermediate states in a more explicit way, i.e. by extending the decomposition in Eq. (1.14) to, for instance [3]

We will, however, not pursue this matter in the present Letter, and leave the discussion of such an extension and of the corresponding multi-channel analysis it requires for future work.

Appendix

In this appendix, we gather some information and give the input values used in order to produce the plots shown in Fig. 3. The running of the Wilson coefficients is given in Eqs. (1.9) and (2.3) using the anomalous dimensions γ J,7V given in [14] and γ I, J given in [33,34]. We have, however, restricted ourselves to the contributions from the two current-current operators C ± = C 2 ± C 1 , neglecting those from the QCD penguin operators C I , I = 3, 4, 5, 6. The Wilson coefficients at next-to-leading order are then given by

, (A.1)

, and

where in the HV scheme. These expressions for C ± (ν) and for C 7V (ν) hold at next-to-leading order, and their truncation the lowest order is obtained upon taking all the coefficients S i j equal to zero. The input values we have used are [14] C + (ν 0 ) = 0.688, C -(ν 0 ) = 2.108, C 7V (ν 0 )/α = -0.026 (A.5) at lowest-order. At next-to-leading order they become at next-to-leading order, and with the obvious truncation at lowest order. The QCD scale for three active flavours is taken as Λ = 332 MeV. Finally, the mass of the lowest resonance is set at M = 1 GeV.