Maxime Schmitt
email: max.schmitt@unistra.fr

Philippe Helluy
email: helluy@unistra.fr

Cédric Bastoul
email: cedric.bastoul@unistra.fr

Automatic Adaptive Approximation for Stencil Computations

Keywords: Theory of computation, API languages, Approximate Computing, Compilation, Code Optimization, Application Programming Interface, Stencil Code

Approximate computing is necessary to meet deadlines in some compute-intensive applications like simulation. Building them requires a high level of expertise from the application designers as well as a significant development effort. Some application programming interfaces greatly facilitate their conception but they still heavily rely on the developer's domain-specific knowledge and require many modifications to successfully generate an approximate version of the program. In this paper we present new techniques to semi-automatically discover relevant approximate computing parameters. We believe that superior compiler-user interaction is the key to improved productivity. After pinpointing the region of interest to optimize, the developer is guided by the compiler in making the best implementation choices. Static analysis and runtime monitoring are used to infer approximation parameter values for the application. We evaluated these techniques on multiple application kernels that support approximation and show that with the help of our method, we achieve similar performance as non-assisted, hand-tuned version while requiring minimal intervention from the user.

INTRODUCTION

Adaptive approximation aims at providing targeted, relaxed computations while preserving good precision on demanding computation kernels. Automatic compiler adaptive optimization unlocks the possibility to use this kind of optimization not available to a majority of non-expert developers. Adaptive techniques are particularly useful in the context of applications such as scientific computing and data analysis where time constraints are important and the demand for computing capacity grows faster than the available resources. Because they respect language standards, compilers strictly adheres to the program semantics and are not allowed to use approximate computing transformations. Only a limited set of compiler options relax compliance to standards to enable optimization, e.g. considering that some floating point operators are associative may dramatically improve automatic parallelization. In this paper we present a new annotation-based technique that relaxes the program semantics and lets the compiler automatically optimize the annotated code section using adaptive approximation techniques.

Compute intensive applications may require a considerable execution time to produce a result. Approximation techniques provide a solution to trade result accuracy for execution time. Adaptive methods is a powerful class of numerical analysis techniques which targets the approximation in the dynamic regions where their effect on the output solution is minimal. These methods provide strong approximation guarantees while lowering the size of the data and computation to an acceptable level [START_REF] Bungartz | Sparse grids[END_REF]. Adaptive techniques are well suited for applications with localized disturbances such as numerical simulations, signal or data processing. In this work, we exploit the fact that substantial solution variation, which implies that precise computations should be done in the region where it happens, may be detected at runtime.

Using adaptive techniques requires a high level of expertise and a fair amount of application engineering to be used on purpose and efficiently. To open such optimization to a wider audience, we propose a semi-automatic, compiler based, optimization which generates adaptive version of an existing code base. We target Adaptive Code Refinement (Section 3), an application programming interface that provides annotations to express approximate computing transformations (code alternatives) and monitoring capabilities for adaptive methods.

The automatic generation of adaptive code relies on the extraction of the following features:

• The monitoring of the application's data at runtime to identify the region of adaptive interest (Section 4.1). This allows the compiler's runtime to pinpoint the application regions where approximation can advantageously be applied. • The information about the type of computation performed by the kernel (Section 4.2). This allows the compiler to generate alternative code versions, by using approximation methods to lower the computational complexity with possible deviation of the result of the application. • The granularity of the adaptive grid at which the code alternatives are applied (Section 4.3). This allows the detection of coarse or fine grained disturbance regions and select the size that presents the best trade-off between computation gain and the deviation of the application's output. We use the information gathered by our method to automatically select the best parameters and code alternatives to provide to ACR. In section 5 we demonstrate that our technique can detect useful metrics and use them to generate an adaptive version on a range of representative applications.

In this paper we make the following contributions:

• We propose a novel method to automatically generate adaptive code from target compute intensive kernels. • We identify the information required by a compiler to generate adaptive codes and provide algorithms to extract them. • We present a code transformation to automatically generate approximate version for stencil computation kernels, which is one of the main class, if not the main class, of computation kernels in application which supports approximation. • We provide experimental evidence that our technique is successful at generating adaptive versions which compete with hand-tuned versions.

MOTIVATING EXAMPLE

Existing automatic approximation tools [START_REF] Ansel | Language and compiler support for auto-tuning variable-accuracy algorithms[END_REF][START_REF] Baek | Green: A Framework for Supporting Energy-conscious Programming Using Controlled Approximation[END_REF][START_REF] Vinay K Chippa | Analysis and characterization of inherent application resilience for approximate computing[END_REF]33] are able to search for applicable approximation, but none of them, to the best of our knowledge, supports adaptive techniques. Our work enables automatic adaptive optimization through language annotations and compiler techniques. To illustrate our approach, we consider the main computational kernel of a steady-state heat equation solver [START_REF] Geankoplis | Transport processes and unit operations[END_REF] shown in Figure 1. An example annotation has been added at line five. This annotation states that the compiler is allowed to generate approximate version for the code block that follows.

The compiler runs both static and profiling analyses on the kernel, here the loop nest, and extracts information required to generate an adaptive version. Three main features are extracted to allow the compiler to generate the adaptive versions automatically:

(1) Code alternatives -Generate kernel alternatives with different levels of approximation: for this example, the compiler detects that the kernel uses stencil computation for which it is able to generate different versions with different approximation levels. (2) Adaptive decision -Find a metric to decide at runtime about the convenient level of approximation to use in a given subset of the computation space: in this example, there is only one written array in this code, hence the deviation of its values drives the selection of the alternative kernels (low deviation corresponds to low precision and high deviation corresponds to high precision). (3) Decision granularity -Select the cell size when decomposing the computation space into a regular grid: profiling information allows to select the convenient cell size to get the best precision/performance trade-off.

We use the ACR programming interface (Section 3) to generate the adaptive code from the extracted high level adaptive features. Section 5.1 shows our technique applied on this example, along with the generated ACR annotations.

During the execution of the optimized program, a runtime system selects, for each cell, the best suited alternative with respect to the dynamic data values, i.e., achieving precise computation only where it matters. In this example, we reduced the number of floating point computation by 74.7% and achieve a speedup of 1.81 with a deviation below 10 -4 , see Section 5 for the complete experimental study.

ADAPTIVE CODE REFINEMENT AND EXTENSIONS

Adaptive Code Refinement (ACR) is a recent language and compiler approach proposed by Schmitt et al. which provides a simple interface for developers that want to exploit approximate computing techniques [START_REF] Schmitt | Adaptive Code Refinement: A Compiler Technique and Extensions to Generate Self-Tuning Applications[END_REF][START_REF] Schmitt | Semi-Automatic Generation of Adaptive Codes[END_REF]. Our goal is to demonstrate that our technique, when applied to ACR, leads to an automation of this tool. In this section we present ACR's use case, its user interface, the underlying compiler infrastructure and the code transformation and generation mechanism that produce alternative versions of the source code automatically. We introduce a new annotation and special construct parameters that triggers the automatic discovery of pertinent approximations and ACR parameters.

ACR targets applications with a compute intensive kernel that updates data multiple times, e.g., numerical simulation or video frame encoding. In such applications, a relation exists between the data used by two consecutive kernel calls. ACR exploits this property to reduce the computation intensity locally, where the expensive update is known or expected to to have little effect on the result.

ACR's Application Programming Interface

ACR framework relies on supplementary semantics provided by the developer through code annotations. These annotations can be inserted at hot spots of the application and allows the compiler to apply additional code transformations, related to approximate computing. The set of annotations is implemented for the C language using compiler directives (#pragma acr). The following compiler directives are provided by ACR to annotate the input source code: #pragma acr grid(Size)

The grid construct provides the ability to break down the original problem data space into hyper-squares shaped chunks of the size provided as input. Each of these chunk is called "a cell". A cell represents the atomic element over which transformation and alternative computation selection are performed.

#pragma acr monitor(Array, Folding_function)

The monitor construct specifies the application's relevant data to the compiler for ACR. This data will be used by the decision mechanism, to drive the dynamic selection of the approximate alternative computation. The folding function is a way for the compiler to extract the information from a cell. This function summarizes the data cell into one unique observation. An observation is defined as a natural number whose value represents the level of approximation that a for ( size_t i = 1; i < sizeX -1; ++ i )

7
for ( size_t j = 1; j < sizeY -1; ++ j ) { Figure 1: Heat conduction solver using the Jacobi method. The left part of the figure shows the kernel which updates the temperature and computes the maximum temperature change during an update step. The solver is executed repeatedly until the simulation settles, which happens whenever the maximum temperature delta is less than a given threshold, e.g., in this example the maximum difference is less than 10 -4 degrees Kelvin. An example simulation output is displayed on the right part of this figure. In this simulation the border conditions are set to a constant temperature. cell should withstand. An observation with a value of zero means that no approximation is allowed.

8 tNext [i ][ j ] = 0.25 * 9 (t[i -1][ j ] + t [ i + 1][ j ] + 10 t[i ][ j -1] + t [ i ][ j + 1]);

#pragma acr alternative AltId(Type, TParam)

The alternative construct defines possible alternative computation to be generated by the compiler. The alternative computation is generated from the original kernel by application of specific code transformations. The type of supported transformations are: Parameter modifies the value of a constant parameter with a constant value or another parameter Code allows the developer to provide an alternative code block Zero-compute treats the target code block as if it is empty Interface-compute disables computation that do not contribute to update the data located on an interface between two grid cells Some alternatives take a parameter here represented as TParam, e.g., parameter type accepts a variable identifier or a constant value. The user provides a unique alternative identifier using AltId. This identifier is used within the strategy construct to select the conditions required to use this alternative.

#pragma acr strategy Type (SParams, AltId)

This construct defines the compiler's strategy, i.e., the necessary conditions for which an alternative should be active.

The strategy may be one of the following types:

Static: the alternative will unconditionally be active. The user can choose whenever to apply it on localized portions or the whole data space. Dynamic: the alternative will be active depending on runtime observation of the monitored data defined with the monitor construct. This strategy applies at the cell granularity. A cell will use the alternative if the folding function returns the same observation as the one defined by this strategy.

In addition to existing ACR pragmas, we propose the following new Automatic Stencil Approximation (ASA) constructs that allows the compiler to automatically generate ACR annotations with pertinent alternatives and parameters: #pragma asa auto (AllowedPrecisionLoss)

This construct replaces all the previously mentioned constructs. It asks the compiler to find the best possible alternatives and parameters for the target kernel.

AllowedPrecisionLoss specifies the maximum deviation for which approximate computation may be used.

#pragma asa interactive (AllowedPrecisionLoss)

This construct is similar to auto but lets the user select between a range of compiler-proposed alternative and parameters.

These new constructs rely on compiler facilities, here ACR, and additional information to generate the approximate versions with little intervention of the user.

ACR's Compiler Infrastructure

With the additional semantic information provided by the user annotations, ACR can generate a statically or dynamically optimized program. The static version is generated at compile time and the dynamic version relies on just-in-time compilation to generate a specialized version of the kernel. ACR allows for a fixed Cartesian grid overlay to delimit the cells (cf. the grid annotation in Section 3.1). The dynamic version allows the compiler to optimize a kernel where dependencies forbid the static transformation or when the user instructs the compiler to select the grid size dynamically.

ACR relies on the polyhedral representation of programs to apply transformations and generate an optimized version of the code. Figure 2 discloses code examples generated from the static and dynamic ACR code generators. During the execution of the program, the monitoring extracts the observation of each cell to decide which alternative should be used. The dynamic runtime uses this information to build a kernel with low control overhead following the same execution order as the source program (Figure 2b). The static version uses a guard statement and data collected during execution to switch between the different alternatives for each cell. Hence, the execution order is modified and a cell is visited entirely before continuing to the next one (Figure 2a).

ACR's dynamic runtime uses a set of threads, each with a specific role. It consists of: a monitor thread that is responsible to extract for each cell, from raw data, the alternative to use. A polyhedral code generation thread that utilizes the information from the monitor thread to generate an optimized kernel with low control overhead. 1 A compilation thread that takes the output of the code generation stage and compiles it to a binary format. Finally, a coordinator thread which is responsible to organize the other threads, schedule their execution and inject the optimized version to be used by the application once it is ready.

APPROXIMATION EXTRACTION

Using approximate computing either requires deep knowledge of the user in this field or, if automated, ask the user to delegate the transformation to the compiler. Each of these two techniques have their pros and cons. An expert developer will select the best approximation technique suitable for his problem. On the other hand, an automatic method may discover potential approximation targets using binary error injection, using genetic/machine learning algorithms or using pattern matching to find a suitable approximation strategy to be applied to the original algorithm [START_REF] Ansel | Language and compiler support for auto-tuning variable-accuracy algorithms[END_REF][START_REF] Vinay K Chippa | Analysis and characterization of inherent application resilience for approximate computing[END_REF][START_REF] Samadi | Paraprox: Pattern-based Approximation for Data Parallel Applications[END_REF]. The goal is to lower the program's computation footprint while keeping the program's output deviation below a user-defined threshold. The deviation is defined as the l ∞ norm of the difference between two states corresponding to the same data, i.e., max(|state istate ′ i |). In this section we provide a new method that automatically extracts adaptive approximation information from a targeted portion of a program. Our method also allows interaction between the user and the compiler to select the best alternatives. Automatic stencil approximation (ASA) relies on the polyhedral model to analyse the code. The polyhedral model abstracts a set of statement having static control (SCoP) [START_REF] Bastoul | Putting Polyhedral Loop Transformations to Work[END_REF]. This includes the class of program with loops, conditional statements and data accesses which depend on affine expression of constant parameters and outer-loop variables.

i j i 1 j 1 i 2 j 2 i 3 j 3 i 4 j 4 i j i j Alt A 0 Alt A 1 Cell contour for (i : {i1, i3}) for (j : {j1, j3}) if (cellUseA 0 (i , j )) for (i ′ : {i, i + 1}) for (j ′ : {j, j + 1}) A 0 else for (i ′ : {i, i + 1}) for (j ′ : {j, j + 1}) A 1 (a) Static for (i : {i1, i2}) for (j : {j1, j4}) A 0 for (i : {i3, i4}) for (j : {j1, j2}) A 0 for (j : {j3, j4}) A 1 (b) Dynamic

Precision Level Discovery

The first semantic information required by the compiler to use approximate techniques is related to the alternative selection mechanism, i.e., in which conditions do we allow a less precise versions of the algorithm to be used. We suppose that the approximate alternative code versions can be ordered according to their level of approximation, e.g.,

V 1 ≺ V 2 ≺ V 3 means that V 1 is more precise than V 2 that
is more precise than V 3 , representing a total order. There will be two or more versions, the original code represented as V 1 and its alternatives (Section 4.2). In this section we present an algorithm to find a function to select the alternative versions at runtime. The input data of the decision function is the application's runtime data, representing the state of the problem that the user has to solve. The compiler has to detect which data has the convenient characteristics to be able to generate a function to select an alternative globally or on a localized portion of the application's data domain. We propose a set of characteristics extracted from our study of typical applications from the field of approximate computing.

The first characteristic which usually provides relevant precision information from the application data is the deviation between the updates. Mathematically, this characteristic represents the first derivative of a function [START_REF] Courant | Introduction to calculus and analysis I[END_REF]. Whenever this rate is measured to be low, it is a good sign that we could use a less precise update function in this area of the application's data domain. Furthermore, this effect is usually local, and in numerical simulations, singularities spread and do not appear in the middle of uniform regions, hence, missing a significant event is unlikely. If the effect is localized, we consider clusters of data as a single entity and apply the alternative to each updates of this cell (more in Section 4.3).

To discover this first characteristic the compiler has to look for the following patterns in the source code:

Swapping arrays They are frequently used in physic simulation and image processing where the transformation uses the data from the previous transformation as input. This characteristic is statically analyzed if the program is simple enough or dynamically asserted by running the application on datasets and monitoring the array addresses accessed and recognize a "n cycle" pattern. Write-only arrays They are usually the sign of a multi-stage update algorithm, e.g., any advection-diffusion equation which updates the physical variables in multiple steps/functions or store statistic information of the application's state. This is addressed with a static analysis of the kernel on written-only data arrays.

A second characteristic, lie in the application's own statistic gathering, e.g., multimedia codec gathers information about previous frames to achieve a better compression ratio, and numerical simulations may use a physical variable or error rate as a simulation stopping condition. Statistics are usually gathered globally and have smaller size than the application's dataset. The application may gather the statistics in a dedicated function which could be marked by the user in the same fashion as the kernel to optimize. Hence, the compiler can extract correlations between the data updated by the kernel and the statistics gathered by the application. In this paper we only consider statistic update located inside the kernel to optimize. Our statistics gathering follows the following pattern:

Accumulators They can be used to store statistic information and are most likely used by the application to take algorithmic decisions. Such accumulators are of particular interest for our method as they may store the information we are interested in. Accumulator are many-to-one information gathering.

The compiler may allocate memory to save this information at a finer granularity if needed.

With these information, the compiler builds the most suitable precision discovery function by monitoring the output of the different functions at runtime. It may ask the user for advice depending on the annotation the user selected, i.e., either interactive or automatic.

In addition to explicitly selecting the target computation kernel, the user inputs the level of deviation under which it becomes acceptable to use approximation, whereas the compiler generates N alternatives. In that case, the compiler needs to affect a deviation level, related to the alternative approximation level, for which each alternative will be selected. Assuming that lower deviation allows for even more aggressive approximation, we can affect decreasing deviation values to more approximate alternatives. The alternative deviation values may be generated following a negative linear or exponential slope. We rely on an exponential slope to increase the density of alternative in the neighborhood of the user's defined deviation threshold. Data:

Input array I [M][N ] Result: Output array O[M -2][N -2] for i ← 1 to M -1 do for j ← 1 to N -1 do O[i][j] ← I [i][j]⊕I [i-1][j]⊕I [i+1][j]⊕I [i][j-1]⊕I [i][j+1];
2 def Five points stencil with two iterations merged:

/* operation • distributive over ⊕ */ Data: Input array I [M][N ] Result: Output array O[M -2][N -2] for i ← 1 to M -1 by 2 do for j ← 1 to N -1 do O[i][j] ← ( 5 • I [i][j] ) ⊕ 2 • ( I [i + 1][j] ⊕ I [i][j -1] ⊕ I [i][j + 1] ⊕ I [i -1][j -1] ⊕ I [i + 1][j + 1] ⊕ I [i -1][j + 1] ⊕ I [i -1][j] ⊕ I [i + 1][j -1] ) ⊕ I [i -2][j] ⊕ I [i + 2][j] ⊕ I [i][j -2] ⊕ I [i][j + 2];

Automatic Alternative Stencil Generation

Alternatives are different versions of the application's kernel which uses approximation techniques to reduce the overall kernel calculation and data access intensity. Such versions suffer a precision loss compared to the original algorithm. In this section we present a technique to identify pertinent approximate computing transformations, generate the alternative versions and classify them by level of approximation.

Stencil code are a class of iterative kernel that uses nearestneighbor computations on a grid or graph data structure, i.e., the update of a point on the grid depends on the adjacent neighbor values according to a fixed pattern. Stencils are extremely commonly found in scientific computing, machine learning and image processing applications. E.g., Algorithm 1 shows a five point stencil kernel on a 2D Cartesian grid. The compiler can rely on a static data access analysis in order to automatically detect the most common form of stencils. In our case, we are interested in data read and write accesses for each statement of the kernel to detect stencil-like computations. We provide Algorithm 2 to identify the stencils. This algorithm builds for each written data at a given iteration, the set of read data that contributes to the computation of the written data (both inside a statement or transitively through several statements). It considers each write statement in the lexicographical order and selects the one with the most reads to the same array. It outputs the set of written + contributing data, with maximum cardinality.

We propose the following adaptive stencil alternatives which reduces the stencil computation and data accesses:

Skipping Skipping some stencil computation entirely (e.g. skip it every N iteration of the computing loop) when the precision function has monitored a low deviation. Narrowing Reduce the size of a stencil to reduce both its calculation and memory load. Border Activation Only activate the stencil computation at the cell borders, skipping the computation in the center of cells while allowing information to pass through the interfaces with their neighbors. Cells having a neighbor with an alternative more precise than theirs will use border activation to capture a possible incoming singularity.

Stencil narrowing may rely on the fusion of multiple stencil steps to grow the size of the halo. For example, the stencil in Algorithm 1 top has a halo of distance one. It is possible to increase the size of this stencil halo by merging two stencil loops, resulting in the bottom stencil in Algorithm 1. A larger halo offers more possibilities to build an approximate version. The halo could be reduced by redistributing the stencil weights of the outside most point of the halo into the remaining weights (narrowing). Another method consists on keeping only one stencil update every N steps (skipping).

The analysis of deviation between updates allows adaptive methods to target the expensive computations in part of the domain where significant updates are expected [START_REF] Marsha | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Huang | Adaptive moving mesh methods[END_REF]. Hence, our algorithm relies on this feature to apply the previously mentioned alternatives.

Stencil skipping allows the compiler to generate a version that reduces the computation intensity, hence the precision, by a factor of N . Stencil narrowing allows to generate a wider range of approximate versions, not available with skipping only. For example, it is possible to remove P points of a stencil with the lowest absolute factors to lower the computation/precision. If the number of points removed is below what stencil skipping deletes from a level N to N -1, it corresponds to fully generated intermediate approximate versions. Stencil merging combined with these techniques allows to generate enough alternatives for an adaptive purpose. In the current state of our implementation choosing a high number of alternative leads to a high overhead, consequently generating only a limited number of versions gives the best results.

Granularity Selection

At this point, the compiler has decided upon one precision feature and multiple associated alternatives. It has yet to decide at which granularity to take the decision. A fine grid will have a better mapping to the application's problem but will have a higher managing overhead whereas a grid that is too coarse may miss approximation opportunities. In this section we present an algorithm to select a pertinent granularity based on an empirical study of a set of applications representative of the approximate computing field.

State of the art adaptive techniques use modular grids, i.e., the grid shape is updated during execution to focus the computation where it is needed. A modular grid creates special cases at the interface of two cell with different precisions. These cases need to be handled by the compiler and complicates the automatic alternative and grid generation. In this work we only consider Cartesian, static grid. We discovered that the application's behaviour to the grid size can be transposed between distinct input problems and data sizes.

Figure 4 shows the three main recognizable patterns seen while monitoring different applications. All of them show a local minimum but with a different distribution of performance. Applications that show one of these three behaviours will present the same behaviour with different data sizes or problem statements (see Section 5). Hence, we propose to match these patterns to categorize the application among one of these three forms and compute the ratio between the maximum grid size and the value read at the local minimum. Our empirical study suggests that this ratio can be reused with different data sizes to scale the grid size accordingly. Figure 3 shows that such algorithm is precise enough to approach the best value for bigger dataset based on data collection on the small dataset. Hence, the compiler has an estimate of the best grid size for any dataset size and the curve indicating the slope direction.

With the current static-grid back end, we cannot guarantee the precision level of the output with respect to any input data. The guarantee we provide is that no approximated version is used when the monitored deviation is greater than the one selected by the user. The user must provide input data to the optimized application within the same granularity scaling than the one used during profiling. The user should run a new automatic profiling step for different input value scales. Ongoing work aims at removing this requirement using dynamic grids and mathematical transformations.

User-Compiler Interaction

Our ASA compiler takes as input a file containing the source code with a kernel annotated with an ASA annotation presented in Section 3.1. Its output is a new file where the kernel has been annotated with the relevant ACR pragmas with relevant parameters, and where new functions corresponding to new code alternatives have been added if necessary. A user wishing to interact with the compiler, for any of the three features presented in the sections 4.1 4.2 and 4.3, can investigate the generated file to see which approximations have been built. The compiler accepts user ACR annotations in addition to the ASA annotation and they always prevails over the compiler decisions. This allows an advanced user to provide specific knowledge to the compiler. The user can also access the data generated by the compiler, as the monitoring information from the applications (see Figure 4) and output deviation. The interaction with the optimizing compiler allows the user to have a high-level feedback about the monitoring data, the approximation strategy and the adaptive grid size.

EXPERIMENTAL STUDY

In this section we present the experimental results achieved by our method on application compute intensive kernel on the set of benchmarks used to evaluate ACR [START_REF] Schmitt | Adaptive Code Refinement: A Compiler Technique and Extensions to Generate Self-Tuning Applications[END_REF] with the addition of the heat equation simulation. This benchmark set includes typical codes that are resilient to approximation. We used a source-to-source approach and targeted the annotations presented in section 3. Because of the runtime overhead of managing many alternatives, we restricted our study to the generation of three alternatives (precise, approximated and border activation). This comes without loss of generality, but we may expect even better results using future ACR compiler versions with smaller runtime overhead. The experimental setup consists of an Intel Xeon E5-2620 with 16 gigabytes of ram. We use the ACR development version2 with GCC 8.2.

Heat Equation

This application, introduced in Section 2, is a solver for the steady state heat equation. It is used for material applications and physics simulations [START_REF] Geankoplis | Transport processes and unit operations[END_REF]. The main computational kernel of this application is shown in Figure 1, left. This application is a good candidate for optimization using approximate computing techniques. By profiling the application, a developer may know that the kernel is the biggest contention part of the program. Hence, this developer may add the annotation present at line five: #pragma asa auto (1e-4). This annotation states that the following kernel can be optimized by the compiler to take advantage of adaptive techniques. In the following of this section, we explain the compiler process to extract the information necessary to generate the approximate computing version.

Heat Equation Precision Level Discovery The compiler generates then executes an instrumented version of the program on a small user-provided dataset to gather information about the array accesses at runtime. In this program, it happens that the kernel's written array switches between two addresses. 3 The compiler can automatically generate the comparison function as there are two arrays to compare. The developer already expressed it's intention to allow approximation when deviation is below 10 -4 . Therefore, the compiler generates the following ACR annotation to select the precision level: #pragma acr monitor(tNext[ i ][ j ], min, diffT ) Where the min function selects the inside cell alternative of the most precise version and diffT is the pre-processing function that embodies the computation of the deviation between the value of arrays tNext and t. The monitor compares the deviation to 10 -4 and returns zero if the original code has to be executed or one if approximation is acceptable, allowing for one alternative along the original code.

Heat Equation Alternative Generation Algorithm 2 detects the stencil pattern. This stencil being too narrow for further optimization, the compiler modifies the kernel to merge two consecutive stencil computations to increase the width of the stencil, creating a loop where the iteration count is defined by a new variable kernel_iter (the value 2 corresponds to the original kernel). A stencil skipping alternative can be created by modifying the number of iterations of the kernel: #pragma acr alternative alt0 (parameter, kernel_iter =2) #pragma acr alternative alt1 (parameter, kernel_iter =1) ACR requires additional annotations to link the alternatives to the value obtained by the monitoring. Therefore, the following annotations link the monitoring folded value zero to the original computation and the monitoring folded value one to the early-terminated number of kernel iterations: #pragma acr strategy dynamic(0, alt0 ) #pragma acr strategy dynamic (1, alt1 ) Heat Equation Grid Selection The compiler tests the generated alternative on a small user-provided dataset to assess the viability of this alternative. The compiler requires a metric to evaluate the level of deviation of the result of the alternative. The preferred metric for ACR is the relative difference or error with respect to the original application's output [START_REF] Törnvist | How Should Relative Changes Be Measured?[END_REF]. The compiler uses this metric to compare the result of the alternative version against the original application and requires that the difference be less than five percent in total. The application is evaluated w.r.t. multiple grid sizes to obtain the data plotted in Figure 4a. Using this curve, it is possible to assess the ratio corresponding to the local minimum that will be used for other datasets, here 24 500 which will be multiplied to the maximum between sizeX and sizeY to compute the new grid size.

Heat Equation Automatic Method Evaluation Table 1 shows the speedup obtained with our algorithm and the best value obtained by a search over all possible grid values. Results show both the benefits to exploit adaptive technique for this benchmark and that the automated process achieves a performance close to the best hand tuned version (see Section 5.6).

Eulerian Fluid Simulation

Figure 5 shows a kernel presented by Stam [42] which solves the Navier-Stokes equation using an iterative solver. This solver must be fast and the precision accurate enough to be human imperceptible as it is targeted for real time video applications or games. for ( unsigned j = 1; j < N -1; ++ j ) Fluid Simulation Precision Level Discovery In this kernel, the compiler has assessed that the arrays x and x0 are swapped between each call to the kernel. Hence, it generates a function that will compare the two arrays to compute the update function derivative: #pragma acr monitor(x[ i ][ j ], min, diffX )

5 x [ i ][ j ] = 6 ( x0 [ i ][ j ] + 7 a * ( x [i -1][ j ] + 8 x [ i +1][ j ] + 9 x [ i ][ j -1] + 10 x [ i ][ j +1]) 11 ) / c ;
Fluid Simulation Alternative Generation The compiler recognizes a stencil in the two inner dimensions of the kernel. This stencil repeats P times during the execution. Therefore, the compiler does not have to widen the stencil size but only has to reduce the number of stencil computation to generate an alternative: #pragma acr alternative alt0 (parameter, k=P) #pragma acr alternative alt1 (parameter, k=1) #pragma acr strategy dynamic(0, alt0 ) #pragma acr strategy dynamic (1, alt1 ) Fluid Simulation Grid Selection Figure 4b shows the application's behaviour to the grid size. The best ratio extracted from the curve is 60 500 × max(M, N ).

Cellular Automaton

Game of life is cellular automaton which is ruled by three easy rules:

• A cell becomes alive if there is at least three alive cells in its direct surrounding at the previous time step. • A cell stays alive if it had two or three alive neighbours at the previous time step. • Otherwise a cell is considered dead at the next step. This automaton has recurring patterns and empty areas that creates potential for optimization [START_REF] Gosper | Exploiting regularities in large cellular spaces[END_REF]. It is possible to create an adaptive version that does little computations in the empty zones. The main kernel of our benchmark implementation is shown in Figure 6.

Cellular Automaton Precision Level Discovery The compiler detects that the written array is switching between two addresses and will use it to compute the derivative. #pragma acr monitor( current_grid [ i ][ j ], min, diff ) Cellular Automaton Alternative Generation To generate a non-approximate stencil version which allows the simulation cell states to spread while disabling the empty regions, the user decides to add the interface_compute alternative which activates the neighboring ACR cells in the locations where the simulation cells are active. The user can help the compiler to select an alternative by adding a back-end annotation, here an ACR alternative shown in Figure 6 lines two and three. The compiler can then generate the link between the monitoring and the added alternative: #pragma acr alternative alt1 (zero_compute) #pragma acr strategy dynamic(2, alt1 ) #pragma acr strategy dynamic (1, alt_interface ) Cellular Automaton Grid Selection This application follows the same pattern as the heat solver Figure 4a. The computed best grid ratio over data size is 275 8000 × max(nb_r ow , nb_col ).

Finite Difference Time Domain (FDTD)

Finite-difference time-domain implements a solver for the Maxwell-Faraday equations of electrodynamics differential equations [START_REF] Ardavan | MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method[END_REF]. The kernel Figure 7 shows the update of the value of the electric plane E as a function of the magnetic field H . This application's singularity resides in its two-phases algorithm, the electric field update followed by the magnetic field update.

FDTD Precision Level Discovery This kernel updates two arrays at the same time, Ex and Ey. The array pointers are not modified during the execution of the program. The compiler arbitrarily chose to monitor Ex (in such case, it is likely that the two updates are correlated, minimizing the impact of that choice: there is actually a correlation which can be statically analysed in this benchmark). The ACR annotation generated is: #pragma acr monitor(Ex[i ][ j ], min, diff ) FDTD Alternative Generation For this application, the stencil merging technique could not help because the arrays are not swapped and the array Hz is not written inside the kernel. Hence, the only remaining possible alternative optimizations are stencil skipping and narrowing: #pragma acr alternative alt1 (zero_compute) #pragma acr strategy dynamic (1, alt1 ) FDTD Grid Selection Figure 4c shows the application's runtime with respect to the grid size. For this application, the maximum grid size is 3000 and the best grid size is located in the plateau between 200 and 300 with 240 3000 × max(I, J )

K-means clustering

K-means is used for data characterization, e.g., image processing, and machine learning [START_REF] Meng | Best-effort parallel execution framework for recognition and mining applications[END_REF]. Its purpose is to place N observations into B buckets where the observations being the closest will appear in the same bucket or centroid. The main kernel of the application is shown in Figure 8. Firstly, the algorithm places the observations into the closest bucket using a comparison function. Secondly, the barycenter of the bucket is updated with the observations belonging to it. The algorithm caries on until all the observations settles.

K-means Precision Level Discovery

The only available array to monitor is point_centr_map: #pragma acr monitor(point_centr_map[pos], min, diff ) K-means Alternative Generation Having access to the whole loop, the compiler uses the stencil merging technique to reduce the iteration count by half in the approximated regions by applying stencil skipping: #pragma acr alternative alt0 (parameter, kernel_iter =2) #pragma acr alternative alt1 (parameter, kernel_iter =1) #pragma acr strategy dynamic(0, alt0 ) #pragma acr strategy dynamic (1, alt1 ) K-means Grid Selection We used the K-mean algorithm to characterize images of different sizes, ranging from 640 × 425 to 4288 × 2848. The best ratio extracted for this application is 90 272000 × points.

Performance Evaluation

In this section we compare the application's performance against the original code version and a hand-tuned version using the ACR annotations. Table 1 shows the performance and deviation of the output results of our method against user search of the parameters. We can see that, while the performance of AMA is lower than what a specialist can extract with the ACR annotations, our technique achieves sensible performance gains often close to hand-tuned versions. The deviation of the output is similar and sometimes better with with the hand optimized versions and well below five percent, besides one over-optimization with K-mean where the automatic method reached nine percent deviation while the hand tuned version was selected with a bigger grid size which allowed for less approximation.

ACR uses a multi-threaded infrastructure to generate and compile the adaptive versions in parallel to the program execution. This requires resources which may not be exploited due to Amdahl's law [START_REF] Gustafson | Reevaluating Amdahl's Law[END_REF], or be better exploited with AMA. Furthermore, adaptive methods are complementary to parallelization techniques, and information from the adaptive grid may be used for a better resource scattering.

STATE OF THE ART

Relaxing Semantics Relaxing the program's semantics provides new optimization potential for the compiler. Precimonious is a tuning assistant which uses this strategy to selectively lower the precision of floating point arithmetics (e.g., double to float) to optimize memory traffic and computation time [START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF]. HELIX-UP chooses to ignore dependencies and synchronizations to increase parallelism in applications [START_REF] Campanoni | HELIX-UP: Relaxing Program Semantics to Unleash Parallelization[END_REF]. SAGE and Meng's et al. work allows for better parallelism on GPU architectures by targeting GPU specific approximate optimizations and using a best effort computing model respectively [START_REF] Meng | Best-effort parallel execution framework for recognition and mining applications[END_REF][START_REF] Samadi | SAGE: Self-tuning Approximation for Graphics Engines[END_REF]. Loop perforation technique allows the compiler to skip loop iterations to gain in performance [START_REF] Sidiroglou-Douskos | Managing performance vs. accuracy trade-offs with loop perforation[END_REF]. Task-based workflows can be approximated by skipping certain task probabilistically [START_REF] Rinard | Probabilistic Accuracy Bounds for Fault-tolerant Computations That Discard Tasks[END_REF]. Rhahimi et al. explored memoization techniques for error-tolerant applications [START_REF] Michie | Memo" functions and machine learning[END_REF][START_REF] Rahimi | Spatial memoization: Concurrent instruction reuse to correct timing errors in simd architectures[END_REF]. Chippa et al. propose an automatic method to characterize the resilience of application to approximation [START_REF] Vinay K Chippa | Analysis and characterization of inherent application resilience for approximate computing[END_REF] a neural network based optimizer that runs on a neural accelerator [START_REF] Hadi Esmaeilzadeh | Neural Acceleration for General-Purpose Approximate Programs[END_REF]. Specialized approximate hardware instructions and memory location can complement software methods to achieve efficient approximate computing [START_REF] Vinay | Scalable Effort Hardware Design: Exploiting Algorithmic Resilience for Energy Efficiency[END_REF][START_REF] Hadi Esmaeilzadeh | Architecture Support for Disciplined Approximate Programming[END_REF][START_REF] Fang | SoftPCM: Enhancing Energy Efficiency and Lifetime of Phase Change Memory in Video Applications via Approximate Write[END_REF][START_REF] San | Load Value Approximation[END_REF][START_REF] Misailovic | Chisel: Reliability-and Accuracy-aware Optimization of Approximate Computational Kernels[END_REF][START_REF] Sampson | Approximate Storage in Solid-State Memories[END_REF].

Compared to existing techniques, our solution monitors the data at runtime and takes local choices. Hence, while previous methods used the same alternative on the whole application's data domain, ours makes sure the approximate versions are active on portions of the computation space where they may be beneficial.

Automatic-Tuning Automatic tuning enables the compiler at automatically searching for the best parameters that allow for approximate computing to be efficient. PetaBricks is a programming language and compiler which uses empirical approach to select the best kernel among a pre-set of algorithms [START_REF] Ansel | PetaBricks: A Language and Compiler for Algorithmic Choice[END_REF]. PetaBricks's compiler has been augmented to automatically generate approximate versions using a genetic algorithm which targets a minimum deviation of the output [START_REF] Ansel | Language and compiler support for auto-tuning variable-accuracy algorithms[END_REF][START_REF] Ding | Autotuning algorithmic choice for input sensitivity[END_REF]. Sculptor has extended loop perforation to target specific instruction and apply perforation at a finer grain with feedback [START_REF] Li | Sculptor: Flexible Approximation with Selective Dynamic Loop Perforation[END_REF]. Paraprox uses pattern recognition to generate approximate versions of kernels [START_REF] Samadi | Paraprox: Pattern-based Approximation for Data Parallel Applications[END_REF]. ASAC constructs a model to minimize a n-dimensional problem of selecting the values for approximate variables while keeping an acceptable quality of the result [33]. ACCEPT proposes a programming model that includes the approximate annotation in the language's grammar [START_REF] Sampson | Accept: A programmer-guided compiler framework for practical approximate computing[END_REF]. Green provides a set of code annotations that allow developers to use a range of approximation strategies with bounded error settings [START_REF] Baek | Green: A Framework for Supporting Energy-conscious Programming Using Controlled Approximation[END_REF].

Our technique relies on the same information extraction methods as some of the previously mentioned works. We improved the previous works by presenting a novel high-level approximation optimization applied to stencil-based computation which generates a set of approximated versions with various performance/precision tradeoffs. Our solution is also the first one, to the best of our knowledge, to automatically apply adaptive approximation techniques to a compute intensive kernel.

Adaptive Methods Adaptive mesh refinement is a mathematical method which can be applied to solve partial differential equations, mainly for physics applications [START_REF] Marsha | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Marsha | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF]. Historically, the adaptive mesh method was tightly coupled with the application that used the method, but more recently packages such as PARAMESH [START_REF] Macneice | PARAMESH: A parallel adaptive mesh refinement community toolkit[END_REF], HAMR [START_REF] Neeman | HAMR: The Hierarchical Adaptive Mesh Refinement System[END_REF], DAGH [START_REF] Parashar | Systems Engineering for High Performance Computing Software: The HDDA/DAGH Infrastructure for Implementation of Parallel Structured Adaptive Mesh[END_REF] and AMRClaw [START_REF] Berger | Adaptive Mesh Refinement using Wave-Propagation Algorithms for Hyperbolic Systems[END_REF][START_REF] Kyle T Mandli | Clawpack: building an open source ecosystem for solving hyperbolic PDEs[END_REF] allows developers to build new application using efficient adaptive mesh refinement implementations. Porting existing applications to these packages is time consuming and error prone. ACR provides a set of easy to use, optional annotations that can be positioned in front of compute expensive kernel [START_REF] Schmitt | Adaptive Code Refinement: A Compiler Technique and Extensions to Generate Self-Tuning Applications[END_REF][START_REF] Schmitt | Semi-Automatic Generation of Adaptive Codes[END_REF]. The annotations instruct the compiler to generate the adaptive version of the code, however they still require a significant domain-specific expertise from the programmer to set them in a pertinent way.

In our work, we propose a set of techniques to automatically set ACR annotations with convenient parameters, which builds the first automatic compiler optimization technique based on adaptive approximation methods.

CONCLUSION

With one simple annotation, our method allows for the first time developers to delegate to the compiler the time consuming task of optimizing an application's kernel using adaptive approximation techniques. We presented three important features that the compiler has to extract from the source code to generate a pertinent adaptive version. Precision level discovery, allows the compiler to choose a precision level from the application's dynamic data. This information is critical to find where the application needs the most precision and where more or less aggressive approximations can be used. The alternative generation requires the compiler to determine what kind of computation is performed by the kernel to apply appropriate code transformations. We used pattern matching to target stencil computation and proposed to merge multiple stencil steps to unleash more optimization opportunities. The third information, the granularity, looks for the best performance/precision tradeoff. With a grid too thin, the adaptive overhead will be too high and with a grid too coarse the perturbation will not be captured by the adaptive grid.

We evaluated our method on a set of representative applications resilient to inexact computations. We showed that we can extract enough information from these kernels and their profiling to generate an adaptive optimized version automatically. We also showed that the compiler can be instructed to use a specific alternative when the user sees it fit. Finally, we provided experimental evidence that the generated adaptive versions perform better than the original version while maintaining a good quality of the result.

Ongoing work includes investigation of adaptive techniques on more computational patterns. The automatic alternative generation could, in the future, use several alternatives simultaneously and evaluate the effect that an alternative has on the application's quality of the output, to e.g., statistically provide an upper bound on the deviation of the result.

5 #

 5 pragma asa auto (1 e -4) 6

Figure 2 :

 2 Figure2: Example of static and dynamic code generated on a four-cell example. Three of these cells require the alternative A 0 and one requires the alternative A 1 . The static version, which is generated at compile time, visits a cell entirely before going to another. A guard is needed to select the alternative before entering a cell. The dynamic version generates a specific code at runtime, which follows the application original schedule shown on the top right.

Algorithm 1 :

 1 Example 2D stencil algorithm 1 def Five points stencil:

Figure 3 :

 3 Figure 3: Heat solver time in function of the grid size, execution with multiple dataset sizes. The grid size that performs the best can be interpolated using a linear regression.

Figure 4 :

 4 Figure4: Running time of three optimized applications as a function of the grid size. Each application is representative of a class of programs sharing a similar behaviour when increasing the grid size. In a given class, we could observe that the behaviour is the same independently of the dataset size or initial conditions of the simulation.

1 #

 1 pragma asa auto (1 e -2) 2 for ( unsigned k = 0 ; k < P ; ++ k ) { 3 for ( unsigned i = 1; i < M -1; ++ i )

4

 4

13 }Figure 5 :

 135 Figure 5: A Gauss-Seidel iterative solver of linear equations.

1 # 25 }Figure 6 :

 1256 Figure 6: Game of Life algorithm to update the state of the cell grid from one generation to the next one.

1 # 5 ( 7 ( 8 } 9 }Figure 7 :

 157897 Figure 7: Kernel of a FDTD update loop of the electric and magnetic field in a 2D space.

3 # 5 unsigned 23 }

 3523 pragma asa auto (0) 4 for ( size_t pos = 0; pos < points ; ++ pos ) { while (! has_converged );

Figure 8 :

 8 Figure 8: K-Means core algorithm where the observations points are placed into the clusters and the center of the cluster is updated.

  Algorithm 2: Search the biggest stencil in a statement set Data: Ordered statement set S Data: For each S i ∈ S, read access set S read For each S i ∈ S, write access set S write Biддest read ) 19 statementWrite(S, N): returns the set of statement in S which writes at the location N. 20 groupeBy(P, S): returns a set of set grouped by P 21 maxCardinal(S): returns the largest cardinal set in S

			i
		Result: Tuple of written and read location of a stencil
	1 Biggest read ← ∅;
	2 Biggest write ← ∅;
	3 n ← ∥S ∥;
	4 while n > 0 do
		/* Statement accesses + transitive accesses	*/
	5	N read ← S read n	;
	6	repeat
	7	N read prev ← N read ;
	8	foreach N i ∈ N read prev do
	9	foreach W i ∈ statementWrite(S, N i ) do
	10	N read ← N read ∪ W read i	;
	11	until N read = N read prev ;
		/* Check for stencil pattern for each arrays	*/

i Data: 12 G ← groupBy(array, N read ); 13 BestLocalMatch ← maxCardinal(G); 14 if ∥Biддest read ∥ < ∥BestLocalMatch∥ then 15 Biддest read ← BestLocalMatch; 16 Biддest write ← S write n ; 17 n ← n -1; 18 return (Biддest write ,

Table 1 :

 1 . Approximate computing application can also take advantage of specialized hardware. Esmaeilzadeh et al. used AMA and hand tuned ACR version performance and deviation of the output results.

	Application	ASA Speedup /	Hand ACR Speedup
		Output Deviation Output Deviation
	Heat	1.25	Training set
	Heat 1000	1.68 / < 0.001%	1.76 / < 0.001%
	Heat 2000	1.81 / < 0.001%	2.03 / < 0.001%
	Fluid	1.1 / 0%	Training set
	Fluid flame	1.17 / 0.8%	1.27 / 1.51%
	Fluid vortex	1.16 / 0.01%	1.19 / 0.01%
	FDTD	1.12 / 0%	Training set
	FDTD 500	1.16 / 0%	1.22 / 0%
	FDTD 3000	1.24 / 0.1%	1.4 / 0.5%
	GOL	1.44 / 0%	Training set
	GOL big	2.09 / 0%	2.2 / 0%
	K-means	2.14 / 3.58%	Training set
	K-means med	1.51 / 0.7%	2.50 / 2.62%
	K-means big	2.18 / 9.03%	1.57 / 4.48%

Note that this does not involve, but is complementary to, polyhedral optimizations techniques such as parallelization or data locality optimization.

http://gauvain.u-strasbg.fr/%7Eschmitt/acr

The corresponding statement is flagged using a specific annotation in our compiler prototype, but it may be discovered automatically as this is a well known pattern in such codes.


E Falgarone

M A Zwaan

B Godard

E Bergin

R J Ivison

P M Andreani

F Bournaud

R S Bussmann

D Elbaz

A Omont

I Oteo

F Walter

Large turbulent reservoirs of

Starburst galaxies at the peak of cosmic star formation 1 are among the most extreme starforming engines in the universe, producing stars over ∼ 100 Myr 2 . The star formation rates of these galaxies, which exceed 100 M per year, require large reservoirs of cold molecular gas [START_REF] Carilli | Cool Gas in High-Redshift Galaxies[END_REF] to be delivered to their cores, despite strong feedback from stars or active galactic nuclei [START_REF] Veilleux | Galactic winds[END_REF][START_REF] Schreiber | The SINS/zC-SINF survey of z ∼2 galaxy kinematics: evidence for powerful active galactic nucleus-driven nuclear outflows in massive star-forming galaxies[END_REF] . Starburst galaxies are therefore ideal targets to unravel the critical interplay between this feedback and the growth of a galaxy [START_REF] Somerville | Physical Models of Galaxy Formation in a Cosmological Framework[END_REF] . The methylidyne cation, CH + , is a most useful molecule for such studies because it cannot form in cold gas without supra-thermal energy input, so its presence highlights dissipation of mechanical energy [START_REF] Flower | C-type shocks in the interstellar medium: profiles of CH + and CH absorption lines[END_REF][START_REF] Lesaffre | Low-velocity shocks: signatures of turbulent dissipation in diffuse irradiated gas[END_REF][START_REF] Godard | Chemical probes of turbulence in the diffuse medium: the TDR model[END_REF] or strong UV irradiation [START_REF] Godard | A complete model of CH + rotational excitation including radiative and chemical pumping processes[END_REF][START_REF] Falgarone | Strong CH + J = 1-0 emission and absorption in DR21[END_REF] . Here, we report the detection of CH + (J =1-0) emission and absorption lines in the spectra of six lensed starburst galaxies [START_REF] Eales | The Herschel ATLAS[END_REF][START_REF] Negrello | The Detection of a Population of Submillimeter-Bright, Strongly Lensed Galaxies[END_REF][START_REF] Bussmann | Gravitational Lens Models Based on Submillimeter Array Imaging of Herschel-selected Strongly Lensed Sub-millimeter Galaxies at z > 1.5[END_REF][START_REF] Swinbank | Intense star formation within resolved compact regions in a galaxy at z = 2.3[END_REF] at redshifts z ∼ 2.5. This line has such a high critical density for excitation that it is emitted only in very dense (> 10 5 cm -3 ) gas, and is absorbed in low-density gas. We find that the CH + emission lines, which are broader than 1000 km s -1 , originate in dense shock waves powered by hot galactic winds. The CH + absorption lines reveal highly turbulent reservoirs of cool (T ∼ 100 K), low-density gas, extending far outside (> 10 kpc) the starburst cores (radii < 1 kpc). We show that the galactic winds sustain turbulence in the 10 kpc-scale environments of the starburst cores, processing these environments into multi-phase, gravitationally bound reservoirs. However, the mass outflow rates are found to be insufficient to balance the star formation rates. Another mass input is therefore required for these reservoirs, which could be provided by on-going mergers [START_REF] Engel | Most submillimeter galaxies are major mergers[END_REF] or cold stream accretion [START_REF] Dekel | Cold streams in early massive hot haloes as the main mode of galaxy formation[END_REF][START_REF] Narayanan | The formation of submillimetre-bright galaxies from gas infall over a billion years[END_REF] . Our results suggest that galactic feedback, coupled jointly to turbulence and gravity, extends the starburst phase instead of quenching it.

Using the Atacama Large Millimeter Array (ALMA), we have detected CH + (J=1-0) lines in high-redshift galaxies: our sample (Table 1) consists of six particularly bright -by virtue of gravitational lensing -submillimeter-selected galaxies (SMGs), with specific star-formation rates > ∼ 3× those of galaxies on the main sequence [START_REF] Schreiber | The Herschel view of the dominant mode of galaxy growth from z = 4 to the present day[END_REF] at z ∼ 2. The SMGs are extremely compact, with half-propagating at v sh ∼ 40 km s -1 in dense pre-shock gas, are those most efficient at producing CH + (see Methods). The CH + column densities of the shocked gas are high enough to emit lines that are much brighter than those detected, allowing for substantial beam dilution. The CH + spectra provide considerable insight into the physics of these extreme star-forming systems. Although the sample size is still small, the concomitance of the broad emission and deep absorption lines of CH + suggests that they trace coeval processes, i.e. the dense shocks with high velocity dispersion and the dilute turbulent reservoirs are tied together. Yet, the disparity between the emission and absorption linewidths is striking. The emission linewidths (1100-1400 km s -1 ) in the three sources, the Cosmic Eyelash, G09v1.40 and SDP17b, where they are well determined, are very similar whereas their star formation rates span a factor of ∼5. This finding suggests that these lines trace the fast thermal expansion of hot winds at similar temperatures T in the three sources. For outflow velocities v out ∼ √ 3c s , where c s ∼ 500 km s -1 ( T 4×10 7 K ) 1/2 is the sound velocity of the hot wind, the range of observed values, v out = 0.7 ∆v em ∼ 700-10 3 km s -1 , provides the range T ∼ 3-6×10 7 K for the hot wind temperature. This high-velocity outflowing hot wind drives large scale shocks in the dense outskirts of the starburst galaxies. Kinetic energy cascades within these shocks down to the velocity of the shocks that are able to form CH + , in a process reminiscent to what is observed in the 35 kpc-long intergalactic shock in the group of galaxies known as Stephan's Quintet [START_REF] Appleton | Shock-enhanced C+ Emission and the Detection of H 2 O from the Stephan's Quintet Group-wide Shock Using Herschel[END_REF] . The outflowing wind eventually escapes from each galaxy, generating turbulence in its environment. It is this transition to turbulence that causes the outflow momentum to change from mostly outward to random and allows the gas to be re-captured gravitationally. This is why the absorption linewidth in the Eyelash (which is 10 × less massive than SDP17b) is smaller than in SDP17b. By preventing a fraction of the outflow from escaping the large-scale galaxy potential well, turbulence makes this fraction available for further star formation [START_REF] Bouché | The impact of cold gas accretion above a mass floor on galaxy scaling relations[END_REF] . Turbulence therefore mitigates the negative feedback of galactic winds on star formation, drawing out the starburst phase. The CH + molecules that are seen in absorption are observed where they form, i.e. in the turbulent dissipation sites of a ∼ 100K dilute molecular phase; this suggests that the hot outflows also induce phase transitions in the galaxy environments that could be the giant haloes of neutral hydrogen at 10 4 K detected in Lyα 29 at high redshift. Importantly, the mass outflow rates that are required to sustain the kinetic luminosity L TR of the turbulent reservoirs closely follow the star formation rates, SFR (see Methods):

Ṁout = 0.05 -0.1 η SFR ( v out 800 km s -1 ) -2 ,
η < 1 being the unknown efficiency of the energy transfer. This is so because the ratios L TR /SFR are observed to be the same within a factor of ∼2 for all the sources (Tables 1 and2). If mass loss occurs at v out ∼ 700-10 3 km s -1 , as suggested by the broad emission lines, then the mass outflow rates are probably lower than the star formation rates. Therefore, given their masses, the turbulent reservoirs that are drained by star formation cannot be fed by outflows only. They must be replenished over the lifetime of the starburst by other inflows, from tidal streams of on-going major mergers or from cold stream accretion. The diffuse gas reservoirs result from the mixing of the gas ejected by the starburst-driven outflows with gas present in the wider galaxy environment. By highlighting turbulent dissipation in cool molecular gas, CH + shows that mechanical energy, fed by gravity and galactic feedback, is largely stored in turbulence and eventually dissipated at low temperature rather than radiated away by hot gas. This turbulent energy storage is observed at two different stages of the feedback process: in the high-velocity dispersion shocked dense gas close to the starburst cores, and very far from them, in the 10kpc-scale massive reservoirs of lowdensity gas. Turbulence appears to be a key process in the triggering and subsequent regulation of star formation.

Figure 1: ALMA rest-frame 360-µm continuum images of our lensed targets. Solid contours are in steps of 3σ × 2 n , n = 0, 1, 2..., and dashed contours at -3σ, where σ is the r.m.s. noise level (Table 1). The synthesised beam is shown in the bottom right corner. In the Planck cosmology [START_REF]Planck 2013 results[END_REF] Figure 2: ALMA CH + J = 1-0 spectra of our lensed targets. Continuum-subtracted CH + J = 1-0 spectra, integrated over the solid angle defined in the caption of Table 1. The pink curves show the final fits; the blue curves show the contributions of the emission (one Gaussian) and absorption (one or two Gaussians). The large velocity overlap of the absorption and emission in NAv1.56 and NAv1.144 makes the fits only indicative (dotted lines). The velocity scale is with respect to the accurate redshifts of Extended Data Table 1. The continuum was removed with a linear baseline to perform the fits.

Table 1: Properties of the lensed SMGs and ALMA CH + line and continuum fluxes. a The IAU full names that provide the galactic coordinates of the sources, as well as the full precision redshifts used are given in Extended Table 1. b Magnification and half-light size inferred from lens modelling of our ALMA data at rest-frame 360 µm (see Methods), except for G09v1.40 and the Cosmic Eyelash where we adopt values from the literature. The uncertainty on the half-light sizes is ∼15%. c Intrinsic L FIR integrated over 40-120 µm [START_REF] Engel | Most submillimeter galaxies are major mergers[END_REF] , unlensed with our lens models. d Star formation rates assuming a Salpeter initial mass function [START_REF] Engel | Most submillimeter galaxies are major mergers[END_REF] and no contribution from active galactic nuclei (AGN). These values are intended to illustrate only the differences among the galaxies. e Observed peak CH + line emission flux integrated over a solid angle defined by the contour level 0.05× the peak continuum level. f Size in the source plane of the brightest star-forming region (source X) [START_REF] Swinbank | The Interstellar Medium in Distant Star-forming Galaxies: Turbulent Pressure, Fragmentation, and Cloud Scaling Relations in a Dense Gas Disk at z = 2.3[END_REF] .

Name a z a µ b log L FIR c r SMG b SFR d S 1.2 mm σ 1.2 mm S CH + e L kpc M
Table 2: Characteristics of the CH + lineprofiles and inferred properties of the turbulent reservoirs a Parameters inferred from fits made simultaneously on the emission and absorption lines: ∆v em and ∆v abs are the FWHM, v off and v em are the offset velocity of the absorption Gaussian and the emission line centroid velocity with respect to the systemic velocity of the galaxy as measured in CO, and τ 0 the absorption peak optical depth (see Methods). b Average CH + column density in a solid angle subtended by the background source of radius r SMG (see Methods). c Radius, mass and kinetic luminosity of the turbulent reservoirs computed for t TR = 50 Myr (see Methods). d In these sources, the emission and absorption velocity coverages broadly overlap and the fits have been performed imposing v em = 0. In NAv1.56, the convergence of the fit is poor. In NAv1.144, the strong residual emission is so narrow that the fit parameters are non-physical and have not been considered.

Name v em a ∆v em a ∆v abs a v off a τ 0 a N (CH + ) b r TR c M TR c L TR c km s -1 km s -1 km s -1 km s -1

Methods

1 -ALMA observations The six targets were observed with ALMA in Cycle 2 as part of project 2013.1.00164.S. The CH + (1-0) line (ν rest = 835.08 GHz) at redshifts z ∼ 2.0-2.6 is shifted into ALMA band 6 at frequencies between 230 and 270 GHz. The ALMA baseline correlator was configured to use three basebands in the upper sideband, so as to cover a maximum contiguous velocity range of approximately 5000 km s -1 centred on the CH + line. The low resolution Time Division Mode was used, which, using spectral averaging, provides 240 channels of each 7.812 MHz per baseband. This corresponds to a spectral resolution of 19 km s -1 after online Hanning smoothing. Between 32 and 34 antennas were used for these observations. A standard observing strategy was used, where a phase reference source was observed every 8 minutes. The total on source time for the six sources combined was 195 minutes. Absolute flux calibration was obtained through the observation of solar system objects and quasars with well-known flux values. The data calibration followed standard routines and was carried out using the CASA 4.4 (common astronomy software applications) package. Continuum images were made for each field using a 'robust' weighting scheme with a Briggs parameter of 0.5. All targets were sufficiently bright to perform one round of phase-only and one round of amplitude and phase self-calibration, which significantly improved the image quality and noise level. These self-calibrated data sets were subsequently used to construct data cubes centred around the redshifted CH + frequency, with a channel spacing of 50 km s -1 . The typical spatial resolution of the final data cubes is 0.5 arcsec, and the r.m.s. noise level is in the range 0.25 to 0.4 5 mJy per channel, depending on the data cube. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00164.S.

-Targets

The Eyelash: This three-image source strongly lensed by a foreground cluster has a combined amplification of 37.5 [START_REF] Swinbank | The Interstellar Medium in Distant Star-forming Galaxies: Turbulent Pressure, Fragmentation, and Cloud Scaling Relations in a Dense Gas Disk at z = 2.3[END_REF] . High-resolution interferometry with the Sub-Millimeter Array (SMA) resolved the rest-frame 260 µm continuum emission in two sets of four clumps.

The clumps have physical (source-plane) sizes of 100-300 pc in a region ∼ 2 kpc within the background galaxy. The total gas mass inferred from CO is M gas ∼ 4 × 10 10 M , the stellar mass is M = 3 × 10 10 M 24 , giving a gas mass fraction f gas = M gas /(M + M gas ) ∼ 60%. The CO spectral line energy distribution suggests that this is a multiple component system with different conditions [START_REF] Danielson | The properties of the interstellar medium within a star-forming galaxy at z=2.3[END_REF] . Low density gas at low temperature is required to model the low-lying transitions of CO. The OH absorption line [START_REF] George | Herschel reveals a molecular outflow in a z = 2.3 ULIRG[END_REF] is interpreted as an outflow originating in the largest star-forming region of the galaxy.

The H-ATLAS sources: The redshift determination in G09v1.40 from the H 2 O lines 32 , z=2.093, is different from the centred value [START_REF] Engel | Most submillimeter galaxies are major mergers[END_REF] of the CH + emission line. 880 µm SMA observation of the dust continuum shows two close components with a separation of ∼ 1 arcsec along the east and west direction. The lensing model [START_REF] Engel | Most submillimeter galaxies are major mergers[END_REF] gives µ = 15.3±3.5. Both the 10m-Keck telescope and the Spitzer space telescope images, in the near-IR and mid-IR respectively, agree about a stellar component magnification of µ = 11.4+0.9/-1 [START_REF] Calanog | Lens Models of Herschel-selected Galaxies from High-resolution Near-IR Observations[END_REF] . For all the others, the results of the lensing model applied to our ALMA data (Table 1) differ by ∼ 20% from those inferred from the SMA data [START_REF] Engel | Most submillimeter galaxies are major mergers[END_REF] , except for SDP130 where the ALMA-inferred magnification is ∼ 4× that from SMA. NAv1.56 is a possible merger at a projected separation of 29 kpc (3.5 arcsec). In SDP17, the total gas mass inferred from CO is M gas = 5.9 × 10 10 M , the stellar mass is M = 2.4 × 10 11 M 25 , so that the gas mass fraction f gas ∼ 20%.

The IAU names of the targets, redshifts and CO and H 2 O lines characteristics are given in Extended Data Table 1.

3 -Exploitation of the CH + lineprofiles CH + has such a highly endoenergic formation route via C + + H 2 (E form ∼ 0.4 eV) that a specific chemistry activated by supra-thermal processes is mandatory [START_REF] Flower | C-type shocks in the interstellar medium: profiles of CH + and CH absorption lines[END_REF][START_REF] Lesaffre | Low-velocity shocks: signatures of turbulent dissipation in diffuse irradiated gas[END_REF][START_REF] Godard | Chemical probes of turbulence in the diffuse medium: the TDR model[END_REF] . Its fast destruction due to collisions with H and H 2 makes its lifetime extremely short, even in dilute gas. Being a light hydride with high dipole moment, the critical density of its J = 1-0 transition is ∼ 10 7 cm -3 , about 10 4 × larger than that of CO J = 1-0, so in low-density gas, most molecules are in the ground state providing large J =1-0 line opacities. In absorption, CH + is therefore a tracer of dissipation of mechanical energy in lowdensity gas (n H < 10 2 cm -3 ). In contrast, in emission, CH + J= 1-0 line is a tracer of dense (n H > 10 5 cm -3 ) ultaviolet(UV)-illuminated gas [START_REF] Eales | The Herschel ATLAS[END_REF] , such as UV-irradiated shocks because, unlike most molecules, the abundance of CH + is enhanced in intense UV fields [START_REF] Godard | A complete model of CH + rotational excitation including radiative and chemical pumping processes[END_REF][START_REF] Falgarone | Strong CH + J = 1-0 emission and absorption in DR21[END_REF] . All our targets are lensed starburst galaxies, with most likely different magnification for the dust continuum emission, the gas seen in absorption and that seen in emission. The angular and spectral resolution of our ALMA data allow some information on the respective distributions of the absorption and the emission in velocity (Extended Data Fig. 1) and of the continuum dust emission and CH + absorption and emission in space (Extended Data Figs. 2 and3). To interpret these data, modelling of differential lensing is warranted. However, the optical depth estimates of the absorption lines are not affected by differential lensing because the background dust emission and the absorbed emission experience the same lensing. Fits have been performed simultaneously on the emission and absorption lines, and the uncertainties on the parameters inferred (Table 2) have been evaluated with the bootstrap method to take into account the correlation between these parameters, in addition to the data noise level. The broad overlap between the absorption and emission velocity ranges makes the uncertainty on the absorption linewidth larger than 100% for NAv1.56. For NAv1.144, no solution with broad emission is found by the fitting procedure, although it is physically plausible. The emission line luminosities, computed as

L CH + = 1.04 × 10 -3 S CH + ∆v em ν rest D 2 L /(1 + z)
where D L is the luminosity distance expressed in Mpc (Extended Data Table 1), S CH + ∆v em is the line-integrated flux in Jy km s -1 and ν rest = 835.07895 GHz, are in the range 2-7×10 8 L . The absorption line optical depth, computed at the velocity v 0 of the absorption peak, is derived from these fits:

τ 0 = -ln(S abs (v 0 )/[S cont + S em (v 0 )],
where S cont , S abs and S em are the continuum, absorption and emission fluxes respectively. The CH + column densities are the average values in a solid angle subtended by the area, πr 2 SMG , of the background source assumed to be uniformly covered by the absorbing screen:

N (CH + ) = 3 × 10 12 cm -2 τ 0 (∆v abs /1 km s -1 )
where ∆v abs is the FWHM of the absorption line. The formal errors on N (CH + ) essentially follow the signal-to-noise ratio of the data. The velocity resolution is coarse and we cannot rule out the possibility that portions of the line might be saturated. The estimated values of N (CH + ) could therefore be lower limits, so in Table 2 we provide them without a compiled error.

-

The turbulent framework Since its discovery in the diffuse interstellar medium (ISM) [START_REF] Douglas | CH + in interstellar space and in the laboratory[END_REF] , the puzzle raised by the large CH + abundances had focused on our Solar Neighbourhood because the line used was an electronic transition (λ ∼ 0.4µm) detected in absorption against nearby stars. The rest frequency of CH + (J =1-0) falling in a domain of high atmospheric opacity, the first observations of CH + (1-0) were achieved with the Herschel satellite. The line was detected in absorption against the dust continuum of remote Galactic star-forming regions with Herschel/HIFI [START_REF] Falgarone | Strong CH + J = 1-0 emission and absorption in DR21[END_REF][START_REF] Falgarone | CH + (1-0) and 13 CH + (1-0) absorption lines in the direction of massive star-forming regions[END_REF][START_REF] Godard | Comparative study of CH + and SH + absorption lines observed towards distant star-forming regions[END_REF] and in absorption or emission with Herschel/SPIRE [START_REF] Naylor | First detection of the methylidyne cation (CH + ) fundamental rotational line with the Herschel/SPIRE FTS[END_REF] . The Herschel/HIFI data reveal the large opacity of the CH + (1-0) line in absorption in the Milky Way and show that the CH + absorption profiles differ considerably from those of other lines [START_REF] Godard | Comparative study of CH + and SH + absorption lines observed towards distant star-forming regions[END_REF] . The Herschel/SPIRE detections in nearby starbursts [START_REF] Van Der Werf | Black hole accretion and star formation as drivers of gas excitation and chemistry in Markarian 231[END_REF][START_REF] Spinoglio | Submillimeter Line Spectrum of the Seyfert Galaxy NGC 1068 from the Herschel-SPIRE Fourier Transform Spectrometer[END_REF][START_REF] Rangwala | Observations of Arp 220 Using Herschel-SPIRE: An Unprecedented View of the Molecular Gas in an Extreme Star Formation Environment[END_REF] lack spectral resolution to definitely disentangle emission and absorption but Arp220 has a strong absorption line while Mrk231 has the line in emission. A recent detection of the visible lines of CH + against a supernova in M82 reveals a highly turbulent component [START_REF] Ritchey | Diffuse Atomic and Molecular Gas in the Interstellar Medium of M82 toward SN[END_REF] . In the Milky Way, the CH + abundances observed in the diffuse molecular gas were found to scale with the turbulent energy transfer rate per unit volume, . The models of turbulent dissipation regions provide an expression of this scaling [START_REF] Godard | Chemical probes of turbulence in the diffuse medium: the TDR model[END_REF] . In a turbulent cascade, the dissipation rate per unit volume is on average equal to the energy transfer rate per unit volume in the cascade:

l = 1 2 ρ l v l 2
/τ l where τ l = l/v l is the dynamical time of scale l, and v l and ρ l are the mean velocity and mean mass density at that scale, respectively. It is this quantity,

l = 1 2 ρ l v l 3 /l,
that is observed to be uniform on average, albeit with a huge scatter, in the inner Galaxy [START_REF] Hennebelle | Turbulent molecular clouds[END_REF] and is found to be an invariant in numerical simulations of compressible turbulence [START_REF] Kritsuk | A supersonic turbulence origin of Larson's laws[END_REF] . This quantity is shortened into "turbulent energy flux" in what follows because it is a flux in the wavenumberspace of the turbulent cascade. The fraction of the turbulent energy flux required to sustain the observed number of CH + molecules in the volume πr 2 bg L subtended by the background continuum source (radius r bg ) and the length of the line of sight L is α = ĖCH + / Ėturb , with

ĖCH + = N (CH + )πr 2 bg E form /t
and Ėturb = × πr 2 bg L. The turbulent energy flux is a large scale average and the CH + lifetime is t = 1yr f -1 H 2 (n H /50 cm -3 ) -1 , where f H 2 = 2n(H 2 )/n H is the H 2 fraction and n H and n(H 2 ) are the proton and H 2 densities, respectively. This is why the CH + opacity of absorption lines that provides N (CH + ) is a direct measure of the turbulent energy flux in the medium. Introducing N H,abs = n H L, the total column density of gas of average density n H sampled by the absorption, the observed abundance is X(CH + ) = N (CH + )/N H,abs , and therefore α = X(CH + ) × n H / × E form /t. In the Solar Neighbourhood 9 , the average values are X(CH + ) = 7.6 × 10 -9 , n H ∼ 1 cm -3 , and ∼ 10 -25 erg cm -3 s -1 , so that α = 7.7 × 10 -4 f H 2 n 50 , where n 50 = n H /(50 cm -3 ) is the local gas density. Galactic observations of [CII] absorption and CI emission along the same background sources where CH + absorption is observed show that this gas is indeed the cold neutral medium (CNM) [START_REF] Gerin | C II] absorption and emission in the diffuse interstellar medium across the Galactic plane[END_REF] . The average H 2 fraction in this low-density molecular gas is found to be f H 2 = 0.4 from CH and HF absorption lines along the same sightlines [START_REF] Godard | Comparative study of CH + and SH + absorption lines observed towards distant star-forming regions[END_REF] , with a density in the range n H = 30-50 cm -3 . The filling factor of the CNM in which CH + forms, f v = n H /n H , is much less than unity.

-Determination of the turbulent reservoir properties

We follow the same line of argument to study the turbulent gas around the SMGs. The turbulent approach is supported by the fact that the narrowest [broadest] absorption components have the smallest [largest] velocity offsets (v off , Table 2). We also rely on a fundamental property of turbulence -a pencil beam across a turbulent volume samples all the velocities because of the velocity correlations -to relate the size of the turbulent reservoir to the mean turbulent velocity provided by the absorption linewidth. We assume a negligible contribution of rotation to the absorption linewidths (see Extended Data Fig. 1). The unknown CH + abundance is therefore an average value over the volume sampled by the absorption, in the same way as in the Milky Way. We again express that a fraction α of the turbulent energy flux

Ėturb = 1 2 M abs v 3 abs /r T R
in the gas mass sampled by the absorption, M abs = µ p m H N H,abs πr 2 SMG , is sufficient to sustain the number of CH + molecules observed in absorption against the background continuum source of radius r SMG , ĖCH + = N (CH + ) πr 2 SMG E form /t. We thus obtain the relation between the radius of the turbulent reservoir and the CH + abundance:

r TR X(CH + ) = α µ p m H 2 t E form v 3 abs
where the mean turbulent velocity is v abs = 0.7∆v abs and µ p is the mean mass per particle. CH + abundances in the range (0.8-5.2)×10 -8 follow from the estimated values of r TR , and hydrogen column densities in the range 10 23 -9×10 21 cm -2 are inferred. The mass of the turbulent reservoir is then computed by assuming mass conservation across the different shells, regardless of the global inward/outward or turbulent motions, so that the average density of the reservoir decreases with radius as n H ∝ r -2 . As a result, the hydrogen column density at impact parameter b decreases as N H (b) ∝ b -1 and the mass is M TR = 2π µ p m H r TR r SMG N H,abs .

Introducing ∆v 400 = ∆v abs /(400 km s -1 ) and N 14 = N (CH + )/(10 14 cm -2 ) for the observed quantities, so that N 14 /∆v 400 = 12τ 0 , X 8 = X(CH + )/10 -8 for the inferred abundances, α 3 = α/10 A mean mass per particle µ p = 1.27 is adopted. These expressions help distinguish the dependence of the inferred sizes and masses on the observed quantities (∆v abs and N 14 ) and free parameters (f H 2 , n 50 and α 3 ). The local gas density n 50 and the H 2 fraction introduce some degeneracy. However, the range of possible values is limited. From the Galactic studies, f H 2 = 0.4 and the gas densities of the CH + -rich regions are in the range 30 -50 cm -3 . All of the determinations are therefore within factors of a few. The CH + abundances are close to the Galactic values, as would be expected if the CH + production rate were comparable to that in the Galaxy and the metallicities were simlar to solar values, as is observed 45 in star-forming galaxies at z ∼ 2.

The kinetic luminosity is equal to the energy input rate needed to balance the dissipation of the kinetic energy E turb = 1 2 M TR v 2 abs in a dynamical time t TR = r TR /v abs or, with the appropriate units, L TR = 2.6 × 10 9 L M TR,10 ∆v 2 400 /t 50 , where M TR,10 = M TR /(10 10 M ). The mass outflow rate at velocity v out that is required to balance turbulent dissipation, L TR = 1 2 η Ṁout v 2 out (assuming an efficiency η < 1 of the energy transfer) is

Ṁout ∼ 19 η M yr -1 L 9 /v 2 800
where L 9 = L TR /(10 9 L ) and v 800 = v out /(800 km s -1 ).

The turbulent reservoirs are more than 10 times larger the SMGs. Such kpc-scale haloes of lowdensity gas, comprising a neutral diffuse molecular phase, have been invoked to interpret the low-lying OH absorption lines in Arp 220 [START_REF] Gonzalez-Alfonso | The far-infrared spectrum of Arp 220[END_REF] and NaI D absorption lines [START_REF] Martin | Mapping Large-Scale Gaseous Outflows in Ultraluminous Galaxies with Keck II ESI Spectra: Variations in Outflow Velocity with Galactic Mass[END_REF] in ULIRGs, as well as CO(1-0) lines broader than CO(3-2) lines in distant SMGs [START_REF] Ivison | Tracing the molecular gas in distant submillimetre galaxies via CO(1-0) imaging with the Expanded Very Large Array[END_REF] . Large kpc-scale haloes of atomic hydrogen are being discovered in the Lyα line around high-redshift galaxies [START_REF] Borisova | Ubiquitous Giant Lyα Nebulae around the Brightest Quasars at z ∼ 3.5 Revealed with MUSE[END_REF][START_REF] Wisotzki | Extended Lyman α haloes around individual highredshift galaxies revealed by MUSE[END_REF] . Last, the gas of the turbulent reservoirs is multi-phasic because its average density n TR = N H,abs /r TR ∼ 0.2-6 cm -3 , is less than that of the CH + -rich gas at n H ∼30-50 cm -3 , that fills only a low (< 10%) fraction of the volume.

-Models of irradiated shocks

A medium shocked at high velocity emits UV photons that interact with the pre-shock gas [START_REF] Hollenbach | Molecule formation and infrared emission in fast interstellar shocks. III. Results for J shocks in molecular clouds[END_REF] . Depending on the shock velocity, this radiative precursor heats, dissociates (for v sh > 90 km s -1 ) or even ionizes (for v sh > 120 km s -1 ) the upstream medium, hence preventing the formation of molecular species. Therefore, while CH + emission has very broad line profiles, it is necessarily formed in shocks of much lower velocity. The extremely broad CH + emission lines trace the large velocity dispersion of lower-velocity shocks, as observed in the kiloparsec-scale shock of the Stephan's Quintet [START_REF] Appleton | Shock-enhanced C+ Emission and the Detection of H 2 O from the Stephan's Quintet Group-wide Shock Using Herschel[END_REF][START_REF] Guillard | Turbulent Molecular Gas and Star Formation in the Shocked Intergalactic Medium of Stephan's Quintet[END_REF] . The production of CH + in low-velocity interstellar shocks has been studied using the Paris-Durham shock code, a state-of-the-art model that was recently improved to take into account the basic physics of photo-dominated regions [START_REF] Lesaffre | Low-velocity shocks: signatures of turbulent dissipation in diffuse irradiated gas[END_REF] . In these starbursts, UV-illumination comes from star formation with a possible contribution from the thermal emission of the optically thick accretion disks of AGNs. Using observations of 22 AGNs 51 , we estimate that the luminosity of a typical AGN with redshift z > 1 is comparable to that due to high-mass star formation, taking dilution within the starburst into account. We ran the code for a range of values of the preshock density (10 3 -10 6 cm -3 ), the UV radiation field (1 -10 4 × the solar neighbourhood value), and the shock velocity (10 -80 km s -1 ), assuming a magnetic field of 2[n H /1 cm -3 ] 1/2 µG, and a visible extinction (i.e. shielding from UV photons) of 1 magnitude. We find that the critical velocity, above which stationary C-type shocks cannot propagate, strongly depends on the density and the radiation field. Below a critical value (v sh ∼ 40 km s -1 i.e. C-type shocks), the large abundance of H 2 and the ion-neutral decoupling are sufficient to produce column densities of CH + as large as about 10 14 cm -2 in the shock front, within timescales smaller than the expansion timescale of the gas. Above this critical value (i.e. J-type shocks), the rapid and self-sustained collisional dissociation of H 2 prevents the immediate formation of CH + in the post-shock where H 2 reforms and the temperature is still high enough.

Column densities of CH + as large as those of C-type shocks are obtained for a kinetic energy input 5 × larger, over distances 100 × larger and for timescales larger than the expansion timescale of the post-shock gas. We conclude that J-type shocks are far less efficient than Ctype shocks to produce CH + . Large CH + column densities, up to N (CH + ) ∼ 10 14 cm -2 and spread over several km s -1 per shock, are formed in v sh = 20 km s -1 C-type shocks propagating in dense pre-shock gas (n H ∼ 10 4 cm -3 ). Lines as bright as 1 mJy/beam, estimated in the Large Velocity Gradient approximation [START_REF] Van Der Tak | A computer program for fast non-LTE analysis of interstellar line spectra[END_REF] , are obtained over a broad range of shock densities ≥ 10 5 cm -3 and gas temperatures ≥ 100K. This result does not greatly depend on the gas density. The thickness of each shock is very small, about a few hundred astronomical units. The number of shocks in the beam is difficult to estimate because it depends on the actual geometry of the shock ensemble and its distribution in velocity. Corrected for lensing with the magnifications estimated for the continuum emission, the intrinsic CH + peak brightnesses are about 0.1 mJy/beam. Hence, millions of low-velocity shocks are required in the beam to reproduce the observed CH + (1-0) brightness spread smoothly over ∼ 1,000 km s -1 . The CH + (1-0) line opacity of individual shocks is less than unity; therefore there is no limitation on the emitted line intensity due to radiative transfer. The lack of detection of CO(1-0) lines as broad as the CH + (1-0) lines in those of our targets (except NAv1.56) where CO(1-0) has been observed (Extended Data Table 1) is due to the much larger emissivity of the CH + (1-0) transition compared to that of CO(1-0). The Einstein coefficients for spontaneous emission are A 10 = 6.4×10 -3 s -1 for CH + (1-0) [START_REF] Eales | The Herschel ATLAS[END_REF] and A 10 = 7.4×10 -8 s -1 for CO(1-0). For an optically thin line, an upper limit of the i → j line integrated intensity is provided by N i A ij hν ij , so the CO column density in the J = 1 level has to offset the huge ratio A 10 ν 10 (CH + )/A 10 ν 10 (CO) ∼ 1.6 × 10 6 . Therefore as long as the column density of CO in the shocks does not exceed that of CH + by a factor of about 10 7 , which is the case in the lowvelocity C-type shocks, the CO(1-0) line remains very weak and probably below the detection level.

  Figure1: ALMA rest-frame 360-µm continuum images of our lensed targets. Solid contours are in steps of 3σ × 2 n , n = 0, 1, 2..., and dashed contours at -3σ, where σ is the r.m.s. noise level (Table1). The synthesised beam is shown in the bottom right corner. In the Planck cosmology 30 , 1 arcsec corresponds to 8.4 kpc at z =2.3 (uncorrected for lensing magnification).

  10 14 cm -2 kpc 10 10 M 10 9 L

	Eyelash	30±53 1301±112 232±40 143±11	1.2	8.3	8	1.4	1.2
	G09v1.40 28±34 1124±87	361±24 111 ±7	0.70 7.4	12	1.1	2.3
	SDP17b	-21±96 1394±157 490±53 -315± 28 0.33 4.7	17	0.9	3.4
	NAv1.56	0 d	> 814	579 d	-249	0.29 4.9	20	0.8	4.5
	NAv1.144 0 d	> 217	799 d	-227	0.25 6.	-	-	-
	SDP130	85±21 248±88	-	-	-

  -3 and t 50 = t TR /(50 Myr) for the most probable age of the turbulent reservoirs, with t TR ∼ t SMG /2 ∼ 50 Myr, the radius and mass of the turbulent reservoirs are r TR = 14.0 kpc × ∆v 400 t 50 and M TR = 3.8 × 10 10 M × τ 0 t 2 50 r SMG,kpc (f H 2 n 50 /α 3 ) and the CH + abundance is X 8 = 2.5 ∆v 2 400 α 3 /(t 50 f H 2 n 50 ).

Acknowledgements: ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. NRAO is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. R.J.I. acknowledges support from ERC in the form of the Advanced Investigator Programme, 321302, COS-MICISM. EF and BG acknowledge support from the national CNRS program Physique et Chimie du Milieu Interstellaire (PCMI).

Data Availability Statement: This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00164.S. The datasets generated and/or analysed during the current study are available in the ALMA archive (http://almascience.eso.org/aq/?project code=2013.1.00164.S) and are also available from the corresponding author upon reasonable request.

Author contributions: EF, EB, FB and DE conceived the initial scientific argument and wrote the ALMA proposal with BG, MAZ, PMA, AO and RSB. MAZ reduced the ALMA data. BG and EF analyzed the spectra. BG provided the results of the shock models. RJI, IOG and FW were invited to join the team at a later stage to provide the results of the lens models (IOG) and contribute to a year-long debate on the data interpretation. EF wrote the paper with contributions from all the authors.

Author informations: The authors declare that they have no competing financial interests. Correspondence and requests for materials should be addressed to E.F. (email: edith.falgarone@ens.fr). a Luminosity distances computed for H 0 = 67.3 km s -1 Mpc -1 , Ω M = 0.315 and Ω Λ = 0.685 [START_REF]Planck 2013 results[END_REF] . b H-ATLAS sources [START_REF] Eales | The Herschel ATLAS[END_REF] , redshift determinations 14 c The CO(1-0) results on the H-ATLAS sources are from GBT/Zspec observations [START_REF] Harris | Blind Detections of CO J = 1-0 in 11 H-ATLAS Galaxies at z = 2.1-3.5 with the GBT/Zpectrometer[END_REF]

Extended Data Table 1: Additional properties of the lensed SMGs

20

Extended Data Figure 1: Position-velocity diagrams of CH + emission and absorption along selected cuts across the sources. The cuts are made along East-West directions, except for the Eyelash where it is made along the long-axis of the lensed images and SDP17b where it is along a NE-SW direction. CH + emission appears in white (blue contours) and absorption in black (green contours). The first contour level and steps are 2σ. A velocity gradient is seen in the Eyelash absorption, twice smaller than that detected in CO [START_REF] Swinbank | The Interstellar Medium in Distant Star-forming Galaxies: Turbulent Pressure, Fragmentation, and Cloud Scaling Relations in a Dense Gas Disk at z = 2.3[END_REF] .

Extended Data Figure 2: CH + emission and absorption overlaid on dust continuum emission for the Eyelash, SDP17b and G09v1.40. The CH + line-integrated emission (blue contours) and absorption (red contours), with contour levels in steps of 2σ, is overlaid on continuum emission (grey scale). All of the images are lensed, so the differences between the distribution of dust continuum and CH + line emission are affected by differential lensing.