Bastien Le Gloannec
email: bastien.le-gloannec@univ-orleans.fr

The 4-way deterministic Periodic Domino Problem is undecidable

Keywords: wang tiles, deterministic tilesets, domino problem, undecidability

The most fundamental undecidable question on tilings is the Domino Problem that asks whether a Wang tileset tiles the discrete plane. Lukkarila proved in 2009 that it remains undecidable when restricting the input to the class of 4-way deterministic tilesets. Due to the existence of aperiodic tilesets, the most natural distinct variant of this problem is the Periodic Domino Problem which asks whether a Wang tileset admits a periodic tiling of the plane. This problem is also undecidable. Jeandel recently discovered a new and elegant proof for this result. Inspired by this new proof technique and some ingredients from Lukkarila's construction, we prove that it remains undecidable when restricted to 4-way deterministic tilesets.

Wang tiles are unit square tiles with a color on each edge. Wang tiles and their purely combinatorial and syntactical formalism of discrete tilings were introduced by Wang [START_REF] Wang | Proving theorems by Pattern Recognition II[END_REF] in relation to the decision question of the ∀∃∀ fragment of the first-order logic. The most common associated decision problem is the Domino Problem which consists in deciding whether a Wang tileset (given as input) tiles the discrete plane. This problem was proved undecidable by Berger [START_REF] Berger | The Undecidability of the Domino Problem[END_REF]. A simpler construction was later introduced by Robinson [START_REF] Robinson | Undecidability and nonperiodicity for tilings of the plane[END_REF]. Several recent more synthetic approaches to this problem also exist in the literature [START_REF] Kari | The tiling problem revisited[END_REF][START_REF] Ollinger | Two-by-two substitution systems and the undecidability of the domino problem[END_REF][START_REF] Durand | Fixed point and aperiodic tilings[END_REF].

The undecidability of the Domino problem is tightly related to the existence of aperiodic tilesets, i.e. tilesets that tile the plane but never in a periodic way. Indeed, a classical variant of the Domino Problem is the Periodic Domino Problem that asks whether a Wang tileset given as input admits a periodic tiling. As the Domino Problem is co-recursively enumerable while its periodic variant is recursively enumerable, both problems would be equal and decidable in a world without aperiodicity. The Periodic Domino Problem was proved undecidable by Gurevich and Koryakov [7]. This result is obtained by a slight modification of Berger's construction [START_REF] Berger | The Undecidability of the Domino Problem[END_REF] for the Domino Problem. Similarly, Robinson's construction from [START_REF] Robinson | Undecidability and nonperiodicity for tilings of the plane[END_REF] can be adapted to prove this result as detailed in [START_REF] Allauzen | Tiling Problems, in: The Classical Decision Problem[END_REF] for instance. It can also easily be seen that the undecidability of the Periodic Domino Problem is equivalent to the recusive inseparability of the aperiodic tilesets from the periodic ones (tilesets that admit a periodic tiling). A new and in many ways simpler proof technique was recently introduced by Jeandel [START_REF] Jeandel | The periodic domino problem revisited[END_REF]. This method does not rely anymore on a previous proof of the undecidability of the Domino Problem, but simply takes as input any arbitrary aperiodic tileset.

A Wang tileset is 4-way deterministic if any of its tilings is uniquely determined and can be locally reconstructed from the tiles from any path crossing each line and column of the discrete plane. Deterministic tilesets were originally introduced by Kari [START_REF] Kari | The nilpotency problem of one-dimensional cellular automata[END_REF] to prove the undecidability of the nilpotency problem of one-dimensional cellular automata by reduction of the Domino Problem with input restricted to one-way deterministic tilesets. This problem had actually already been considered and proved undecidable in a somehow forgotten proof of the undecidability of the Domino Problem by Aanderaa and Lewis [START_REF] Aanderaa | Linear Sampling and the ∀∃∀ Case of the Decision Problem[END_REF]. A two-way deterministic aperiodic tileset also already existed in the literature: the 16 Wang tiles derived from the geometrical aperiodic set of Ammann A2 [START_REF] Grunbaum | Tilings and Patterns[END_REF]. A 4-way deterministic aperiodic tileset was first built in 1999 by Kari and Papasoglu [START_REF] Kari | Deterministic aperiodic tile sets[END_REF] and an isomorphic tileset was recently discovered starting from a different construction in [START_REF] Le Gloannec | Substitutions and strongly deterministic tilesets[END_REF]. The Domino Problem with input restricted to the class of 4-way deterministic

t W t S t E t N Figure 1: A Wang tile t
tilesets was proved undecidable ten years later by Lukkarila [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF]. Although the Periodic Domino Problem was proved undecidable when restricted to one-way deterministic tilesets by Mazoyer and Rapaport [START_REF] Mazoyer | Global fixed point attractors of circular cellular automata and periodic tilings of the plane: Undecidability results[END_REF] to derive the undecidability of the nilpotency of one-dimensional cellular automata over periodic configurations, the far more restrictive setting of 4-way deterministic tilesets remained untreated.

We prove in this article that the 4-way deterministic Periodic Domino Problem is undecidable. A quite natural strategy to obtain such a result would be to determinize a construction from an existing proof of undecidability in the non-deterministic case. Some constructions are however intrinsically non-deterministic. Moreover, determinism is a low-level, i.e. syntactical, property of the tileset, hence the determinization process requires a comprehensive knowledge and understanding of the construction at the level of the tiles and a careful treatment of non-determinism issues. Instead, we are inspired here by the general idea of the proof technique from Jeandel [START_REF] Jeandel | The periodic domino problem revisited[END_REF], which allows a rather neat high-level presentation and prevents us from a lot of purely technical annoying considerations. Our proof also builds upon the 4-way deterministic Turing machine simulation from Lukkarila [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF].

Wang tiles and determinism

A Wang tile is a unit square tile with a color on each edge. Formally, it is a quadruple t = (t W , t S , t E , t N ) ∈ C 4 where C is a finite set of colors, as depicted in figure 1. A tileset is a finite set of Wang tiles. Given a tileset τ , a tiling of the discrete plane is a map c : Z 2 → τ that associates to each point of the plane a tile of τ such that adjacent tiles share the same color on their common edge: for all (x, y) ∈ Z 2 , t(x, y) E = t(x + 1, y) W and t(x, y) N = t(x, y + 1) S .

It is often convenient to describe Wang tiles by drawing simple patterns on them, instead of explicitely writing color letters. In that case, the color matching condition for tilings is represented by the continuity of pattern lines when crossing edges.

A tiling c is periodic if there exists a periodicity vector p ∈ Z 2 such that for all x ∈ Z 2 , c(x + p) = c(x). A tiling is bi-periodic if it admits two linearly independent periodicity vectors, or equivalently if it admits both a horizontal and a vertical period. The following lemma (folklore) links these notions regarding sets of tilings.

Lemma 1 (folklore). The set of tilings by a tileset τ contains a periodic tiling if and only if it contains a bi-periodic tiling.

A tileset is aperiodic if the set of tilings by τ is non-empty and does not contain any periodic tiling. A tileset τ is ne-deterministic if each tile t ∈ τ is uniquely determined by its pair of colors (t W , t S ). Equivalently, for any pair of tiles (t W , t S ) ∈ τ 2 , there is at most one tile t ∈ τ that is simultaneously compatible to its west side with t W and to its south side with t S . {nw, sw, se}-determinism are defined symmetrically. A tileset is 4-way deterministic if it is simultaneously deterministic in the four directions. In this article, we will essentially consider 4-way deterministic tilesets. As we will often make use of this without mentionning it, observe that it clearly follows from the definitions that any subset of a 4-way deterministic tileset remains 4-way deterministic.

The earliest and most common decision problem in the literature on tilings is the following.

Problem 1 (Domino Problem, [START_REF] Wang | Proving theorems by Pattern Recognition II[END_REF]). Given a tileset τ , decide whether there exists a tiling by τ .

This problem was proved undecidable in [START_REF] Berger | The Undecidability of the Domino Problem[END_REF]. The periodic variant of the Domino Problem is the following. Problem 2 (Periodic Domino Problem). Given a tileset τ , decide whether there exists a periodic tiling by τ .

This problem was introduced and proved undecidable in [START_REF] Gurevich | Remarks on Berger's paper on the domino problem[END_REF]. A new proof technique was recently discovered in [START_REF] Jeandel | The periodic domino problem revisited[END_REF].

In this article, we are interested in the 4-way deterministic restrictions of these problems where the input tileset is furthermore required to be 4-way deterministic. The 4-way deterministic Domino Problem was proved undecidable in [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF]. We prove in the following of this article that the 4-way deterministic Periodic Domino Problem is also undecidable.

Some transformations on tilesets

In this section, we introduce some simple operations or transformations on tilesets that preserve the determinism and will be useful in our construction.

Product. Given two tilesets τ 1 ⊆ C 4 1 and τ 2 ⊆ C 4 2 , the product tileset τ = τ 1 × τ 2 is formally defined as the Wang tileset τ = {((w 1 , w 2 ), (s 1 , s 2 ), (e 1 , e 2 ), (n 1 , n 2 )), (w i , s i , e i , n i ) ∈ τ i , 1 ≤ i ≤ 2} over the set of colors C 1 × C 2 .
Although the set τ is not strictly speaking the cartesian product of sets τ 1 and τ 2 , it is convenient to interpret it as such, seeing a product tileset as a two-layered tileset whose tiles hold a tile of τ 1 on the first layer and a tile of τ 2 on the second layer with a local matching condition requiring that the matching conditions of both layers are verified. It is not difficult to see that if τ 1 and τ 2 are 4-way deterministic, then τ 1 × τ 2 is also 4-way deterministic.

Disjoint mirror. Given a tileset τ over a set of colors C, the horizontal disjoint mirror tileset of τ is the tileset τ h obtained by duplicating all colors and swapping west and east colors of tiles from τ . More formally, we create a disjoint copy C h of C whose colors are denoted a h for each a ∈ C and we define τ h = {(e h , s h , w h , n h ) ∈ C 4 h , (w, s, e, n) ∈ τ )}. We symmetrically define the vertical disjoint mirror τ a v over a disjoint copy C v of C (whose colors are denoted a v for each a ∈ C) by

τ v = {(w v , n v , e v , s v ) ∈ C 4 v , (w, s, e, n) ∈ τ )}.
We finally define the central disjoint mirror τ c over a disjoint copy C c of C (whose colors are denoted a c for each a ∈ C) by τ c = {(e c , n c , w c , s c ) ∈ C 4 v , (w, s, e, n) ∈ τ )}. This last transformation corresponds to the combination of the two previous ones (hence to a central symmetry on the tiles). These three definitions are illustrated on figure 2.

It is clear that if a tileset is 4-way deterministic, any of its disjoint mirrors is also 4-way deterministic (even though some determinism directions are interchanged in the transformation process). As their sets of colors are disjoint, any of the (disjoint) unions of a tileset with some of its disjoint mirrors, in particular τ ⊔ τ h ⊔ τ v ⊔ τ c , is also 4-way deterministic.

Grouping. Given a tileset τ , the 2 × 2 grouping of τ is the tileset τ 2×2 over the set of colors τ 2 defined by

τ 2×2 = (ac, cd, bd, ab) ∈ (τ 2 ) 4 , a b c d is a valid 2 × 2

pattern by τ

The tileset τ 2×2 simply is a coding in Wang tiles of 2 × 2 patterns by τ with the local matching condition that two horizontally (resp. vertically) adjacent patterns overlap on one column (resp. row). Considering a tile of τ 2×2 , we will use the notations , , and to identify the four underlying tiles of τ . Observe that if τ is 4-way deterministic, then its 2 × 2 grouping τ 2×2 is also 4-way deterministic.

The deterministic reduction

In this section, we prove the following result.

Theorem 1. The Periodic Domino Problem restricted to 4-way deterministic tilesets is undecidable.

We proceed by reduction of the Halting problem of Turing machines on the empty tape. Both problems are recursively enumerable. Given a Turing machine M , we build a 4-way deterministic tileset that always admits non-periodic tilings (some of them will actually not even contain Turing computations) but admits periodic tilings if and only if the computation of M halts from the empty tape.

The built tileset will be a product of several layers, each of which will be 4-way deterministic.

First layer: Aperiodic layer with vertical mirror lines. We consider here a 4-way deterministic aperiodic tileset τ a on the set of colors C. As already mentioned, such tilesets exist, some are built in [START_REF] Kari | Deterministic aperiodic tile sets[END_REF][START_REF] Le Gloannec | Substitutions and strongly deterministic tilesets[END_REF]. Let us stress that we do not need here to know anything about the way they are built or the structure of their tilings, the only required properties are 4-way determinism and aperiodicity. Let us consider the horizontal disjoint mirror τ a h of tileset τ a over a disjoint copy C h of C (whose colors are denoted a h for each a ∈ C). It is clear by construction that τ a h is aperiodic. We also add a set of vertical mirror line tiles by introducing a new color # / ∈ C ∪ C h and defining M v = {(a, #, a h , #) and (a h , #, a, #), a ∈ C} as depicted in figure 3(a) (where the color # is painted black). We then define the tileset τ 1 as the disjoint union of the three previously considered tilesets τ 1 = τ a ⊔ τ a h ⊔ M v . A typical tiling by τ 1 is represented on figure 3(b). Due to the addition of mirror lines, this set is no longer aperiodic.

Lemma 2. τ 1 satisfies the following properties:

1. τ 1 is 4-way deterministic.

Every horizontally periodic tiling by τ 1 contains a vertical mirror line. Hence, by periodicity, it contains infinitely many such lines, the distance between two consecutive lines being bounded.

Proof. For the point 1, as τ a is 4-way deterministic, τ a h also has this property. As their sets of colors are disjoint, τ a ⊔ τ a h is 4-way deterministic. The tileset M v is also 4-way deterministic since the color on each side of a mirror line is uniquely determined by the color on the other side. The color # introduced in M v being new, τ 1 is in turn 4-way deterministic.

For 2, if such a tiling does not contain any vertical mirror line, then it must be a tiling by one of the aperiodic tilesets τ a or τ a h . It is hence non-periodic, a contradiction.

τ 2 = τ 1 × τ v .
Lemma 3. τ 2 satisfies the following properties:

1. τ 2 is 4-way deterministic.

2. Every bi-periodic tiling by τ 2 contains a (irregular) grid (whose vertical and horizontal lines are to be read on the first and second layer respectively) with bounded rectangle size.

Proof. The point 1 is clear as τ 2 is a product of 4-way deterministic tilesets. The point 2 is a direct consequence of the argument of lemma 2 (point 2) applied both horizontally on τ 1 and vertically on τ v .

Third layer: Turing computation. We consider here the 4-way deterministic tileset introduced by Lukkarila in [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF] (section 5, theorem 6) to simulate any given Turing machine in the tilings. This tileset is used to prove the undecidability of the 4-way deterministic Domino Problem with a seed tile, that asks whether a 4-way deterministic input tileset containing an identified seed tile admits a tiling of the plane in which the seed tile appears. It is generally easier to deal with that problem than the general Domino Problem as the seed tile can for instance be used to properly initialize an empty Turing machine tape with exactly one machine head. This tileset is a subset of a product of several layers, one of them (the first layer in [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF], section 5.2) being a classical simulation of a Turing machine (whose head always moves at each transition) in Wang tiles using the tiles represented on figure 4(a) (where q i is the initial state and B the blank letter). Such a tileset is obviously not itself 4-way deterministic, one of the difficulties overcame in [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF] was to determinize it by coupling it with a particular structure on other layers. The built tileset draws the space-time diagram of the machine in the tilings, as depicted in figure 4(b). The lines of the tiling realize the transitions, each line of the diagram of the machine can be read at the edge between two lines of the tiling.

Consider a Turing machine M and the associated 4-way deterministic simulation tileset τ M with unique seed tile t 0 (which simulation layer contains the indicated seed tile of figure 4(a)). We denote as F ⊆ τ M the set of tiles in τ M whose Turing machine simulation layer transmits a final state of M from its west or east color to its north color. Tiles of F mark the end of the simulated Turing computation. Note that a tiling of the whole plane cannot contain both t 0 and a tile of F : Assuming that transitions from a final state are undefined in the Turing machine (which is the usual convention), there is no transition tile in the set with a final state as south color, hence the line directly on top of the line containing a tile of F cannot be tiled without error. Let us consider the three disjoint mirrors of τ M , over the set of colors C, denoted as τ M h , τ M v , τ M c , over the respective sets of colors C h , C v , C c (disjoint copies of C). We assemble these four tilesets using a parity grid tileset G M whose tiles are represented on figure 5(a). Their purpose is to draw a grid with four types of rectangles alternating depending on the horizontal and vertical parities, as depicted in figure 5(b). The following constraint on G M will be crucial for our final result and cannot be read on the figure: Basically any color of τ M (and its mirrors) can appear along grid lines of τ M except for horizontal colors containing a machine state that cannot appear along vertical grid lines. The purpose of this alphabetical restriction is to forbid a Turing machine head from leaving its rectangle and complementarily to forbid any other Turing machine head from appearing unexpectedly in a rectangle. As the grid lines of G M play the role of mirror lines for the colors of the four mirror tilesets of τ M , G M is easily seen 4-way deterministic. We define

τ ′ M = τ M ⊔ τ M h ⊔ τ M v ⊔ τ M c ⊔ G M .
Because the colors associated to the grid lines of G M do not appear in any of the four mirror tilesets of τ M , it is not difficult to see that τ ′ M is 4-way deterministic. The tilings by τ ′ M are simply irregular parity grids (with possibly infinite rectangles) with patterns by each of the four mirror versions of τ M in each of the four types of rectangles, as depicted in figure 5(b). Patterns are locally symmetrized along grid lines.

Let us assemble τ ′ M with our previous construction τ 2 by defining a three-layered tileset τ 3 as a subset of τ 2 × τ ′ M (as τ 2 has two layers, we consider here the tiles of τ ′ M as the third layer) with the following restriction on tiles: A vertical (resp. horizontal) mirror line appears on the first (resp. second) layer (tiles of τ 1 , resp. τ v ) if and only if it is coupled with a vertical (resp. horizontal) grid line tile or a cross of G M on the third layer (tiles of τ ′ M ). Note that as a consequence, a vertical mirror line on the first layer is coupled with a horizontal mirror line on the second layer if and only if the third layer holds a cross of the grid. Moreover, a tile of τ M ⊔ τ M h ⊔ τ M v ⊔ τ M c appears on the third layer if and only if it is coupled with tiles from the aperiodic sets on the first and second layers. Lemma 4. τ 3 verifies the following properties: Proof. τ 3 is 4-way deterministic as a subset (restriction on tiles) of the product of the 4-way deterministic tilesets τ 2 and τ ′ M . By lemma 3 (point 2), the third layer of any bi-periodic tiling by τ 3 contains a grid of finite rectangles. By construction of τ ′ M , each of the four types of rectangles contains a finite pattern by its associated mirror tileset of τ M . Note however that any patterns can appear provided that they are correctly locally symmetrized along grid lines (and that they do not contain a machine head in their horizontal border colors due to the restriction on G M ). grid tiles G M , no other machine head can appear in the rectangle and the original machine head cannot leave the pattern. Hence the local rules of the tileset τ M (whose computation layer is depicted in figure 4) enforce the following rows to contain exactly the successive configurations of the Turing machine tape. As the head cannot leave without tiling error, the pattern width is necessarily sufficiently large for the head not to leave the pattern at the right border. Because of the additional requirements on the machine, its computation can be simulated without the head leaving the pattern at the left border. Hence the computation of the machine starting on the empty tape is correctly simulated in the pattern. As a final state is required in the top-right corner, the computation is required to halt. Thus the simulated machine M halts from the empty tape, which is the requested result. Moreover, the height of the pattern is exactly the halting time of the machine and its width is exactly the distance between the starting cell and the cell where the head halts, which exactly corresponds to the space used by the computation (due to the additional requirements on the machine). Thus, as considering patterns by the three mirror tilesets of τ M would lead to symmetrical descriptions, the third layer of a bi-periodic tiling must be of the kind of the figure 7.

1. τ 3 is 4-way deterministic. t 0 t c 0 t v 0 t h 0 (a) cross A ∈ F c ∈ F ∈ F h ∈ F v (b) cross D
This concludes the reduction and the proof of theorem 1.

Figure 2 :

 2 Figure 2: A tile by τ and the corresponding tiles in τ h , τ v and τ c

tiling by τ 1 Figure 3 :

 13 Figure 3: Description of tileset τ 1

Figure 4 :

 4 Figure 4: Computation layer of τ M

Figure 5 :

 5 Figure 5: Description of tileset τ ′ M

Figure 6 : 4 2.

 64 Figure 6: Alphabetical restrictions on the third layer of the grouped tileset τ 4

Both steps could actually have been done simultaneously combining directly the four mirrors of τa, we however believe that a two-step presentation is lighter and clearer.

It only remains to enforce correctly initialized and terminating Turing computations in the rectangles of the third layer of τ 3 . We denote as t h 0 , t v 0 , t c 0 and F h , F v , F c the mirror versions of the seed tile t 0 ∈ τ M and the final set F ⊂ τ M in the mirror tilesets τ h M , τ v M , τ c M respectively. Let us define τ 4 as a subset of the 2 × 2 grouping of τ 3 with the following restriction on the third layer of the 2 × 2 grouped tiles:

• if the tile is the cross A, then the tile must be t 0 ;

• if the tile is the cross D, then the tile must be in F ;

• if the tile is the cross A, then the tile must be t h 0 ; • if the tile is the cross D, then the tile must be in F h ;

• if the tile is the cross A, then the tile must be be t v 0 ; • if the tile is the cross D, then the tile must be in F v ;

• if the tile is the cross A, then the tile must be be t c 0 ; • if the tile is the cross D, then the tile must be in F c .

These constraints are summarized on the figure 6.

Lemma 5. τ 4 is 4-way deterministic.

Proof. Clear as τ 4 is a subset of τ 2×2 3 which is 4-way deterministic as τ 3 is.

For the following result, we furthermore require the Turing machine M to verify the following properties:

• during the computation, the head never moves to a position situated to the left of its starting position;

• when the computation halts (if it halts), the position of the head is exactly the rightmost position that has been reached during the computation. These requirements can be made without loss of generality as any Turing machine is simulated by a machine satisfying these constraints (which can be recursively obtained from M ).

The tileset τ M always admits tilings that do not contain Turing computation (all blank on their computation layer), hence τ 4 always admits non-periodic tilings that do not embed any Turing computations on their third layer. Proposition 1. τ 4 admits periodic tilings if and only if the Turing machine M halts from the empty tape.

Proof. If the machine halts, then it is easy to describe a bi-periodic tiling by τ 4 . Indeed, there exists a rectangle pattern P by τ M that fully describes the space-time diagram of the computation of M from the empty tape, with the tile t 0 in the bottom-left corner and a tile of F in the top-right corner (due to the previous additional requirements on M ). Then this pattern can be embedded in a rectangle of the grid of the third layer of τ 4 and the three mirror patterns by τ h M , τ v M and τ c M can be embedded in the three neighboring rectangles (east, north and north-east) of the grid. The obtained pattern of τ 4 can be repeated bi-periodically. Such a tiling is represented on figure 7.

Conversely, if τ 4 admits periodic tilings, then it admits bi-periodic tilings by lemma 1. By lemma 4 (point 2), the third layer of such a tiling draws a grid whose rectangles contain finite patterns by τ M and its three mirror copies. Let us consider the case of such a pattern by τ M (one can also observe that the three other cases would be symmetrical). Due to the restrictions imposed on the grouping, the tile at the bottomleft corner of this pattern is the seed tile t 0 . This tile enforces the first (bottom) line of the pattern to be an empty Turing machine tape with a (unique) head on the first cell. Recall that, due to the restrictions on the


H Courtois
email: herve.courtois@neel.cnrs.fr

F W J Hekking

H Q Nguyen

C B Winkelmann

Electronic coolers based on superconducting tunnel junctions: fundamentals and applications

Keywords: 74, 78, Na, 74, 45, +c

Thermo-electric transport at the nano-scale is a rapidly developing topic, in particular in superconductor-based hybrid devices. In this review paper, we first discuss the fundamental principles of electronic cooling in mesoscopic superconducting hybrid structures, the related limitations and applications. We review recent work performed in Grenoble on the effects of Andreev reflection, photonic heat transport, phonon cooling, as well as on an innovative fabrication technique for powerful coolers.

Introduction

Let us consider a NIS tunnel junction made of a Normal metal and a Superconductor, in contact through an Insulator barrier. The superconductor's electronic density of states features an energy gap of width ∆. Tunneling through the junction is possible for electrons with an energy E (compared to the Fermi level E F ) such that |E| > ∆. Only elastic processes are considered here. The charge current through such a junction biased at a voltage V is [START_REF] Tinkham | Introduction to Superconductivity[END_REF]:

I = 1 eR N ∞ -∞ n S (E)[f N (E -eV ) -f S (E)]dE (1) = 1 eR N ∞ 0 n S (E)[f N (E -eV ) -f N (E + eV )]dE, (2)
where R N is the normal state resistance, f N the electron energy distribution in the normal metal and n S (E) the normalized BCS density of states in the superconductor. No Dynes parameter smearing the superconductor's density of states is taken into account. Expression [START_REF] Tinkham | Introduction to Superconductivity[END_REF] above is the usual one, while expression (2) can be obtained by a symmetry argument. The latter shows that the current does not depend on the superconductor temperature, only on the value of the gap. When using it as a thermometer [START_REF] Nahum | [END_REF], a NIS junction is biased at a small and constant current I th . The current can be adjusted so that the related heat flow (see below) can be safely neglected. The voltage drop V th in the thermometer junction pair is then measured and compared to its value calibrated against the bath temperature T bath . This provides a measure of the electronic temperature T e with a sub-mK resolution. Let us note that the calibration is realized at equilibrium, with the superconductor and the normal metal being at the bath temperature. In contrast, practical experiments are usually conducted in quasi-equilibrium conditions where, to a first approximation, only the normal metal temperature changes. If one considers temperatures above about half the superconductor critical temperature T c /2, the related decrease of the superconductor gap ∆ has to be taken into account. The calibration then requires the full calculation of the voltage V th (T e ) using Eq. [START_REF] Nahum | [END_REF].

The same tunnel current can cool down or heat up the normal metal electron population [3,4]. The heat current through a NIS junction reads:

Q0 N = 1 e 2 R N ∞ -∞ (E -eV )n S (E)[f N (E -eV ) -f S (E)]dE, (3)
where f S is the energy distribution function in the superconductor. Equation 3 differs from Eq. 1 by an additional energy term in the integral, which also changes the symmetry compared to the charge current case: Q0 N (-V ) = Q0 N (V ) while I(-V ) = -I(V ). Equation 3 describes both the cooling and heating regimes. With a voltage bias V smaller than the gap ∆/e, the tunnel current is selectively made out of high-energy electrons (or holes) from the tails of the distribution function. This cools the electronic population of the normal metal, which means that Q 0 N is positive. Every tunnel event extracts a heat of about k B T while a heat ∆ is dissipated in the system, so that the cooling efficiency Q0 N /IV is about k B T e /∆. Above the gap (|V | > ∆/e), the heat flow Q 0 N becomes negative (the metal is heated). At large bias |V | ≫ ∆/e, the usual Joule heating regime is reached, Q 0 N is equal to -IV /2. In every case, the full Joule power IV is deposited in the device, so that a power QS = IV + Q0

N is transferred to the superconductor. In the cases studied here, we assume a quasi-equilibrium situation: the normal metal electrons and phonons follow thermal energy distribution functions at respective temperatures T e and T ph . For electrons, this is justified because the inelastic scattering time (of the order of the phase-coherence time) is much shorter than the mean residence time estimated to about 100 ns for an island of typical dimension 50 nm × 1 µm 2 [5]. The thermalization of an isolated normal metal electron population to the thermal bath occurs through electron-phonon coupling. We will make use of the standard expression:

P e-ph (T e , T ph ) = ΣU (T 5 e -T 5 ph ), (4)
of the electron-phonon coupling power, where Σ is a material-dependent constant and U is the metal volume. Several experiments have demonstrated the validity of this expression in mesoscopic (micron-sized) devices [6,7], although it is derived for a bulk material. In what follows, we take for Cu the established value Σ = 2 nW.µm -3 .K -5 .

Although electronic cooling can be obtained in a single NIS junction where the normal metal is isolated from its electrical contact pad by the way of a thermally-insulating NS contact [3], the symmetric SINIS configuration with two NIS junctions acting as coolers brings the advantage of simplicity and a doubled cooling power [4]. Thanks to the current singularity at the gap edge, the SINIS geometry is rather immune to a small junctions' resistance asymmetry [8].

The above considerations form the basis for electronic cooling. This topic has been discussed in recent reviews [9,10]. Starting from a bath temperature of about 300 mK, an electronic temperature reduction by a factor of about 3 is routinely achieved in aluminum-copper hybrid devices that are biased at a voltage just below the gap of superconducting Al. The cooling power is in the 10 pW range for a tunnel junction area of about a square micron. Microcoolers based on large lithographic junctions patterned on a membrane have demonstrated the cooling a bulk material [11] or a detector [12]. Exciting new developments towards larger cooling powers [13] and cooling of macroscopic objects [14] are emblematic of the present activity in the field.

In practice, electronic cooling in a normal metal is limited by several mechanisms of fundamental interest, see Fig. 1. Electron-phonon scattering provides the main channel for thermalization of the electronic bath with respect to the substrate considered to be a thermal bath. The question of phonon cooling is therefore of primary importance. In recent years, we have also studied how other channels of thermalization contribute. In the following, we will discuss how Andreev current or photonic heat channel alter the electronic cooling in a SINIS device. We will not discuss quasi-particle relaxation in the superconducting electrodes, although it is a major and long-standing issue [15][16][17].

Methods

We study SINIS samples based on a normal metal island symmetrically coupled to two superconducting leads through tunnel barriers, see Fig. 2a. In general, they are prepared by electronic beam lithography on a resist bilayer, followed by two-angle evaporation in a high vacuum chamber [18]. We focus on aluminium as a superconductor because of its unique high-quality oxide. As a thin film, its critical temperature T c is about 1.3 K. The normal metal is usually copper, with a diffusion constant D of about 100 cm 2 /s. The total normal-state resistance R n is in the range 2-3 kΩ for a tunnel junction area of about 1 µm 2 . In addition to the two cooling junctions, one can add superconducting tunnel probes in order to probe the normal metal electronic temperature. In one case, we have built two such complete devices one on top of the other, in order to study the phonon bath behavior.

Low temperature measurements were performed in a 3 He cryostat or a dilution cryostat. Filtering was provided by π-filters at room temperature and lossy micro-coaxial lines thermalized at the cold plate. Four-wire d.c. transport measurements were performed using home-made electronics combining independent current bias sources, one of them providing the voltage reference, the other being floating. The differential conductance dI/dV (V ) is obtained by numerical differentiation. We have taken special care to determine most accurately the sub-gap conductance of our current-biased samples down to a level of about 10 -4 of the normal state conductance. Plots are usually a combination of multiple curves covering different measurement ranges, extending over several decades of current values. Semi-log plots are useful to provide a good understanding of the charge and heat transport. In such a plot, the sub-gap differential conductance of a NIS junction at constant temperature follows a linear behavior with the voltage drop. This comes from the fact that the Fermi distribution function of the energy can be approximated to a decaying exponential at energies well above or below the Fermi level.

In order to understand some of the experimental observations, we have investigated theoretically various aspects of heat and charge transport in NIS tunnel junctions. Since the devices measured experimentally are in the diffusive limit, the appropriate theoretical framework to describe thermoelectric transport is based on quasi-classical Keldysh-Usadel equations [START_REF] Larkin | Nonequilibrium Superconductivity[END_REF][START_REF] Belzig | [END_REF]. Specifically, using this formalism, we have calculated the two-particle current due to Andreev reflection [21]. This is relevant for the experiments discussed in Section 3 below. We also considered the effect of inelastic relaxation in the superconductor lead [22]. In the absence of strong relaxation, both the electric current and the cooling power at low voltages are suppressed. We could attribute this suppression to the effect of back-tunneling of non-equilibrium quasiparticles from the superconductor into the normal metal, a topic relevant for the experiments described in Section 5 below.

Heating induced by the Andreev current

Operation and limitations of NIS electronic coolers in the low-temperature regime (T ≪ T c ) have been questioned in the past [23]. A residual smearing of the superconductor density of states has been proposed as the origin of the inferior performance with respect to theoretical expectations.

At low energy |E| < ∆, charge transfer at a NIS interface occurs through Andreev reflection [24,25]. In the normal metal, an electron (a hole) impinging on the superconducting interface is reflected as a hole (an electron), enabling the transfer of a Cooper pair into (out of) the superconductor. For a NIS tunnel junction of intermediate or low transparency, the Andreev reflection probability is predicted to be vanishingly small. Taking into account the quasi-particles confinement in the vicinity of the interface, this is no longer true. This confinement can be induced by disorder or the presence of a second barrier in the normal metal. A single quasiparticle then experiences several collisions with the interface [26,27] so that the actual Andreev reflection transmission coefficient corresponds to the coherent addition of many individual transmission probabilities. Therefore, the Andreev sub-gap current can significantly exceed the ballistic case prediction [28].

Fig. 2b displays on a logarithmic scale the differential conductance of a SINIS device for a series of cryostat temperatures [29]. We observe in the high temperature (T > 200 mK) data an upward curvature which constitutes a clear signature of electronic cooling. In the low temperature regime (T < 200 mK), clearly a different characteristic is obtained with a differential conductance peak at zero bias. The zero-bias differential conductance increases while the temperature is lowered below about 200 mK. We ascribe the low bias differential conductance peak to an Andreev current, i.e. a double particle tunnel current created by Andreev reflections at the NIS junctions.

In order to calculate the Andreev current, we take into account the finite gap ∆ and the disorder both in the normal metal and in the superconductor. We consider the 1D regime where the coherence length of an Andreev pair in the normal metal L E = hD/E [30] is much larger than the junction dimension. Usadel equations formalism can also be used and provide the same result [21]. In the probe junction, the electronic temperature can be considered as constant and very close to the cryostat temperature. In this case, the isothermal calculation fits nicely the experiment (not shown). In the cooler junction, the non-linear behavior in a semi-log scale (Fig. 2c) shows that electronic cooling (and heating) has to be taken into account.

The work performed by the current source feeding the circuit with a (Andreev) current I A generates a Joule heating I A V (> 0) that is deposited in the normal metal. The full heat balance equation for the normal metal electrons can then be written as:

Q0

N + P e-ph -

I A V = 0. ( 5
)
Here the factor 2 accounts for the presence of two tunnel junctions. In the low-temperature limit considered here, the normal metal phonons can be considered as thermalized at the bath temperature: T ph = T bath . With this complete heat balance equation taken into account, we calculate the currentvoltage characteristic at every cryostat temperature. The agreement is very good at every temperature, see Fig. 2c, and covers 4 orders of magnitude for the current. At temperatures above about 300 mK, phonon cooling introduces significant corrections, see below.

The numerical solution of the heat balance equations also provides the electron temperature at every bias. Fig. 2d shows the calculated dependence of the normal metal central island electron temperature with the voltage across the cooler. At very low temperature, the electron temperature first increases with the bias due to Andreev current heating, before decreasing due to the single quasi-particle tunnel current-based cooling effect. Although the Andreev current is a small effect in a such junction if one considers the charge current, this is no longer true if one considers the heat current. Compared to a Joule power IV , the Andreev current contributes fully to heating while the tunnel current cools with a moderate efficiency k B T e /∆ of about 5 % at a 100 mK electron temperature.

Photonic heat channel

In metallic systems, heat conduction can be achieved by electrons, phonons and also photons. The photonic channel [31] was recently revealed experimentally [32,33] at very low temperature in devices including superconducting transmission lines. With a good matching between the source and the drain, the thermal conductance of a superconducting transmission line is equal to the thermal conductance quantum [34]:

K Q = k 2
B T π/6h. The photonic channel for heat transfer can in principle couple metallic systems that are galvanically isolated, e.g. through a capacitor, or weakly coupled, e.g. through a tunnel junction. This effect can be beneficial in some cases, but also detrimental when one wants to maintain two electronic populations at different quasi-equilibrium temperatures. This is precisely the case encountered in electronic coolers.

We have investigated theoretically the photonic heat transfer through a general reactive impedance, i.e. a linear coupling circuit that contains a capacitor, an inductance, a resonant circuit or a transmission line [35]. Only the case of a capacitor will be discussed here, but the same approach can be used in all cases. We follow a simple circuit approach [36], valid at low temperatures when the relevant photons have wavelengths larger than the size of the typical circuit element. The metallic parts can then be treated as lumped elements characterized by an electrical impedance.

Let us consider a circuit of two pure resistors 1 and 2 of resistance R coupled capacitively through an impedance Z c = 1/iωC, see Fig. 3b inset. With a small temperature difference and using the dimensionless frequency

x = hω 2k B T , (6)
the exchanged power P 1,2 and the thermal conductance K C are given by:

P 1,2 = K C ∆T = T k 2 B πh ∞ 0 dxT C (x) x 2 sinh 2 (x) ∆T. ( 7
)
Fig. 3 (a) Case of a capacitive coupling. Normalized spectrum of the thermal noise power density x 2 / sinh 2 x (black dots), the photon transmission coefficient TC(x) = x 2 /(x 2 +α 2 C ) (dotted blue line) and of photonic heat (thin red full line) as a function of the dimensionless frequency, for values of the parameter αC = 0.03, 1 and 30 from top to bottom. The frequency is plotted in units of the thermal frequency 2kBT /h. We consider the case of perfect resistance matching R1 = R2 = R. (b) Dependence of the photonic power through a capacitive coupling impedance on the parameter αC for different values of the relative temperature difference ∆T /T , in units of KQ∆T , the maximum photonic power in the case of a linear response.

taking into account the frequency-dependent photon transmission coefficient:

T C (x) = x 2 x 2 + α 2 C . (8)
The cross-over frequency is determined by the parameter

α C = h 4RCk B T ; ( 9
)
it separates a low-frequency regime x ≪ α C where the capacitor is opaque and T C (x) ∼ x 2 /α 2 C ≪ 1 from a high-frequency regime x ≫ α C where the capacitor is transparent and T C (x) ∼ 1.

Fig. 3a displays the photon transmission coefficient T C . The limit α C ≪ 1 means that the capacitance is large, i.e. it has a negligible impedance over most of the thermal spectrum. The transparency T C is then equal to unity and one recovers K C = K Q . In the limit α C ≫ 1, the photonic signal is strongly suppressed by the RC filter composed of the series capacitance and the receiver resistance, leading to

K C ≪ K Q .
The total photonic power, integrated over the full frequency range, is shown in Fig. 3b as a function of the parameter α C . We compare both linear and non-linear response, changing the values of the relative temperature difference ∆T /T . The total power is maximal for small α C ; it decays as 1/α 2 C when α C is large. A cross-over between the linear regime P ∝ ∆T and the non-linear regime occurs at ∆T /T ≈ 1. When the temperature difference is large, the photonic thermal conductance is larger than its quantum because of the broader frequency range of the emitted photons. Nevertheless, the thermal conductance remain of the order of K Q , i.e. 1 pW/K at 1 K, and thus little contributes in usual situations.

Phonon cooling

As mentioned above, the electron-phonon coupling provides the main channel for coupling of the cooled electron population to the thermal bath of the cryostat. The question arises of the actual temperature of the phonon bath in the normal metal. At a temperature T , phonons have a thermal wavelength of the order of hc/k B T , where c is the material-dependent sound velocity. This dominant wavelength amounts to about 200 nm in Cu at 1 K, which is the order of magnitude of an usual device dimensions. It is thus generally assumed that phonons in a mesoscopic quantum device are strongly mixed with the substrate phonons and are thus thermalized at the bath temperature [9].

Nevertheless, phonon cooling is at the heart of the possibility of cooling a bulk detector [37] or a quantum device [38] supported on a membrane cooled by superconducting tunnel junctions. The same situation also holds for suspended metallic beams [39,40]. Recently, measurements of the electronphonon coupling strength in a thin metallic film at T ≈ 0.1-0.3 K demonstrated that it is nearly completely substrate-insensitive [41]. This supports the idea of an independent phonon population in the metallic thin film. While this phonon bath could exhibit specific properties due to its reduced dimensionality [42,43], only small deviations from bulk material properties were observed in suspended devices [44].

We have devised a dual-level sample with two independent but stacked SINIS devices including electron thermometers, see Fig. 4a. The experiment [45] consists in current-biasing one of the two level's cooler junction pair while monitoring simultaneously the related voltage drop V as well as the two levels' electronic temperatures T e,b and T e,t , where "b" stand for bottom and "t" for top. As no power is directly injected in the unbiased electronic bath, its temperature is equal to the included metal's phonon temperature. For voltages below 2∆/e applied to the bottom level, we observe the expected electronic cooling, see Fig. 4d. At voltages V above 2∆/e, we observe a hotelectron regime: the temperature T e,b increases and goes well above the bath temperature of 432 mK.

Remarkably, when the bottom level electronic temperature decreases, the top electronic temperature T e,t also diminishes with a variation ∆T e,t reaching a maximum of -2.0 mK, see Fig. 4e. As the operation of the electronic cooler is dissipative as a whole, i.e. heat is dissipated in the chip, this observation cannot be related to an improper thermalization of the chip or of electrical leads. The observed cooling of the top level therefore demonstrates phonon cooling in the normal conductors of the experiment.

The cooling/heating power was calculated by using Eq. 3, the measured electronic temperature, and a value ∆ = 214 µeV . As a fraction of the hot quasi-particles injected in the superconductor tunnels back in the normal metal, we have described this as a correction to the power proportional to N , while all other data are plotted as a function of the corrected power QN = Q0

N + αQS. The dash-dotted and the dotted lines corresponds to the point where the power Q0

N or QN respectively change sign. The full lines are fits calculated using a single set of parameters. The black dotted lines indicate the calculated phonon temperature of the island that is cooled or heated. Bath temperature is 281 mK.

QS [46]: the net power then writes QN = Q0

N + α QS . The fit parameter value α = 0.087 is comparable to what appears in the literature [46]. The electronic temperature as a function of the net power QN absolute value then follows a single curve when one goes through the maximum cooling point, see Fig. 4g.

We consider a thermal model, see Fig. 4c, assuming two distinct phonon populations at quasi-equilibrium at temperatures T ph,b and T ph,t in the bottom and top metallic islands respectively. We assume that two neighbouring phonon populations (here x and y) are coupled through a Kapitza heat flow of the form Kx-y = k xy a xy [T 4 ph,x -T 4 ph,y ] (10) [9], where a xy holds for the contact area between the two considered populations and k xy is an interface materials-dependent parameter. We have chosen to take a common value for the substrate-bottom and bottom-top Kapitza parameters. The fit-derived value k bt = k bs = 45 pW.µm -2 .K -4 compares well to values from the literature for a Cu-Si interface [47]. As for the topsubstrate coupling, considering the contact area a ts to be the area of the tunnel junctions connected to the island, one obtains a Kapitza coefficient of 920 pW.µm -2 .K -4 , which is much larger than anticipated. Thus the heat transfer occurs presumably also along the continuous Si layer separating the two levels, which explains the relatively modest amplitude of the observed phonon cooling.

From the thermal balance relations, one can calculate the phonon temperature variation in the cooled or heated metal, see dotted lines in Fig. 4g. The phonon population temperature decoupling is significant in the temperature range 0.3 -1 K, consistently with previous estimates [48]. At lower temperature, the electron-phonon bottle-neck makes the phonon temperature tend to the bath temperature. At a bath temperature below 100 mK, phonon cooling becomes negligable.

Our experimental study thus demonstrates the existence of an independent phonon bath in a quantum device. The thermal couplings are well described with the usual laws for electron-phonon coupling and Kapitza resistances. This new understanding has significant outcomes in the analysis of quantum nano-electronic devices thermal behavior.

Towards large electronic cooling power

Optimizing a SINIS device for electron cooling can be achieved by increasing the heat current and/or isolating better the cooled electron bath. As the heat current is proportional to the tunnel barriers' conductance, reducing the barriers' thickness is the first option. Nevertheless, this can lead to the appearance of two-particle Andreev reflection processes at low energy, which deposit heat in the normal metal [29]. The obvious alternative is to increase the junction area at a fixed transparency [49]. We present a method to fabricate large area SINIS devices of high quality and with a suspended normal metal [50]. The method is based on a pre-deposited multilayer of metals, which can be prepared at the highest quality. The upper normal part is suspended in the first lithography, which keeps it isolated from the substrate. The second lithography defines the junction area with any geometry of interest.

The fabrication starts with depositing a Al/AlO x /Cu multilayer on an oxidized silicon substrate. A first deep ultra-violet lithography is used to define the overall device geometry. The central part is a series of adjacent holes of diameter 2 µm and with a side-to-side separation of 2 µm, see Fig. 5b. The copper layer is etched away over the open areas using either Ion Beam Etching (IBE) or wet etching. The aluminum is then etched through the same resist mask, using a weak base. The etching time is controlled as to completely remove aluminum from the circular region within a horizontal distance of about 2 µm starting from the hole side. The line of adjacent holes visible in Fig. 5b therefore creates a continuous gap in the Al film, bridged only by a stripe of freely hanging Cu. The area of the NIS junctions is defined by a second lithography. Through the open areas, trenches are etched into the copper layer only, using one of the two methods cited above. These trenches allow to isolate a copper island, which forms the central normal metal part of the SINIS device.

Fig. 5d,e shows the current-voltage characteristic and the numericallyderived differential conductance of a typical sample at various cryostat temperatures. We have measured the electronic temperature as a function of the cooler voltage bias by using the two attached smaller junctions, see Fig. 5b, as an electron thermometer. At a bath temperature T bath of 300 mK and at the optimum bias point, the measured electronic temperature T e reaches a minimum of 240 mK.

Removing the contact between the substrate and the cooled metal by suspending the latter is quite promising for electronic refrigeration applications as it can significantly improve cooling of electrons and phonons. Our approach also has the advantage that as fabrication starts with preparing the multilayer, the wafer can be baked in ultra-high vacuum environment, which is an essential ingredient for obtaining pinhole-free large-area NIS junctions. This process has recently been improved by including a quasi-particle drain in good contact with the superconducting electrodes. A nanoWatt cooling power [13] as well as a base temperature down to 30 mK have been demonstrated [51]. This high performance can be understood only by assuming some phonon cooling in the suspended metal.

Conclusion

As a summary, we have reviewed a series of experimental studies combined with theoretical analyses of different phenomena at work in SINIS electronic coolers. Andreev current heating, coupling to phonons, phonon cooling, photonic coupling are all important limiting factors and either can dominate, depending on the practical case. The field is in constant progress, for example with the demonstration of improved performance [13] and suitability for the cooling of macroscopic objects [14].

Fig. 1

 1 Fig.1Thermal diagram for electronic cooling: a heat Q0N is extracted from a normal metal N by running a current through a SINIS device, generating also a voltage V. The current is made of a single-particule tunnel current IT and an Andreev current IA. Also shown are the different contributions discussed in this work that affect the actual cooling.

Fig. 2

 2 Fig. 2 (a) Scanning electron microscope (SEM) micrograph of a typical cooler sample made of a normal metal Cu electrode (light grey) connected to two superconducting Al reservoirs (dark grey) through tunnel junctions. One additional probe junction connected to one Al reservoir is visible at the bottom. (b) Normalized differential conductance as a function of the voltage and at the cryostat temperatures of 90 (red curve), 230 (blue), 330 (green), 440 mK (brown) and 600 mK (red). (c) Current-voltage characteristic of the S-I-N-I-S cooler junction as a function the voltage at cryostat temperature of 90 (red line), 230 (blue) and 330 mK (green) together with best-fit calculated curves. Dot-dashed line: model including the cooling by the tunnel current, with parameters 2∆ = 0.43 meV, K.A = 144 W.K -4 . Dashed line: model including in addition the Andreev current, but not the related heating, with the parameters D = 80 cm 2 .s -1 , Lϕ =1.5 µm. Dotted lines: full model taking into account the Andreev current and the related heating. Compared to theoretical calculation, the Andreev current was multiplied by 1.55. (d) Dependence of the calculated electronic temperature with the voltage as derived from the fit of the experimental results and for a series of cryostat temperatures: 90 (red symbols), 140 (green) and 230 mK (blue).

Fig. 4

 4 Fig. 4 (a) SEM image of the device. Each of the two normal islands (colorized in red or in green) is inserted between two sets of superconducting junctions. One junction pair is used as a cooler (or heater) and the other one is used as thermometer. (b) Schematic side-view of the set-up. The top (t) and bottom (b) levels are galvanically isolated from each other by a 40 nm thick layer of Si. (c) Heat transfer model. Electrons and phonons of the bottom (top) island at a respective temperature T e,b and T ph,b (Te,t and T ph,t ) exchange a heat power Qb e-ph ( Qt e-ph ). Phonons of each island are coupled together via the Kapitza power Kb-t , and to the bath phonons via Kb-s and Kt-s. (d) Electronic temperatures T e,b and Te,t = T ph,t as a function of the voltage V across the bottom cooler/heater junctions at a bath temperature T bath = 432 mK. (e) Temperature variations ∆T e,b and ∆Te,t = ∆T ph,t based on the same data, the latter temperature change being amplified by a factor 30. (f) Maximum of ∆T e,b (circles) and ∆Te,t (squares) as a function of the bath temperature. (g) Electronic temperatures T e,b (green circles) and Te,t (red diamonds) as a function of the absolute value of the power injected in the bottom island. Grey squares show T e,b plotted as a function of the raw power Q0N , while all other data are plotted as a function of the corrected power QN = Q0N + αQS. The dash-dotted and the dotted lines corresponds to the point where the power Q0N or QN respectively change sign. The full lines are fits calculated using a single set of parameters. The black dotted lines indicate the calculated phonon temperature of the island that is cooled or heated. Bath temperature is 281 mK.

Fig. 5

 5 Fig. 5 (a) From top: fabrication starts with an Al/AlOx/Cu multilayer, on which a photoresist is patterned with contact pads and holes. Then, Cu and Al are successively etched, leaving a suspended membrane of Cu along the line of adjacent holes. A second lithography and etch define the Cu central island. (b) Optical microscope image showing regions by decreasing brightness: bare Al, Cu on Al, suspended Cu and substrate. On the top, two thermometer junctions are added. (c) Colorized scanning electron micrograph of a sample cut using Focused Ion Beam, showing the Cu layer suspended over the holes region. The thickness of Al and Cu is 400 and 100 nm, respectively. (d) Current-voltage characteristic and (e) differential conductance of a sample at different cryostat temperatures.

Acknowledgements This work was funded by the EU FRP7 low temperature infrastructure MICROKELVIN and by Institut universitaire de France. Samples were fabricated at Nanofab platform -CNRS. This work has been led in collaboration with S. Rajauria, L. M. A. Pascal, A. Fay, B. Pannetier, A. Vasenko, T. Crozes and T. Fournier. We acknowledge fruitful discussions along the years with J. P. Pekola, F. Giazotto and M. Meschke.


Ion Candel
email: ion.candel@gipsa-lab.grenoble-inp.fr

Cornel Ioana
email: cornel.ioana@gipsa-lab.grenoble-inp.fr

Bertrand Reeb
email: bertrand.reeb@edf.fr

ROBUST SPARSE REPRESENTATION FOR ADAPTIVE SENSING OF TURBULENT PHENOMENA

Keywords: sparse representation, wide band waveforms, underwater acoustic propagation, turbulence, decomposition

In this paper, we propose a method for turbulence characterization using sparse representation of channel's impulse response.

We consider the case of moving vortices created naturally or artificially that do not conserve their physical properties when observed at two distinct positions in space. Existing amplitude-based techniques fail to provide an accurate representation when the physical properties of the dynamic turbulence are altered. A two stages approach is proposed in this paper. The first one deals with the design of robust waveforms for sensing of turbulent phenomena. The second stage consists in sparsely representing the decomposition of the turbulence's impulse response, based on a physically driven decomposition basis. Tests carried out in a reduced scale experimental facility show, on real data, the efficiency of the turbulence tracking. We compare several types of signals and show that wide band signals are best suited for the application, achieving a high resolution combined to excellent results in terms of robustness.

INTRODUCTION

The dynamics of aquatic phenomena can be understood as the propagation of vortices produced either naturally (natural turbulence embedded in water flow, marine life such as fish and mammals) or artificially (underwater obstacles, submarines, vessels). These vortices are detected and analyzed using pairs of acoustic transducers. If a vortex intersects the path between two pairs of active transducers, the signatures observed on the received signal's amplitudes are quantifiable.

However, specific phenomena associated to the context (attenuation, dispersion, obstacles, etc.) alter the signature when the turbulence travels to another pair of transducers and thus the similarity between the two signatures can be weak.

An alternative is to decompose the impulse responses corresponding to the two transducers pairs on basis of elementary functions. In order to adapt the basis to the context of our application, two concepts are studied in this paper.

First, we propose to design a decomposition basis inspired by the physical characteristics of the vortex "seen' by the impulse response of two propagating acoustic paths.

The second concept is the one of the sparse representation [START_REF] Candés | Decoding by Linear Programming[END_REF]: the analyzed phenomena can be recovered only using some of the decomposition coefficients that are relevant with respect of the minimization criterion, thus obtaining a sparse representation of the phenomena. At this level, we show the interest of the L 1 norm.

The article is structured as follows: section 2 presents an analysis of the best waveform to carry turbulence information, section 3 presents the theoretical concepts used in our work: waveform decomposition and sparse representation. Section 4 presents the main results and remarks of our work and section 5 presents the concluding remarks and future work.

TURBULENT FLOW INTERACTION WITH DIFFERENT WAVEFORMS

The analysis of turbulent phenomena is of major importance in the field of hydraulic systems or monitoring of natural environments. Conventional techniques rely on short sine pulses that are transmitted into the turbulent environment.

The turbulence modifies the amplitudes of pulses in a random way and the analysis of the amplitude modification, when analyzed in bi-static configurations, conducts to the estimation of the turbulence's dynamics [START_REF] Lampa | Turbine Flow Measurement in Intakes with the Acoustic Scintillation Flow Meter[END_REF]. However, these techniques quickly reach their limitations because they provide little information regarding the dynamics of the turbulent flow. It is therefore necessary to construct more advanced waveforms to better understand the evolution of underwater turbulence.

The first research direction consists of adaptive waveform design that will be used for robust sensing of the turbulence. As a reference, typical to existing techniques, the short pulses are used as emitted signal. Two types of signals are then investigated for their potential to deliver a most accurate description of the turbulence dynamics: frequency shift keying (F.S.K.) signals and linear frequency modulations (CHIRPS). These signals are illustrated in the figure 1. In order to illustrate the motivation of our approach, the next example shows the evolution of impulse responses obtained via a matched filtering procedure, corresponding to the three types of signals, when the sensing waveforms interact with a moving vortex. Figure 2 shows that due to the short duration of the pulse and thus low level of transmitted energy, the turbulence signature is reduced and therefore this type of signal is not suited for our application. Little information can be extracted, as most methods rely on the time delay between the emitted and received signal and amplitude changes. In addition to that, the signals can be easily corrupted by unwanted interference (equipment noise, vibrations, echoes, etc.) and therefore the range of applications is limited. In figure 3 however, the results are improved, but a relevant signature is still unavailable. This is due to the short duration of the F.S.K. steps. Initially developed for coding binary information, the F.S.K. signals' capability of encoding turbulence information is limited. In addition, the on/off commutation of transducers in a short time (successive pulses that are injected into the transducers in short amount of time), such as the pulses in the first 120µs in the second subplot of figure one can cause more problems. The acoustic transducers do not support very short energy keying times (of the order of tens of µs) and received signals may be distorted. Also due to the short duration between these steps reflection might occur which render the technique less robust to interference.

Figure 4 illustrates the case of CHIRP signals. A more precise description of the turbulent phenomena is provided as expected. The reason is that the wide band signal's duration is large enough to capture the passage of a turbulent front even in situations where the transmission has a high attenuation (small detail in figure 4). This represents an increase in the sensitivity of the method.

Consider the emission of two wide band signals from figure 1 using two acoustic transducers, s e1 (t) and s e2 (t), in the form of two linear frequency modulations, described by the following equation:

2 e1 e2 0 t s ( t ) s ( t ) A exp j t , 2                       (2.1)
where 0 + γt is the linear frequency variation law. The received signals s r1 (t) and s r2 (t) for each transducer can be expressed, ignoring the noise, as:

r1 r1 e1 r2 r2 e2 s ( t ) A ( t ) s ( t ) s ( t ) A ( t ) s ( t )         (2.2)
where τ is the propagation time of the vortex between the emission and reception transducers. The envelopes A r1 and A r2 in the equations (2.2) contain the information about the turbulence evolution. We calculate the impulse response of the system for each transmission by correlating the received signal with the emitted one [START_REF] Bell | Information Theory and Radar Waveform Design[END_REF]:

r1 r1 e1 e1 r2 r2 e2 e2 h ( t ) A ( t ) s ( t ) s ( t )d h ( t ) A ( t ) s ( t ) s ( t )d                 (2.3)
Figure 5. The signature of turbulence on the envelope of the impulse response.

DECOMPOSITION BY SPARSE REPRESENTATION

The general layout of the analyzed phenomena is illustrated in figure 6. Wide band signals are transmitted using two pairs of sensors placed on the Upstream and Downstream sides. The reason for a second pair of sensor resides in the fact that turbulence dynamics changes with distance. The sources of these changes must be identified (natural energetic decay, artificial and natural sources) or reduced in order to recover almost identical signatures of the same turbulence at two different paths separated by a certain distance.

In the absence of turbulence (such as a vortex), there are no changes in the shape of the impulse responses of the channels defined by the two pairs of transducers. The passage of turbulence between the transducers introduces modifications on the amplitudes of the impulse responses, calculated for each transmission, a fact that can be seen as signatures of the turbulence's source.

However, the vortex can change the characteristics while moving from one pair of sensors to the other, the similarity of the two signatures is low and can be considered as belonging to two different phenomena. The simple observation of the impulse responses is not enough for turbulence dynamics estimation.

Therefore, the turbulence dynamics must be estimated by analyzing the impulse responses using their decomposition based on a physicallydriven dictionary. That is, the decomposition attempts to highlight the common elements that will allow us to identify the turbulence. However, the extraction of this signature must take into account the fact that due to occurring reflection of echoes, a secondary much smaller peak appears. This representation of the signature (amplitude data) cannot reveal accurate amount of information regarding the physical properties of the vortex. Therefore, we propose the analysis via a decomposition using an appropriate dictionary function constructed from the waveform received by the first pair of sensors (upstream acoustic path)

having an amplitude variation similar to the one introduced by the vortex (see figure 7).

The basic principle of the decomposition is illustrated, using real data, in figure 7.

Step one of the sparse representation consists in constructing the elements of the initial dictionary which is the same for both h r1 (t) and h r2 (t), with a zero mean and a rapid decrease. The dictionary starts from an initial waveform 0 , which is a modulated sine wave centered on the resonance frequency of the transducers as illustrated in equation (2.4):

0 0 ( t ) A( t ) exp( j t ),     (2.4)
where 0 corresponds to the resonance frequency of the transducers. The term A(t) represents the amplitude modulation produced by the dynamic vortex over the impulse responses. The choice of this term is based on the physics of the vortex:

from one path to another, the envelopes are similar, but due to the heterogeneity induced by the vortex the shape of the envelope suffers a contraction in time (figure 8).

Figure 8. The contraction of the impulse response envelope due to the dynamic vortex.

Thus, the choice of A(t) is made from the envelopes of the impulse responses that are affected the most by the turbulence, in order to provide a physically adapted dictionary for the decomposition. As the turbulence progresses through water, it keeps on enlarging due to its inertial character [START_REF] Jimenez | The Contributions of A.N. Kolmogorov to the Theory of Turbulence[END_REF] and its energy decreases with distance.

Having defined the initial waveform we construct a set of functions that are shifted in time and scale, according to the continuous wavelet methodology [START_REF] Farge | Wavelet Transforms and Their Applications to Turbulence[END_REF], thus creating a dictionary of elementary functions, at different scales, shifted in time:

s,u 0 1 t u ( t ) , s 0.01 n; n 1...N s s           (2.5)
where N is the maximum decomposition resolution, s and u are, respectively, the scale and shift parameters of the initial waveform 0 . The envelopes of the impulse responses from the Upstream side will be more contracted in time than the ones from the Downstream side. The reason for selecting the initial waveform among the impulse responses where the envelope's contraction is the highest resides in the decomposition process. The impulse responses will match the initial waveform 0 at different scales highlighted by the s parameter in (2.5). The representation that we wish to obtain for h r1 (t) and h r2 (t) can be written, using the resulting coefficients, as follows [START_REF] Baraniuk | Model-Based Compressive Sensing[END_REF] :

N 1 r1 s,u s,u i1 N 2 r2 s,u s,u i1 h ( t ) C ( t ) h ( t ) C ( t )         (2.6)
where and are the N resulting coefficients issued from the scalar product between the signal and each element of the s,u (t) dictionary (basis) for the two impulse responses. The and coefficients represent the projection of the impulse function h r (t) on the s,u (t) set of functions described in (2.4):

i s,u r,i s,u C h ,   (2.7)
The second step of the sparse optimization consists in minimizing the number of coefficients obtained from the decomposition. This is defined as a problem of minimization [START_REF] Lampa | Turbine Flow Measurement in Intakes with the Acoustic Scintillation Flow Meter[END_REF] with solutions derived from convex optimization [START_REF] Boyd | Complex Optimization[END_REF].

In the case of the turbulence estimation, a possible solution to the problem can be formulated as follows [START_REF] Chen | Atomic Decomposition by Basis Pursuit[END_REF]: quantify the impact of each coefficient -basis pair between the two estimated turbulence signatures, in our context, and minimize the number of coefficients needed to construct two signatures ̃ and ̃ with minimum errors.

Based on (2.6), we calculate a residual for each one of the N coefficients (C i,j s,u ) and its corresponding dictionary element:

P i i, j p ri s,u s,u j1 R ( t ) h C , P, j 1...N ; i 1,2         (2.8)
where R P j (t) is the residual calculated for the P th coefficients of the decomposition.

In order to quantify the impact of each residual with respect to the decomposition and minimize the number of coefficients, we can use the L 1 or L 2 norm for each of the N residuals and select only the K coefficients corresponding to the number of coefficients for which the norm is minimum [START_REF] Candés | Stale Signal Recovery from Incomplete and Inaccurate Measurements[END_REF]. For each transmission, we obtain a residual matrix R P containing N residuals calculated form (2.8). The optimum number of coefficients was fund by calculating a L 1 or L 2 norm for each line of the matrix (a residue):

R R

i P 1 1 i P 2 2 K min ( t ) u sin g the L norm K min ( t ) u sin g the L norm       (2.9)
We choose to use the L 1 norm over the L 2 norm because the former leads to fewer coefficients than the latter as proved by the figure 9. In this figure, the plots for the L 1 and L 2 norms are zoomed in around their minima points.

The minimum of the L 1 norm correspond to the lowest number of coefficients used to represent sparsely the two signatures ̃ and ̃ :

P ~1 r1 s,u s,u i1 P ~2 r2 s,u s,u i1 h ( t ) C ( t ) h ( t ) C ( t )         (2.10)
Since the basis of the decomposition is the same for both ̃ and ̃ , the two signatures will have a very strong likelihood, as it will be shown in the next sections.

RESULTS

The tests were carried out in our reduced scale experimental facility (figure 10). For our experiment, we used two pairs of 1 MHz transducers placed on the outer walls of the tank. The length of the acoustic path was 1 meter at a depth of 40 cm. Wide band signals were generated with a linear frequency modulation between 800 kHz and 1.2 MHz, and were downloaded into a signal generator for transmission. The repetition rate of the signals was 1msec and the duration of the wide band signal is of 200µs in order to avoid the overlap of echoes over the received signals. First, the pairs of transducers had to be set apart by a certain distance because of the occurrence of crosstalk caused by the wide beam angle of each transducer. The scale factor of s (decomposition basis parameter) was set to 0.01 due to resolution considerations (the number of elements in the dictionary) and the initial waveform, was selected among the zones extracted from the impulse response where the turbulence had the strongest effect, thus creating a physically driven basis. Each received signal yielded decomposition and a number of K coefficients was calculated using the L 1 norm. It was no surprise that except for the moments where the effects of turbulence were strongest, the number of coefficients used for decomposition remained constant.

We compare our results in two ways: first, we show a significant improvement with a typical representation technique consisting in calculating the maxima values of the impulse response corresponding to each transmission as illustrated in figure 11. This technique is used extensively in certain flow metering applications presented in [START_REF] Candel | Using Wide Band Signals for Obstacle Path Correction in Acoustic Scintillation Flow Meters[END_REF] and [START_REF] Reeb | Investigation of the Performance of Acoustic Scintillation Flow Meter when Turbulence Levels are Low[END_REF] and is sensible to interference induced by the measurement conditions. Secondly, we use the results presented in figure 13 to calculate the average flow velocity between the acoustic paths.

The accuracy of estimating the average velocity depends on the precision of estimating the time delay between the two acoustic paths. This is done by computing the cross correlation of the two signatures (Upstream and Downstream), as illustrated in figure 14 (a zoomed region displaying the main peaks of the two cross correlations).

In the second subplot of figure 13, we can observe an enhanced similarity between the signatures for the two acoustic paths, which describes an initial dynamic and random vortex. This result is possible only because of the physically driven dictionary that was used in the decomposition.

The second value (the real value calculated after the decomposition, highlighted on the red dotted trace with the vertical black line) corresponds to the correct orientation of the flow and a more pertinent 1.47m/s flow velocity is computed.

Using the average flow velocities from multiple levels over a square section output the values of flow on that section. Accurate flow values calculated with existing methods are biased by unfavorable conditions: low turbulence levels [START_REF] Reeb | Investigation of the Performance of Acoustic Scintillation Flow Meter when Turbulence Levels are Low[END_REF] and transducer spacing [START_REF] Marko | Negative Bias in ASFM Discharge Measurement in Short Intakes -Transducer Spacing[END_REF]. Figure 14 shows that, by comparing with the initial estimation in figure 13, the sparse representation manages to provide a correct estimation of a turbulence translated in space. The purpose of sparse representation, as shown in [START_REF] Candés | Stale Signal Recovery from Incomplete and Inaccurate Measurements[END_REF], was to recover a signal using minimum number of coefficients from incomplete and altered measurements.

CONCLUSIONS

In this paper, we show that sparse representation and adaptive waveform design techniques can be successfully used to recover a complete representation of turbulence from incomplete and contaminated measurements (sparse representation). We highlight the importance of adapting the content of the waveforms by studying the physics of the phenomena under measurement. This is due to the failure of a general accepted waveform to provide relevant and accurate information. Our method takes into account the turbulent nature of flow and the changes suffered by a waveform as it interacts with turbulence.

A small scale experiment consisted in generating turbulence was carried out and wide band signals were used to highlight the passage of the turbulence at two separated acoustic paths, all in the least favorable conditions.

Results obtained with our technique prove the theoretical concepts of sparse representation and adaptive waveform using large band signals: more information can be extracted, leading to more accurate results due to the robustness of the processing method.

Our work will focus in the future on combining adaptive waveform techniques with sparse representations aiming at finding a signal that is appropriate to underwater turbulence estimation using an orthonormal decomposition basis.

Figure 1 .

 1 Figure 1. Different types of signals used to characterize turbulent phenomena.

Figure 2 .

 2 Figure 2. Impulse responsespulse (conventional techniques).

Figure 3 .

 3 Figure 3. Impulse responses -FSK signals.

Figure 4 .

 4 Figure 4. Impulse responses -CHIRP signals.

Figure 6 .

 6 Figure 6. Vortex travelling with the flow of water

Figure 9 .

 9 Figure 9. Variation of the L 1 and the L 2 norms.

Figure 10 .

 10 Figure 10. Reduced scale experiment test bench.

Figure 11 .

 11 Figure 11. Turbulence signature for Upstream path.

Figure 12 .

 12 Figure 12. Turbulence signature for Downstream path.

Figure 13 .

 13 Figure 13. Turbulence representation before and after the sparse decomposition.

Figure 14 .

 14 Figure 14. Flow velocity calculation using the cross correlation between the acoustic paths.