Anna Maria Vegni
email: amvegni@uniroma3.it

Valeria Loscrí
email: valeria.loscri@inria.fr

Alessandro Neri
email: neri@uniroma3.it

Marco Leo
email: marco.leo@uniroma3.it

A Bayesian Packet Sharing Approach for Noisy IoT Scenarios

Keywords: IoT-based, Belief Propagation, Cloud. Subnetwork N Subnetwork 2 Subnetwork Cloud-based Architecture IoT Device

Cloud computing and Internet of Things (IoT) represent two different technologies that are massively being adopted in our daily life, playing a fundamental role in the future Internet. One important challenge that need to be handled is the enormous amount of data generated by sensing devices, that make the control of sending useless data very important. In order to face with this challenge, there is a increasing interest about predictive approaches to avoid to send high spatio-temporal correlated data. Belief Propagation (BP) algorithm is a method of performing approximate inference on arbitrary graphical models that is becoming increasingly popular in the context of IoT. By exploiting BP, we can derive effective methods to drastically reduce the number of transmitted messages, while keeping high the data throughput in the global information system.

In this paper, we propose a BP approach in a hierarchical architecture with simple nodes, gateways and data centers. We evaluate the error bounding and propose a corrective mechanism to keep a certain quality of the global information in the architecture considered.

I. INTRODUCTION

Recently, with the evolution of the Internet and related technologies, there has been an evolution of a new emerging paradigm, namely the Internet of Things (IoT) [START_REF] Atzori | The internet of things: A survey[END_REF]. In IoT scenarios, a large number of devices -and more in general objects-are seamlessly connected to each another for information sharing through the Internet. All these devices connected to the IoT may be of heterogeneous types with respect to their operational mode, and communication technologies. As one of the main strengths behind the IoT paradigm, it is the high impact on several aspects of everyday-life, from working to the domestic fields. As an instance, domotics [START_REF] Miori | Anticipating health hazards through an ontology-based, iot domotic environment[END_REF], e-health [START_REF] Hu | On the application of the internet of things in the field of medical and health care[END_REF], and smart cities [START_REF] Bonino | Almanac: Internet of things for smart cities[END_REF], [START_REF] Jiong Jin | An information framework for creating a smart city through internet of things[END_REF] are main application scenarios where the IoT paradigm is expected to play a leading role in the next future.

From the above considerations, and due to the huge amount of heterogeneous devices, information sharing among IoT devices is one of the biggest challenges. Classic Internet approaches need to be revised to address the complex requirements imposed by IoT. This asks for the development of intelligent algorithms for routing [START_REF] Alvi | Energy efficient green routing protocol for internet of multimedia things[END_REF], information sharing security [START_REF] Giuliano | Security access protocols in iot networks with heterogenous non-ip terminals[END_REF], [START_REF] Cervantes | Detection of sinkhole attacks for supporting secure routing on 6lowpan for internet of things[END_REF], novel network paradigms [START_REF] Jiong Jin | Network architecture and qos issues in the internet of things for a smart city[END_REF], [START_REF] Zhijing Qin | A software defined networking architecture for the internet-of-things[END_REF], [START_REF] Xingmei | Research on the basic characteristics, the key technologies, the network architecture and security problems of the internet of things[END_REF], new services [START_REF] Li | Study on the architecture of platform for internet of things service based on edsoa[END_REF], [START_REF] Tnjes | A test-driven approach for life cycle management of internet of things enabled services[END_REF], and advanced techniques for data fusion [START_REF] Bijarbooneh | Cloudassisted data fusion and sensor selection for internet-of-things[END_REF].

A few related works have addressed the issue of forwarding data among IoT devices, by modelling the IoT network as a Bayesian network [START_REF] Bijarbooneh | Cloudassisted data fusion and sensor selection for internet-of-things[END_REF], [START_REF] Kumar | Bayesian coalition game for the internet of things: an ambient intelligence-based evaluation[END_REF], [START_REF] Liu | Toward reliable data analysis for internet of things by bayesian dynamic modeling and computation[END_REF]. Under this hypothesis, Bijarbooneh et al. [START_REF] Bijarbooneh | Cloudassisted data fusion and sensor selection for internet-of-things[END_REF] present an adaptive sensing belief propagation algorithm, where each node updates its belief about the environment status by incorporating its local measurement with the beliefs of its neighboring nodes and the belief obtained in the past. Finally, a cloud-based network architecture allows reducing energy consumption for data fusion and storage tasks through a carefull selection of the set of IoT nodes involved in the distributed estimate of the environment status.

In this paper, we address the connectivity issue among IoT heterogeneous devices for data sharing, under the hypothesis of Bayesian network i.e. . We assume each IoT device represents a node in the IoT network with some sensing and processing capabilities. Moreover, these devices may be located at different places across the globe. They are connected to the Internet, althgough the rate of data transfer, and the supported security level may be different. Each node needs to get information about its local environment, in order to perform some task and/or to provide this information to a higher decision level. As an instance, an IoT node deployed in a domestic network may need information about current and future usage of some limited resources, like energy or communication bandwidth, to orchestrate their comsumption with the help of the opther IoT devices controlling specific appliances. Data exchange among devices allows a single node to increase its own knowledge of global information and optimize the scheduling of the tasks that it has to accomplish with the usage of a shared resource.

In this paper, we present a data sharing approach based on Pearl's Belief Propagation (BP) algorithm in the IoT context with a cloud-based architecture. BP is an iterative technique mainly used for solving inference problems. In the IoT context, the belief of a device (e.g., a sensor node) is the data measurement. The BP infers the measurements of other neighboring devices, especially in cases where the data is missing. Moreover, the BP technique allows to correct the errors that can occur in the data propagation. At each run, the BP provides to the devices both spatial and temporal cooperation. Indeed, in BP-based approaches, each sensor node determines its belief by incorporating its local measurement with the beliefs of its neighboring nodes, as well as its beliefs obtained in the past run. Then, at each run, we assume a node is able to (i) reduce the own distortion level (i.e., estimation error on global information), and (ii) provide an update of the This paper is organized as follows. Section II describes the reference cloud-based architecture for the IoT scenario. Due to the heterogeneity of the IoT devices, we assume a multinetwork scenario with a plethora di interconnected devices. In Section III we present our technique based on BP algorithm, for data sharing in a noisy IoT scenario. The presence of errors along the message reconstruction phase occurring at each receiver node can be mitigate through our BP's algorithm. Finally, conclusions are drawn at the end of this paper.

II. NETWORK MODEL

In our architecture, we consider different entities with specific computational and communication capabilities and functionalities. As illustrated in Figure 1 the IoT network model may comprise several sub-networks (i.e., from 1 to N ), associated with different applications. Indeed, due to the huge amount of IoT devices, we assume each sub-network is composed by IoT devices connected to each others for data sharing, and a gateway that interacts with the "external world". The role of the gateways consists into relaying the messages to the cloud, which is responsible of typical cloud-based services (e.g., data fusion, storage, etc.).

Each device performs sensing and processing activities. Our network model supports multi-hop routing, and the gateways collect data and forward them to the cloud. To support high scalability of the architecture, gateways implement publish-subscribe message passing mechanisms based on message brokers. Depending on the device capabilities, publishsubscribe mechanisms can be supported even at IoT node level. In this scenario, we assume that each IoT node can be in two states, either idle or active. A node is in a idle state, when it disconnects its radio, and it cannot send and receive data. On the other hand, a node is in active state when it can perform sensing activities, and can send and receive messages.

Moreover presence of errors in the connectivity links among IoT devices has to be accounted for. This can affect communication reliability, and also cause packet losses. Packet retransmission are then necessary to overcome packet errors In order to keep low the effective number of retransmissions, each active sensor implements an error conceiling strategy based on a BP message passing algorithm with a twofold objective i.e., (i) to recover missing data through the BP algorithm, and (ii) to reconstruct "incomplete" or "corrupted" messages.

The proposed architecture has been adopted by the authors in the design of the proof of concept of DAHMS and LOGON projects, both funded by the Italian Ministery of Economic Development in the framework "New Technologies for the Made in Italy", carried out at the Radiolabs research center labs. DAHMS (Distributed Architecture Home Modular Multifunctional Systems) scope is the improvement of the quality of life and the degree of self-sufficiency of chronically ill and elderly and disabled persons through the integration of Home Automation and remote heathcare functionalities. LogOn (Logistic Open Network) concerns the goods logistic in historical art cities with high level of tourism economy. The Secure Mediation GateWay (SMGW) represents the conjunction element among the devices within the Control Room and on-board (vehicular) devices.

The on-board devices are able to send events and receive commands from a set of control interfaces accessible both from the Control Room, the vehicles, and also from mobile devices (Android or iOS).

The information messages sent from each device follow a publish/subscribe framework able to forward MQTT (MQ Telemetry Transport) messages. MQTT is a messaging protocol working on the top of TCP/IP, that has been projected for specific situations where low impact and limited bandwidth are required. All the devices are equipped with MQTT and will be in charge to forward informations about events or commands. Through the federated system of SMGW, the informations published via MQTT will be available in a seamless way from each intraSMGW domain, as described in the architectural scheme of SMGW. Within each domain, the devices -Raspberry or Arduino equipped-are in charge of exchanging event/commands and messages with the SMGW. Also, the SMGW exports each message in a secure way towards all the other SMGWs.

III. A BP TECHNIQUE FOR ERROR CORRECTION

In this section, we investigate the proposed BP technique in IoT noisy scenarios. Our aim is to reduce the message errors through an iterative algorithm that corrects and updates the received data at each run. With regards to Figure 1, each device is initiated in active mode, and transmits messages to its own neighbors. A message is related to local data measurements, sampled at a fixed time step. The global information, related to a given sub-network, is then obtained from the contributions coming from each IoT device within the subnetwork.

Let us assume that the distributed sensing system consists of N IoT nodes, interconnected in various ways. Each IoT node collects a set of data provided by several sensors. Our scope is then to estimate the state X of the sensed environment starting from the sets {D i } of data collected by the individual nodes related to non overlapping regions {R i }. Here X is modeled as a dynamical Random Field. More in detail we focus our attention to the case in which the dynamical (i.e. temporal) behaviour of the system can be described by a linear model, while the spatial behavior is described by a Markov Random Field.

We incidentally recall that given a finite rectangular lattice

L = {(i, j), 1 ≤ i ≤ N 1 , 1 ≤ j ≤ N 2 }
, a neighborhood system associated with L is, by definition, a collection of subsets η = {η ij } with the property that each subset η ij , namely the neighborhood of i, j, is such that:

• (i, j) ∈ η ij , • if (k, l) ∈ η ij , then (i, j) ∈ η kl , ∀(i, j) ∈ L
It follows that a random field X is said to be a Markov Random Field w.r.t. (L, η) if and only if:

P (X(i, j)/X(k, l), (k, l) ∈ L -{(i, j)}) = = P (X(i, j)/X(k, l), (k, l) ∈ η ij -{(i, j)}) , ∀(i, j) ∈ L.
Then, in the distributed estimation scheme we can take advantage of the Hammersley-Clifford Theorem, stating that the joint distribution of a Markov Random Field w.r.t. (L, η) is of the form

P X (x) = 1 Z exp {-U (x)} , (1)
where Z is a normalizing constant and

U (x) = ∀clique c V c (x), (2)
is the energy function, and V c (x) is the potential associated with clique c ∈ C.

The only constraint on the clique potential V c (x) is that it depends only on the restriction of x to C. Nevertheless, here we focus our attention on those Markov Random Fields for which the potential function consists only of a set of singleton potentials, defined on single variables, and on a set of pairwise potentials, defined on pairs of variable.

To derive the distributed state estimation model, here we resort to the unified representation for both Bayesian Networks and Random Markov Fields, constituted by the Factor Graphs. Factor Graphs use factor nodes to describe the factorization property of the joint distribution, as the one stated by the Hammersley-Clifford Theorem.

At this aim, for sake of compactness of the notation, without loss of generality, we assume that the sensor data D i and D j provided by the sensors respectively connected to the i-th and j-th node cover two non-overlaping areas, and that the overall state space X is the Cartesian product of the state subspaces X i associated to the environmental variables related to the areas covered by the individual nodes.

By associating each node i of a sub-network, with a random variable X i that represents the local information, and by considering a set of edges E, we can write the joint distribution as:

P X (x) = i ψ i (x i ) (i,j)∈E ψ ij (x i , x j ), (3)
where the function ψ ij () represents the message exchange among node i and j. In practice, p(x i ) represents the marginal distribution of i-th node, and the BP allows the computation of the marginal distribution at each node i.

From rate-distortion theory, given a one-dimensional random variable X (X) is the representation of X, so that

X ∈ 1, 2, ..., 2 nR , ( 4
)
where R are the bits needed for the representation of X.

Then, the distortion function is a mapping

d : X × X → R + , (5)
from the set of source alphabet pairs X into the set of nonnegative real numbers. It measures the cost of representing symbol x by x. Each node has the aim to minimize the own distortion level. A distortion measure is said to be bounded if the maximum value of the distortion is finite i.e.,

d max := max x∈X , x∈X d (x, x) < ∞. (6)
For the distortion, we can assume a squared-error definition i.e.,

d (x, x) = (x -x) T (x -x), (7)
from which we can derive the distortion between sequences x n and xn as

d (x n , xn ) = 1 n n i=1 d (x i , xi ). (8)
It follows that the distortion associated with a (2 nR , n) code is defined as:

D = E [d (X n , g n (f n (X n )))] , (9)
where f n : X n → 1, 2, ..., 2 nR , and

g n : 1, 2, ..., 2 nR → X n .
Finally, we can derive the information rate distortion function R(D) for a source X with distortion measure d(x, x) as:

R (D) = min I X; X , (10)
where I(X; X) is the mutual information. Notice that Eq. ( 10) is subject to the following constraint

p ( x| x) : (x,x) p (x) p ( x| x) d (x, x) ≤ D, (11)
that is, the minimization of the mutual information is over all conditional distribution p ( x| x) for which the jointly distribution p(x, x) satisfies the expected distortion constraint.

From ( 9), the expectation value respect to the probability distribution on X is as follows:

D = x n p (x n ) d (x n , g n (f n (x n ))). (12)
Now, in order to solve previous equation, we need the estimation of x i . This can be provided through the BP algorithm. This is a message passing algorithm for the calculation of a posteriori probabilities of nodes of a loop-free Factor Graph, given a priori probabilities and observations. As known, the BP algorithm is a graphical model to represent conditional independence relations of large numbers of random variables.

Since the BP algorithm is a message-passing technique between nodes, and represents an update to the outgoing message from i-th node to j-th neighboring node, we can state that the message from i-th to j-th node related to the local information x i is proportional to

m ji (x i ) ∝ ψ ji (x j , x i )ψ j (x j ) u∈Γj \i m uj (x j )dx j , ( 13
)
where the incoming messages from previous iteration are represented by m uj . This equation represents the message update operation that is performed in the BP's algorithm.

Notice that the BP is capable to compute the exact marginalization in the case of tree-structured graphical models, and this means that [START_REF] Miori | Anticipating health hazards through an ontology-based, iot domotic environment[END_REF] converges in a finite number of iterations, limited to a superior bound, that is the length of the longest path in the graph.

The BP algorithm starts with a "belief updating" phase, where the a posteriori probabilities of the random variable

x i u 2 y 1 y 2 u 1 π(u 2 ) π(y 1 ) π(u 1 ) λ(u 1 ) π(y 2 ) λ(y 2 ) λ(y 1 ) λ(u 2 ) u n π(u n ) λ(u n ) y m λ(y m ) π(y m )
Figure 4: Architecture of the IoT Bayesian network as a graph, for the computation of the BP algorithm.

x i associated to the i-th node, i.e. BEL(x i ), is computed through the information about the evidence coming from the neighboring nodes i.e., BEL(x i ) = αµ(x i ), where µ(x i ) represents the double contribution from "child" and "parent" nodes w.r.t. the i-th node i.e.,

µ (x i ) = λ (x i ) π (x i ) , (14)
with λ (x i ) = j λ xj (x i ) , (15)
π (x i ) = u1,...,un

P ( i | u 1 , ..., u n ) k π xi (x k ), (16)
where j and k are the indexes for the child and parent nodes, respectively. Figure 4 depicts a schematic of an IoT Factor Graph, assumed as a graph with parent and child nodes with respect to the x i node. The computation of the a posteriori probability of the node x i , given all evidence except for the information coming from the j-th child node, is obtained through the parent-to-child message for the child node whose information is excluded. The message from the i-th parent node x i to the j-th child node y j is denoted as π yj (x i ), whose expression is:

π yj (x i ) = α m =j λ ym (x i ) u1,...,un P ( x i | u 1 , ..., u n ) k π xi (u k ) . (17)
Finally, the computation of the conditional probability of the evidence coming from the children of x i given different possible values for the random variable corresponding to the i-th node is obtained through the message exchange from the j-th child node to the k-th parent node as:

λ xi (u i ) = β xi λ (x i ) u k ;k =i P ( x i | u 1 , ..., u n ) k =i π xi (u k ) . (18)
In order to discuss the effects of propagation errors introduced to the BP messages, we can consider multiplicative error functions, which describe the difference between a true message m ji (x i ) and its estimation i.e.,

m ji (x i ) = m ji (x i ) • e ji (x i ), (19)
where m ji represents the message from j-th to the i-th node, and e ji is the error function associated to the message m ji . Then, for the message error propagating from node j to node i, we can apply a particular functional measure d(e ji ) defined as:

d (e ji ) = max xi,xj e ji (x i ) e ji (x j ) , ( 20
)
where e is the error function related to the information x i and x j received from i-th and j-th node, respectively. Then, we have that m ji (x) = mji (x)∀x, if and only if, log d(e ji ) = 0.

To facilitate our analysis, we can split (13) into two parts i.e., (i) message products, and (ii) message convolution. In the first part, we obtain:

M ji (x j ) ∝ ψ j (x j ) u∈Γj \i m uj (x j ), (21)
where as usual, the proportionality constant is chosen to normalize M . We show the message error metric is additive, i.e. that the errors in each incoming message add in their impact on M . The second operation is the message convolution, and we obtain:

m ji (x i ) ∝ ψ ji (x j , x i ) ψ j (x j ) M ji (x j ) dx j , (22)
where M is a normalized message or product of messages. Notice that the log of d(e ji ) is additive, since for several incoming messages muj we have:

subject to certain constraints, such as positivity of the messages and potentials. From (24), we obtain two bounds:

d(e ji ) 2 ≤ d(E ji ) 2 , d(e ji ) 2 ≤ d(ψ ji ) 4 . (25)
IV. CONCLUSIONS

In this paper, we addressed an heterogenous IoT networks scenario, with a plethora of devices for sensing applications. The issue of data sharing and message correction in a multihop IoT environment has been investigated through a Bayesian approach. Specifically, a BP algorithm for message correction, and information update has been presented in its infancy. Future works will address the assessment of the proposed algorithm in an extended simulated scenario.

ACKNOWLEDGMENT

This work has been partially supported by DAHMS project.

Figure 1 :

 1 Figure 1: IoT network model with several subnetworks. The main entities are: IoT devices, equipped with sensing instrumentation, that selectively send data to the gateway and implement the BP algorithm, the gateways, and the cloud.

Figure 2 :Figure 3 :

 23 Figure 2: DAHMS's proof of concept architecture

Figure 2 and

 2 Figure 2 and Figure 3 depict the set of technologies and communication systems that are part of the DAHMS and LOGON proofs of concept.

= d mji m ji 2 .

 2 log d (e uj ). (23) Finally, we can derive a minimum rate of contraction on the errors. Let us consider the message from j-th to i-th node. The error measure d(e ji ) is given by d(e ji ) 2


Stéphane Louvet

Dr Ulf Niebling

Dr Mouham Tanimou

Model Sharing to leverage customer cooperation in the ECU software development

Keywords: Engine Control, embedded software, Model Based Design, automatic code generation, Software Sharing, Model Sharing, Simulink®, AUTOSAR

Motivation

Software Sharing in automotive embedded software development has continuously grown over the last decades and is still getting more importance and attention. Main advantages for Software Sharing, as seen by the OEM, are acceleration of development time and cycle as well as a full flexibility to customize the final product (with introduction of its own software). Although classic software sharing has shown its capabilities in practice, it shows some limitations and need to be extended in order to cover additional use cases. The concept of model sharing can help to address those use cases and the following paper outline five level of model sharing to enhance customer cooperation.

Introduction

An increasing number of car manufacturers develop their own functionalities or their own software platform in order to communize functions across several control units supplied by different Tier 1 companies. This happens in general as stub-software (already compiled) to be integrated in the final Software. This process can however be enhanced to increase development efficiency. In order to achieve such goals, OEMs need the ability to test functions in a very early stage of development and continuously along the V-Cycle development up to the final software. Model Based Design is a mean to do this and the tool Simulink® from MathWorks is becoming a quasi-standard in the automotive industry for this task. Modeling-tools like Simulink® allow OEMs to verify their function in an early phase using virtual and/or rapid prototyping. The so verified model together with the achieved data specification can then be exchanged with suppliers as executable specification; this is one face of the medal of what is called "Model Sharing" at Bosch. The other face consists of provision of development environment as well as use of common modeling guidelines and libraries. Intent for cooperation in this manner always leads first to following questions: What kind of Model Sharing partnership can be established? What are the minimum requirements to set up for Model Sharing cooperation? What are the technical barriers and what are the solutions proposed by Bosch and applied in development ? This paper describes five levels of "model sharing" cooperation between OEMs and suppliers that Bosch has identified, from "Use of ECU virtualization environment" up to "Joint Development". It represents a framework for any contract that has to be established to clearly define fields of responsibilities and guarantees for intellectual properties of each party. Some of the results as well as benefits achieved with these cooperation forms in pilot projects will also be presented.

Model-Sharing

As part of the Model Based Development method, the Model-Sharing is a way of exchanging models and model artefacts in order to ease Software development and cooperation between two or several partners. It is also a way to allow rapid prototyping as well as simulation. Basically, motivations for Model Sharing implementation in embedded SW development process are:

• the OEM focuses on physical modeling, simulation and rapid prototyping. The OEM can also, via model exchanges, take benefit from supplier physical modeling and simulation know-how. • the supplier is in charge of the code generation and its industrialization Simulink® and/or ASCET models can be used as executable specifications for SW development. Using a common tool (with it a common language, and standardized files ) for modeling allows to reduce the amount of to be exchanged documents, and thus enhance the understanding of OEM requirements by the suppliers. This can significantly lead to reduction of development effort and development time as well as delays in the development. Typical use cases for Model Sharing are described in the sections below.

Use case 1: OEM describes the expected functionalities in customer specific functions developed by the supplier

• the desired changes can be implemented directly in the model

Architecture

In order to cope with the complexity of the Software development for Engine Control, a strong and highly flexible Software architecture is required. The answer of Bosch to this challenge is the creation of the Software architecture VeMotionSAR TM that strengthens the cooperation with the OEMs and that has been published in Internet (http://www.bosch-vemotionsar.com). One of the main benefits of this architecture is the modularity of the function description and the definition of the interfaces. It permits to describe individually the various specifications of the corresponding development partners' functionalities in a highly modular architecture domain according to their resources / tools.

Cooperation levels

Since each customer has its own dedicated development process, it has therefore specific requirements for development cooperation with suppliers. On the other hand, suppliers need to classify the requests of OEM to find a common approach. Following this line of thought a five-level classification for model Sharing between Tier1 and OEM has been performed at Bosch: from "Use of ECU virtualization environment" up to "Joint Development". This build a base for discussions with OEMs before starting a new project and it provides a framework for development methodology and contracting. The cooperation usually contains aspects of different types of contracts like service or licensing contract. However, these aspects are typically covered by one agreement describing the complete scope of the cooperation. To build a virtual test bench, virtual engine ECU-SW (engine electronic control unit -software) is needed. PC-executable ECU-functions in a powerful simulation environment are then needed.

Current engine ECU-SW (for Bosch systems) has been developed over many years (through several generations) and shows thereby a high level of maturity. However, this ECU-SW originates from different source artifacts with specifications in various formats (models in different languages and maturity, pdf documents, word documents …) that are not always runnable in common behavior modeling tools (e.g. Matlab/Simulink) environments. During generation process of executables for the target "MEDC control device" these different specifications are dissolved and transferred in machinestranslatable code.

Figure 3 -Workflow for generating executable code in BOSCH environment

For virtual applications, PC-Executables will be generated only for the SW-architecture area ASW (Application SoftWare). However, certain functional characters from the BSW (Base SoftWare) -area (e.g. diagnosis …) have to be made available in order to better represent the behavior of the real control device in the simulation environment. This approach enables to use existing MEDC functions in an easy and flexible manner on a desktop PC as well as to develop own functions or pre-calibrate these functions and/or to verify these functions in combination with Bosch-functions. Furthermore, this environment can be used (a conglomerate of functions) to check the quality of the used plant model. Such a development environment can be realized using the tool INTECRIO from ETAS.

Figure 4 -Integration and configuration platform INTECRIO

The tool INTECRIO (from ETAS) is an integrated environment for prototyping of electronic vehicle functions. It allows users to integrate functional models and code from various sources (e.g. ASCET, MATLAB ® /Simulink ® , C-code and atomic AUTOSAR software components) for different purposes:

• Prototyping in real time on dedicated prototyping devices with interfaces to the real hardware, e.g., to ECU devices and bus signals in the real vehicle (so-called rapid prototyping).

• Prototyping independent of real time on a Windows

® PC using MEDC function models, as well as plant models of and/or recorded signal profiles (so-called virtual prototyping).

For these prototyping purposes, INTECRIO integrates functions with the real time operating system RTA-OSEK which is also used for example on BOSCH production ECUs; thus a real ECU close behavior can be achieved. With this approach much more significant (or meaningful) results can be achieved than this would be possible with a pure model simulation.

Typically results achieved here are prior to the V-Cycle development. This means that this level of cooperation always lead to different additional cooperation levels (to be described in the next sections) that can then be applied for production development.

For Bosch, a next step to enlarge virtual prototyping possibilities will be to introduce tool and features which allows access to more BSW services as today (EEP, ECU coordination states, Com…) Characteristics for contracting on this cooperation level for Model Based Development Bosch to deliver Bosch model as executables object for INTECRIO Behaviour model (e.g. as Simulink and/or ASCET models) or source code are not scope of the delivery to the OEM for this cooperation level The OEM can simulate a large number of Bosch functions together with its own developed functions. The Model sharing contract describes also rights to use for deliverables. All warranty and liability is excluded as the cooperation takes place in the prototyping phase. Ownership of all intellectual property remains with the owner of the models.

Level 1: Processing of OEM Model

This level of model based cooperation is the so called "Commissionedd Development" within the development process at Bosch; the OEM delivers a (simulatable) model that has to be processed by the Tier1 supplier to generate production code without any functional change on the model. Alternatively, the OEM can also deliver a specification (e.g. PDF-document) and the supplier creates a model according to this specification. This model will then be delivered to the OEM for further development.

Using same modeling tool by the OEM and the Tier1 supplier helps to avoid converting the model to another modeling environment. Since model corrections (shall) can only be performed by the OEM, the Tier1 supplier has to wait for a corrected model in case of model errors. And this is a significant drawback that needs to be considered here, especially during the project preparation phase. Use of the updated model to generate the production code is only possible if model preparation can be fully automated. Otherwise differences between previous and updated model have to be identified and changes implemented on Tier1 side in the model (previous or updated one) where the code changes are minimum.

Level 2: OEM Model Enhancement and Processing

This model sharing cooperation level is also considered as part of the development process "Commissioned Development" at Bosch. The OEM delivers a model and the Tier1 supplier can introduce non-functional necessary modifications in the model for several purposes, such as:

• optimisation for production code generation • model corrections in case of mistakes found during model implementation • specific interfaces e.g. for diagnostic In comparison with the cooperation Level 1, there is no need to wait for a correction model from the OEM. However, the OEM shall agree on proposed modifications, for example via a requirement review document. In addition to cooperation Level 1, The Tier1 supplier adds value here to the OEM model. Furthermore the development efficiency is getting enhanced, since the modified model by the supplier is then delivered to the OEM and thus modifications done by the supplier are then available to the OEM and so risks of model discrepancies are reduced since the corrections are done at a single site.

Level 3: Level 3a: Use of OEM Models in Cooperation

The OEM owns a model and delivers it to the supplier who is charge of the code generation for industrialisation. Different to cooperation Level 2, the supplier can introduce functional modifications in the OEM model. The kind of Know-How added by the supplier can be for example specific calculations mechanisms such as special filters, or functionalities related to the OBD (= "On Board Diagnostic"). The modified model is delivered back to the OEM. Naming conventions for labels (e.g. variables and/or calibration parameters) have to be applied according to the OEM rules.. In general, for all cooperation levels, a CIL (Customer Interface Layer) is needed to combine OEM functions with Bosch functions. However, if the OEM uses interfaces from the published architecture VeMotionSAR TM , then there is no need for such a Customer Interface Layer.

Figure 7 -OEM model augmented with BOSCH-IP

Contracting for Model Based Development

The contract between Bosch and the OEM is dealing with close development cooperation for selected topics. The Model Sharing contract describes the rights to use for transferred OEM models and for modified OEM models including Bosch Intellectual Properties. The regulation of rights to use in serial products is typically excluded.

Level 3b: Use of Bosch Models in Cooperation

The supplier provides its own models to the OEM who can then enhance the model by its own functionalities and Know-How. The OEM delivers back the modified model to the supplier who is then in charge of industrialisation and production code generation. The final responsibility for the model (physical Behaviour) remains with the supplier who should review the OEM functional modification for approval. The OEM shall respect the supplier naming convention for the labels like variables or calibration parameters. The development time of OEM specific functions can be reduced drastically since there is no need for the OEM to deliver any specification or schematic that would have to be analysed by the suppliers and then interpreted in the model. Since the model remains property of the supplier, the OEM is not authorised to deliver the model to a 3rd party (e.g. another supplier).

Contracting for Model Based Development

The contract between Bosch and the OEM is dealing with close development cooperation for selected topics. The Model Sharing contract describes the rights to use for transferred Bosch models and for modified Bosch models including OEM Intellectual Properties. All warranty and liability for the use of the models incl. their modifications is excluded. The regulation of rights to use in serial products is typically excluded. It can be typically used in parallel with the Level 0.

Level 4 : Joint development

The joint development deals with a common ownership of the models. The base model can be created at the start of the partnership or can be provided by one of the party. The level of contribution from the supplier and the OEM is comparable and the created Intellectual Properties are shared. It leads to accelerated development: a unique model can be continuously exchanged between the supplier and the OEM or modified by one party with the contribution of the other party. As a constraint, it is necessary to define clear responsibilities between the supplier and the OEM. Regarding Verification and Validation, we can say that for Level 1 and 2, Tier1 is only responsible of Verification, i.e assure that the code generated respect agreed standards (MISRA, ISO…) and that back to back MIL/SIL tests are accepted by OEM. Hence, we can difference introduces with code implemented in fixed point.

For Level 3a, 3b and level4, the functional validation remains under the responsibility of the model owner and most of the time defined and executed by the model owner.

Results

In this section, we will present some results achieved in cooperation at some level presented above

Use case with Level 0

With a German car maker several use cases based on different roles (early calibration, function verification, plant models verification) have been performed.

The plant model used here is an EGT (Exhaust Gas Aftertreatment) (see [1]) model designed in GT-Suite. In the left picture the model predicts very well the transient behavior whereas on the right side, the accuracy of the accumulative tailpipe emissions prediction can be seen to be satisfactory.

Results of this virtual test setup has already been published with [1] MTZ in September 2013

Use case with Level 1

Commissioned development were done in the past at Bosch just by using ASCET or manual code; i.e. since Simulink models couldn't be directly processed, the model itself had to be transformed to ASCET or manually coded. Introduction of support of Simulink toolchain for code generation at Bosch has considerably improved efficiency within "Commissioned Development" projects for those cases where the OEM is able to provide Simulink models to the supplier. In a recent project performed within the scope of this cooperation model, the complete production code for the Application Software Layer has been automatically generated out of the shared models and with respect to the deadlines for a usual ECU development cycle. The code has been generated based on Simulink models and data dictionaries delivered by the OEM. Since the OEM did not expect any functional changes to be performed on his model by Bosch, only small adaptations have been done to make the model ready for code generation. The major challenges were to ensure the functional behavior of model and code, whereas the OEM was not providing any test cases. The approach used in this project was to perform unit tests by doing back-to-back tests with the tool TPT from the company Piketek. A set of stimuli is partly autogenerated with a tool (called GAIO) and permits to compare:

• MIL (Model in The Loop) test with the original non-implemented OEM model

• HIL (Hardware in The Loop) test in Labcar bench with the compiled code in an ECU. The stimuli are injected into the functions via bypasses generated with the tool EHOOKS from ETAS. Approaches at Bosch to ease the Model Sharing with Levels 1 to 4

The Model Sharing approach in a Model Based Development environment might face many technical barriers that can reduce the expected benefits. Based on past and current experiences, several ways have been identified at Bosch to ease the introduction of the Model Sharing in a new project.

Toolchain : An evident pre-requisite to be able to introduce Model Sharing is to use the same tool at OEM and supplier sides. The experience shows that conversion tools (for example to convert a Simulink® model into ASCET) require a huge investment for a result that is not always satisfying and may require a lot of manual rework after conversion. Using the same tool version avoids saving a model into a former version and thus prevents the risk of regressions or incompatibility. It is difficult for a supplier to be able to support different versions of the same tool because in most of the case the tool are customized with specific add-on, especially for the code generation.

Data dictionary : Providing a data dictionary as a model artifact helps a lot the supplier in charge of the code industrialization. The data dictionary can be processed with scripts that reduce the workload and the risk of copy errors. In case of a Simulink® toolchain, it is the main input to generate the Workspace elements. Perl scripts interfaced with Matlab scripts and a Matlab GUI (Graphic User Interface) can be used to import automatically the OEM data dictionary into the Matlab® Workspace. To ease the import process, an XML file (eXtensible Markup Language) is generally preferred to an Excel sheet. The data dictionary can be automatically generated by the OEM if all the labels are stored in a database such as ADD from Visu-IT! Architecture : Before development of group of functions, it is mandatory to define prior the architecture.

A cooperation with the supplier helps to fit the OEM architecture with the supplier one and can drastically reduce the workload to adapt the OEM functions into the supplier's Software (less adapter mappings, no need to cut some supplier functions). The experience shows that too large functions may lead to validation difficulties in case of long calculations chains. In order to ensure an Intellectual Property protection of a part of a function (with Levels 3 and 4) it is possible to hide some subsystems with protected Model Reference mechanism in case of a Simulink model. A preferred solution is to split the functions, it is easier to handle from a contracting and licensing point of view.

MDX / AUTOSAR standard : the standards MDX or AUTOSAR ease the integration of the generated code of a 3 rd party into the supplier Software. It avoids the integrating suppliers developing an integration toolchain for each OEM. The support of the AUTOSAR standard by the code generators is a strong requirement of the ECU suppliers and is a worldwide challenge for the coming years.

Service Library : The use of OEM specific Library routines requires an additional development effort on supplier side, especially for the adaptation is the Library element to code generation. Bosch is using almost exclusively AUTOSAR R4.x service routines with its new Simulink® toolchain and strongly recommends its customers to use this Library. An important aspect is that the AUTOSAR R4.x service routines do not use any AUTOSAR RTE (Real Time Environment) mechanism, it is fully compatible with other standards such as MDX. Despite the fact the AUTOSAR R4.x Library is not available with MathWorks toolbox, Bosch has developed its own set of Simulink blocks and is able to deliver it to its customers after contracting and licensing agreement.

Functions verification and validation : With Model Sharing cooperation there is a good potential to reduce the tests efforts. One interesting way is to reuse the stimuli from the OEM that have been generated during the modeling or Rapid Prototyping phase. Since the supplier is not responsible for the functionalities, it avoids writing test cases. Thanks to the code coverage tool from Simulink® it is possible to know the coverage level (for example with the MC/DC method) of the set of stimuli provided by the OEM and complete the tests by additional test cases in order to achieve a sufficient coverage level that has been agreed with the OEM before starting the project. Back to back tests between the OEM model in floating point and the implemented model in fixed point help the Software developers to detect the calculation errors due to the fixed point implementation (for example : loss of precision)

Summary

Built up of a complete model sharing tool chain (from the scratch) at car manufacturers site requires a clear and stringent roadmap since the investment is considerable and return on invest come along with a long term strategy and sufficient volume of models shared. Model Sharing is currently applied at Bosch for several serial projects and is leading to clear improvements in development efficiency.

Before starting a project, a contract has to be established to clearly define fields of responsibilities and guarantees for intellectual properties of each party. Beside the legal aspects, organizational, software architecture and technical aspects have also to be clarified in order to ease the cooperation.

Bosch is adapting the Process, Method and Tools to car manufacturers having their own approaches. Use of Standards (e.g. ASAM (MDX, CDF …) AUTOSAR (service libraries, methodology …) as well as standardized interfaces specifications (e.g. as published with AUTOSAR R4.x) can significantly ease exchange of models. A flexible tool chain is hence created, since a "one fits all" approach for all customers is not possible.

Bibliography

Figure 1 -

 1 Figure 1 -Different cooperation levels for model sharingThe picture below illustrates those cooperation levels in a contracting and process flow view

Figure 2 -

 2 Figure 2 -Workflow for Model Sharing and contracting Level 0: Virtual Prototyping Since several years, it can be observed that trend for ECU software engineering is moving towards PC simulation in order to 1) Save costs by increased development productivity and 2) Manage complexity and apply holistic engineering approaches. 3) Allow different development speed between HW, mechanics and SW And OEM are requesting from Tier1 supplier a stringent support of development cycle oriented usecases at OEM in a virtual manner e.g. for • ECU function development and testing from idea to production, • Optimization and co-simulation (GT-Power, Matlab) and • Pre-calibration and testing of ECU functions

Figure 5 -

 5 Figure 5 -Workflow for simple processing of OEM model

Figure 6 -

 6 Figure 6 -Workflow for OEM Model enhancement and processing

Figure 8 -

 8 Figure 8 -Joint development: Contribution of both partners with IP Contracting for Model Based Development The contract between Bosch and the OEM is dealing with a common model development for selected topics.The Joint Development contract regulates the rights to use for generated Intellectual Properties including licensing and possible exclusivity phases.

Figure 9 -

 9 Figure 9 -Cumulated AdBlue mass dosed in WHTC Figure 10 -Cumulated NOx emissions by mass in WHTC

Figure 11 -

 11 Figure 11 -Verification workflow for OEM model processing For a specific case where the complete development has been done using Simulink (Matlab 2012b,) and Autosar 3.1 SW-methodology for a diesel air management, following steps have been performed: • Implementation of OEM Simulink Model to fixed point. • code generation of implemented model for Autosar 3. • Performing of Model In the Loop (MIL) tests on basemodel as well as on implemented model using TPT; stimuli were given by the OEM and additional stimuli have been created at Bosch to perform code coverage higher than 75% • Preparation of S-function (for Software in the Loop testing) from the generated code to be used on the ECU • Performing of Software In the Loop tests using the S-function and TPT in Simulink mode (MILlike, but using the compiled code instead of the Simulink model).

Figure 12 -

 12 Figure 12 -Detailed development and verification workflow for OEM model processing

Figure 13 -

 13 Figure 13 -back to back MIL floating point, SIL fixed point

  • Simulation/Rapid Prototyping can be carried out w/o waiting for an "official" SW release Use case 2: OEM develops its owns model and the supplier is in charge of the "so-called" industrialization of the code • This corresponds to Commissioned Development in Bosch process development • The code implementation is simplified if both parties are working with the same modeling environment (examples : Simulink®, TargetLink®, ASCET) • Functionalities can be tested with Rapid/virtual Prototyping before delivery to the supplier • corrected (enhanced) model can be delivered back to the OEM: This helps to avoid additional traceability documents and minimize risk of forgotten improvements in own (original) model

  : 1. Dr. M. Tanimou, Dr. M. Münzenmay, K. Zimmermann, Robert Bosch GmbH; Dr. B. Lumpp, Dr. E. Trapel, E. Bouillon, Dr. M. McMackin, MAN Truck & Bus AG: Software-in-the-Loop durch Co-Simulation von EDC-Modell und GT-Modell; Mainz, 14 th MTZ-Conference -virtual powertrain creation, 2013


H Kurkjian

Y Castin

A Sinatra

Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover

Keywords: numbers: 03.75.Kk, 67.85.Lm, 47.37.+q

Concavity of the collective excitation branch of a Fermi gas in the BEC-BCS crossover

H Kurkjian, Yvan Castin, A Sinatra

I. INTRODUCTION

Cold atomic gases offer a broad flexibility of the microscopic parameters in exploring the many-body problem. In particular, in spin-1/2 Fermi gases, the interaction strength can be adjusted experimentally using Feshbach resonances without inducing strong three-body losses. This degree of freedom, unique among Fermi systems, allowed cold-atom experiments [1][2][3][4][5][6][7][8][9] to study the crossover between a superfluid of Cooper pairs in the socalled Bardeen-Cooper-Schrieffer (BCS) regime of interaction and a superfluid of tightly bound, almost bosonic, dimers in the Bose-Einstein Condensation (BEC) regime. Another advantage of cold atomic gases is the simple theoretical description of the interactions that the cold and dilute regime in which they occur allows for. For a Fermi gas in two spin states ↑ and ↓, one can show that the only significant interactions at low temperature and weak density occur between opposite spins fermions in the s wave and can be fully characterized by a single parameter called the scattering length and denoted by a. In this propitious theoretical framework, entirely analytical studies of experimentally accessible properties of the gas are possible. This article is one of those.

At zero temperature, the three-dimensional spatially homogeneous unpolarized Fermi gas is fully paired and its excitation spectrum consists of two branches: a fermionic branch of excitation of the internal degrees of freedom of the ↑↓ pairs of fermions and a bosonic branch of excitation of their center-of-mass motion, which has a phononic behavior in the long wavelength limit. The latter branch is sometimes said to be collective since it involves a large number of the fermionic modes of internal excitations of the pairs. The fermionic branch is described to lowest order by the BCS theory. To tackle the bosonic branch several approaches have been proposed: The Random Phase Approximation (RPA) of Anderson [10], a Gaussian approximation of the action in a path integral framework [11,12], a Green's function approach associated with a diagrammatic approximation [13], and a linearization of the time dependent BCS variational equations [14]. Re-markably, those theories all lead to the same approximate spectrum of bosonic excitations, which they describe by the same implicit equation.

The concavity of this spectrum has been studied in the weak-coupling BCS limit [11] and, in a qualitative way, in the rest of the BEC-BCS crossover [13,15,16]. A complete quantitative study is thus missing, a gap that this article intends to bridge. In particular, we obtain the spectrum analytically up to order 5 in the wave number q of the center of mass of the pairs. This allows us to conclude on the concavity of the branch of excitation in a neighborhood of q = 0 over the whole BEC-BCS crossover.

Several physically relevant problems can be addressed after our study. First, the processes that dominate the collective mode damping at low temperature can be identified. If the branch is convex over a neighborhood of q = 0 then the Landau-Beliaev 2 phonons ↔ 1 phonon interaction processes [17,18] dominate while if it is concave, those processes are forbidden by momentum and energy conservation and the Landau-Khalatnikov 2 phonons ↔ 2 phonons processes [19] take over. At low temperature the contribution of the gapped fermionic branch to the collective mode damping is exponentially small [20]. Second, the quantitative knowledge of the concavity parameter γ is required in order to predict the phonon damping rate due to the 2 phonons ↔ 1 phonon processes in the convex case beyond the quantum hydrodynamics approximation [16], or due to the 2 phonons ↔ 2 phonons processes in the concave case where the effective interaction predicted by quantum hydrodynamics involves virtual 2 phonons ↔ 1 phonon processes and depends on γ. Finally, the knowledge of γ gives access to the phase diffusion coefficient of the condensate of pairs, a quantity responsible for an intrinsic and fundamental limit to the coherence time of the gas [14].

II. THE RPA EQUATION ON THE EXCITATION SPECTRUM

The RPA equation yielding implicitly the energy ω q of the collective excitations as a function of their wave vector q is I ++ (ω q , q)I --(ω q , q) = 2 ω 2 q [I +-(ω q , q)] 2 .

(1)

The collective nature of the bosonic modes is visible in the quantities I σσ ′ which are integrals on the relative wave vector k of the pairs and depend on ω 2 q via the denominator of their integrand [13]:

I ++ (ω, q) = d 3 k (ǫ k+q/2 + ǫ k-q/2 )(U k+q/2 U k-q/2 + V k+q/2 V k-q/2 ) 2 2 ω 2 -(ǫ k+q/2 + ǫ k-q/2 ) 2 + 1 2ǫ k , (2)
I --(ω, q) = d 3 k (ǫ k+q/2 + ǫ k-q/2 )(U k+q/2 U k-q/2 -V k+q/2 V k-q/2 ) 2 2 ω 2 -(ǫ k+q/2 + ǫ k-q/2 ) 2 + 1 2ǫ k , (3)
I +-(ω, q) = d 3 k (U k+q/2 U k-q/2 + V k+q/2 V k-q/2 )(U k+q/2 U k-q/2 -V k+q/2 V k-q/2 ) 2 ω 2 -(ǫ k+q/2 + ǫ k-q/2 ) 2 . ( 4
)
We have introduced the amplitudes U k and V k and the eigenenergies ǫ k of the fermionic eigenmodes of the BCS theory [21]:

ǫ k = 2 k 2 2m -µ 2 + ∆ 2 , (5)
U k = 1 2 1 + 2 k 2 2m -µ ǫ k , (6)
V k = 1 2 1 - 2 k 2 2m -µ ǫ k , ( 7
)
where m is the mass of a fermion. The two natural parameters of the BCS theory, with which we express the energy ω q , are the chemical potential µ, identical for the two spin states, and ∆, the gap in the BCS spectrum of the fermionic excitations when µ is positive. If needed, they can be replaced by the scattering length a and the total density ρ of the gas by inverting the two following relations [START_REF] Castin | Ultra-cold Fermi Gases[END_REF][START_REF] Kurkjian | [END_REF]:

m 4π 2 a = d 3 k (2π) 3 m 2 k 2 - 1 2ǫ k , (8)
ρ = d 3 k (2π) 3 2|V k | 2 . ( 9
)
In practice, the Fermi wave number k F defined by ρ = k 3 F /(6π 2 ) is often used instead of the density.

III. GLOBAL NUMERICAL STUDY OF THE CONCAVITY

From a numerical solution of the dispersion equation (1) we obtain the dispersion relation q → ω q over its existence domain. We show an example in Figure 1 in the unitary limit 1/(k F a) = 0, where ∆/µ ≃ 1.162 according to the BCS theory. Rather than q → ω q itself, we plot as a black solid line the function q → ω q -cq, where c is the speed of sound and q → cq is the linear part of the spectrum. The concavity properties of this function are the same as those of q → ω q , but they are more visible graphically at low q, because they are not masked by the dominant linear part, which anyway plays no role in the selection of the damping processes mentioned at the end of the Introduction. In the figure, it is apparent that the dispersion relation exhibits an inflection point at q/k µ ≃ 0.795, where k µ = (2mµ) 1/2 / , separating a low-q interval over which the dispersion relation is convex and a high-q interval over which it is concave. As a consequence, according to the RPA, the leading damping processes of the collective excitations of the unitary gas at low temperature are the Beliaev-Landau 2 phonons ↔ 1 phonon processes. Figure 2 synthesizes our numerical results on the concavity of the bosonic branch for all values of ∆/µ and q/k µ . Several domains can be identified depending on the values of ∆/µ, or equivalently 1/(k F a). (i) When the scattering length is negative a < 0, that is for ∆/µ < 1.162 (lower dotted line), the existence domain of the solution to Eq.( 1) as a function of q is compact and simply connected [13], hence of the form [0, q sup ]. The dispersion relation is entirely concave for ∆/µ < 0.869 while for 0.869 < ∆/µ < 1.162 it is first convex at small q and then concave. Between those two zones it goes through an inflection point, whose position q infl (∆/µ) we compute analytically in the small q limit [see the black dashed line and Eq.( 29) of section IV]. (ii) On the other side of the resonance (a > 0) and up to ∆/µ = 1.729 (upper dotted line), the existence domain of the solution to Eq.( 1) splits up into two connected components a) -1 = 0, which is here ∆/µ = 1.162, the dispersion relation of the collective excitations is plotted as a black solid line after subtraction of its phononic part cq to better reveal its concavity properties.

The wave number q is in units of kµ = (2mµ) 1/2 / . The linear part q → cq is shown for comparison (steep violet dashed straight line), as well as the cubic (red dashed upper line) and the quintic (blue dashed lower line) approximations of ωq -cq, see Eq. ( 10). The second order derivative d 2 ωq/dq 2 , plotted as a gray solid line, vanishes at the inflection point q infl ≃ 0.795kµ marked by the vertical dotted line. To the left (right) of this point the dispersion relation q → ωq is convex (concave).

[0, q sup ] and [q inf , +∞[ [13]. While the branch is always convex in the second component [q inf , +∞[, the first one [0, q sup ] exhibits interesting variations: From small q to high q, the branch is convex and then concave for 1.162 < ∆/µ < 1.22, convex, concave and then convex again for 1.22 < ∆/µ < 1.710 and finally entirely convex for 1.710 < ∆/µ < 1.729. (iii) When ∆/µ is greater than 1.729, or negative, the two components of the existence domain merge and a solution to Eq.( 1) exists for all q [13]. The branch is then entirely convex.

All these numerical values are predicted by the RPA or the BCS theory. They are therefore approximate. Up to now the only value that can be compared to experiments is that of ∆/µ at the unitary limit: From the measured values ∆ ≃ 0.44 2 k 2 F /(2m) [24] and µ ≃ 0.376 2 k 2 F /(2m) [8] we get ∆/µ ≃ 1.17, which is remarkably close to the BCS-theory prediction ∆/µ ≃ 1.162. One must also keep in mind that the RPA spectrum results from a linearized treatment of the pair-field quantum fluctuations, which neglects the interactions among the bosonic quasiparticles. In reality, these interactions will shift the eigenenergies ω q . They will also give rise, even at zero temperature, to an imaginary part in ω q , corresponding to a finite lifetime of the excitations, provided that the concavity of the dispersion relation allows for resonant 1 phonon → 2 phonons Beliaev processes [16]. FIG. 2: (Color online) Local concavity of the bosonic branch q → ωq depending on ∆/µ and on the wave number q in units of kµ = (2mµ) 1/2 / . The values of 1/(kF a) corresponding to those of ∆/µ are given on the right vertical axis. Red points show that the branch is locally concave and orange points that the branch is locally convex. The thick solid line is the border of the existence domain of the solutions of the dispersion equation ( 1). The black dashed line is the low-q analytical prediction (29) of the boundary between the red zone and the orange zone, that is of the locus of the inflection points of the curve q → ωq, indicated by a thin solid line.

The ordinate ∆/µ = 0.869 of the point where this line meets the q = 0 axis, and the ordinate ∆/µ = 1.710 above which the concavity zone disappears are indicated by arrows. For the values of ∆/µ or 1/(kF a) in between the two dotted lines indicated by arrows, the q-existence domain of the solution of equation ( 1) is not simply connected. The lower dotted line indicates the unitary limit, where the scattering length diverges |a| → +∞.

IV. ANALYTICAL STUDY OF THE CONCAVITY IN THE LONG WAVELENGTH LIMIT

The dispersion relation can be obtained analytically in the long wavelength limit q → 0. To this end we expand the eigenenergy of the collective mode up to order 5 in q:

ω q = q→0 cq 1 + γ 8 q mc 2 + η 16
q mc 4 + O q mc 6 . ( 10
)
To lowest order, as for any superfluid system, the energy is phononic with a sound velocity given by the hydrodynamic expression

mc 2 = ρ ∂µ ∂ρ a , (11)
where the derivative is taken for a fixed scattering length a, as indicated by the notation. When applied to the approximate equation of state (9), the hydrodynamic expression (11) gives the RPA sound velocity, as shown in reference [13] by an expansion of the solution ω q of equation (10) up to first order in q. We give here an explicit expression in the form of a rational fraction

mc 2 µ = 2(xy + 1) 3(y 2 + 1) , (12)
in terms of the variables

x = ∆ µ and y = ∂∆ ∂µ a . (13)
In turn, the y variable is written as a function of x,

y = d 3 k 2 k 2 /(2m)-µ ǫ 3 k d 3 k ∆ ǫ 3 k = +∞ 0 du u 2 (u 2 -x -1 ) [(u 2 -x -1 ) 2 +1] 3/2 +∞ 0 du u 2 [(u 2 -x -1 ) 2 +1] 3/2 (14)
by taking the derivative of equation ( 8) with respect to µ at fixed scattering length a and by expressing the wave vectors in units of k ∆ = (2m∆) 1/2 / to form the dimensionless integration variable u. The integrals over u in the numerator and in the denominator of the right-hand side of equation ( 14) may be expressed in terms of complete elliptic integrals of the first and second kinds [11]. At the unitary limit one has y = x, since ∆ and µ are proportional due to scale invariance.

To obtain the expression (12) of the reduced sound velocity, we take the derivative of the equation of state (9) with respect to µ at fixed a and we express all resulting integrals as functions of x and y using ( 14) and the relation

d 3 k (2π) 3 ∆ 3 ǫ 3 k ρ = +∞ 0 du u 2 [(u 2 -x -1 ) 2 +1] 3/2 +∞ 0 du u 2 1 - u 2 -x -1 [(u 2 -x -1 ) 2 +1] 1/2 = 3x/2 1 + xy , (15)
which may be derived by integrating by parts the integral over u in the denominator (u 2 is the function to be integrated).

To obtain the dimensionless coefficients γ and η of the terms of ω q of higher order in q in equation ( 10), we cannot rely on known thermodynamic expressions and we must face the laborious expansion of equation ( 1) in powers of q. Still the result can be put into a rational form in terms of x and y:

γ = 4 i=0 P i (y)x i 135x 2 (x 2 + 1) (y 2 + 1) 3 , (16) η
= 8 i=0 Q i (y)x i 1020600 (y 2 + 1) 6 x 4 (x 2 + 1) 2 . ( 17
)
The P i (y) and Q i (y) polynomials that appear as coefficients of x i in the numerators of γ and η are given by P 0 (y) -4 13y 4 + 16y 2 + 8 , P 1 (y) = 4y 13y 4 + 41y 2 + 8 , P 2 (y) = 50y 6 -21y 4 -252y 2 -61, (18) P 3 (y) = 2y y 4 + 32y 2 + 71 , P 4 (y) = 35y 6 + 56y 4 -13y 2 -54, and by Q 0 (y) = 16 7745y 8 + 19528y 6 + 20304y 4 + 8384y 2 + 1088 , Q 1 (y) = 32y 2857y 8 + 67y 6 -3186y 4 -7920y 2 -2624 , Q 2 (y) = -8 12882y 10 + 28061y 8 -26936y 6 + 7221y 4 -24496y 2 -5232 , Q 3 (y) = -8y 8456y 10 -9859y 8 + 9977y 6 + 145295y 4 + 3523y 2 + 23720 , Q 4 (y) = -17500y 12 -247996y 10 -1249743y 8 -1341332y 6 + 337202y 4 -694392y 2 + 18321, (19) Q 5 (y) = -4y 25564y 10 + 36027y 8 -66984y 6 + 92206y 4 + 387932y 2 -56121 , Q 6 (y) = -2 12250y 12 + 115637y 10 + 558246y 8 + 1071518y 6 + 589478y 4 -248499y 2 + 53082 , Q 7 (y) = -4y 12957y 10 + 33764y 8 -41904y 6 -173106y 4 -96189y 2 + 53406 , Q 8 (y) = -8575y 12 -44544y 10 -149742y 8 -360644y 6 -477615y 4 -270756y 2 -20412.

Our analytical expressions ( 16) and ( 17) result from a Taylor expansion of the integrals I ++ (ω q , q), I --(ω q , q) and I +-(ω q , q) after replacement of ω q with the expansion (10). At each order, we reuse the results of the lower orders, that is the value (12) of c to obtain γ, and then those of c and γ (16) to obtain η. We encounter integrals involving in the denominator high powers of ǫ k (or of [(u 2 -x -1 ) 2 + 1] 1/2 after the k-to-u change of variable). They can be evaluated by repeated integration by parts, as explained in Appendix A.

We have plotted in Figure 3 the coefficients γ and η as functions of the parameter 1/(k F a) (which we have preferred here to ∆/µ). Let us briefly review their asymptotic behaviors in the BEC 1/(k F a) → +∞ and BCS 1/(k F a) → -∞ limits, and their values in some specific relevant regimes. FIG. 3: Dimensionless coefficients (a) γ and (c) η of the terms q 3 and q 5 in the RPA dispersion relation of the bosonic excitations for q → 0, see equation (10), as functions of the interaction parameter 1/(kF a). In the BEC limit, they have a finite limit corresponding to the Bogoliubov dispersion relation, see equations (21,[START_REF] Castin | Ultra-cold Fermi Gases[END_REF]. In the BCS limit, they diverge, and it is more appropriate to express the wave number q in units of ∆/( c), which amounts to considering the quantities (b) γ∆ 2 /(m 2 c 4 ) and (d) η∆ 4 /(m 4 c 8 ) which have a finite limit. In between these two limits, γ and η vanish and change sign for 1/(kF a) ≃ -0.144 and 1/(kF a) ≃ -0.389 respectively, and they are weakly positive (γ ≃ 0.0838) or weakly negative (η ≃ -0.0354) at the unitary limit. The value of η is relevant mainly at the point where γ vanishes. At this point, η ≃ -0.0428 so the RPA dispersion relation is concave over a neighborhood of q = 0.

BEC limit

In the BEC limit k F a → 0 + , the system is equivalent to a weakly interacting gas of bosons of mass 2m, with a chemical potential

µ B = 2µ -E dim , (20)
where E dim = -2 /(ma 2 ) is the internal energy of a dimer [25,26]. The dispersion relation of the bosonic excitations is then known to be convex and to take the Bogoliubov form at chemical potential µ B when q = o(1/a) [13]:

ω Bog q = 2 q 2 4m 2 q 2 4m + 2µ B 1/2 , (21)
in which case the sound velocity is given by 2mc 2 = µ B .

The coefficients γ and η are thus expected to have the following limits:

γ → kF a→0 + 1 4 and η → kF a→0 + - 1 128 . (22)
This is confirmed by equations ( 16) and (17) taken at the BEC limit, that is for

x = O(k F a) 3/2 → 0 -[22]
and, as shown by equation ( 14) after the change of variable u = u ′ /|x| 1/2 , for y ∼ -4/x.

BCS limit

When k F a → 0 -, the lower border of the twofermionic-excitation continuum (at fixed total wave number q < q sup ) becomes exponentially weak and forces the bosonic excitation branch, which cannot enter into this continuum, to bend downward [13,15]. This bending takes place over a wave number range q c such that cq c = ∆, that is such that the leading term in the expansion (10) is of the order of ∆. This means that the collective modes are affected by the pairs internal structure when the mode wavelength becomes comparable to the pairs size in real space ≈ 2 k F /(m∆) [11], a quantity that is indeed of order 1/q c since c ∝ k F /m in the BCS limit. This qualitatively explains why the dispersion relation is a concave function of q for low q and low ∆/µ in Figure 2, and why γ < 0 in the BCS limit.

More quantitatively, we expect that the normalized energy ω q /∆ becomes a universal function of q/q c = cq/∆ when k F a → 0 -, in which case all the terms in square brackets in equation ( 10) are of the same order of magnitude for q = q c , that is |γ( q c /mc) 2 | ≈ 1, |η( q c /mc) 4 | ≈ 1. This is indeed what we find by taking the limits x → 0 and y → 0 in equations ( 16) and (17):

γ ∼ kF a→0 - - 8 15
mc 2 ∆ 2 and η ∼ kF a→0 - 136 1575 mc 2 ∆ 4 , (23)
the first result reproducing that of reference [11].

Crossover region

In the crossover region between BEC and BCS, γ is an increasing function of 1/(k F a). It vanishes and changes sign for the value x 0 of x = ∆/µ given by x 0 ≃ 0.868567. (24) This value corresponds to 1/(k F a) ≃ -0.144292, in agreement with the numerical result of Figure 2 and with reference [16]. The RPA prediction is then that the dispersion relation of a unitary gas is convex close to q = 0:

γ ≃ (kF a) -1 =0 0.083769. ( 25
)
The coefficient η changes sign for a value x 1 of ∆/µ given by

x 1 ≃ 0.566411 (26)
corresponding to 1/(k F a) ≃ -0.389027. It is negative both at unitarity η ≃

(kF a) -1 =0 -0.035416, (27)
and at the point x 0 where γ = 0:

η(x 0 )≃ -0.042794. (28)
At that very point the sign of η is important as it determines the concavity of the dispersion relation close to q = 0.

Locus of the inflection points

The coefficients γ and η allow us to find analytically the border between the orange and the red zones of Figure 2 for small q, that is the ensemble of points with coordinates (q infl /k µ , ∆/µ) such that the second derivative d 2 ω q /dq 2 is zero. Using the expansion (10) for ω q and expanding the coefficients γ(x) and η(x) around x = x 0 , to order one and order zero in x -x 0 respectively, one obtains the equation

q 2 infl k 2 µ ∼ x→x + 0 - 3γ ′ (x 0 ) 10η(x 0 ) mc 2 µ (x -x 0 ) ≃ 2.015858(x -x 0 ) (29
) plotted as a black dashed curve in Figure 2, which reaches the axis q = 0 with a horizontal tangent.

In contrast, the border between the red and orange zones in Figure 2 reaches the border of the existence domain of the collective excitation branch with an oblique tangent. This is due to the fact that the third derivative of q → ω q is nonzero at the contact point q = q sup contrarily to what happens at q = 0.

V. CONCLUSION

We have considered a spatially homogeneous unpolarized gas of spin-1/2 fermions at zero temperature, and we have obtained analytically the spectrum ω q of the bosonic excitation branch predicted by the RPA up to order 5 included in the wave vector q close to q = 0. The coefficients of the obtained expansion are rational fractions of two variables ∆/µ and (∂∆/∂µ) a , where the second variable can be analytically related to the first one using the BCS equation of state. This allows us to show analytically that the dispersion relation q → ω q is concave close to q = 0 when 1/(k F a) is between -∞ and a value close to -0.144, a point where the first correction to the linear dispersion relation is of order q 5 with a slightly negative coefficient. For -0.144 < 1/(k F a) < 0.157 the branch is convex close to q = 0 but becomes concave when q increases, and it remains so for q increasing up to the maximal possible value q sup if 1/(k F a) < 0.022, while it becomes convex again in the opposite case 0.022 < 1/(k F a) < 0.157. Beyond 1/(k F a) = 0.157 the bosonic branch is convex over its whole existence domain.

A straightforward application of our quintic approximation (10) for the spectrum is to determine if a lowq collective excitation of the Fermi gas can decay via a Beliaev process, that is into two collective excitations of wave vectors q 1 and q 2 = qq 1 . Energy conservation allows such a process if ω q > Ω inf q where Ω inf q = inf q1 (ω q1 + ω q-q1 ) is the lower border of the twoexcitation continuum at fixed total wave vector q. If ∆/µ is away from the critical value x 0 ≃ 0.869 where the coefficient γ of q 3 vanishes, the dispersion relation is either entirely convex or entirely concave at low q, and the Beliaev decay is respectively allowed or forbidden. If ∆/µ is close to the critical value x 0 , the dispersion relation has an inflection point at low q. We then apply the analysis of reference [15] to Eq. ( 10) and we find Ω inf q = min(2ω q/2 , ω q ) [START_REF]Since the dispersion relation q → ωq is here an increasing function, only the cases (i) and ([END_REF]. The Beliaev decay is thus allowed if q mc 2 < 8γ 5|η| (30) to leading order in γ, that is in ∆/µ -x 0 . At low but nonzero temperature, there exist additional decay mechanisms: (i) The Landau mechanism q+q 1 → q 2 is forbidden whenever the Beliaev one is since it can be viewed as an inverse Beliaev mechanism with an initial wave vector q 2 of modulus > q; (ii) the higher-order Landau-Khalatnikov decay process q + q 1 → q 2 + q 3 is always allowed, but it is subleading when the Beliaev or Landau processes are present.

Our results on the concavity of the dispersion relation close to q = 0 can be tested experimentally in a gas of cold atoms trapped in a flat bottom potential [27]. This can be done either (i) indirectly by measuring dissipative effects such as the damping of collective excitations at low temperature, or dispersive effects such as the spreading of a wave packet of sound waves created by a laser pulse [9,28], or (ii) directly by accessing the dispersion relation at low temperature via the dynamic structure factor of the Fermi gas through Bragg excitation at a selected wave vector q [29-31] [33]. Supplemented by kinetic equations for the collective mode occupation numbers, our results open the way to an analytical determination of the phase diffusion coefficient at low temperature, hence to the intrinsic limit to the coherence time of the condensate of pairs in a finite-size Fermi gas [14].

FIG. 1 :

 1 FIG.1:(Color online) At unitarity (kF a) -1 = 0, which is here ∆/µ = 1.162, the dispersion relation of the collective excitations is plotted as a black solid line after subtraction of its phononic part cq to better reveal its concavity properties. The wave number q is in units of kµ = (2mµ) 1/2 / . The linear part q → cq is shown for comparison (steep violet dashed straight line), as well as the cubic (red dashed upper line) and the quintic (blue dashed lower line) approximations of ωq -cq, see Eq.(10). The second order derivative d 2 ωq/dq 2 , plotted as a gray solid line, vanishes at the inflection point q infl ≃ 0.795kµ marked by the vertical dotted line. To the left (right) of this point the dispersion relation q → ωq is convex (concave).

Appendix A: Expressing integrals in terms of the variables x and y

In the expansion of I ++ (ω q , q), I --(ω q , q) and I +-(ω q , q) at low q, and after a rescaling of the wave vectors by k ∆ as in equation ( 14), k = uk ∆ , we encounter integrals of the form

and the total density ρ is given by equation (9). The integrals giving I n,p and J n,p are convergent for n-p ≥ 2 and for n -p ≥ 3, respectively. Integrals that depend on the direction of k can be expressed in the forms (A1) and (A2) after angular integration:

where f (k) is an arbitrary function of the modulus of k.

Let us first establish the four recurrence relations

holding under the conditions 1 ≤ p ≤ n -4 for the first relation, 1 ≤ p ≤ n -3 for the second one, 1 ≤ p ≤ n -2 for the third one and 1 ≤ p ≤ n -3 for the last one. In order to derive the relation (A6), we integrate by parts the integral

selecting u → u 2p+1 ξ u as the function to be differentiated. In order to derive the relation (A7), we integrate by parts the integral defining J n,p in equation (A2), selecting u → u 2p+1 as the function to be differentiated.

In both cases, we note that the function u → uξ u /ǫ n u admits the primitive u → -[2(n -2)ǫ n-2 u ] -1 . Finally, we simply write u 2p+2 = u 2p (ξ u + x -1 ) in the integrand of (A1) in order to obtain (A8), and we write

in the integrand of (A2) in order to obtain (A9). This procedure generalizes that of reference [11].

We now show by induction using the relations (A6,A7,A8,A9) that I n,p and J n,p can be expressed as functions of I 3,0 and J 3,0 , for all odd n ≥ 3 and for all positive p within the existence domain of the integrals. Let n be odd and ≥ 3 and assume that we know all the I n,p , 0 ≤ p ≤ n -2, and all the J n,p , 0 ≤ p ≤ n -3. Then (i) I n+2,1 and J n+2,1 can be deduced using (A6) and (A7), (ii) using (A9) and (A8) we obtain a Cramer system for I n+2,0 and J n+2,0 :

which we solve, (iii) we use (A8) and (A9) to access the values of I n+2,p and J n+2,p for p ≥ 2. We set the induction basis at n = 3, by expressing I 3,1 as a function of I 3,0 and J 3,0 due to (A8).

Finally we relate I 3,0 and J 3,0 to x and y thanks to the relations ( 14) and (15) of the main text, which take the form y = J 3,0 /I 3,0 and I 3,0 = 3x/[2(1 + xy)] with the notation of this appendix.


Rémy Mensire

Keyvan Piroird

Elise Lorenceau

Capillary imbibition of aqueous foams by miscible and non-miscible liquids

published or not. The documents may come

Capillary imbibition of aqueous foams by miscible and nonmiscible liquids

INTRODUCTION

The earliest studies on impregnation of a single vertical capillary tube from an infinite reservoir date back to the 1700's with the pioneer works of Jurin, Newton and later on, Laplace. Subsequent studies have shown that capillary rise exhibits several dynamical regimes, such as the famous Lucas-Washburn t 1/2 scaling behavior for early times [1,2].

In such a model geometry, the radius of the tube is the only length scale of the system and both sets the amplitude of the capillary pressure drop and the viscous dissipation. In homogeneous porous media with more complex geometries such as packed beds of granular beads, dry soils or paper, various dynamics have been experimentally reported. The imbibition dynamics of porous media characterized by a unique length scale, is logically described by the t 1/2 law [2][3][4][5]. Yet, other dynamics also occur in more complex media, such as deformable or heterogeneous media, where a unique length scale can not be defined unequivocally. The capillary rise of water in column of glass beads follows a t 1/4 regime for late times due to micronic heterogeneities of the glass surface [4,6,7] and deformable sponges wick with an unexplained t 1/4 dynamics at late times [7]. A similar t 1/4 dynamics is also observed in heterogeneous media with gradients of pore size [8] and the imbibition between flexible boundaries, which results from an interplay between capillary, gravity and elasticity of the boundaries, exhibits different regimes [9].

In this context, dry aqueous foams, which provide perfectly smooth boundary conditions for the liquid, are interesting examples of complex deformable porous media. In those porous systems, the liquid is distributed between the thin films, the slender liquid-filled channels, called Plateau borders and the nodes as depicted in Figure 1-b) andc). When a dry liquid foam is put into contact at its bottom with its constitutive foaming solution, the liquid spontaneously rises within the liquid network following a t 1/2 law as experimentally observed in microgravity conditions [10][11][12]. On Earth, theoretical predictions concerning semi-infinite foams with an initial liquid fraction equal to zero report that the capillary rise continues indefinitely despite gravity, following a t 1/3 law [13] that recalls capillary rise observed in corners [14][15][16]. However, to the best of our knowledge, this scaling law has never been observed experimentally.

In this letter, we therefore consider the imbibition dynamics of dry aqueous foams. We use both miscible liquid or an immiscible oil. Indeed, the remarkable imbibition ability of a liquid foam also holds when it is put into contact with an immiscible liquid, provided that it is not a defoaming agent for the foam, as recently shown at the scale of a single Plateau border or in confined geometries [17][18][19].

The foaming solution we use is described in [17,20,21] and consists in a mixture of different components (cocoamidopropyl betaine, sodiumlaurylethyl sulfate, myristic acid and glycerol). It provides rigid boundary conditions of the air/foaming interfaces due to the presence of the fatty myristic acid [20,21]. The air/foaming solution surface tension, viscosity and density of this mixture are γ aw = 23.7 ± 1 mN/m, η F = 1.4 mPa.s and ρ F = 1026 kg/m 3 . The liquids of viscosity and density η and ρ, are either miscible or immiscible with the foaming solution. Aqueous miscible solutions with viscosity η ranging from 1.4 to 100 mPa.s, are obtained by adding glycerol to the foaming solution. To maintain the architecture of the foam unaltered during the imbibition of immiscible fluids, we use two different oils which do not exhibit any defoaming activity. This requires that the oil droplets trapped within the foaming liquid must neither penetrate the interface nor spontaneously spread at the air/foaming solution interface. In other words, the oil droplets in the foaming solution do not wet the air/foaming solution interfaces: the immiscible liquid is non-wetting. Thus, E = γ aw + γ ow -γ ao and S = γ aw -γ ow -γ ao , the entry and spreading coefficients defined using γ aw , γ ow and γ ao the air/foaming solution, oil/foaming solution and oil/air surface tensions, must both be negative [20,[START_REF] Garrett | Defoaming: theory and industrial applications[END_REF]. With our foaming solution and organic oils, such as olive oil or sunflower oil, both E and S exhibit negative values, as shown in Table I. The viscosity η and the density ρ of the two oils are also displayed in Table I.

To generate the foam, a slow flow of gaseous nitrogen is blown through a syringe needle into a 27 cm-high rectangular column with a 6 cm-wide square base filled with the foaming solution. We use different syringe needles to obtain monodisperse foams with an average radius R of 0.5 mm, 1.8 mm or 2.5 mm. Besides, the 5-cm high upper part of the column is removable. Thus, we can extract at the top of the column a foam sample for which the volume Ω and the mass m are known and measured by using a precise weighing scale. We deduce the liquid fraction of the sample from = m/Ωρ and typically 0.7.10 -4 < < 1.2.10 -3 .

With this set-up, we can therefore vary independently and R and thus obtain curvature Characteristics of the oils used for imbibition with the foaming solution described in [20]. γ aw , γ ow and γ ao are the air-foaming solution, oil-foaming solution and oil-air surface tensions in mN/m, η is the dynamic oil viscosity in mPa.s, ρ is the density in kg.m -3 . E and S are the entry and spreading coefficient in mN/m.

radii of the Plateau border, noted r F , between 15 and 85 µm (see Figure 1-c)). Indeed, in the limit where < 10 -2 , most of the liquid is comprised within the Plateau borders and [START_REF] Cantat | Foams : Structure and Dynamics[END_REF][START_REF] Weaire | The Physics of Foams[END_REF].

r F = √ 3 √ R from volume conservation
Given that the contrast of optical indexes between the foam and the liquid is small or null, we add a small quantity of fluorescent dyes to the liquid at a concentration of 10 -2 g/g (fluorescein from Sigma-Aldrich for aqueous solutions and Yellow Black from Rohm and Haas for organic oils). Those fluorescent markers, which are trapped in the liquid, provide a means of visualizing the swollen part of the foam [19]. With this set-up and under illumination with excitation at 488 nm, the liquid is luminous while the aqueous foam -only constituted of air and foaming solution -is not visible. The fluorescence intensity is followed with a camera recording at typically 10 frames per second.

Figure 1-d) illustrates such an experiment where a dry dark foam is put into contact with a large 20 µL fluorescent oil drop sitting on a solid surface. We first observe a lateral spreading of the oil drop, squeezed between the solid surface and the foam. Yet, after ten seconds, oil wicks into the foam, revealing the foam architecture and forming a complex oilladen foam structure. After 5 minutes, oil has risen in the foam up to h = 10 mm. The foam film withstands without any problem the oil invasion: we barely observe any film break-up during the whole process. This confirms that the imbibition only occurs within the Plateau borders and the nodes and does not interfere with the fragile structure of the foam films as previously observed at the scale of a single Plateau border [17]. Yet, the structure of this oil-laden foam is not frozen: elementary topological changes known as T1 are frequently observed due to liquid redistribution [19]. Our fluorescent set-up allows us to follow h, the upper boundary of the imbibition front in the foam as defined in the fifth frame of Figure 1-d). Figure 2 shows h as a function of time t, when a 20 µL olive oil drop is put into contact with foams at various liquid fractions. As classically observed in porous media, the early times dynamics is fast. Moreover, here, the drier the foam, the faster the early-time dynamics, a feature which underlines the importance of the liquid fraction of the foam prior to imbibition. When plotted in log-log scale, those curves do not exhibit a well-defined power law as can be seen on the inset of Figure 2, where two straight lines corresponding to the t 1/2 early and t 1/4 late-time dynamics are drawn. We also emphasize that for the driest foam, the front velocity is small: typically, dynamics with the diffusive coefficient given in [12], while the dashed line illustrates the t 1/4 dynamics.

ḣ ∼ 0.2 mm/s at t = 10 s. For a centimeter length scale, the Reynolds number is around Re = ρ ḣL/η ∼ 10 -8 << 1. Note that the imbibition of glycerol or sunflower oil exhibits the same trends.

The rigid boundary conditions provided by the foaming solution and the small value of Re suggest that the fluid motion in the foam network corresponds to a Stokes flow where the capillary driving force is balanced with the viscous dissipation. In the limit where γ ow /γ aw < 1, recent experimental works at the scale of a single horizontal Plateau border have revealed the following features [17,[START_REF] Cohen | [END_REF]. On one hand, the driving capillary pressure that sustains the imbibition in the Plateau border is set by γ aw /r F , where r F is the curvature radius of the Plateau border prior to imbibition. On the other hand, the bulk viscous dissipation only occurs in the oil phase and scales as η ż/r 2 (z, t), where r 2 (z, t) is proportional to the swollen Plateau border cross section. In reference [17], the experiments were either driven at constant flow rate or constant volume of oil, thus r(z, t) was deduced from mass conservation [26], while here imbibition occurs from a large reservoir. Thus, to gain quantitative understanding of the spatial and temporal evolution of r, we take advantage of the fluorescence intensity of our images. In Figure 3, we report the spatial evolution of the fluorescence intensity of a single node located at x = 1 mm and z = 7.9 mm for different times. The oil front reaches the height z = 7.9 mm at t = 270 s, thus the signal prior to t = 270 s corresponds to noise fluctuation. At t = 274 s, a peak in the fluorescence intensity is observed. The intensity of this peak grows until reaching a steady shape observed 100 s after the front has gone through the node, as can be seen in the inset of Figure 3 where a saturation of the peak's width is observed after t = 400 s. This suggests that for the foam below the oil front, the distribution of oil quickly reaches a steady-state profile. This steadystate profile must match the equilibrium profile given by a balance between capillarity and gravity [START_REF] Cantat | Foams : Structure and Dynamics[END_REF][START_REF] Weaire | The Physics of Foams[END_REF] . At late time and for z << h, r does not depend anymore on t and:

r(z, t) = r(z) ∼ R 1 + αz (1)
where α = ρgR γaw+γow ∼ ρgR γaw , has the dimension of an inversed length. This classic profile illustrates the complex multi-scale geometry of foam under gravity: for z < 1/α, the curvature is only set by the bubble size r ∼ R, while for z > 1/α, r ∼ R/αz. To describe the imbibition dynamics within the aqueous foam, we consider the flow throughout the slender Plateau borders below the front [START_REF] Cantat | Foams : Structure and Dynamics[END_REF][START_REF] Weaire | The Physics of Foams[END_REF] and use the Hagen-Poiseuille law:

v = - K c δr(z) 2 η dP dz (2)
where r(z) is given by Eq. 1 and v, η, K c represent the average velocity of the liquid, the liquid dynamic viscosity and the permeability coefficient for infinite slender channels, K c = 1/50. δ = 0.16 is a numerical factor linking r and A, the cross-section of a Plateau border so that A = δr 2 [27][28][29]. Then, we write mass conservation, which yields v(z) = v(h)r 2 (h)/r 2 (z) = ḣr 2 (h)/r 2 (z). Integrating Eq. 2 between 0 < z < h with the boundary conditions p(z = 0) = p o and p(h) = p o -γ aw /r F and using Eq. 1 gives:

ḣ (1 + αh) 3 - 1 (1 + αh) 2 = αδγ aw R 2 10ηr F (3)
Eq. 3 is then integrated with the condition h(t = 0) = 0:

(1 + αh) 4 + 4 1 + αh -5 = 2δ 5 √ 3
α 2 Rγ aw √ η t (4)
At late times or for αh >> 1, Eq. 4 predicts an evolution of the rising front in t 1/4 , while for αh << 1, a second-order limited expansion of Eq. 4 gives h 2 ∼ Rγ aw t/ √ η, thus suggesting an early time dynamics of imbibition in t 1/2 . These regimes intersect for h ∼ 1/α, the typical height for which the curvature radius of the foam changes from R to R/αh. Moreover, Eq. 4 suggests a collapse of the data using the dimensionless variables αh and t * = α 2 Rγaw √ η t, which is in very good agreement with our experimental results obtained for gives [13]:

h ∼ γ 2 aw ηρg t 1/3 (5)
The scatter in the data makes it difficult to distinguish between the different temporal regimes of equations 4 and 5. Yet, contrary to equation 4, equation 5 neither captures the influence of nor R that is observed in the experiments.

This can be understood recalling that Equation 5 has been derived i) solving the Partial Differential Equation (PDE) ruling the evolution of r ( z, t) ii) assuming a foam with an initial liquid fraction equal to zero and iii) seeking for self-similar solutions of the PDE.

This procedure is particularly useful to describe the asymptotic behavior of solutions in the limit where these no longer depend on the detail of the initial and boundary conditions.

However, in our experiment, the foam sample is only 5-cm high, not infinitely dry and the time-scale is limited. Thus, the initial and boundary conditions can not be neglected. We therefore never observe any convergence of the dynamics towards the t 1/3 self-similar regime.

We highlight the ability of dry liquid foams to displace and drive upward miscible or immiscible liquids on Earth. This remarkable feature is of considerable interest in various industries such as soil remediation, detergency, shampoo industry, enhanced oil recovery and dismantling of nuclear power plants, where foams are used for their low density, non-newtonian rheology and the low quantity of liquid waste they generate. Considering steady-state solutions and rigid interfaces in the framework of porous media, we provide a physical picture for this imbibition dynamics, which captures well our experimental observations concerning the bubble radius R, the initial liquid fraction and the long-term dynamics. Yet, integrating this physical picture in the framework of the foam drainage equation, which has proven successful for both mobile and rigid interfaces remains to be done.

  Organic phase γ aw γ ao γ ow η ρ E S Sunflower oil 23.7 32.5 4.5 55 0.92 -4

FIG. 1

 1 FIG. 1. a) Experimental set-up. b) Geometrical elements of foams: Nodes, Plateau borders and foam films. c) Slender Plateau border with a radius of curvature r F . d) Capillary rise of oil in a dry aqueous foam. The olive oil drop appears white thanks to a fluorescent dye and invades the network of Plateau borders. Time interval between each frame is 67 s. The scale bar is 3 mm.

FIG. 2 .

 2 FIG. 2. Height of the rising front for olive oil as a function of time for different initial liquid fractions. Inset: log-log scale of the data corresponding to = 0.7.10 -4 . The plain line corresponds to the t 1/2

FIG. 3 .

 3 FIG. 3. Typical evolution of the fluorescence intensity of a single node during sunflower oil imbibition as a function of the node's vertical dimension. Initially, the node is located at x = 1 mm and z = 7.9 mm . The intensity of fluorescence is in arbitrary units and the different curves correspond to fluorescence intensity measurements taken at different times. The black arrow indicates the evolution of time. Inset: Width of the peak of intensity (defined at half of the maximum intensity) as a function of time. The red arrow highlights the time at which the oil front has reached the node.

FIG. 4 .

 4 FIG. 4. Normalized imbibition front αh as a function of normalized time t * in log-log plot. In the legend, AS1, AS10 and AS100 respectively stand for Aqueous Solution with η =1.4,10 or 100 mPa.s, while OO and SO stands for Olive Oil and Sunflower Oil. The black line corresponds to (1 + αh) 4 + 4 1+αh -5 = b 2δ 5 √ 3 t *

√ 3 t * , with a coefficient b deduced from the experiments of b ∼ 0.1.Foam drainage theory provides different scaling laws for foam imbibition, which can be confronted to our results. In particular, the diffusive regime observed under microgravity conditions[10,11] is in very good agreement with our data as can be seen in the inset of figure2, thus underlining the negligible role of gravity when the foam is wet at the bottom. At late times, analytical and numerical calculations for the imbibition of infinitely dry foams on Earth

We thank O. Pitois, M. Adler, A. Delbos, F. Rouyer, A.L. Biance, I. Cantat, Y. Peysson, F. Chevoir and E. Reyssat for fruitful discussions. Samples of surfactants were kindly provided by M. Pepin from Stepan Europe and R. Roth from Evonik. We gratefully acknowledge financial support from Agence Nationale de la Recherche (ANR-11-JS09-012-WOLF).