Senior Member, IEEE Maria Pia Fanti
email: mariapia.fanti@poliba.it

Member, IEEE Gregory Faraut
email: gregory.faraut@lurpa.ens-cachan.fr

Member, IEEE Jean-Jacques Lesage
email: jean-jacques.lesage@lurpa.ens-cachan.fr

Michele Roccotelli
email: michele.roccotelli@poliba.it

D M P Akopian

J.-J Lesage

An Integrated Framework for Binary Sensor Placement and Inhabitants Location Tracking

Keywords: AAL), binary sensors, integer linear programming (ILP), smart home. I

This correspondence paper deals with the sensor placement optimization problem in the context of indoor multiple inhabitants location tracking to solve ambient assisted living problems. Binary sensors, like passive infrared (PIR) sensors, are used to guaranty specific coverage requirements and allow privacy respecting. Moreover, within real home environments, different kinds of obstacles (like walls, high furniture, etc.) can affect the detection capacity of PIR sensors. This paper proposes an integrated framework devoted to optimize the placement of sensors and PIR sensors in smart homes by taking into account physical topologies and coverage precision constraints. An integer linear programming problem is formalized and a case study illustrates the applicability of the proposed approach and the scalability of the optimization method.

Index Terms-Ambient assisted living (

of health problems or to improve the environment comfort [START_REF] Danancher | Online location tracking of a single inhabitant based on a state estimator[END_REF], [START_REF] Danancher | Model-based location tracking of an a priori unknown number of inhabitants in smart homes[END_REF]. For these purposes, it is relevant to optimize the number and placement of sensors in smart homes in order to improve the location tracking performance [START_REF] Feng | Optimal design of infrared motion sensing system using divide-and-conquer based genetic algorithm[END_REF]- [START_REF] Vlasenko | The smart-condo: Optimizing sensor placement for indoor localization[END_REF] as well as minimizing the cost of instrumentation.

In this correspondence paper, an integrated framework devoted to provide the designers with the best solution for the binary sensor placement is proposed. The system receives in input, the data about the possible sensor choices and the apartment topology, and optimizes the placement of the binary sensors like passive infrared (PIR) sensors, by considering different aspects: 1) the obstacles, like walls, that can limit the sensors detection; 2) a topological classification of the zones in which the tracking has to be performed; and 3) coverage and precision constraints of the location. In particular, the dwelling is divided into independent zones on the basis of the presence of walls and, for each independent zone, an independent sensor placement optimization problem is formalized as an integer linear programming (ILP) problem. Moreover, the application of the presented method to a case study, including different scenarios-in terms of precision and coverage of the location-is presented and discussed.

This paper is organized as follows. Section II discusses the related works. Section III presents the structure and the components of the integrated system. Moreover, Section IV describes the phases of the optimization strategy for sensors placement and Section V shows the effectiveness of the proposed method by considering a real flat topology. Finally, Section VI summarizes the conclusions and gives some perspectives for further works.

II. RELATED WORKS

In the literature related to AAL problems, some studies propose several technological solutions to realize different applications mainly related to activities of daily living (ADLs) recognition and location tracking. Rashidi and Mihailidis [START_REF] Rashidi | A survey on ambient-assisted living tools for older adults[END_REF] summarized the state-of-the-art of AAL technologies, tools, and techniques. They also explore successful case studies and deployed systems and discuss the present and future challenges. Magherini et al. [START_REF] Magherini | Using temporal logic and model checking in automated recognition of human activities for ambient-assisted living[END_REF] described a system that uses propositional temporal logic and model-checking to support automated real-time recognition of ADLs within a smart environment. The effectiveness of the proposed approach is shown by applying it to a smart kitchen. Moreover, Andò et al. [START_REF] Andò | An event polarized paradigm for ADL detection in AAL context[END_REF] explored the possibility to detect ADLs, with a particular focus on fall detection and present a smartphone-based ADL detector, which adopts an advanced classification paradigm to distinguish between different kinds of falls. Nag and Mukhopadhyay [START_REF] Nag | Occupancy detection at smart home using real-time dynamic thresholding of flexiforce sensor[END_REF] dealt with the problem of correctly recognizing the output of flexi-force sensors that may vary depending on the time usage. These sensors are typically used to monitor the occupancy or nonoccupancy status of chairs, sofa, and bed. Nag and Mukhopadhyay [START_REF] Nag | Occupancy detection at smart home using real-time dynamic thresholding of flexiforce sensor[END_REF] proposed a real-time method based on dynamic thresholding that allows differentiating the two situations even in case of long-time usage. Furthermore, Zhou et al. [START_REF] Zhou | A case-driven ambient intelligence system for elderly in-home assistance applications[END_REF] presented a case-driven ambience intelligence (C-AmI) system that aims to manage elderly people ADLs at home. The C-AmI system architecture integrates wireless sensor network, activity recognition, case-based reasoning, elderly in-home assistance customized knowledge within a coherent framework. It supports the sensing, predicting, and reasoning of assistive actions by exploiting the a priori knowledge about solved cases to meet the aged people's needs in a domestic environment.

Other related works in the AAL context refer specifically to the indoor inhabitants location tracking and different technologies and methodologies are used to locate the position of humans within a smart environment. Some authors develop tracking systems based on invasive sensors like wearable and/or cameras devices [START_REF] Mehdi | A survey of human location estimation in a home environment[END_REF]- [START_REF] Koch | Indoor localisation of humans, objects, and mobile robots with RFID infrastructure[END_REF]. In this context, Tabar et al. [START_REF] Tabar | Smart home care network using sensor fusion and distributed vision-based reasoning[END_REF] proposed a wearable identification badge, to localize a person in a flat. The user badge is equipped with two accelerometers for measuring sudden movements representing the fall of an elderly person and the approximate location of the user is computed by evaluating the signal strength measurement in the home environment. Furthermore, a distributed vision-based system for managing the accident in the assisted living environment has been developed in [START_REF] Aghajan | Distributed vision-based accident management for assisted living[END_REF]. The posture of people is detected by merging features from different camera images in order to reduce false alarms. In [START_REF] Schindler | A wearable interface for topological mapping and localization in indoor environment[END_REF], a wearable ear-mounted device for localizing an inhabitant in the home environment has been developed. It is composed of an infrared proximity sensor and a dual axis accelerometer. When the person walks wearing the device, the accelerometer detects the footsteps while sometimes the proximity sensor fails to recognize the doorways. The authors remark that sometimes the proximity sensor fail to recognize the doorways accurately and, as a result, the human localization is not precise.

Wearable sensor systems can localize inhabitant positions with good accuracy but several drawbacks are derived from their use. Elderly people can forget to follow the instructions and do not like to carry special devices with them and there are also privacy issues. Moreover, many of these systems require complex infrastructure to be installed at home for inhabitant localization.

Some authors deal with the localization of inhabitants by using mobile robotic systems [START_REF] Mehdi | A survey of human location estimation in a home environment[END_REF], [START_REF] Koch | Indoor localisation of humans, objects, and mobile robots with RFID infrastructure[END_REF], [START_REF] Volkhardt | Playing hide and seek with a mobile companion robot[END_REF]. Such systems are invasive and the main results of the presented methodologies are derived from simulation studies while the real life experiments are still limited.

In order to recognize the inhabitant activities in real time, by using cheap sensors and without violating the privacy, PIR motion sensors can be used. In this context, Yang et al. [START_REF] Yang | Distributed multi-human location algorithm using naive bayes classifier for a binary pyroelectric infrared sensor tracking system[END_REF] presented a distributed human location procedure for a binary pyroelectric infrared sensor tracking system and propose a two-level regional location, with static and dynamic partitioning. Crandall and Cook [START_REF] Crandall | Tracking systems for multiple smart home residents[END_REF] used multiple PIR sensors and two location tracking algorithms. Crandall and Cook [START_REF] Crandall | Tracking systems for multiple smart home residents[END_REF] remarked that it is needed to install a dense grid of sensors for the full coverage of the environment but anyway in some situations the human activity is not monitored correctly by one or more sensors resulting in less accurate results.

Hence, it appears that the study of the optimal design of infrared motion sensing systems for human motion localization is an open problem. In particular, Feng et al. [START_REF] Feng | Optimal design of infrared motion sensing system using divide-and-conquer based genetic algorithm[END_REF] aimed at finding the optimal number and placement of PIR sensors in order to improve the localization performance. A multiobjective, mixed-integer-discretecontinuous optimization problem is presented and solved by using a divide and conquer-based genetic algorithm method. Moreover, Emmons and Kamangar [START_REF] Emmons | Evaluating the optimal placement of binary sensors[END_REF] evaluated the optimal placement of binary sensors by providing several definitions for sensor coverage and spatial probability and by focusing on 1-D sensor arrays. A metric for measuring the overall utility of a given set of choices is presented and includes the precision evaluation of measurement and coverage, and the associated costs. In addition, Zhao et al. [START_REF] Zhao | Optimal deployment of mobile sensors for target tracking in 2D and 3D spaces[END_REF] faced the problem of finding the optimal sensor placement that can minimize the target tracking uncertainty by studying a control strategy that provides a flexible solution to autonomous optimal sensor deployment in 2-D and 3-D spaces.

In the previously mentioned works, different methods to solve the problem of sensors placement in a given space (1-D, 2-D, or 3-D) for location tracking purposes are developed but they do not solve the problem of multiple inhabitants location tracking in smart homes by considering real environments. In this context, Gaddam et al. [START_REF] Gaddam | Selection and optimization of wireless sensors in a smart digital home for the elderly[END_REF] presented a framework to care elderly people in smart home based on wireless sensors. The sensor selection architecture is performed by an evaluation module that chooses the appropriate sensor nodes on the basis of performance requirements for observability, fault detection, fault robustness, and cost. However, Gaddam et al. [START_REF] Gaddam | Selection and optimization of wireless sensors in a smart digital home for the elderly[END_REF] do not show the applicability of the method on a real home environment. On the other hand, Vlasenko et al. [START_REF] Vlasenko | The smart-condo: Optimizing sensor placement for indoor localization[END_REF] proposed a hardware/software platform that aims to help inhabitants in performing several daily home activities. The goal is to recognize the general activities in real time, by respecting the privacy and developing a method for planning the optimal deployment of binary PIR motion sensors. Vlasenko et al. [START_REF] Vlasenko | The smart-condo: Optimizing sensor placement for indoor localization[END_REF] optimized the sensors placement on the basis of the zones visitation frequency by considering only one inhabitant. In [START_REF] Danancher | A discrete event model for multiple inhabitants location tracking[END_REF]- [START_REF] Danancher | Model-based location tracking of an a priori unknown number of inhabitants in smart homes[END_REF], a method that aims at the systematic construction of a discrete event system (DES) model for real-time location tracking of one or multiple inhabitants into a smart home is proposed. In addition, a location tracking algorithm as well as a procedure to evaluate the relevance of a given instrumentation for location tracking, are presented. Several criteria are defined to evaluate a priori the performance of the chosen instrumentation and to help the expert for manually improving it. In particular, one of these criteria measures the inability of the sensors network to detect the presence of inhabitants in certain zones of the flat, defined as unlocationable zones. On the basis of the defined criteria, the performance of the location tracking algorithm is evaluated by an expert that can decide if a modification of sensors in terms of type, number, and placement is necessary.

In order to overcome the problem of manually modifying the sensors network by means of an expert, this paper starts from the contributions in [START_REF] Danancher | A discrete event model for multiple inhabitants location tracking[END_REF]- [START_REF] Danancher | Model-based location tracking of an a priori unknown number of inhabitants in smart homes[END_REF] and provides an automatic procedure to optimize the placement of the sensors by guaranteeing the respect of the coverage requirements.

III. STRUCTURE OF THE INTEGRATED FRAMEWORK

FOR SENSOR PLACEMENT In this section, we describe the proposed methodology devoted to optimize the sensor placement in a generic flat topology. We present an integrated system that can be used by the designers to select the best solution for the placement of PIR sensors in order to satisfy suitable constraints and optimize a defined objective function. In particular, we are not interested in determining the identity of the inhabitants and the kinds of activities they are doing, but we aim to track the location of inhabitants in the flat with a specified accuracy. For this purpose, we are interested in optimizing the sensing system placement to reach the desired location tracking performances.

The integrated system is defined by considering the following assumptions: 1) binary sensors, like PIR sensors, are used because they are low cost, and ensure the privacy respect and 2) the living space is discretized by means of bidimensional grids because the placement of sensors and the localization of inhabitants in 2-D spaces are required for AAL purposes [START_REF] Hiromori | Sensor placement optimization method for people tracking[END_REF]- [START_REF] Younis | Strategies and techniques for node placement in wireless sensor networks: A survey[END_REF].

Moreover, the proposed integrated system structure is constituted by two main components: 1) the data component, which includes In particular, the model component is the core of the integrated system and is specified by the following three phase procedure.

1) Classification Phase: The zones of the dwelling that are limited by walls or obstacles (in the most cases the zones are the rooms) are classified by means of a strict order inclusion relation, which is graphically represented by a topological tree. This procedure determines a set of closed zones: we define closed zone a part of the flat where the optimization can be separately applied. 2) Discretization Phase: For each closed zone, two discretization scales are defined according to the required accuracy for sensor placement and location tracking: a) a first scale defines a grid for sensor placement and b) a second scale defines a grid for zone coverage. 3) Optimization Phase: An independent optimization problem is defined for each closed zone. Fig. 1 gives an overview of the optimization framework structure. The main components of the system and the relationships among them are pointed out. The outputs are the solutions provided by the optimization problems that give the optimal choices of type, position, and number of the PIR sensors for each closed zone.

IV. DATA AND MODEL COMPONENTS OF

THE INTEGRATED FRAMEWORK This section specifies the data and the model components of the presented system and explains in detail the proposed procedure to optimize the sensor placement.

A. Data Component

PIR sensors belong to a set of types with cardinality S, where each element h ∈ H denotes a motion detector type. Moreover, each sensor type h is characterized by the following attributes. A detection angle da h that can be da h = 180 • or da h = 360 • , a detection radius dr h ∈ N + ; a cost c h ∈ R + . In particular, the PIR sensors can be placed in different zones on the basis of their types. For example, sensors with da h = 180 • can be positioned on the walls, and sensors with da h = 360 • can be placed on ceilings. As an example, Fig. 2(a) shows the geometric definition of the coverage zone of a PIR sensor type with da h = 180 • and dr h = R. The presence of a wall into this zone constitutes an obstacle, which generates an uncovered zone. Fig. 2(b) provides a discretization of the same coverage zone for a given coverage scale. Moreover, the closets, the doors and the windows can also be considered as obstacles for the motion detectors within the domestic environment.

B. Classification Phase

During the classification phase the topology of the flat is considered in order to define and classify the set of closed zones by means of a topological tree expressing an inclusion relationship between zones.

Here, we define as closed a zone whose space is limited by walls or materials, which constitute obstacles to the signals emitted by PIR sensors.

Each node of the tree represents a closed zone that may or not include smaller closed zones. The whole flat can be decomposed in closed zones and at the end of the decomposition the leaves of the tree have to be considered for the independent optimization application.

The procedure to build the tree is as follows.

1) The root node is the flat and is placed at the level 1, i.e., the highest level.

2) The main zones Zc i of the flat limited by the main walls are singled out and form the nodes of the tree second level. 3) Each zone Zc i is inspected in order to determine other zones Zc i,j limited by cupboards, furniture, partitions, etc. 4) The procedure goes to an end when no other area is highlighted. By adopting the proposed zone classification, an independent sensor placement optimization problem can be defined for each leaf of the tree with different granularity for each level and each zone. The following example clarifies the proposed classification phase.

Example 1: Let us consider the flat shown in Fig. 3 and composed by the following rooms. An open space including the living room (5.5 × 5.2 m), the kitchen and the dining room (3.6 × 4.2 m); a laundry room and a pantry room; the bathroom; the master suite zone including a master suite (4.2 × 5.2 m), a bathroom and a dressing room; a second bedroom (4 × 3.4 m); a garage (7 × 6.5 m).

Considering the topology of the flat described in Fig. 3, internal walls of the main rooms define a set of closed zones as represented in Fig. 4. Starting from level 1 (the flat) the four main closed zones belonging to level 2 are singled out: the zone Zc 1 , the master suite zone Zc 2 , the second bedroom Zc 3 , and the garage Zc 4 (see Fig. 5).

Moreover, the following closed zones of level 3 are included in zone Zc 1 : the main open space that is represented by the node Zc 1.1 ; the second bathroom that is represented by the node Zc 1.2 ; the laundry and pantry rooms zone that is identified by the node Zc 1.3 . In addition, the following closed zones of level 3 are included in zone Zc 3 : the master suite (node Zc 3.1 ); the dressing room (node Zc 3.2 ); and the bathroom (node Zc 3.3 ).

Four closed zones compose level 4: the laundry room (Zc 1.3.1 ) and the pantry room (Zc 1.3.2 ) that are included in Zc 1.3 ; the bath

C. Discretization Phase

In order to model the optimization problem by an ILP formulation, it is necessary to discretize the flat area. Two independent discretization grid scales are defined for each zone: 1) the grid scale g c divides the zones on the basis of the user requirements in term of location accuracy and 2) the grid scale g p defines the granularity of the sensors placement and has to be fixed by considering the technology of the sensors, especially, the detection radius. Note that it is possible to define different grid scales, for both g c and g p , for each of the 2-D: for the sake of the simplicity in the rest of this paper the same scales are considered for the 2-D. cretizing the considered zone within a 2-D grid scale g p that is suitable for defining the precision of sensor placement. 4) P h = { j ∈ M| j = 1, . . . , P h }: P h ⊂ M, set of grid points where sensors of type h can be placed. 5) P = ∪ h∈H P h : Set of grid points (with cardinality P) where sensors can be placed. Note that the discretization can be obtained by standard techniques [START_REF] Zhao | Optimal deployment of mobile sensors for target tracking in 2D and 3D spaces[END_REF].

D. Optimization Problem Formulation

In this section, we formalize the sensor placement optimization problem for each closed zone by an ILP problem. In the proposed optimization framework, the objective is to optimize the placement of the sensors, by minimizing costs and satisfying the coverage requirements.

The grid coverage is denoted by the following binary variables, given as input data for the optimization problem:

g ijh = ⎧ ⎨ ⎩ 1, if sensor of type h that is placed at point j covers point i 0, otherwise ∀h ∈ H, j ∈ P h , i ∈ N /O.
The values of g ijh are determined in Algorithm 1.

It is evident that g ijh = 0, ∀ i ∈ O with h ∈ H and j ∈ P h . The procedure considers for each sensor of type h a possible position j ∈ P h (step 2). Then, each point i ∈ N /O of the grid is selected (step 3) in order to evaluate the distance ij . If it holds ij ≤ dr h and the segment ij does not intercept any obstacle, then the considered sensor can cover point i (step 4) and it holds g ijh = 1.

Moreover, the decision variables are the binary variables x jh defined as follows:

x jh = 1, if sensor of type h is placed at point j 0, otherwise ∀h ∈ H, j ∈ P h , i ∈ N /O.
The goal of the problem is minimizing the costs of the sensors that have to cover the grid points i ∈ N /O according to specific requirements and by excluding the obstacles. To this purpose, we introduce the following parameters.

a i ∈ N + if i ∈ N /O
, where a i represents the minimum number of sensors that have to cover the grid point i that is not an obstacle.

a i = 0 if i ∈ O, i.e.
, i is an obstacle. Now, the ILP problem is formalized as follows:

min P j=1 S h=1 c h x jh (1a) subject to P j=1 S h=1 g ijh x jh ≥ a i , ∀i ∈ N /O. (1b)
Equation ( 1a) is the objective function that minimizes the sum of the sensor costs by considering the grid points j ∈ P available for sensor placements. Constraints (1b) impose that the minimum number of sensors that have to cover the grid point i ∈ N /O is equal to a i . The ILP problem (1a), (1b) is known in the related literature as the set multicover problem that is NP-hard [START_REF] Chekuri | On the set multicover problem in geometric settings[END_REF]. Even if for the dimension of the considered problems the exact solution is obtained in reasonable time, several algorithms are available in the literature and allow solving large problems [START_REF] Hua | Exact algorithms for set multicover and multiset multicover problems[END_REF].

V. CASE STUDY

In this section, the application of the optimization method to the flat of Example 1 and shown in Figs. 2 and4 is discussed by considering four scenarios where different closed zones are considered. Each scenario represents a particular situation that is analyzed to show the powerful of the proposed method by considering suitable zones under specified constraints.

1) In scenario 1, Zc 1.1 is considered to show the covering of a zone by properly placing the sensors in presence of permanent walls. 2) In scenario 2, Zc 2 is considered to show how to cover a zone in presence of open windows and open doors that can be obstacles for the sensor detection. 3) In scenario 3, the ability of the methodology to take into account the sensor failures is studied by considering zone Zc 1.1 . In the considered case study, the database of sensors includes five types of motion detectors: H = {h = 1, . . . , 5}. Table I reports, for each sensor type the corresponding attributes: note that sensors of type h = 1, 2, 3 can be positioned on the wall and sensors of type h = 4, 5 can be positioned on the ceiling. Moreover, two 2-D grids are defined for each zone and identified through their scales that are reported in Table II. A grid of scale g p is determined for the sensor placement and a grid of scale g c defines the coverage precision.

In such a way, we show that it is possible to choose the values of the grid scales on the basis of the accuracy of placement and coverage precision requirements. In the considered scenarios, the discretization of the grids is determined as follows.

Let us consider a rectangular zone (A × B) and a grid scale of size g. The resulting grid set N = {i = 1, . . . , N} is composed by N = ( A/g + 1) * ( B/g + 1) points and each point i corresponding to the grid position (x i , y i ) is determined by the following formula:

i = y i /g × ( B/g + 1)+x i /g + 1, for i = 1, . . . , N (2)
where α denotes the floor function of α, i.e., the largest integer not greater than the real number α.

Example 2: Consider the rectangular zone (A × B) = (2×1) m 2 depicted in Fig. 6 where a grid scale of g = 1 m is applied. The resulting grid is composed by 6 points N = {i = 1, . . . , 6} and each point corresponds to the grid position (x i , y i ) determined according to [START_REF] Danancher | Online location tracking of a single inhabitant based on a state estimator[END_REF] and shown in Fig. 6.

In the following sections, we solve the ILP problem (1a), (1b) by considering the three selected scenarios. The ILP problems are solved by using a standard solver, i.e., GNU linear programming kit [31] running on an Intel-Core i5-2450M, CPU at 3.40 GHz, with 6 GB RAM. Moreover, the solution of the problems is obtained within less than 30 min.

A. Scenario 1

In scenario 1, the closed zone Zc 1.1 is considered. The objective is to optimize the placement of sensors by imposing that at least one sensor covers each zone not representing obstacles, i.e., a i =

1 if i ∈ N /O and a i = 0 if i ∈ O.
Assuming that the walls and the closets are obstacles and using the grid scales reported in Table II, the zone Zc 1.1 is discretized as

B. Scenario 2

In scenario 2, the closed zone Zc 2 is considered and includes two cases. In case 2a, the door and window are supposed to be closed and are no obstacles; in case 2b, the door and the window may be closed or opened, i.e., they have to be considered as possible obstacles. Also in this scenario, we impose that at least one sensor covers each zone not representing obstacles, i.e., a i = 1 if i ∈ (N /O) and a i = 0 if i ∈ O. The grids reported in Table I result in the following discretization of zone Zc 2 : N = M = 378, P = 308.

In case 2a, the ILP problem provides the following set of optimal solutions, composed of sensors of type h = 1, very close together and exhibiting the same value of the objective function: S(2a) = {(x s , y s ) = (10 + s, 0) for s = 0, . . . , 4}. Furthermore, the ILP problem solution of case 2b provides a set of 20 optimal solutions and each solution is composed of two sensors: 1) one sensor of type h = 1 and 2) one sensor of type h = 4. The set of the solutions is the following (sensors of types 1 and 2 are placed, respectively, at (x sa , y sa ) and (x sb , y sb ), for s = 1, . . . , 20): In particular, Fig. 9 shows the first solution of set S(2b).

S(2b) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (x 1a , y 1a ) = (

C. Scenario 3

Scenario 3 considers the closed zone Zc 1.1 and imposes that at least two sensors cover each grid point i not belonging to an obstacle, i.e., a i = 2 if i ∈ N /O and a i = 0 if i ∈ O. Using the grids shown in Table I, the zone is discretized as follows: N = M = 483, P = 270.

The ILP solution returns only one optimal solution that includes four sensors: one sensor of type h = 1 is placed at the grid coordinates (x, y) = [START_REF] Hiromori | Sensor placement optimization method for people tracking[END_REF][START_REF] Andò | An event polarized paradigm for ADL detection in AAL context[END_REF] and three sensors of type h = 4 are placed at the grid coordinates (x, y) = (7, 16), (x, y) = (13, 7), and (x, y) = (13, 12), respectively.

The placements of the sensors are shown in Fig. 10. In this case, it is possible to perform a location tracking even in case of a sensor failure.

VI. CONCLUSION

This correspondence paper presents an integrated framework to optimize the sensor placement in a multiple inhabitant locationtracking by using binary sensors. The integrated system is composed of a database describing the flat and the sensors and a model component that provides the optimal sensor placements.

To this aim, the flat is divided into a set of independent closed zones and the proposed method optimizes the placement of sensors within each zone by minimizing the sensor costs. The proposed method allows an accurate placement of sensors and a complete coverage of each zone by discretizing the space in bidimensional grids with different scales and by considering the problem of obstacles for motion detectors. On the basis of a real smart home topology, three different optimization scenarios are presented and discussed to demonstrate the effectiveness of the optimization method that may provide different alternative optimal solutions that the user can choose. The solutions presented in the considered case study show the benefit of the proposed method: it is possible to cover the flat with the minimum number of sensors by guaranteeing full coverage and addressing sensor failures.

Future research will consider other kinds of sensors in order to improve the accuracy of the coverage and the performance of the location tracking.

Fig. 1 .

 1 Fig. 1. Main components of the integrated framework.

Fig. 2 .

 2 Fig. 2. Motion detectors coverage including obstacles. (a) Coverage zone. (b) Coverage points on the grid.

Fig. 3 .

 3 Fig. 3. Flat topology.

Fig. 4 .

 4 Fig. 4. Closed zones definition.

Fig. 5 .Algorithm 1

 51 Fig. 5. Topological tree.Algorithm 1 Grid Coverage DeterminationStep0 h = 1 Step 1 Set N s = N /O, Ps = P h Step 2 Select a point j ∈ Ps, set Ps = Ps -{j} Step 3 Select a point i ∈ N s, set N s = N s -{i}Step 4 Determine the Euclidean distance between i and j denoted by ij : if ij ≤ dr h and segment ij does not intercept any obstacle then set g ijh = 1, otherwise set g ijh = 0. Step 5 If N s = ∅ then go to Step 3 Step 6 If Ps = ∅ then go to Step 2 Step 7 h = h + 1. If h ≤ S then go to Step 1 Step 8 End

Fig. 6 .

 6 Fig. 6. Two-dimensional grid description of a simple zone (2 m × 1 m).

Fig. 7 .

 7 Fig. 7. ILP solution for Zc 1.1 . (a) Zone covered by a sensor of type h = 1; (b) Zone covered by a sensor of type h = 4.

Fig. 8 .

 8 Fig. 8. ILP solutions for Zc 2 in case 2a. (a) Zone of optimal solutions. (b) First solution (x 1 , y 1 ) = (10, 0), of the set S(2a).

Fig. 9 .

 9 Fig. 9. ILP solution of case 2b. (a) Sensor of type h = 1 in (x 1a , y 1a ) = (10, 0). (b) Sensor of type h = 4 in (x 1b , y 1b ) = (6, 6).

Fig. 8 (

 8 Fig. 8(a) shows the zone of the possible locations of the sensors and Fig. 8(b) shows the first solution of set S(2a).

Fig. 10 .

 10 Fig. 10. ILP solution for Zc 1.1 . (a) Zone covered by the sensor of type h = 1. (b)-(d) Zone covered by the three sensors of type h = 4.

  10, 0) and (x 1b , y 1b ) = (6, 6) (x 2a , y 2a ) = (10, 0) and (x 2b , y 2b ) = (6, 7) (x 3a , y 3a ) = (20, 6) and (x 3b , y 3b ) = (6, 8) (x 4a , y 4a ) = (12, 0) and (x 4b , y 4b ) = (3, 2) (x 5a , y 5a ) = (13, 0) and (x 5b , y 5b ) = (4, 3) (x 6a , y 6a ) = (14, 0) and (x 6b , y 6b ) = (5, 4) (x 7a , y 7a ) = (15, 0) and (x 7b , y 7b ) = (5, 3) (x 8a , y 8a ) = (16, 0) and (x 8b , y 8b ) = (3, 1) (x 9a , y 9a ) = (17, 17) and (x 9b , y 9b ) = (3, 3) (x 10a , y 10a ) = (20, 10) and (x 10b , y 10b ) = (3, 4) (x 11a , y 11a ) = (20, 3) and (x 11b , y 11b ) = (3, 5) (x 12a , y 12a ) = (20, 15) and (x 12b , y 12b ) = (4, 4) (x 13a , y 13a ) = (20, 14) and (x 13b , y 13b ) = (4, 5) (x 14a , y 14a ) = (20, 12) and (x 14a , y 14a ) = (5, 7) (x 15a , y 15a ) = (20, 7) and (x 15b , y 15b ) = (3, 6) (x 16a , y 16a ) = (16, 17) and (x 16b , y 16b ) = (4, 7) (x 17a , y 17a ) = (20, 13) and (x 17b , y 17b ) = (4, 6) (x 18a , y 18a ) = (20, 8) and (x 18b , y 18b ) = (5, 5) (x 19a , y 19a ) = (20, 9) and (x 19b , y 19b ) = (5, 6) (x 20a , y 20a ) = (20, 1) and (x 20b , y 20b ) = (4, 2)⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ .

This work was supported by the Italian Minister Ministero dell'Istruzione, dell'Università e della Ricerca (MIUR) in the "Smart Cities and Communities and Social Innovation" section, in the framework of project "ASMARA-Applications post Directive 2010/65 in Italian port realities of the Suite MIELE to support the Authority to optimize the inteRoperability in the intermodAlity of city-port flows" SCN_00529. This paper was recom-


J Rech

T Jonckheere

T Martin

B Douçot

D Feinberg

R Mélin

A proposal for the observation of nonlocal multipair production: the biSQUID

Keywords: numbers: 74.50.+r, 74.78.Na 85.25.Dq 74.45.+c, 73.63.Kv

We propose an all-superconducting three-terminal setup consisting in a carbon nanotube (or semiconducting nanowire) contacted to three superconducting leads. The resulting device, referred to as a "biSQUID", is made of four quantum dots arranged in two loops of different surface area. We show how this biSQUID can prove a useful tool to probe nonlocal quantum phenomena in an interferometry setup. We study the measured critical current as a function of the applied magnetic field, which shows peaks in its Fourier spectrum, providing clear signatures of multipair Josephson processes. The device does not require any specific fine-tuning as these features are observed for a wide range of microscopic parameters -albeit with a non-trivial dependence. Competing effects which may play a significant role in actual experimental realizations are also explored.

I. INTRODUCTION

Multiterminal setups offer a great way to explore nonlocal quantum effects as well as entanglement in condensed matter devices, attracting both experimental and theoretical attention. Recently, these efforts have focused on all-superconducting hybrid structures, involving quantum dots or metallic islands connected to multiple superconducting leads. [1][2][3][4][5] In particular, the prediction of nonlocal quartet production, 6 that is the emission of spatially correlated pairs of Cooper pairs, opens new perspectives in the realization of electronic entanglers.

Motivated by this, a recent work 7 considered an allsuperconducting bijunction consisting of a central superconductor coupled via gate-controllable quantum dots to two lateral voltage-biased superconductors. It uncovered the presence of a coherent transport mechanism away from equilibrium, manifesting as multipair phasecoherent Josephson resonances in the current, appearing on top of the usual local dissipative transport of quasiparticles. Similarly, experimental work performed on a three-terminal voltage-biased Josephson junction involving a central T-shaped metallic region revealed features in the electronic subgap transport consistent with the production of nonlocal quartets 8 though an unequivocal signature is still lacking at the moment. In this context, alternative ways of detecting such nonlocal multipair processes are highly desirable, in particular in the coherent dissipationless regime, a route that could be provided by an interferometric setup, such as a SQUID.

A superconductor quantum interference device, or SQUID, consists in a superconducting loop, defining two paths each interrupted by a Josephson junction. As a direct consequence of phase coherence, the Cooper pairs flowing along these two paths interfere, in a way that is controlled by the magnetic flux through the loop. In a re- cent achievement of molecular electronics, a carbon nanotube SQUID has been realized experimentally. 9 There, the two constitutive Josephson junctions are both made of a nanotube quantum dot, allowing to tune their transparency with the help of an external gate voltage, therefore providing a new generation of versatile sensors. Interestingly, similar carbon nanotube (or nanowire) devices have been used as Cooper pair beam-splitters, [10][11][12][13][14] a source of entangled electron pairs [15][16][17] where the quasi one-dimensional nanostructure is contacted to a central superconducting source and two metallic drains. These rely on a specific nonlocal process, referred to as CAR for crossed Andreev reflection, [18][19][20][21][22] which amounts to separating the constituents of a Cooper pair into spatially distinct contacts (provided they are within a distance set by the coherence length).

It is therefore only natural to devise a setup bringing together the interferometric properties of the SQUID with the nonlocal aspects associated with CAR processes. Such a setup could be constructed from a carbon nanotube with readily accessible technology, as illustrated by the many examples of multiterminal nanotube-based devices now available. [23][24][25][26] In this paper, we propose such a device, dubbed a biSQUID, consisting of two twinned SQUIDs, where the superconducting contact which emits the multiple pairs is common to both loops in order to reveal multipair processes through interferometry. The two loops realize a three-terminal structure made of a single carbon nanotube (or nanowire) contacted to three superconducting leads (see Fig. 1), delimiting four quantum dots which are controlled by external gate voltages. The device is ideally operated by fixing the total current flowing through, which we express in terms of the magnetic fluxes piercing the two loops.

Recently, a two-terminal SQUID geometry involving two loops was realized experimentally using superconductor -normal metal -superconductor junctions. 27 While somewhat similar in spirit to the present proposal, this setup was specifically designed for sensing applications rather than to explore nonlocal effects. It shows, however, that not only the actual realization of these systems is within our grasp but also that such superconducting nanodevices constitute an active field of research.

The outline of the paper is as follows. In section II, we introduce the setup and present a first simple phenomenological approach. We then derive the expression for the critical current in section III, starting from a microscopic description of the setup. Section IV is devoted to our results, where we show and analyze the dependence of the critical current on the external magnetic field, and comment on the robustness of the observed features. Finally, in section V, we explore the competing effects which might spoil the expected signatures of pure CAR processes, before concluding in section VI.

II. SETUP AND PHENOMENOLOGY

The biSQUID setup is presented in Fig. 2. It basically amounts to twinning two nanotube-based SQUIDs by a common central electrode whose width is smaller than the coherence length, in order to support nonlocal Andreev scattering processes. 29 A carbon nanotube (or nanowire) is contacted with three superconducting electrodes, referred to as S a , S 0 and S b . This defines four quantum dots, labeled a1, a2, b1 and b2, controllable via nearby electrostatic gates. The area of the two resulting SQUID loops (denoted A and B respectively) are chosen different, so as to ensure that electrons feel different magnetic fluxes depending on the loop they flow through.

In order to make the upcoming discussion as clear as possible, it is important to properly define the various su-perconducting phases involved in this setup. A phase ϕ s is attributed to the central superconducting finger, while the other two superconducting electrodes are characterized respectively by a phase ϕ a and ϕ b . Furthermore, the magnetic fluxes through the two loops A and B are defined as Φ A = BS A and Φ B = BS B , where S j is the surface area enclosed by loop j.

It follows from flux quantization that the phase differences seen by each quantum dot Josephson junction, are given by

δϕ a1 = ϕ s -ϕ a (1) δϕ b1 = ϕ s -ϕ b ( 2
)
δϕ a2 = ϕ s -ϕ a -2π Φ A Φ 0 ( 3
)
δϕ b2 = ϕ s -ϕ b + 2π Φ B Φ 0 , (4)
where Φ 0 = h/(2e) is the flux quantum. The biSQUID is controlled by fixing the total current I = I a + I b measured from the common output, which depends on all of these four phase differences, or alternatively on the superconducting phases ϕ a and ϕ b as well as the magnetic fluxes through the loops Φ A and Φ B (from this point on, we set ϕ s = 0 as the reference for the superconducting phases).

A key quantity of interest for such a superconducting interferometer is the flux-dependent critical current, i.e. the maximum dissipationless current which can flow through the device, defined here as

I c (Φ A , Φ B ) = Max ϕa,ϕ b |I (δϕ a1 , δϕ b1 , δϕ a2 , δϕ b2 )| . (5)
Out of the four quantum-dot based junctions, only a1 and b1 are coupled by nonlocal effects, leading to signatures in the critical current involving specific combinations of the external fluxes, which can be revealed through a simple phenomenological approach.

In the limit of low transparency, we can perform a simple perturbative treatment in the Cooper pair tunneling. Processes involving a single pair contribute to the total current I measured from S 0 in the standard Josephson form of a sinusoidal current-phase relationship for each of the four junctions, leading to

I 1P = I J [sin δϕ a1 + sin δϕ b1 + sin δϕ a2 + sin δϕ b2 ] . (6)
Processes involving two pairs all lead to second-order harmonics of the Josephson current. However, one must distinguish between two different contributions. Local processes, amounting to two pairs crossing any given junction, involve twice the phase difference seen by each quantum dot:

I 2P,local = I JJ [sin (2δϕ a1 ) + sin (2δϕ b1 ) + sin (2δϕ a2 ) + sin (2δϕ b2 )] . (7)
In addition to these, the present device allows the possibility for nonlocal processes, and one thus needs to consider the so-called pair cotunneling and quartet supercurrents. 6 The latter corresponds to the splitting of two correlated pairs from S 0 into S a and S b , which leads to a Josephson-like current-phase relationship involving both the phase difference between S 0 and S a , and the one between S 0 and S b , namely

I 2P,quartet = I Q sin (δϕ a1 + δϕ b1 ) . (8)
Due to the exchange process intrinsic to quartet emission, we expect I Q to be negative. 7 This nonlocal process relies on crossed Andreev reflection through S 0 , which naturally coexists with normal transmission through the central superconducting electrode without electron-hole conversion, the so-called elastic cotunneling. It follows that the quartet process comes with a similar partner corresponding to the exchange of a pair from S a to S b via double elastic cotunneling through S 0 , therefore contributing to the supercurrent as

I 2P,pair cotunneling = I PC sin (δϕ a1 -δϕ b1 ) , (9)
where the various current scales introduced above satisfy

|I PC |, |I Q |, |I JJ | I J
, as a consequence of the low transparency of the junctions.

Combining Eqs. ( 1) though (9), one obtains the following expression for the total critical current, up to second order in the pair tunneling

I c (Φ A , Φ B ) = Max ϕa,ϕ b 2I J sin ϕ a + π ΦA cos π ΦA + 2I J sin ϕ b -π ΦB cos π ΦB + 2I JJ sin 2ϕ a + 2π ΦA cos 2π ΦA + 2I JJ sin 2ϕ b -2π ΦB cos 2π ΦB + I Q sin (ϕ a + ϕ b ) + I PC sin (ϕ a -ϕ b ) (10)
where Φj = Φ j /Φ 0 is the magnetic flux through loop j, in units of the flux quantum. Apart from small regions near integer values of ΦA,B (where the single pair Josephson current is near suppression), the critical current can readily be obtained without any further calculation in most of the (Φ A , Φ B )-plane, and takes the form

I c (Φ A , Φ B ) =2I J cos π ΦA + cos π ΦB + |I Q | sin π ΦA -π ΦB + |I PC | sin π ΦA + π ΦB . (11)
There, the first term corresponds to the critical current in the absence of a quartet supercurrent. The two loops are decoupled in this case and I c (Φ A , Φ B ) splits into two independent contributions corresponding to the simultaneous maximization of I a and I b . Interestingly, quartet

A B S a S b S 0 Ia I b I a2 a1 b1 b2 δϕa2 δϕa1 δϕ b1 δϕ b2
Figure 2. The biSQUID setup as modeled here. Three superconducting terminals and four quantum dots (originating from a single nanotube or nanowire) define a two-loop system with enclosed areas of different sizes. Each dot Josephson junction sees a phase difference labeled δϕα combining the phases of each superconducting electrodes as well as the enclosed magnetic fluxes. The arrows near each junction set the conventions for the flowing currents and phase differences.

and pair cotunneling processes are responsible for an extra contribution to the critical current with a very specific flux dependence, leading to a macroscopic manifestation of nonlocal effects. In particular, this simple calculation points out that a measurement of the quartet current is possible by detecting the flux periodicity of the biSQUID critical current.

III. MICROSCOPIC THEORY

This section is devoted to a microscopic calculation of the current through the biSQUID. Our goal is to justify the expression (11) obtained from our simple phenomenological approach, to generalize it to more transparent junctions, and to show that the energy levels in the dots can be chosen such as to promote multi-pair transport, and in particular reveal the presence of a quartet resonance.

A. Hamiltonian

The model Hamiltonian for the biSQUID setup is expressed as the sum of three contributions

H = j=a,0,b H S,j + α=a1,a2, b1,b2 H D,α + H T . (12)
Here H S,j is the Hamiltonian associated with the superconducting lead S j (j = a, 0, b), which is given by the following compact form

H S,j = k Ψ † jk (ξ k σ z + ∆ j σ x ) Ψ jk , (13)
where ξ k = k 2 2m -µ and ∆ j is the superconducting gap of lead j. Each quantum dot α (α = a1, a2, b1, b2) is modeled by a single non-interacting level, with energy α , described by the Hamiltonian

H D,α = α d † α σ z d α . (14)
In both these expressions, we used Pauli matrices acting in Nambu space, and introduced Nambu spinors for the lead and dot electrons, defined respectively as

Ψ jk = ψ j,k,↑ ψ † j,-k,↓ d α = d α↑ d † α↓ . ( 15
)
Introducing a tunneling amplitude t jα between lead j and dot α, and performing a gauge transformation to incorporate the superconducting phases ϕ j in the tunneling term, one has for the tunneling part of the Hamiltonian

H T = k Ψ † jk T jkα d α + H.c. ( 16
)
where, in absence of external magnetic field, we have T jkα = t jα e -ikrjα σ z e iσzϕj /2 , r jα being the position of lead j where tunneling to/from dot α occurs. The resulting Hamiltonian is fully quadratic in terms of both the leads and the dot electrons. It is therefore convenient to integrate out the leads degrees of freedom and derive an effective theory involving only the dot electrons.

This is achieved through standard techniques, starting from the partition function

Z = D ψ, ψ, d, d e -S[ ψ,ψ, d,d] , (17)
with action

S = β 0 dτ j,k Ψjk (τ ) (D τ + ξ k σ z + ∆ j σ x ) Ψ jk (τ ) + β 0 dτ α dα (τ ) (D τ + α σ z ) d α (τ ) + β 0 dτ j,α,k Ψjk (τ )T jkα d α (τ ) + dα (τ )T * jkα Ψ jk (τ ) , (18)
where D τ is defined in Nambu space as

D τ = → ∂ τ 0 0 - ← ∂ τ .
Carrying out the integration over the ( ψ, ψ) Grassmann fields, one is left with an effective action of the form

S eff = β 0 dτ α dα (τ ) (D τ + α σ z ) d α (τ ) + β 0 dτ dτ α,γ dα (τ )Σ αγ (τ -τ )d γ (τ ) (19)
where we introduced the tunneling self-energy

Σ αγ (τ ) = j,k T * jkα G j,k (τ )T jkγ (20)
which depends on the leads electrons Green's function, defined in Matsubara frequency space as

G j,k (iω n ) = (iω n 1 -ξ k σ z -∆ j σ x ) -1 . ( 21
)
It follows that the effective field theory, quadratic in ( d, d), can be described uniquely in terms of the Matsubara Green's function Ĝ(iω n ) for the dot electrons, which takes the form of a 8×8 matrix in Nambu-dot space given by Ĝ

-1 (iω n ) =    g a2,a (iω n , Φ A ) 0 0 0 0 g a1,a (iω n , 0) f (iω n ) 0 0 f (iω n ) g b1,b (iω n , 0) 0 0 0 0 g b2,b (iω n , -Φ B )    (22)
where

g α,j (iω n , Φ) =     iω n 1 + πν(0) √ ω 2 n +∆ 2 (t 2 0,α + t 2 j,α ) -α -∆e iΦ (t 2 0,α + t 2 j,α e iδϕα ) πν(0) √ ω 2 n +∆ 2 -∆e -iΦ (t 2 0,α + t 2 j,α e -iδϕα ) πν(0) √ ω 2 n +∆ 2 iω n 1 + πν(0) √ ω 2 n +∆ 2 (t 2 0,α + t 2 j,α ) + α     (23)
f (iω n ) = πν(0)t 0,a1 t 0,b1 e -R/ξ(ωn)   iωn √ ω 2 n +∆ 2 cos(k F R) -sin(k F R) - ∆ √ ω 2 n +∆ 2 cos(k F R) - ∆ √ ω 2 n +∆ 2 cos(k F R) iωn √ ω 2 n +∆ 2 cos(k F R) + sin(k F R)   (24)
where we introduced the width R of the central super-conducting lead, and assumed all superconducting elec-trodes to have the same gap energy ∆ j = ∆, and density of states at the Fermi level ν(0). The energy-dependent coherence length is defined as

ξ(ω n ) = ξ 0 ∆/ ω 2 n + ∆ 2 .

B. Current

The current through a given quantum dot Josephson junction α, can be readily expressed as the time derivative of the number of electrons out of the superconducting reservoir, and as such can be related to the phase derivative of the tunneling Hamiltonian. In particular, the average current through dot α reads

I α = 2e ∂H T ∂δϕ α = 2e ∂F ∂δϕ α . ( 25
)
where F = -k B T log Z is the free energy. Keeping in mind that the effective theory contains all the relevant phase-dependent degrees of freedom, this can be further simplified as

I α = - 2e β ∂ log Z eff ∂δϕ α = - 2e β ∂ ∂δϕ α n log det Ĝ-1 (iω n ) (26)
where we performed explicitly the integration over the Grassmann fields ( d, d). The current through any given junction can thus be obtained from Eq. ( 22), after performing a Matsubara frequency summation. The total current through the device is readily obtained by summing up the contribution coming from each junction. The width R of the central superconducting lead, which corresponds to the separation between tunneling points from dot a1 to S 0 and dot b1 to S 0 , enters the above expression for the current in two important ways. First, the total current contains terms which decrease exponentially with R. These correspond to nonlocal contributions which expectedly vanish if the separation between tunneling points exceeds the superconducting coherence length, therefore setting a typical order of magnitude for R. Second, some of these terms are also oscillating on a scale set by the Fermi wavelength λ F . Since λ F is typically several orders of magnitude smaller than the superconducting coherence length, these contributions are rapidly oscillating as a function of R. It is thus natural to average over these rapid oscillations, assuming that the separation between tunneling points is susceptible to fluctuate slightly on a scale given by a few Fermi wavelengths. Defining R = R 0 + r, where R 0 is a fraction of the coherence length ξ 0 and r is of the order of the Fermi wavelength, we introduce the r-averaged total current as Considering the maximization with respect to the superconducting phases, the critical current through the device is now expressed in terms of the microscopic parameters as

ĪTot = 1 2λ F λ F -λ F dr -2e β n α ∂ log det Ĝ-1 (iω n ) ∂δϕ α . (27)
I c (Φ A , Φ B ) = Max ϕa,ϕ b
ĪTot (δϕ a1 , δϕ b1 , δϕ a2 , δϕ b2 ) , (28) where one needs to replace all phase differences δϕ α by their expression, Eqs. ( 1)-( 4), prior to evaluating the maximum.

C. Low transparency expansion

It is possible at this stage to provide a firmer basis for our earlier phenomenological treatment. Indeed, starting from Eq. ( 26), one can perform a perturbative expansion of the total current I coming out of S 0 in powers of the tunneling amplitude.

Focusing on the simpler case of equivalent junctions, i.e. assuming constant tunneling amplitude t 0 , and considering all dot energies to be the same (up to a possible sign, α = ± ), one can show analytically that the total current takes the form ĪTot =I J [sin δϕ a1 + sin δϕ b1 + sin δϕ a2 + sin δϕ b2 ]

+ I JJ [sin (2δϕ a1 ) + sin (2δϕ b1 ) + sin (2δϕ a2 ) + sin (2δϕ b2 )] + I Q sin (δϕ a1 + δϕ b1 ) , (29)
when expanding up to forth order in the tunneling rate Γ = 2πν(0)t 2 0 , which amounts to taking into account processes involving up to two Cooper pairs. The various contributions are expressed in terms of the microscopic parameters introduced in Sec. III A as

I J = e β n ∆ 2 (∆ 2 + ω 2 n )( 2 + ω 2 n ) Γ 2 e Γ 2 2
(30)

I JJ = e β n -∆ 4 Γ 4 4(∆ 2 + ω 2 n ) 2 ( 2 + ω 2 n ) 2 - e Γ 4
16 3 (31)

I Q = e β n -∆ 4 Γ 4 2(∆ 2 + ω 2 n ) 2 ( 2 + ω 2 n ) 2 cos 2 (k F R)e -2R/ξ - e Γ 4 16 3 e -2R0/ξ0 (32)
where we only kept the leading order contribution to I J (discarding the third and forth order corrections) and provided simplified forms valid in the large-gap, lowtemperature limit ∆ Γ, β -1 . At this level of approximation, the quartet component is of the same order of magnitude as the Josephson second harmonics (up to the nonlocal prefactor in R 0 ) and does not depend on the specific arrangement of the dots energy levels (I Q is unchanged whether a1 = b1 or a1 = -b1 ). These properties are specific to the present perturbative treatment and generally not expected to remain valid once one includes contributions from all orders in tunneling. The result of Eq. (32) also reveals that the amplitude of the quartet current is negative. This π-type behavior has been observed in Ref. 7 for multipair dc resonances in an out-of-equilibrium voltage-biased bijunction. It can be attributed to the internal structure of a Cooper pair via the antisymmetry of its wavefunction.

Interestingly, there is no trace of the pair cotunneling current in this derivation. This is an artifact of the low-order expansion, combined with the symmetry of the limit considered here. Indeed, a pair cotunneling current of the order of I Q × ∆ 2 does appear, arising from both S a and S b , but these two terms end up contributing to the total current with an opposite sign. This calculation not only validates the phenomenological form proposed in Sec. II, but also justifies that

I J |I JJ |, |I Q |, |I PC |,
as assumed in our discussion of the critical current.

IV. CRITICAL CURRENT

We now make use of the general expression (28) of the critical current of the device in order to explore its evolution for a broad range of parameters, going beyond the low transparency regime.

As argued from our simple phenomenological treatment presented in Sec. II, we expect quartets (and more generally all multipair processes) to lead to specific signatures in the critical current of the biSQUID, identifiable through their periodicity in the magnetic flux. In what follows, we thus focus on the Fourier spectrum of the critical current.

A. Signatures in Fourier space

The critical current can be written in a very general way in terms of its harmonics in Φ A and Φ B , namely

I c (Φ A , Φ B ) = +∞ n=-∞ +∞ m=-∞ I n,m e i2nπ Φ A Φ 0 e i2mπ Φ B Φ 0 . (33)
While we could study independently the variations of the critical current with respect to Φ A and Φ B , it makes more sense at this stage to introduce a new set of variables. Indeed, in practice, the magnetic flux through the device is provided by an external homogeneous magnetic field, which affects both Φ A and Φ B in a correlated way. The magnetic fluxes through loops A and B are therefore proportional to one another and only differ as a result of the different surface area enclosed by each loop, so that one can write

Φ A,B = Φ (1 ∓ η) , (34)
where we introduced the average magnetic flux Φ and an asymmetry parameter η defined as

Φ = Φ A + Φ B 2 (35) η = S B -S A S B + S A . (36)
In terms of these new variables, the expansion of the critical current I c takes the form

I c (Φ) = +∞ n=-∞ +∞ m=-∞ I n,m e 2iπ Φ Φ 0 [n(1-η)+m(1+η)] . (37)
The most convenient way to probe for these harmonics is to introduce the Fourier transform with respect to the average phase accumulated around the loops ϕ = 2πΦ/Φ 0 , which is defined as

Ĩc (N ) = 1 Φ 0 dΦ I c (Φ) e -2iπN Φ Φ 0 = +∞ n=-∞ +∞ m=-∞ I n,m δ (n(1 -η) + m(1 + η) -N ) . ( 38
)
It follows that the Fourier-transformed critical current Ĩc (N ) is given by a set of peaks, located at N = m + n + (m -n)η, whose height is related to the amplitude of the (n, m) harmonics. Identifying which of these features are present in the spectrum provides a direct and clear indication of what processes are at play in the biSQUID setup. Since the critical current is real and even in Φ, so is its Fourier transform. The Fourier spectrum is thus even in N , and we focus only on positive N . As a first step, we consider the simplifying approximation of identical junctions. This gives us an opportunity to discuss the typical signatures observed and dress up a physical picture, without obscuring it by dealing with too many parameters.

( - 2 , 2 ) ( - 3 , 3 ) ( 1 , 0 ) ( 0 , 1 ) ( 2 , 0 ) ( 3 , - 1 ) ( 1 , 1 ) ( 0 , 2 ) ( - 1 ,
3 ) ( 3 , 0 ) ( 2 , 1 ) ( 1 , 2 ) ( 0 , 3 ) ( 3 , 1 ) ( 1 , 3 )
Figure 4. (Left) Relative weight of the absolute-valued Fourier components (in N space) of the critical current as a function of Γ, for a symmetric arrangement of the energy levels α = = 0.1∆ (and η = 0.1, β = 100/∆, e -R 0 /ξ 0 = 0.9). (Right) Same relative weight in the specific case of high transparency, Γ = 0.8∆. There we also restored the sign of the various Fourier amplitudes and identified those in terms of a pair of integers (n, m) corresponding to the periodicity nΦA + mΦB of the critical current.

B. Identical junctions

We treat here the case of identical junctions, characterized by the same dot-lead tunneling rate Γ. For simplicity, we assume that the energy levels α of the dots can only take two values, either + or -, and distinguish various scenarios, based on the possible arrangements of these levels.

We focus on the low-temperature regime (β = 100/∆) and set the attenuation factor e -R0/ξ0 = 0.9, compatible with current experimental realizations.

Non-resonant case

We consider first the non-resonant case, for which the energy levels of the dots are non-zero, i.e. detuned with respect to the chemical potential in the superconducting leads. The results are summarized in Fig. 4.

In particular, we focused on two situations: (a) a symmetric arrangement of the energy levels of the dots involved in the two loops (all chosen equal, α = ) and (b) an antisymmetric arrangement of the energy levels between the two loops (where a1 = a2 = while b1 = b2 = -). Surprisingly, the obtained Fourier profiles are strictly identical for the symmetric and antisymmetric cases, signaling that the r-averaged critical current is insensitive to the arrangement of the dots energy levels, even beyond the low-transparency regime for which we could establish this property semi-analytically. The energy level symmetry has no influence on the Fourier spectrum, which one might interpret as a consequence of some kind of induced particle-hole symmetry, due to the proximity effect from the neighboring superconducting electrodes the dots are in contact with. In particular this means that the quantum dots no longer play the role of an energy filter, as they typically do in hybrid Normal metal -Superconductor -Normal metal structures, such as the Cooper-pair beam splitter, 10 where it allows to favor crossed Andreev reflection over pair cotunneling processes.

In the low transparency regime, the nonlocal effects are not strong enough to be readily identified and one expects the critical current to be insensitive to the arrangement of the energy levels. Indeed, in both cases, the Fourier profile is largely dominated by the usual Josephson current associated to two independent loops. The Fourier spectrum is composed of two sets of peaks, located at multiples of 1-η and 1+η respectively. These correspond to harmonics in Φ A and Φ B (of order (n, 0) and (0, m) respectively) and can be readily associated with the usual Josephson current flowing through loops A and B in absence of any nonlocal coupling between the two. Indeed, in the case of independent loops, one expects the critical current to be of the form 2I J j cos π Φj . This, in turn, leads to a characteristic evolution of the Fourier profile, which we could compare to the one obtained in the present low-transparency regime. In Fig. 3, we show that the weights associated with the different harmonics follow the predicted behavior, confirming that the critical current is dominated by the usual local harmonic Josephson component. Quartets

( - 1 , 1 ) ( - 2 , 2 ) ( - 3 , 3 ) ( - 4 , 4 ) ( 1 , 0 ) ( 0 , 1 ) ( 2 , 0 ) ( 3 , - 1 ) ( 1 , 1 ) ( 0 , 2 ) ( - 1 ,
3 ) ( 2 , 1 ) ( 1 , 2 ) ( 2 ,
2 )

Figure 5. (Left) Relative weight of the absolute-valued Fourier components (in N space) of the critical current as a function of Γ, for a resonant arrangement of the energy levels α = = 0 (and η = 0.1, β = 100/∆, e -R 0 /ξ 0 = 0.9). (Right) Same relative weight in the specific case of high transparency, Γ = 0.8∆. There we also restored the sign of the various Fourier amplitudes and identified those in terms of a pair of integers (n, m) corresponding to the periodicity nΦA + mΦB of the critical current.

As the transparency is increased, new structures appear in addition to the ones already observed in multiples of 1 ± η (see Fig. 4). Observing components which involve the fluxes associated with each loop (i.e. with both n and m nonzero) signals the presence of nonlocal processes contributing to the critical current. They are relatively faint for low values of Γ, but increase as one goes into the high transparency regime. For the highest value of the tunneling amplitude considered here (Γ = 0.8∆) these nonlocal contributions represent an appreciable total weight of about 30% of the signal.

Out of these new features, the most pronounced one -a peak in N = 2η -corresponds to a component of the r-averaged critical current with a periodicity (n = -1, m = 1) in the fluxes, i.e. periodic in Φ A -Φ B . Since it depends on both fluxes, this contribution can only arise from a process involving both the exchange of pairs between S 0 and S a as well as between S 0 and S b . Moreover, the periodicity in Φ A -Φ B requires such an exchange to be correlated, and to occur in the same direction, that is from S 0 to S a,b or the other way around. This is precisely the microscopic definition of the so-called quartet process, an exchange of two correlated pairs from the central superconducting electrode S 0 to S a and S b . The quartet contribution increases rapidly with Γ, and becomes substantial in the limit Γ , representing close to 10% of relative weight, much more than any other components (apart from the leading harmonics in Φ A and Φ B ). The physical nature of this structure is further confirmed by analyzing the sign of its amplitude (Fig. 4, right), which turns out to be negative, as expected for the quartet contribution. 30 In the non-resonant case, the biSQUID thus seems a promising candidate to observe signatures of the quartet current, by monitoring the periodicity of the critical current in the total flux, i.e. as a function of the external magnetic field.

Note that other, much weaker signatures are also present with a periodicity in both fluxes. These could be attributed to either anharmonicities of the quartet current (which would be expected for such high transparency) or to higher order multipair contributions.

It seems that the nonlocal effects are more pronounced when Γ > , which naturally takes us to the regime where this condition is most strongly fulfilled: the resonant case.

Resonant case

We now turn to the resonant case, where all energy levels are chosen equal to the chemical potential in the superconducting leads (set to 0 by convention). The results are summarized in Fig. 5.

The low transparency regime is again largely dominated by the local harmonic Josephson current , and absent of any nonlocal signatures. As a result, while the global prefactor I J might differ, the Fourier spectrum of the critical current is identical to the one obtained in the non-resonant case.

Contributions associated with nonlocal processes appear rapidly, and already represent about 20% of the signal for tunneling amplitudes as low as Γ = 0.1∆. As the transparency is further increased, the nonlocal components get stronger, totaling about half of the signal once one reaches Γ ∼ ∆.

As in the non-resonant case, there is a proliferation of signatures corresponding to various periodicities in the fluxes Φ A and Φ B (identified through the pair of integers (n, m)). The strongest of these nonlocal components is again observed at N = 2η and corresponds to the quartet current. This contribution grows much more rapidly as a function of Γ, compared to the non-resonant case, and ends up being twice as large for a comparable value of the tunneling rate. This is consistent with the results obtained in the previous section, as the regime Γ > proved to be the most suitable to observe signatures of multipair processes, it was to be expected that the resonant case would strengthen this trend further.

Along with the main harmonics (1, 0) and (0, 1), and the quartet component (-1, 1), the Fourier-transformed critical current also shows a substantial contribution for N = 2, corresponding to the periodicity Φ A + Φ B (n = 1, m = 1). Like the quartet component, this contribution also increases with the transparency, and while it was only marginally present in the non-resonant case, it becomes significant at resonance. This feature can be attributed to the pair cotunneling process, which amounts to sending a Cooper pair from S a to S b through the central superconducting electrode S 0 . As argued in Ref. 6, the pair cotunneling process goes hand in hand with the quartet process, and it is not surprising that they appear side by side and behave in a somewhat similar way. As opposed to the low transparency expansion presented in Sec. III C where it vanished for symmetry reasons, the pair cotunneling contribution is not negligible here, despite the apparent symmetry of the setup. This reveals the importance of high-order tunneling processes, not accounted for in the low-Γ expansion, which involve the fluxes enclosed in the two loops, thus breaking the symmetry between the two branches of the setup.

The presence of particularly strong quartet signatures makes this resonant case the best candidate to observe such nonlocal contributions to the critical current. We thus continue our study of the resonant case but move on to a more realistic scenario of different junctions in order to determine to what extent our results are robust.

C. Robustness of the results

The results of the previous section were obtained in the simpler case of identical junctions, and one might thus argue that they require some specific fine tuning in order to be observed. While the presence of external gates enables a rather precise control of the energy levels of the four quantum dots, the tunneling rates of the junctions depend on the sample design and cannot be manipulated. In particular, making four identical highly transparent Josephson junctions constitutes an experimental challenge.

We therefore focus now on a more realistic situation, where we set the dots to be resonant with the super- Relative weight of the main harmonic (0, 1) in Φ B (%) Probability of occurrence Figure 6. Probability of occurrence of the main harmonics and in ΦA (Top) and ΦB (Bottom) as a function of the relative weight they represent in the Fourier spectrum. These were obtained by computing the critical current for 1000 different realizations of the setup with resonant dots ( α = 0) and randomly chosen tunneling amplitudes tjα between 0 and ∆. conducting electrodes, but do not make any assumptions concerning the tunneling parameters. To take this even further, we contemplate the possibility of having different tunneling amplitudes t jα for each dot-lead junction. We compute the critical current as a function of the average flux for 1000 different realizations of the setup, each corresponding to a random pick of the eight tunneling amplitudes, chosen following a uniform distribution over the range [0; ∆]. From the resulting Fourier spectra, we monitor the fate of the four most prominent features identified earlier, namely the main harmonics in Φ A and Φ B (corresponding to the features in (1, 0) and (0, 1)), as well as the quartet and pair cotunneling contributions (corresponding to (-1, 1) and (1, 1) respectively). In Figs. 6 and7, we show histograms representing the probability of occurrence of each of these contributions as a function of the relative weight they represent in the Fourier spectrum.

From Fig. 6, one sees that the main harmonics in Φ A and Φ B follow a very similar distribution. These are clearly the leading contributions to the critical current averaging a relative weight of about 30% each, while together they add up to over 50% of the signal in about half the realizations. The results compiled in Fig. 7 suggest that the quartet and the pair cotunneling contributions bear some striking similarities when it comes to their weight distribution, further underlining the deep connection between the two processes. Focusing on the quartet component, one readily sees that for a completely random selection of the tunneling amplitudes, the probability of observing a strong quartet signal is rather weak, as less than 60 out of the 1000 realizations show a relative weight beyond 10%. However, the quartet signature is more often present than not, its relative weight reaching over 3% for more than half the realizations. On average, the contribution to the critical current associated with quartets culminates a little above 4%. To put things into perspective, in the resonant case considered in the previous section, a junction transparency of Γ = 0.08∆ was sufficient to reach the same kind of values, stressing the importance of making the four junctions as identical as possible in practice.

While on average the quartet signature is small, it is still detectable and could be further enhanced by filtering out the main Josephson harmonics in Φ A and Φ B which are easily identifiable. One also has to keep in mind that the situation considered here is among the most unfa-vorable ones, and while practical realizations should be tailored to promote nonlocal multipair processes, even a random realization still has a solid chance of revealing their specific signatures.

V. COMPETING EFFECTS

Various physical phenomena can be responsible for a coupling between degrees of freedom from the two loops, leading to nonlocal effects whose signatures might mask a pure quartet signal. We hereby consider two of the most likely candidates susceptible to appear in our biSQUID device.

A. Mutual inductance

The setup under consideration is comprised of two loops where current can flow, so that one should take into account the geometrical inductance of such a circuit. While there could be both local and mutual inductances at play in the setup, the latter is the most likely to lead to strong nonlocal signatures as it couples the currents of pairs flowing through junctions a1 and b1. We could check that the presence of a self inductance only marginally modifies our results and thus decided, for simplicity, to focus on the effect of the mutual inductance alone.

Due to the mutual geometrical inductance M between loops A and B, the fluxes felt by the electrons circulating in the circuit is no longer set by the external magnetic field alone. Instead, there exists an additional contribution on top of this external flux (hereby labeled Φ ext A,B ) which depends on the current flowing in the nearby loop. Following the conventions introduced in Fig. 2 and Eqs. ( 1)-( 4) for the currents and fluxes, we write

Φ A = Φ ext A -M (I b1 -I b2 ) (39)
Φ B = Φ ext B -M (I a2 -I a1 ) , (40)
where the mutual inductance M is positive. Since some of the junction currents I α depend explicitly on the fluxes, it is obvious from Eqs. (39)-(40) that the total current through the device now has to be computed self-consistently. Performing this self-consistent treatment in addition to the r-averaging and on top of the maximization (in order to extract the critical current) makes the numerical calculation a lot more demanding. In order to estimate the competing effects associated with the mutual inductance, we thus consider the simpler case of two disconnected loops (i.e. without any microscopic nonlocal coupling) and turn on the mutual inductance progressively, monitoring the Fourier-transformed critical current. The results are provided in Fig. 8, where for illustrative purposes, we focused on the resonant highly transparent regime.

Whereas there are no noticeable modifications of the Fourier spectrum for weak values of the mutual inductance, new signatures start appearing beyond 0.1M 0 (where M 0 = 2 /(e 2 ∆)), a threshold which tends to increase as the tunneling rate Γ is reduced (not shown).

These signatures grow rapidly with M and are mainly located at N = 2 and N = 2η (i.e. the (1, 1) and (-1, 1) components respectively). Unlike the general situation considered in the previous section, the most pronounced component involving both fluxes corresponds here to (1, 1) rather than (-1, 1), a characteristic specific to the presence of mutual inductance, which could be used in actual experiments to detect the presence of such effects. Indeed, this can be explained by noticing that the mutual inductance tends to anticorrelate the currents in branches a1 and b1 in order to decrease the total current in the central electrode. On the contrary, the quartet mechanism tends to correlate the currents I a1 and I b1 reinforcing the current in the central electrode. The antagonistic effects of quartets and mutual inductance allows to discriminate both mechanisms. For the largest value of the mutual inductance considered here, this structure at N = 2 even becomes the leading component of the Fourier signal, overcoming the contributions which depend on only one of the two fluxes. While a peak at N = 2η is visible (and as such, could be mistaken with signatures from the quartet process), it stays relatively weak (5% at best) and even decreases back for large values of M .

In the end, although mutual inductance leads to structures in the critical current located at the same position in N space, the behavior of these new features is specific enough to not be confused with quartet and pair cotunneling processes. Still, the mutual inductance needs to be carefully estimated in experimental setups 31 and taken into account when attempting to fit the data.

B. Direct tunneling

Although absent from nanowire-based systems, direct tunneling between quantum dots is known to occur in carbon nanotubes. This generates a new contribution to the bare dot Hamiltonian, of the form

H D,direct = t d d † a1 σ z d b1 + H.c.. ( 41
)
where t d is the direct tunneling amplitude between dots a1 and b1. This term opens a new channel for current to flow between the two loops, competing with the nonlocal exchange of pairs through S 0 .

To get a flavor of the effect of direct tunneling, let us compute the critical current as a function of t d in the simpler case where the nonlocal coupling through the central superconducting electrode has been turned off. The new term that now appears in the total Hamiltonian changes the form of the Matsubara Green's function for the dot electrons by affecting its anomalous part f (iω n ),

M = 0.1M 0 M = 0.4M 0 M = 0.6M 0 M = 0.9M 0 M = 1.3M 0 Figure 8.
Relative weight of the Fourier components of the critical current as a function of the mutual inductance M between loops (expressed in units of M0 = 2 /(e 2 ∆)), in the resonant ( = 0) highly transparent regime (Γ = 0.8∆). All other parameters have been chosen identical to Figs. 4 and5, apart from the attenuation factor, e -R 0 /ξ 0 = 0. which becomes

f (iω n ) = -t d 0 0 t d . (42)
The critical current is then obtained from Eqs. ( 27) and (28) as before. The Fourier-transformed critical current as a function of the direct tunneling t d is presented in Fig. 9 for the resonant highly transparent case. As it turns out, direct tunneling can lead to signatures that are very similar to nonlocal multipair processes. In the resonant case considered here, even for rather small values of t d (say, 0.1∆) new peaks appear in the Fourier spectrum corresponding precisely to the contributions of interest, namely (-1, 1) and (1, 1) (associated earlier to quartet and pair cotunneling processes). Moreover, these features grow rapidly as t d is increased, and for a typical value of only t d = 0.4∆, the resulting Fourier-transformed critical current looks almost identical to the one obtained in Fig. 5 in presence of nonlocal crossed Andreev reflection (same set of harmonics with similar weights, only the ones in the vicinity of N = 4 differ). Indeed, from inspecting the anomalous Green's functions in presence of a coupling through S 0 , Eq. ( 24), and in presence of direct tunneling, Eq. (42), one readily sees that t d roughly assumes the same role as the prefactor Γe -R0/ξ0 /2 (up to complications related to r-averaging). Indeed, direct interdot tunneling offers an alternative channel for the production of correlated pairs, as shown in Refs. 11 and 12. Yet, it is interesting to discriminate this mechanism from the quartet one, for instance by synchronous detection while slowly varying one key parameter of the setup.

t d = 0.2∆ t d = 0.4∆ t d = 0.6∆ t d = 0.8∆ Figure 9
. Relative weight of the Fourier components of the critical current as a function of the direct tunneling amplitude t d , in the resonant ( = 0) highly transparent regime (Γ = 0.8∆). All other parameters have been chosen identical to Figs. 4 and5, apart from the attenuation factor, e -R 0 /ξ 0 = 0.

VI. SUMMARY AND CONCLUSION

To summarize, our main results are as follows: (a) the phenomenological model is sufficient to capture the main ingredients of the physics involved in this system in the low transparency regime;

(b) studying the flux-dependence of the critical current through a full microscopic calculation leads, for identical junctions in the non-resonant regime, to new signatures. These could be attributed to nonlocal multipair processes;

(c) in the case of resonant dots, and for highly transparent junctions, these quartet signatures become particularly strong (of the same order as the Josephson contribution) making this regime the most promising to investigate multipair production;

(d) our results are robust against strong variations of tunneling amplitude throughout the setup. Even in the worst case scenario of 8 very different tunneling parameters, there is a strong probability of observing a quartet signal;

(e) the most likely competing effects either lead to qualitatively different signatures or open a new channel that contributes to the formation of correlated pairs, in addition to the standrad crossed Andreev reflection.

This work thus provides an alternative experimental setup for the observation of nonlocal multipair processes. Unlike previous work on a voltage-biased allsuperconducting bijunction 7 , the current study focuses on an equilibrium situation, in the coherent dissipationless regime. Making use of interferometry, this allows for the observation of specific features associated with quartet processes that do not require to study intricate phase and voltage dependences, and hopefully provides an unequivocal experimental signature of these phenomena.

In conclusion, we proposed a setup consisting of four Josephson junctions, defining two twinned loops in a nanotube (or nanowire) based system connected to three superconducting electrodes, which we dubbed a biSQUID. We presented a phenomenological argument to motivate our analysis then derived a microscopic theory allowing a careful description of the setup in the full parameter space. We showed that by measuring the critical current as a function of the average flux through the loops (i.e. the external magnetic field), the device could reveal signatures associated with nonlocal multipair contributions, in particular the so-called quartet process, thus making them experimentally observable in an equilibrium situation. We also suggested that potential experimental realizations pay specific attention to some competing effects which might interfere with these very signatures.

A natural extension of the present work consists in taking into account local Coulomb interaction on the dots. This constitutes a real challenge, which could be tackled through various approximate treatments, such as perturbative diagrammatic resummation or self-consistent mean-field approaches. 28 Interactions are expected to be particularly relevant when dealing with closed dots, with weak tunneling amplitudes to the leads. We showed, however, that the most interesting regime to observe quartet signatures in the biSQUID setup corresponds to the opposite case of large values of the tunneling rate Γ, for which Coulomb interactions play a more limited role. As long as one stays away from the deep Kondo regime (for which the occupation of the dots is close to 1 at all times), we therefore do not expect any dramatic qualitative change of our results due to Coulomb interactions. 29 For simplicity, the other two electrodes Sa and S b are chosen wide enough so as to avoid nonlocal effects. 30 This negative sign is however unrelated to the expected π-type behavior of the quartet current, as the sign of the Fourier amplitudes of the critical current is not directly related to the sign of the different contributions to the total current. 31 For a gap energy ∆ ∼ 100 µeV, the mutual inductance considered here is between 50 and 500 pH, which should be a reasonable order of magnitude of what is realized in actual setups.

Figure 1 .

 1 Figure 1. Artistic view of the biSQUID setup where a single nanotube or nanowire is contacted to three superconducting electrodes (all at the same potential), with magnetic field biasing.

Figure 3 .

 3 Figure 3. The relative weight of the numerically obtained Fourier components (in N space) of the critical current in the low transparency regime (Γ = 0.01∆), for identical loops (η = 0) at low temperature, β = 100/∆. This is compared to the analytic solution corresponding to the usual harmonic Josephson current (dotted line).

  Relative weight of the pair cotunneling contribution (%) Probability of occurrence

Figure 7 .

 7 Figure 7. Probability of occurrence of the quartet (Top) and pair cotunneling contributions (Bottom) as a function of the relative weight they represent in the Fourier spectrum. These were obtained by computing the critical current for 1000 different realizations of the setup with resonant dots and randomly chosen tunneling amplitudes tjα.

ACKNOWLEDGMENTS

We acknowledge the support of the French National Research Agency, through the project ANR-NanoQuartets (ANR-12-BS1000701). Part of this work has been carried out in the framework of the Labex Archimède ANR-11-LABX-0033 and of the A*MIDEX project ANR-11-IDEX-0001-02 (J.R., T.J. and T.M.). The authors acknowledge seminal discussions with T. Kontos concerning the biSQUID scheme.

Appendix A: Tunneling self-energy

The tunneling self-energy introduced in the text is given, in terms of Matsubara frequency, by

where we substituted in Eq. ( 20) the expression for the lead electrons Green's function, and the tunneling parameters.

Performing the k-integral requires information on the dimensionality of the system. It is believed 12 that the contact between dots and superconducting leads actually occurs within the carbon nanotube, between an electrostatically confined region (the dot) and a proximityinduced superconducting region, which then acts as an effective 1D superconductor. Within this assumption, the integral over momentum reduces to

where we introduced R j,αγ = |r jα -r jγ | as well as

Here ν( ) is the density of states of the superconducting region, which we assumed to be constant close to the Fermi level. Substituting this result back into Eq. (A1), one has, for the tunneling self-energy in Matsubara frequency space

In practice, the only nonlocal terms allowed by the tunneling amplitudes involve coupling between dots a1 and b1 through lead S 0 . Therefore, we only need to introduce one length scale R j,αγ corresponding to R 0;a1,b1 = R 0;b1,a1 = R.