Marie Maurel
email: marie.maurel@veolia.com

Damien Chenu
email: damien.chenu@veolia.com

Pierre Mandel
email: pierre.mandel@veolia.com

Olivier Piller
email: olivier.piller@irstea.fr

Demand of the End Users of the Online Security Management Toolkit

System architecture -Data flows -Data format -End users Objectives Specifications of the data flows to be considered in order to reach the objectives of the online security management toolkit. In particular, the type (frequency, volume, format etc.) of flows should be determined, based on the already existing system (hardware and software) architecture of the end users.

Summary

The SMaRT-Online WDN project aims at developing an online security management toolkit for drinking water distribution systems (DWDS) based on the combination monitoring/modelling. Many modelling software tools for DWDS have been available for more than twenty years. However, very few of them are able to process all the real-time data, that are increasingly being generated by sensors. Handling these data and making sensors smart raises the issues of: (i) transferring/storing the output data they generate; (ii) processing continuously these data in order to provide operators with reliable and valuable information. The SMaRT-Online WDN project shall deal with them.

The objective of this deliverable is to specify the data flows to be considered in order to reach the objectives of the online security management toolkit. In particular, the type (frequency, volume, format etc.) of flows should be determined, based on the already existing system (hardware and software) architecture of the end users.

Here are the main conclusions:

1) A general system architecture is proposed with the following modules (hardware or software tools): Sensors, Supervisory Control And Data Acquisition (SCADA) system, Data treatment, Model (hydraulics and quality), Online calibration, Database, Geographic Information System (GIS), Event detection algorithm, Statistical model, Optimisation algorithm and Source Identification; 2) The new modules to be developed are: Data treatment, Model (hydraulics and quality), online calibration, Database, Event detection algorithm, Statistical model, Optimisation algorithm, Source Identification; 3) Aside from the development of the new modules, the following modules, particular attention should be paid to the input/output data flows of the following elements: GIS and model (hydraulics and quality).

4

List of Figures: .

List of Tables:

Table 1: Main characteristics of the existing elements of the online security management toolkit for the end users . High-frequency workflows ................................................................................................................................................. Medium-frequency workflows ........................................................................................................................................... Low-frequency workflows ................................................................................................................................................ Exceptional workflows ......................................................................................................................................................

Proposal for

Background

While multiple modelling software tools for drinking water distribution systems (DWDS) have been available for more than twenty years, very few of them are able to work with real-time data. However, recent developments in water quality monitoring are quickly popularising the use of sensors, which generate large amounts of real-time data. Handling these data and making sensors smart raises the issues of: (i) transferring/storing the output data they generate; (ii) processing continuously these data in order to provide operators with reliable and valuable information.

Modelling software tools for DWDS have first been developed to reproduce average hydraulic phenomena (pressure and flow). This coincided with the advent of the first computers affordable by research institutions in the early 1960s [START_REF] Ormsbee | The History of Water Distribution Network Analysis: The Computer Age[END_REF]. In the beginning of the 1980s, researchers started modelling water quality in DWDS and proposed several approaches, generally combined with the already existing hydraulic models [START_REF] Graymann | An overview of water quality modeling in a distribution system[END_REF].

Recently, there has been a renewed interest in quality modelling since DWDS may be potential targets for terrorists. A review of recent water quality modelling methods can be found in [START_REF] Ostfeld | A review of modeling water quality in distribution systems[END_REF]. Nowadays, while hydraulic simulation software is in widespread use within the water industry for strategic supply analysis, design of control strategies, network extensions, and maintenance planning [START_REF] Machell | Online modelling of water distribution systems: a UK case study[END_REF], water quality simulation software remains much less used.

Sensors for water quality monitoring in DWDSsmall in size and able to communicate in short distancehave been developed recently, benefiting from the advances in wireless communications and electronics [START_REF] Akyildiz | Wireless mesh networks: a survey[END_REF]. These sensors differ according to the type of parameter they measure: On one hand, low-cost, low-power and multi-functional measure physicochemical parameter such as conductivity, pH, redox potential, diffusion, chlorine concentration, turbidity, UV absorbance (Storey et al., 2010). On the other hand, broadband biosensor systems measure the toxicity of water by monitoring the agility of microorganisms highly sensitive to the presence of toxic ingredients in water [START_REF] Tecon | Bacterial biosensors for measuring availability of environmental pollutants[END_REF]. Due to their high cost, biosensor systems will not be considered within the SMaRT-Online WDN project.

Data mining in sensor networks has lately received considerable attention because of the increasing amounts of data stored in databases, data warehouses or other repository information [START_REF] Han | Data Mining: Concepts and Techniques[END_REF]. Data mining may involve integration of techniques from multidiscipline such as database technology, statistics, machine learning, neural networks, information retrieval etc. [START_REF] Hand | Principles of Data Mining[END_REF]. Aside from these techniques, model and inference considerations, visualisation and online updating (for models, GIS etc.) are topics of particular relevance to the field of DWDS modelling.

Aims of the Online Security Management Toolkit

According to [START_REF] Piller | SMaRT-Online WDN : Online Security Management and Reliability Toolkit for Water Distribution Networks[END_REF] and for online management, the main innovative tasks of the SMaRT-Online WDN Project are listed below:

1) Detection of abnormal events with a binary classifier of high accuracy. How to differentiate from non reported but normal change in the system? 2) Generation of real-time, reliable (i) flow and pressure values, (ii) water quality parameter values of the whole water network; which are physical phenomena to consider?

3) Consideration of water quality measurement for online calibration of the hydraulic model; 4) Semi-automatic aggregation of complex networks; 5) Semi-automatic update and online-calibration of the model; Another important task for the project but to consider in the offline concept is using risk analysis methods for risk assessment and impact evaluation. It will no longer discuss in this report.

System Architecture and Workflow of the Online Security

Management Toolkit

General View

In order to detail the data flows that have to be considered to reach the aims set in Part 2, the general system architecture has been designed and is schematically presented in Figure 1. Various components, existing or not, hardware or software, have been assembled and the main data flows passing from one to another have been represented. This figure does not represent the data entered by the end-users to control the online security management tool (via a GUI for example). This will be presented in Part 4. The components taken into account in the system architecture are listed below. The components to be developed are written in bold characters: 1) Sensors: considered here as completely independent, they are only sending regular data at arelatively to all other data flows of the considered architecturehigh frequency; 2) Supervisory Control And Data Acquisition (SCADA) system: receives and dispatches data possibly to other sites. The effective control of the network (actions on valve, pump and other devices) is not represented as it is beyond the scope of the SMaRT-Online WDN project; 3) Data treatment: Software function, which (i) "cleans up" the raw data retrieved by the SCADA system and (ii) delivers the appropriate (select, recalculate etc.) data for subsequent tasks; 4) Online Model (hydraulics and quality): Software tool, which reproduces the behaviour of the DWDS. As its position in Figure 1 suggests, this component plays a central role and its accuracy will be crucial; slow transient model will be considered for hydraulic prediction in place of extended period simulation model. For the transport of usual physico-chemical products the dispersion phenomenon will be added and non-perfect mixing laws at nodes will apply for T and cross junctions. This should be regularly (e.g. 5 min) updated from some data measurements and the demand should be calibrated online. 5) Database: Structured collection of data that are stored at different frequency rates according to their type (treated measurements from the sensors, events, model settings etc.). As the database may not grow indefinitely, a regular purge should be scheduled.

The corresponding data flow has not been represented on; 6) Geographic Information System (GIS): Software tool already existing mainly used in the proposed architecture as a display tool, possibly offering GUI possibilities; 7) Event detection algorithm: Software function, which determines if an event has occurred in the DWDS. Basically, the algorithm compares treated data series

7

(representing the current state and previous states) issued from the sensors to the results given by a statistical model. Previous data shall also be extracted from the database (previous states on a larger timescale, previous events etc.). 8) Source Identification: this module will use velocity recorded in the database and a reverse transport model to determine the mix and the spread of a detected event (contamination of water quality for example) and its evolution in the future. 9) Statistical metamodel (for water quantity and quality at sensor network): Software function, which reproduces specific behaviours of the DWDS at sensor locations. It shall be designed to be less generic than the model presented before, targeting only the some parameters that may be affected by one of the event to be detected by the management toolkit. It shall be continuously "fed" with fresh data extracted from the database. For example, an Artificial Neural Network (ANN) may be used; 10) Optimisation algorithm: Software function, which finds optimal model settings in the sense of minimising the gap between observations (data issued from sensors) and predictions (data issued from the model) to calibrate the model.

Detailed Workflows

We review in the following section the main characteristics of each workflow (and give some methods that could be used) and detail the data flows (content, type, size, frequency). Also, the data inputs from the operators via GUI will be mentioned.

High-frequency workflows

Figure 2 represents the workflows that have to be routinely applied to every data batch. The frequency shall depend upon the sensors settings and eventually on the frequency of acquisition.

Figure 2: System architecture for high-frequency workflows

Data Treatment

The software tool to be developed should fulfil three tasks: 1) Basic data treatment: a. eliminate gross errors (such as negative concentrations etc.) b. filter (removing high frequency inputs, which may alter the signal) c. interpolate (recreating missing data) 2) Data reconciliation: some data should be used as inputs to the model in order to generate up-to-date hydraulic and water quality predictions every 5 minutes (for the current state and for the coming hours for example). If some data are inconsistent, for example water levels in the tanks and pressures in the DWDS, they should be corrected in some way not to lead to a crash of the model; 3) Possible intermediate calculation if the model needs inputs that differ from treated data; this may not be the case for display purposes in the GIS.

The treated data should be directly readable by the model and the GIS, so the data format must be specified.

Online Model: Hydraulic and quality Hydraulic and water quality models should be fed by treated data at high frequency (for instance every 5 minutes). Water demands in the hydraulic model will be obtained directly by online 9 calibration to match measurements in the system. Some measurements will substitute the predictions while other will be used for demand calibration. The overall defines the online model.

The online hydraulic and quality model may be run in several modes. The first one, or reference mode, represents the current state of the systems and include a look back in the past. A second one, or look ahead mode, will give point and confidence interval predictions for the next hours. It will use the demand forecast (from the statistical metamodel). A last one, or what if mode, will permit testing different control actions and their impacts onto the system.

Statistical metamodel: water quantity and quality

This software function will use some machine learning method to model an updated "normal" behaviour of the DWDS. In particular, the outputs of this model should be particularly impacted by the events sought to be detected (for example, conductivity if saline intrusion should be detected etc.). Data coming from the database should regularly "train" (calibrate) the metamodel; other data may be used for validation purposes. The frequency rate, at which data is brought to the model, may vary, according to the type of model selected, but it should be relatively high if both calibration and validation are desired. The model should continuously be based on a "moving" amount of data (i.e., last current month). This amount should be: (i) relatively large in order not to be influenced by a recent event not being detected yet; (ii) and relatively small to be representative of the current state of the DWDS. A similar software tool will be developed for water demand prognosis purpose.

Medium-frequency workflows

Figure 3 represents the workflows that have to be routinely applied to data batch at a "medium" frequency (every hour, for example). The frequency shall depend upon the frequency of acquisition and the end-users.

10

Event detection algorithm

Based on the comparison of treated data with the predictions from the statistical metamodel, the detection algorithm should decide whether an event has occurred on the DWDS. The amount of treated data transmitted should represent a longer period than the acquisition period in order to be representative of an event. Many algorithms in the field of error detection might be used for that purpose. The statistical metamodel should also deliver results with uncertainty levels for the event detection algorithm to distinguish between natural variations in the data and variations related to an event.

Source Identification Once an event has been detected, a specific module should be launched to evaluate the potential contaminant sources by backtracking in reverse time. Thanks to online modelling and storing of the velocities every 5 min source identification will be more robust and less sensitive to demand uncertainty.

The future impact of the event shall also be simulated by using the statistical metamodel for the demand forecast and the transport model in look ahead mode. All the output data regarding the event, its possible past (source) and future (impact) should be exportable to the GIS (through data tables) and finally stored into the database.

11

Low-frequency workflows Figure 4 represents the workflows that have to be routinely applied to data batch at a "low" frequency (once a month, for example). The frequency shall depend upon the end-users.

Figure 4: System architecture for low-frequency workflows

In case of no event on treated data found by the event detection algorithm, these data shall be used to recalibrate the models (hydraulics and or quality) and to feed the optimisation algorithm.

Optimisation algorithm

Based on observed data considered as "safe" (not having been considered as representative of an event), and with a determined frequency (once a month, for example), the model will be automatically recalibrated. Communication function shall be developed in order to: (i) modify the input files of the model according to the variable values determined by the optimisation algorithm; (ii) get the model output values to calculate the objectives and constraints. Despite its routine character, the optimisation procedure might require some user inputs, such as:  The definition of the optimisation problem (and especially the objective or particular constraints);  The selection of a configuration found optimal by the algorithm (in case of a multi-objective problem, one of the points belonging to the Pareto frontier for example).

Exceptional workflows

Exceptional workflows are manual updates. As such, they will not be detailed in this section. Automatic back-ups to the database should however be scheduled with every manual update.

12 4. Proposal for Software Implementation and Interfaces

Characteristics of existing architecture modules

The answers given by three end users (CUS, BWB and VEDIF) to a questionnaire on the existing elements of the online security management toolkit are summarised in Table 1. Some existing characteristics are shared by all the end users, they have been specified in green cells in Table 1; the characteristics that differ from each other have been specified in red cells.

Characteristics of the data flows of the online security management toolkit

Based on the proposed system architecture presented in Part 3 and on the answers to the questionnaire presented in the previous section, a summary of the characteristics of the data flows within the system architecture is proposed in Table 2. If the characteristic is already that of the existing structure, the cell background is green. Particular attention should be paid to the items written in cells with yellow backgrounds (those related to the GIS and the model).

14 1) A general system architecture is proposed with the following modules (hardware or software tools): Sensors, Supervisory Control And Data Acquisition (SCADA) system, Data treatment, Model (hydraulics and quality), Database, Geographic Information System (GIS), Event detection algorithm, Statistical model (transport and quality), Optimisation algorithm; 2) The new modules to be developed are: Data treatment, Online transport model (it include the online demand calibration), Database, Event detection algorithm, Source identification, Statistical model, Optimisation algorithms for offline calibration; 3) Aside from the development of the new modules, the following modules, particular attention should be paid to the input/output data flows of the following elements: GIS and model (hydraulics and quality).

Further developments include: All data flows related to user inputs should be specified; in particular their interconnection with the data flows presented herein should be specified.

Figure 1 :

 1 Figure 1: General system architecture and workflows .

  Software Implementation and Interfaces .................................................. 4.1. Characteristics of existing architecture modules ...................................................................................... 4.2. Characteristics of the data flows of the online security management toolkit ..............................

Figure 1 :

 1 Figure 1: General system architecture and workflows

Figure 3 :

 3 Figure 3: System architecture for medium-frequency workflows

Table 2 :

 2

Main characteristics of the data flows .....................................................................................................................................

Contents

: 1. Background ................................................................................................................... 2. Aims of the Online Security Management Toolkit ........................................................... 3. System Architecture and Workflow of the Online Security Management Toolkit .............

4.1. General View................................................................................................................................................................. 4.2. Detailed Workflows ...................................................................................................................................................

Table 1 : Main characteristics of the existing elements of the online security management toolkit for the end users End user

 1

	CUS	BWB	Vedif

VERI TZW ENGEES Irstea

13

Conclusion

Here are the main conclusions:


Nicolas Augier

Madalena Chaves

Jean-Luc Gouzé

Control for synchronization of bistable piecewise affine genetic regulatory networks

Keywords: Control of networks, Robust control, Synchronization, Gene networks

In this paper, we study the diffusive coupling of a network of identical bistable switch piecewise-affine systems. We propose a control strategy which synchronizes every bistable switch sub-system of the network towards the same steady state.

INTRODUCTION

Piecewise-affine systems are now a frequently used framework for modeling genetic regulatory networks (see, for instance, de [START_REF] De | Modeling and simulation of genetic regulatory systems: a literature review[END_REF]; [START_REF] Farcot | A mathematical framework for the control of piecewise-affine models of gene networks[END_REF]). This type of systems has been used to model several distinct dynamical behaviors, such as, multistability (see, for instance, [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF]) or oscillatory behaviors (see, for instance, [START_REF] Farcot | Periodic solutions of piecewise affine gene network models with non uniform decay rates: The case of a negative feedback loop[END_REF]). Diffusive coupling of identical sub-systems, where each pair of subsystems are coupled by diffusion, and related synchronisation issues have attracted a large interest for years (see [START_REF] Aminzare | Synchronization of diffusively-connected nonlinear systems: results based on contractions with respect to general norms[END_REF]; [START_REF] Scardovi | Synchronization of interconnected systems with applications to biochemical networks: an input-output approach[END_REF]). In [START_REF] Chaves | Coupling and synchronization of piecewise linear genetic regulatory systems[END_REF], the authors have introduced coupling of piecewise affine systems by diffusion and have studied synchronization issues. The systems are coupled by discrete diffusion, that is, the dynamics is the sum of a reaction term corresponding to the individual bistable dynamics and a diffusion term involving symmetric Laplacian matrices, which define a coupling topology, so that the dynamics of the network can be seen as a space discretization of a system of two coupled reaction-diffusion pde's (see, for instance [START_REF] Pouchol | Phase portrait control for 1D monostable and bistable reactiondiffusion equations[END_REF]; [START_REF] Baccoli | Boundary control of coupled reaction-diffusion processes with constant parameters[END_REF] for results about the control of reaction-diffusion equations, [START_REF] Aminzare | Synchronization of diffusively-connected nonlinear systems: results based on contractions with respect to general norms[END_REF] for synchronisation issues in reaction-diffusion pde's). It has been shown that in the case of the coupling of identical piecewise affine systems, the coupled system may have a much more complex dynamics than the individual one. For instance, the coupling scheme introduces new steady states, that we call in this paper diffusive steady states, which may be locally stable (see [START_REF] Chaves | Coupling and synchronization of piecewise linear genetic regulatory systems[END_REF]). For biological purposes, these diffusive steady states may be undesirable, so it becomes interesting to find a control strategy that prevents convergence to those states. Controlling coupled piecewise affine systems is a challenging task, and to the best of the authors knowledge, very few results have been developed in this direction. We propose to study a network of two dimensional bistable systems, and consider the problem of stabilizing this network towards its steady states. To take biological constraints into account, we consider qualitative measurements, and propose to use a control which depends only on qualitative knowledge of the state variable, as in [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF]. In this paper, we propose a control method which works for every type of coupling topology, and which is robust w.r.t. the coupling strength. In particular, we prove a result of stabilization of parametric ensembles of systems at steady states (see, for instance, [START_REF] Ryan | On simultaneous stabilization by feedback of finitely many oscillators[END_REF], or [START_REF] Chittaro | Asymptotic ensemble stabilizability of the Bloch equation[END_REF], [START_REF] Augier | Adiabatic ensemble control of a continuum of quantum systems[END_REF] for results occuring in different frameworks).

SINGLE SYSTEM DYNAMICS AND CONTROL

In this section, we recall some results from [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF] concerning the control of a single bistable switch system. The two variables represent two proteins which mutually inhibit each other.

Individual dynamics

Consider first the uncontrolled equation ẋ1

= -γ 1 x 1 + k 1 s -(x 2 , θ 2 ) ẋ2 = -γ 2 x 2 + k 2 s -(x 1 , θ 1 ), (1)
where for j ∈ {1, 2}, x j ∈ R, and for

θ ∈ R, s -(•, θ) : R → R is such that s -(x, θ) = 0 if x > θ, and s -(x, θ) = 1 if x < θ. It is assumed that s -(x) ∈ [0, 1] for x = θ.
The positive constants (γ j ) j∈{1,2} , (k j ) j∈{1,2} correspond, respectively, to the degradation and the production rates of each variable. It is a classical fact (see [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF]) that the domain K = [0, k1 γ1 ] × [0, k2 γ2 ] is forward invariant by the dynamics of Equation (1). From now on we consider only solutions evolving in K.

Define the regular domains

B 00 = {x = (x 1 , x 2 ) ∈ R 2 | 0 < x 1 < θ 1 , 0 < x 2 < θ 2 }, B 01 = {x = (x 1 , x 2 ) ∈ R 2 | 0 < x 1 < θ 1 , θ 2 < x 2 < k 2 γ 2 }, B 10 = {x = (x 1 , x 2 ) ∈ R 2 | θ 1 < x 1 < k 1 γ 1 , 0 < x 2 < θ 2 }, B 11 = {x = (x 1 , x 2 ) ∈ R 2 | θ 1 < x 1 < k 1 γ 1 , θ 2 < x 2 < k 2 γ 2 }.
For a piecewise affine system defined on R n whose dynamics restricted to a regular domain B is a linear dynamical system having an asymptotically stable equilibrium, define the focal point on B as this equilibrium point. Notice that this point may belong to B or not. Each region B ij has a focal point, φ ij = (x i , xj ) corresponding to xi = ki γi s -(x j , θ j ). and Equation ( 1) has two locally asymptotically stable steady states, φ 10 = ( k1 γ1 , 0) ∈ B 10 and φ 01 = (0, k2 γ2 ) ∈ B 01 , and an unstable Filippov equilibrium point at (θ 1 , θ 2 ). In addition, there exists a curve, called separatrix, passing through (θ 1 , θ 2 ) and dividing K in two regions (above and below) such that the solutions of Equation ( 1) reach B 01 or B 10 , respectively, in finite time. Moreover, B 10 (respectively, B 10 ) is included in the bassin of attraction of φ 10 (respectively, φ 01 ). We make the following assumptions on the parameters of the system (for more details, see [START_REF] Chaves | Exact control of genetic networks in a qualitative framework: the bistable switch example[END_REF]):

θ j < k j γ j , j ∈ {1, 2}; θ 2 θ 1 > k 2 k 1 γ 1 γ 2 ; θ 2 θ 1 < k 2 k 1 . (H)
2.2 Controlling the individual system to φ 01 (respectively, φ 10 )

Consider now the controlled equation:

ẋ1 = -γ 1 x 1 + uk 1 s -(x 2 , θ 2 ) ẋ2 = -γ 2 x 2 + uk 2 s -(x 1 , θ 1 ) (2)
The control u ≡ u(t, x(t)) is assumed to act on the production rates of each variable. It is assumed to depend only on the domain (B jk ) j,k∈{0,1} to which the solution x(t) of Equation ( 9) at time t belongs, and u has values in a finite set of the form {u min , 1, u max }, where u max ≥ 1 and u min ≥ 0. Note that u changes the location of the focal points φ ij .

Here we present a control strategy, proposed in (Chaves and Gouzé, 2011, Section 5). The separatrix (S u ), which separates B 00 into two regions (S u ) -and (S u ) + is defined, in the coordinates (x 1 , x 2 ) ∈ B 00 , for u ≥ 0, as the curve of equation

x 2 = α(x 1 , u) = k 2 u γ 2 - k 2 u γ 2 -θ 2 k1u γ1 -x 1 k1u γ1 -θ 1 γ 2 γ 1
. ∈ (S umax ) -).

One

Control algorithm:

• First phase: Choose u ≡ u 01 min (respectively, u 10 min ) during a time T > 0 large enough.

• Second phase: Choose u ≡ u max until x enters in B 01

(respectively, B 10 ). • Third phase: Choose u ≡ 1.

During the first phase, every focal point of the system belongs to B 00 , hence the solution x(t) of Equation ( 9) converges towards the point (u min k1 γ1 , u min k2 γ2 ) ∈ B 00 when t → ∞. During the second phase, x(t) reaches B 01 or B 10 in finite time, depending on the choice of u ≡ u 01 min or u ≡ u 10 min . During the third phase, x(t) converges towards φ 01 or φ 10 .

COUPLED SYSTEM

In this section, we consider a network of N ∈ N identical systems whose individual dynamics is given by Equation (1), which are coupled by diffusion, as it has been studied in a slightly different setting in [START_REF] Chaves | Coupling and synchronization of piecewise linear genetic regulatory systems[END_REF]. Then we propose an adaptation of the control algorithm of the individual dynamics presented in Section 2.2 to synchronize the coupled system at the same individual steady state.

Dynamics of the coupled system

Let N ∈ N. Define, for k ∈ {1, 2}, x k = (x k,j ) j∈{1,...,N } ∈ R N , and the N -dimensional vector

q(x k ) = t s -(x k,1 , θ k ), . . . , s -(x k,N , θ k ) . Define for every x = (x k,j ) (k,j)∈{1,2}×{1,...,N } ∈ R 2N the canonical projection π k (x) = x k for k ∈ {1, 2}.
Consider the following equation ẋ1

= -(Γ 1 + L 1 )x 1 + k 1 q(x 2 ) ẋ2 = -(Γ 2 + L 2 )x 2 + k 2 q(x 1 ). (3)
where, (γ j ) j∈{1,2} , (k j ) j∈{1,2} satisfy Assumption (H), Γ j is the N -dimensional diagonal matrix with coefficients γ j .

The matrix L j is a Laplacian N -dimensional symmetric matrix, that is, its coefficients (l j kl ) k,l∈{1,...,N } satisfy

l j kl =      N i=1 a ki , l = k -a kl , l = k,
where a kl ≥ 0 for every k, l ∈ {1, . . . , N }.

From now on, assume that for j ∈ {1, 2}, L j defines a strongly connected graph. Then it follows that Γ j + L j is invertible, -(Γ j + L j ) is Hurwitz and (Γ j + L j ) -1 is positive.

Define the regular domains of Equation (3) as the cartesian products of the domains (B jk ) j,k∈{0,1} defined in Section 2.1. More precisely, for a sequence (j l , k l ) l∈{1,...,N } such that for every l ∈ {1, . . . , N }, (j l , k l ) ∈ {0, 1} 2 , we say that

x ∈ B j1,k1 × • • • × B j N ,k N if for every l ∈ {1, . . . , N }, (x 1,l , x 2,l ) ∈ B j l ,k l .
One can show that the full domain K N , where K is defined in Section 2.1, is forward invariant by the dynamics of Equation (3).

For a regular domain B of K N the focal point corresponding to the domain B is given by

x1 = k 1 (Γ 1 + L 1 ) -1 q(π 2 (B)) x2 = k 2 (Γ 2 + L 2 ) -1 q(π 1 (B)), (4)
where, by abuse of notations, q(π j (B)) is the constant value taken by q in the set π j (B), for every j ∈ {1, 2}. It

follows easily that (x 1 , x2 ) is a steady state when (x 1 , x2 ) ∈ B, that is x1 = k 1 (Γ 1 + L 1 ) -1 q(x 2 ) x2 = k 2 (Γ 2 + L 2 ) -1 q(x 1 ). (5) Define x = Φ 10 ∈ R 2N such that for every k ∈ {1, . . . , N }, (x 1,k , x 2,k ) = φ 10 . Define Φ 01 ∈ R 2N similarly.
Because of the Laplacian structure of (L j ) j∈{1,2} , Φ 10 and Φ 01 , that we call synchronized steady-states are locally stable steady states of Equation ( 3), which are independent from the coupling L 1 and L 2 .

However, as already noticed in [START_REF] Chaves | Coupling and synchronization of piecewise linear genetic regulatory systems[END_REF] for a sligthly different class of systems, the dynamics of Equation ( 3) may converge to other not synchronized steady states, given by Equation ( 5), that we call diffusive steady states.

These diffusive states all belong to K N . In the following we assume that ||L 1 || and ||L 2 || are small enough to guarantee that the diffusive focal points are in the same domains as their respective non-diffusive focal points, that is, for every regular domain B of K N , the focal point given by ( 4) with (L 1 , L 2 ) lies in the same domain as the focal point given by ( 4) with L 1 = L 2 = 0.

Control strategy for the coupled system

This section, where proofs are not given, constitues the main contribution of this paper. We are interested in controlling the dynamics of the coupled system ẋ1 = -(

Γ 1 + L 1 )x 1 + uk 1 q(x 2 ) ẋ2 = -(Γ 2 + L 2 )x 2 + uk 2 q(x 1 ) (6)
towards Φ 01 or Φ 10 . The control u is assumed to act on every sub-system of the network and satisfies the same hypothesis as in Section 2.1. For this, we show that the control algorithm exposed in Section 2.2 can be adapted to this case. Let u 01 min , u 10 min , u max be chosen as in Section 2.2. For a regular domain B of K N the focal point of Equation ( 6) corresponding to the domain B is given by x1

= uk 1 (Γ 1 + L 1 ) -1 q(π 2 (B)) x2 = uk 2 (Γ 2 + L 2 ) -1 q(π 1 (B)). (7)
We have the following lemmas. Lemma 1. Every focal point of Equation ( 6) with u ≡ u 01 min (respectively, u 01 min ) belongs to

B 00 × • • • × B 00 . Lemma 2. Let u ≥ 0. Assume that x 0 ∈ (S u ) + × • • • × (S u ) + (respectively, x 0 ∈ (S u ) -× • • • × (S u ) -). Then the solution x(t) of Equation (6) such that x(0) = x 0 reaches B 01 ו • •×B 01 (respectively, B 10 ו • •×B 10 ) in finite time.

Control algorithm:

First step: synchronization Choose u ≡ u 01 min (respectively, u 01 min ). By Lemma 1, every system of the network reaches the domain B 00 in finite time. We obtain the following. Theorem 3. The solution x of Equation ( 6

Fix

> 0 small enough. Lemma 3 allows to consider T > 0 such that for every initial condition x 0 ∈ K N , the solution of Equation ( 6) such that x(0) = x 0 satisfies ||x(T ) -Φ * || < , that is, every subsystem of the network can be driven to an arbitrary small neighborhood of (u 01 In particular, one can choose > 0 small enough to guarantee x(T ) ∈ (S umax )

+ × • • • × (S umax ) + (respectively, x(T ) ∈ (S umax ) -× • • • × (S umax ) -).
Second step: Choose u ≡ u max . Using Lemma 2, one can prove the following result. Theorem 4. Assume that u is defined as in the two steps above. Then the solution x(t) of Equation ( 6) reaches the regular domain

B 01 ו • •×B 01 (respectively, B 10 ו • •×B 10 ) in finite time.
Third step: Choose u ≡ 1. Theorem 5. Assume that u is defined as in the three steps above. Then the solution x(t) of Equation ( 6) converges to Φ 01 (respectively, Φ 10 ) when t → ∞.

Robustness

Fix a Laplacian matrix L j of a strongly connected graph. Consider the system (3), replacing L j by Lj , where Lj = α j L j , with α j ∈ R is a fixed but unknown parameter in a certain given interval, that is ẋ1

= -(Γ 1 + α 1 L 1 )x 1 + uk 1 q(x 2 ) ẋ2 = -(Γ 2 + α 2 L 2 )x 2 + uk 2 q(x 1 ). (8)
Up to increasing the time T of phase 1 of the control algorithm, by continuity of the flow of Equation ( 8) w.r.t. α = (α 1 , α 2 ), we obtain a control strategy which is uniform w.r.t. α in every compact subset I of [0, d 1 )×[0, d 2 ), where d 1 , d 2 > 0 depend on (L 1 , L 2 ).

We deduce that Equation ( 8) is ensemble stabilizable at Φ 10 and Φ 01 , uniformly w.r.t. α = (α 1 , α 2 ) ∈ I, that is, there exists a control u : [0, +∞) → {u 10 min , 1, u max } (respectively, u : [0, +∞) → {u 01 min , 1, u max }) independent of α, such that the parametric solution x α of Equation ( 8) where, for every α ∈ I, x α (0) = x 0 ∈ K N , with x 0 independent of α, is such that x α (t) converges to Φ 10 (respectively, Φ 01 ) when t → ∞, uniformly w.r.t. α ∈ I.

EXAMPLE: CASE

N = 2 Consider the following piecewise-affine differential system in R 4 , ẋ1,1 = -γ 1 x 1,1 + k 1 us -(x 2,1 , θ 2 ) + a 1 (x 1,2 -x 1,1 ) ẋ2,1 = -γ 2 x 2,1 + k 2 us -(x 1,1 , θ 1 ) + a 2 (x 1,1 -x 1,2 ) ẋ1,2 = -γ 1 x 1,2 + k 1 us -(x 2,2 , θ 2 ) + a 1 (x 1,1 -x 1,2 ) ẋ2,2 = -γ 2 x 2,2 + k 2 us -(x 1,2 , θ 1 ) + a 2 (x 1,2 -x 1,1 ), (9)
where, (γ j ) j∈{1,2} , (k j ) j∈{1,2} are positive constants satisfying Assumption (H), and a 1 ≥ 0 (respectively, a 2 ≥ 0) are the diffusion couplings on the first variables (respectively, the second variables). We have made simulations of such a system with (γ 1 , γ 2 ) = (0.2, 0.8), (k 1 , k 2 ) = (0.5, 1.4), (θ 1 , θ 2 ) = (1, 1.2), a 1 = a 2 = 0.05. For this choice of parameters, the value of u max can be chosen as u max = 1. One can choose u 01 min = 0.1, u 10 min = 0.3, and the time T of the second phase of the control strategy is chosen as T = 20. Consider the initial conditions: (x 1,1 (0), x 2,1 (0)) = (0.5, 1.3), and (x 1,2 (0), x 2,2 (0)) = (0.8, 0.3). Figure 1 illustrates the convergence of x(t) towards a diffusive steady state for the uncontrolled system u ≡ 1, the dashed line (S) = (S umax ) represents the separatrix, as defined in Section 2.1. Figure 2 illustrates the convergence of x(t) towards Φ 01 for the controlled system. Finally, Figure 3 illustrates the convergence of x(t) towards Φ 10 for the controlled system. We notice that after the first phase of the control strategy, that is the phase during which the two systems converge towards Φ * , the two systems (red and green curves) follow very close trajectories.

  ) converges when t → ∞ to the synchronized state Φ * = (x 1 , x2 ) defined by x1,k = u 01 min k1 γ1 , and x2,k = u 01 min k2 γ2 (respectively, x1,k = u 10 min k1 γ1 , and x2,k = u 10 min k2 γ2 ) for every k ∈ {1, . . . , N }.

Fig. 1 .

 1 Fig. 1. Convergence towards a diffusive steady state for the uncontrolled system

This work was supported in part by the ANR project Maximic (ANR-17-CE40-0024-01), Labex SIGNALIFE (ANR-11-LABX-0028-01), and the ANR project ICycle (ANR-16-CE33-0016-01)


I T Davies

C Raufaste

Effect of gravity on the orientation and detachment of cubic particles adsorbed at soap film or liquid interfaces

We investigate the interaction that occurs between a light solid cube falling under gravity and a horizontal soap film that is pinned to a circular ring. We observe in both experiments and quasi-static simulations that the final orientation of a cube that becomes entrapped by a soap film is strongly dependent on the Bond number. A cube is rotated by a soap film into one of three main orientations in a process that is driven by energy minimisation. The likelihood of observing each of these final orientations is shown to depend on the Bond number, and the most energetically favourable orientation depends on the terminal height reached by the cube. We also find a critical value for the Bond number, above which a cube is no longer supported by a soap film and detachment occurs, to be less than one.

Introduction

Aqueous foams are complex two-phase materials consisting of gas bubbles, liquid films and channels [START_REF] Cantat | Foams Structure and Dynamics[END_REF]. The interaction between a soap film and a solid object or between particles and an aqueous foam is important in froth flotation for mineral separation, paper deinking, waste water treatment and soil remediation [START_REF] Prud | Foams in mineral flotation and separation processess[END_REF][START_REF] Shean | A review of froth flotation control[END_REF][START_REF] Feng | Effect of contact angle and contact angle hysteresis on the floatability of spheres at the air-water interface[END_REF][START_REF] Xing | Recent experimental advances for understanding bubble-particle attachment in flotation[END_REF][START_REF] Schellenberger | Detaching microparticles from a liquid surface[END_REF]. The efficiency of the flotation process depends on the interaction that occurs between a foam's films and the particles it needs to carry and those it needs to leave behind. Experiments and a simple scaling argument show that the energy required for a particle to pass through a film scales as γr 2 , the typical excess of energy due to the interface deformation, with γ the surface tension of the gas-liquid interfaces and r the radius of the object [START_REF] Goff | Energy absorption in a bamboo foam[END_REF]. More recently, foams have found applications in microfluidics, where capillary forces tend to dominate and gravitational effects become negligible [START_REF] Huerre | Bubbles and foams in microfluidics[END_REF][START_REF] Whitesides | The origins and the future of microfluidics[END_REF]. The interaction that occurs between a foam's films and a solid object or particle also has many potential applications at this scale, for example in pharmaceuticals and medicine, where controlled transportation of particles, objects or solvents through confined geometries is needed [START_REF] Drenckhan | Rheology of ordered foamson the way to discrete microfluidics[END_REF] or where a liquid filter is needed [START_REF] Stogin | Free-standing liquid membranes as unusual particle separators[END_REF].

In another context, capillary forces are key elements in patterning assemblies of nanoparticles at liquid interfaces [START_REF] Di Leonardo | Very-Long-Range Nature of Capillary Interactions in Liquid Films[END_REF][START_REF] Yao | Near field capillary repulsion/ Soft Matter[END_REF][START_REF] Mcgorty | Colloidal self-assembly at an interface[END_REF][START_REF] Furst | Directing colloidal assembly at fluid interfaces[END_REF][START_REF] Law | Self-Assembly of Two-Dimensional Colloidal Clusters by Tuning the Hydrophobicity, Composition, and Packing Geometry[END_REF], which provide materials with peculiar optical and electronic properties. For instance nanocubes self-assemble into graphenelike honeycomb and hexagonal lattices [START_REF] Evers | Low-Dimensional Semiconductor Superlattices Formed by Geometric Control over Nanocrystal Attachment[END_REF][START_REF] Boneschanscher | Long-range orientation and atomic attachment of nanocrystals in 2D honeycomb superlattices[END_REF][START_REF] Beugeling | Topological states in multi-orbital HgTe honeycomb lattices[END_REF][START_REF] Wang | Interplay between spherical confinement and particle shape on the self-assembly of rounded cubes[END_REF] or the pattern can be tuned as the particle shape changes from a sharp cube to a rounded cube [START_REF] Wang | Interplay between spherical confinement and particle shape on the self-assembly of rounded cubes[END_REF]. The capillary deformation induced by adsorbed cubes depends on the Young's contact angle [START_REF] Morris | Modelling the self orientation of particles in a film[END_REF][START_REF] Sharp | Adsorption of shapeanisotropic and porous particles at the air-water and the decane-water interface studied by the gel trapping technique[END_REF]. The orientation of a single cube at the interface drives the interface deformation and the arrangement. On the example of a 90 • contact angle, a single-adsorbed cube prefers a configuration with one corner between three faces pointing upward, which generates a hexapolar deformation to the interface. This hexapolar deformation drives multiple cubes to self-assemble into hexagonal or graphenelike honeycomb lattices [START_REF] Soligno | Self-assembly of cubes into 2d hexagonal and honeycomb lattices by hexapolar capillary interactions[END_REF]. The orientation of the particles at the interface as a function of the packing, confinement and shape has attracted many studies [START_REF] Sharp | Adsorption of shapeanisotropic and porous particles at the air-water and the decane-water interface studied by the gel trapping technique[END_REF][START_REF] Anjali | Basavaraj Shape-Induced Deformation, Capillary Bridging, and Self-Assembly of Cuboids at the Fluid-Fluid Interface[END_REF][START_REF] Soligno | Self-assembly of cubic colloidal particles at fluid-fluid interfaces by hexapolar capillary interactions[END_REF][START_REF] Song | Control of Orientation, Formation of Ordered Structures, and Self-Sorting of Surface-Functionalized Microcubes at the Air-Water Interface[END_REF][START_REF] Song | Reconfigurable Microcube Assemblies at the Liquid/Air Interface: The Impact of Surface Tension on Orientation and Capillary-Force-Interaction-Driven Assembly[END_REF][START_REF] Anzivino | Chains of cubic colloids at fluid-fluid interfaces[END_REF].

In both contexts, objects are considered small for gravity to be neglected. We can expect this hypothesis to fail as the size or density of the object increases. We therefore aim to provide criteria to account for gravity. We expect the relevant dimensionless quantity for cubic objects interacting with a soap film to be the Bond number Bo = ρgr 2 s /2γ, with ρ the object density, g the standard gravitational acceleration, r s the half side-length and γ the surface tension. (Note that ρ should be replaced by ρ -ρ and 2γ by γ for the interaction at the interface of two liquids of the same density ρ .) We have previously showed that a cube with a Bond number of 1 is reoriented by a soap film to a flat orientation as it descends through a horizontal soap film [START_REF] Davies | Simulating the interaction between a descending super-quadric solid object and a soap film[END_REF]. We now ask what are the conditions for a cubic object to be trapped or to detach from a soap film in terms of the Bond number and its initial orientation? Although we will probe the final orientation of a single cube in terms of its Bond number, this final orientation is crucial to determining the interaction between multiple cubes trapped at an interface.

In this work, we study the interaction between a soap film and a cube both in simulation and experiments. We will consider what happens when we vary the Bond number and the initial orientation of the cube systematically while the contact angle is set to 90 • . We aim to probe how a cube is reoriented by a soap film it falls onto, and when does it become trapped by the soap film. In section 2, we include details of the methodology for our simulations and experiments. The result are presented in section 3 and section 4 includes our conclusions.

Methodology

Simulations

We use the Surface Evolver [START_REF] Brakke | The Surface Evolver[END_REF] to simulate the interaction between a horizontal soap film and a cube with curved corners and edges. The cube is defined by the superquadric equation

(x -x 0 ) λ + (y -y 0 ) λ + (z -z 0 ) λ = r λ s , (1)
where (x 0 , y 0 , z 0 ) denotes its centre of mass and r s is the minimum distance between its centre and its surface. The shape parameter λ, which determines the curvature of the edges and corners of the cube, is kept fixed at λ = 10 throughout. The volume of the superquadric cube defined by equation ( 1) is calculated using the formula

V s = 8r 2 s 3λ 2 Γ 1 λ 3 Γ 3 λ (2)
(as derived in [START_REF] Jaklic | Segmentation and Recovery of Superquadrics, volume 20, chapter Superquadrics and Their Geometric Properties[END_REF]), where the gamma functions Γ (a/λ) = ∞ 0 t a/λ-1 e -t dt are computed numerically. As in previous work [START_REF] Davies | Simulating the interaction between a descending super-quadric solid object and a soap film[END_REF], we set the volume of the cube to that of a sphere with radius 1, that is V s = 4 3 π, so that its sides have a half length of r s ≈ 0.817. The surface area of the cube is A s ≈ 14.24. We vary the Bond number by varying the weight of the cube, choosing values of Bo ∈ [0.1, 1.0].

The cube is initially positioned just above a horizontal soap film. The initial orientation of the cube is prescribed (and varied), and is described by the angles α i between the normal vectors N i that define it and the vertical axis, for i = 1, 2, 3 (see figure 1a). The boundary of the soap film is fixed on a horizontal circular ring of radius r c = 4 (chosen to be large enough so that boundary effects are negligible) and height h = 5.

As in previous work [START_REF] Davies | Simulating the interaction between a descending super-quadric solid object and a soap film[END_REF], when the object is in contact with the soap film, we assume that its motion is slow and over-damped so that inertial effects can be neglected, and a quasi-static model is appropriate. In this case, the soap film is at equilibrium between small increments in the position and orientation of the cube.

The soap film exerts a network force F n on the solid object due to the pull of surface tension, which is calculated geometrically as

F n = 2γ i l i n i . (3)
Here l i denotes the contact length of the triangular facet i of the soap film that is in contact with the object and n i denotes the unit normal vector to the surface at (x i , y i , z i ), the mid-point of edge i (see figure 1b),

z (a) (b)
Figure 1: (a) A rounded cube with side length 2r s is initially positioned just above a horizontal soap film that is pinned at the boundary to a circle with radius r c and height z = h. The orientation of the cube is prescribed and then traced by recording the angles that the normal vectors that define it, N 1 , N 2 and N 3 (obscured by the cube) make with the vertical axis. (b) The film exerts a network force on the cube due to surface tension; for each edge i of the triangulated surface of the soap film that contacts the cube, an outward network force is exerted in the normal direction, n i over its length l i . and 2γ is the surface tension of the film (which we set to 1 in our simulations). Since the surface of the film contacts the object at 90 • , it applies a network torque, τ n :

τ n = 2γ i l i r i × n i , (4)
where r i = (x i -x 0 , y i -y 0 , z i -z 0 ) is the vector that connects the centre of the object with the midpoint of edge i. Therefore the force on the cube is

F = -mg z + F n (5)
where z denotes the unit vector in the positive z direction. At each time step of a simulation, the soap film is equilibrated while the position and orientation of the cube is fixed. The cube is then moved by ε F and rotated it by ε τ where ε is a small parameter that sets the effective time scale of the simulations. We choose ε = 2.5 × 10 -3 throughout this work. The soap film is then re-equilibrated and the process is repeated.

The equilibrium states are found by minimizing the energy of the film, which is proportional to its surface area A, using gradient descent and conjugate gradient iterations. Our energy minimization procedure continues until convergence of the energy E to within a tolerance of 1 × 10 -5 has been achieved, which can take up to 10,000 iterations. The iterations are interspersed with upkeep of the tessellation and checks for soap film detachment. We also apply small perturbations to the surface of the soap film during the equilibration process by jiggling the vertices slightly. In this case, a random displacement is applied to each vertex independently using a Gaussian distribution with a deviation of 0.02 times the mean edge length of the triangular mesh [START_REF] Brakke | The Surface Evolver[END_REF]. This perturbation was found to be robust enough to optimize the process of reaching a minimum energy for the soap film under the given constraints.

Experiments

Soap films are formed from a SLES/CAPB/SLES mixture [START_REF] Golemanov | Surfactant Mixtures for Control of Bubble Surface Mobility in Foam Studies[END_REF]: we mixed 6.6% of sodium lauryl ether sulfate (SLES) and 3.4% of cocamidopropylbetaine (CAPB) in mass in ultrapure water; we then dissolved 0.4% in mass of myristic acid (MAc), by stirring and heating at 60 • C for one hour, and we diluted 20 times this solution with water. The surface tension of the final solution equals 23.8 mN.m -1 . This solution is known for its high dissipative properties [START_REF] Denkov | The role of surfactant type and bubble surface mobility in foam rheology[END_REF] which prevent any significant flow from the film to the meniscus at the contact with the object at the time scale of one experiment.

A second solution was obtained by dissolving tetradecyl trimethyl ammonium bromide (TTAB) into deionizedwater. The concentration in TTAB was 3 g/l for the aqueous solution and the surface tension equals 38 mPa.m [START_REF] Cohen | Inertial Mass Transport and Capillary Hydraulic Jump in a Liquid Foam Microchannel[END_REF][START_REF] Cohen | Drop coalescence and liquid flow in a single Plateau border[END_REF].

In both solutions a small quantity of fluorescein sodium salt (0.5 g/L) was added for better visualization. The soap films are formed and held on a circular wire, 8 cm in diameter.

Cubes are 3D printed in polylactic acid (PLA) with a radius ranging between 0.5 and 2.8 mm. The density of PLA is 1.25 g.cm -3 . This leads to Bond numbers between 0.06 and 1.9. Cubes are initially pre-wetted and are held with tweezers at the contact with the soap film or slightly above it. For each cube, 10 experiments are performed with various initial orientations. Once the cubes are released, the dynamics is recorded with a color highspeed camera at 50-500 frames per second.

Results

In this section, we present results from simulations and experiments. First we investigate the final orientation of a cube that is trapped by a soap film and how this varies with the Bond number. We go on to determine the critical value of the Bond number for which a cube detaches from a soap film and continues its descent under gravity. Throughout this section, we will refer to three main (stable or meta-stable) orientations for a cube in contact with a soap film: (i) a flat orientation, where one of the angles α i is zero and the other two are equal to π/2, (ii) a diagonal orientation in which two of the angles are equal to π/4 and the other is π/2, and (iii) a rotated orientation where all three angles are equal. We will from now on refer to orientations (i)-(iii) as {100}, {110} and {111} respectively.

In simulations, we varied the initial orientation of the cube as well as the Bond number. We set the initial orientation of the cube by rotating it around the x-axis by an angle nπ/10, and then around the y-axis by mπ/10, considering all possible combinations of n, m ∈ {0, 1, .., 5}. We vary the Bond number between 0.1 and 1.

Figure 2 includes snapshots from three simulations that show how the orientation of a cube changes over time as it interacts with a soap film. The initial orientation of the cube is the same for the three simulations shown, while the Bond numbers are Bo = 0.1, 0.375, 0.5. The cube becomes trapped by the soap film in all three cases, but its final orientation varies with the Bond number. The cube settles into the {111}, {110} and {100} orientations in order of increasing Bond numbers. We will try to explain this surprising result later, but let us first investigate whether or not the final orientation of a cube that becomes entrapped by a soap film is independent of its initial orientation.

We investigate how the orientation of a cube that becomes entrapped by a soap film varies from all initial orientations considered in our simulations. Figure 3a shows that the most likely settling orientation for a cube with Bond number Bo = 0.375 is the {111} orientation. In this case, the cube only settles into the {110} or {100} orientations when the initial orientations are very close to these (as shown in figure 2). For a Bond number of Bo = 0.25 (data not shown), we observed that a cube rotates to the {111} orientation unless it is initially very close to the {100} orientation. When Bo = 0.1, a cube always rotated to the {111} orientation, even when it started from the {100} orientation. Figure 3b shows that a cube with Bond number Bo = 0.5 always rotated to the {100} orientation. This is also the most likely final orientation for cubes with a higher Bond number, even those that detach from a soap film, as shown in previous work [START_REF] Davies | Simulating the interaction between a descending super-quadric solid object and a soap film[END_REF]. The final orientation of a cube that interacts with a horizontal soap film versus the Bond number is summarised in full in figure 4.

We find a critical value of the Bond number, denoted Bo c , for which a cube's weight is no longer supported by the soap film and detachment occurs. In simulations, we found that a cube with Bond number 0.6 was trapped by the soap film, while a cube with a Bond number of 0.7 was not supported and detached from the soap film. We therefore extrapolate that the critical value of the Bond number is Bo c = 0.65 ± 0.05 in simulations (as shown in figure 4).

We now seek to validate these results on the final orientation of an entrapped cube as well as the critical value of the Bond number by comparing with experiments.

First we present the results with the SLES/CAPB/SLES solution. In experiments the Bond number is varied between 0.06 and 1.9 and we focus on the range of interest in figure 5. Only two stable configurations are observed. For 0.06 ≤ Bo ≤ 0.25, the {111} orientation is systematically observed, but we noticed that the {110} orientation can appear transiently (Movie 1 in the Supplemental Material). For Bo = 0.38 both the {111} and the {100} orientations were observed, which indicates a significant change around this value in agreement with simulation. For Bo = 0.55, only the {100} orientation is observed, but the cube detaches in some cases. This value of the Bond number corresponds to the entrapment/detachment transition in experiments and the cube detaches systematically for larger Bo (Movie 2 in the Supplemental Material). This value for the critical Bond number is slightly lower than Bo c ≈ 0.65 found in simulations. This lower value can be attributed either to the small amount of liquid that collects on the lower surface of a cube in experiments, thus contributing to an additional weight that is not accounted for by the Bond number, or to the inertia of the object that is neglected in the model.

A second round of experiments were performed with the TTAB solution. Here we observe a flow from the soap film to the meniscus due to capillary suction [START_REF] Sonin | Role of surface elasticity in the drainage of soap films[END_REF]. This flow makes a difficult comparison between experiments and simulations since the later assume an infinitely small meniscus. As an example in figure 6, Bo = 0.24 should lead to the diagonal configuration as expected from the first set of experiment with the SLES/CAPB/SLES solution and from the simulations. This is indeed observed initially, but after a few seconds the growing meniscus triggers a transition toward a flat orientation with the bottom/top side of the meniscus connecting the bottom/top side of the cube respectively. This shows that the notion of final orientation depends on the timescale considered: it depends on the time it takes for the liquid to rearrange inside the soap film and the meniscus. As observed, this dynamics is strongly correlated with the nature of the surfactant solution.

We now return to simulations to discuss the multiple scenarios at play in more detail. To explain why a cube with a given Bond number is rotated to a particular final orientation by a horizontal soap film, we consider how the energy of the soap film varies in our simulations. Figure 7a shows how the energy of the film varies with the height of the descending cube. The data shown here is for Bond numbers ranging between 0.1 and 0.75, with the cube initially oriented in the {100}, {110} or {111} orientations. The energy is given in terms of a dimensionless absorption energy per unit area, E 1 /2γA s , where E 1 = 2γ(A -πr 2 c ) is the difference between the energy of the soap film in its current state and at an unperturbed state. As already noted, a cube with Bo = 0.75 is too heavy to become entrapped by the soap film, but it is also not reoriented by the soap film from any of these three main orientations. This is shown by the three distinct curves labeled A, B and C in figure 7a. A drop in the surface energy occurs when the cube touches and perturbs the soap film for the first time, and the energy continues to decrease to a minimum when the height of the cube is perfectly aligned with the height at which the outer boundary of the soap film is pinned. At this point, the energy of the soap film is at its lowest for the cube in the {111} orientation. Once the cube has passed this point, the surface energy increases as the film becomes stretched, before the cube eventually detaches from the soap film. The film is stretched the most before detachment for the cube in the {100} orientation.

It can be seen in figure 7a that a cube with Bo ≤ 0.5 can be rotated by a soap film from one metastable orientation to a more energetically favourable orientation. This is shown by the jumps from one of the distinct energy curves (labelled A, B and C) to another. These transitions are surprising given the symmetric deformation caused to the soap film by each of the three main orientations. They are triggered by small imperfections in our mesh approximation of the structure. The orientation of an entrapped cube that results in perturbing the soap film so that its energy is minimised is dependent on the height reached by the cube, which in turn depends on the initial orientation of the cube and the Bond number. If the centre of the cube has advanced less than r s past the height h at which the film is pinned at the boundary, then it is almost always rotated by the soap film to the {111} orientation, which is the orientation that achieves the minimum energy for the soap film here. The stability of the {100} and {110} orientations decreases with the Bond number, and as a result a cube with Bo = 0.1 is rotated from both of these orientation into the stable {111} orientation. If a cube advances beyond this height, then it is always rotated by the soap film to the {100} orientation, which results in the most energetically favourable perturbation of the soap film here. It is interesting to note what happens when the settling position of the cube is very close to h -r s , as happens in a couple of simulations where Bo = 0.375. The three energy curves in figure 7a intersect here, thus the three main orientations are as energetically favourable as each other at this point. This explains why a cube with this Bond number can settle on the soap film with the {110} orientation from a similar initial orientation (as seen in figure 2).

In a second set of simulations we probe the stability of the three main final orientations in the trapped regime. At mechanical equilibrium, we have that F = 0 in Eq. ( 5), from which we can write the equilibrium Bond number Bo eq = F n z r 2 s /(2γV s ). Thus measuring the vertical component F n z of the drag force exerted by the film gives a direct measurement of the Bond number of the trapped object. For each (fixed) orientation, the object is artificially displaced through the soap film from first contact to detachment to sweep all the possible combinations of Bond numbers and surface energies. Figure 7b shows how the energy of the soap film varies with Bo eq for the three main cube orientations. This confirms that the orientation {111} is the most energetically favourable when the Bond number is small and that only the {100} orientation is achievable when the Bond number is greater than 0.5 with a critical Bond number around 0.69 that gives the maximal value before detachment. The {110} orientation is never optimal in terms of perturbing the soap film in a way that achieves a minimum energy. Note that for a given Bond number, there exists two branches for orientations {111} and {110} and only the lower one is stable.

Discussion and conclusions

We have observed in both simulations and experiments that the final orientation of a cube that becomes entrapped by a soap film depends on the Bond number, which determines the settling height of the cube. A light cube with a Bond number Bo < 0.375 becomes entrapped by a soap film and is almost always rotated to the {111} orientation. For these Bond numbers, a cube descends less than half its side length past the height of the soap film, and in this range it is the {111} orientation for the cube that perturbs the soap film so that its energy is minimized. A cube with a Bond number Bo ≈ 0.375 is almost always rotated to the {111} orientation as before. However if the cube has an initial orientation that is very close to either the {110} or {100} orientations, then it will rotate into the nearest of these orientations. The {110} orientation was rarely observed for other Bond numbers. This special case can be explained by the fact that the absorption energy of the soap film is very similar for all three cube orientations for Bo = 0.375. For larger Bond numbers, a cube that becomes entrapped in the soap film is always rotated to the {100} orientation, which provides the minimum energy for the soap film when the cube has descended well beyond the height of the soap film. These results should help in patterning capillary assemblies with prescribed orientations. Orientation {111} is quite easily obtained as long as objects are light enough. On the opposite, observing the {100} orientation would require a precise control of the Bond number, between 0.5 and 0.65, while the {110} orientation would be almost impossible to obtain in a controlled way.

We found that the critical Bond number (Bo c ), above which a cube is no longer supported by a soap film, is less than one. This can be explained by a simple force balance argument. From our results, we know that the maximum drag force applied by a soap film on a cube occurs when the cube is in the {100} orientation and the film contacts it over its upper horizontal face, forming a stretched catenoidal neck (as seen in figure 5). The catenoidal shape of the soap film means that the contact with the cube is well approximated by a circle. For a cube with sharp edges and corners, the largest circle that can be inscribed onto this surface has radius r s , and since we assume that the soap film meets the obstacle at 90 • , an upper bound for the network force is given by 4πr s γ. Balancing this out with the weight (mg) of the cube gives an approximation of the critical Bond number of

Bo c = mgr 2 s 2γV s ≈ (4πr s γ)r 2 s 2γV s = 2πr 3 s V s ≈ 0.82.
This approximation can be improved if we also take into account the curved edges and corners of the cubes that we used. This reduces the size of the largest circle that can be inscribed into the horizontal region of the upper surface of the cube. For the superquadric definition that we used for the cube (in equation 1), with λ = 10, a horizontal plane z = z 0 + r s approximates the upper surface of the cube (to within 1% accuracy) over a diameter of 1.5r s . Therefore our estimate for the critical Bond number reduces to Bo c = 1.5πr 3 s /V s ≈ 0.61, which is similar to what we observed in our simulation and experiments. In our simulations, we kept the contact angle between the soap film and the cube fixed at 90 • , as this provided a reasonable approximation for what we observed in experiments as long as the dynamics is quasistatic [START_REF] Raufaste | Deformation of a Free Interface Pierced by a Tilted Cylinder[END_REF][START_REF] Raufaste | Deformation of a Free Interface Pierced by a Tilted Cylinder: Variation of the Contact Angle[END_REF]. This contact angle is clearly an important parameter, and changing it will surely affect the final orientation of entrapped cubes, as already observed in the small Bo limit [START_REF] Morris | Modelling the self orientation of particles in a film[END_REF][START_REF] Sharp | Adsorption of shapeanisotropic and porous particles at the air-water and the decane-water interface studied by the gel trapping technique[END_REF], as well as the value of the critical Bond number Bo c . For instance, we predict that decreasing the contact angle (meaning the liquid on top has a larger affinity with the object than the liquid at the bottom in the case of an interface between two liquids) will increase the value of the critical Bond number (Bo c ), but it is not clear what would become the most dominant final orientation of a cube as it becomes entrapped by a soap film. We did not consider in detail what happens when a cube detaches from a soap film. We observed in experiments that bubble entrainment occurs as it does for a spherical particle [START_REF] Chen | The distortion of a horizontal soap film due to the impact of a falling sphere[END_REF][START_REF] Cox | Bubble entrainment by a sphere falling through a horizontal soap foam[END_REF]. We expect that the size of bubbles generated during detachment will depend on the contact angle as well as the orientation and the sharpness of the edges and corners of the cube. A more detailed investigation of the effect of film thickness is also left as future work. Finally, we note a fundamental difference between cubes at the interfaces of a soap film and at the interface between two immiscible liquids. In the former case the presence of a meniscus can trigger effects that are not accounted for by our model. As emphasized with experiments, a meniscus triggers a flow, with a speed that depends on the nature of the surfactant and on the viscosity of the solution. At a large time scale we can expect the width of the meniscus to be comparable to the size of the cube, and a change of orientation for the cube is possible. Conversely, capillary assemblies at the interface between two immiscible liquids should be stable with respect to this process.

Figure 2 :

 2 Figure2: Snapshots from simulations in which a cube, with an initial orientation that is just off the {110} orientation (with α 1 = π/2, α 2 = π/5 and α 3 = 4π/5), interacts with a soap film for Bond numbers Bo = 0.1, 0.375, 0.5. The snapshots in the left column are taken after 200 time steps, while the centre and right columns are after 500 and 1,000 time steps respectively. The right column represents the final position of the cubes.

Figure 3 :Figure 4 :

 34 Figure3: The orientation of a cube with (a) Bo = 0.375 and (b) Bo = 0.5, given in terms of the angles (α 2 , α 3 ), is tracked over time from all initial orientations considered in our simulations. Note that the pairwise angles (0, π/2), (π/2, π/2) and (π/2, 0) represent the {100} orientation here, (π/4, π/2) and (π/2, π/4) both denote the {110} orientation, and the {111} orientation is given by tan -1 1/ √ 2 , tan -1 1/ √ 2 .

Figure 5 :

 5 Figure 5: Stable or transient orientations for various Bo in the range 0.2 -0.8 found in experiments.

Figure 6 :

 6 Figure 6: Evolution of the orientation due to the growth of the meniscus. TTAB solution, r s = 1.25, Bo = 0.24, 5 s between the first and last images (Movie 3 in Supplemental Material).

Figure 7 :

 7 Figure7: The variation of the dimensionless absorption energy per unit area of the soap film (E 1 /2γA s ) versus (a) the normalised height of the cube, (z 0 -h)/r s , and (b) the equilibrium Bond number Bo eq . In (a), cubes of different Bond numbers descend (from left to right on the height-axis) from the three main orientations {100}, {110} and {111}, labelled A, B and C respectively. The markers denote the settling position of a cube that becomes entrapped by the soap film. In (a) and (b), trajectories obtained for the three main orientations {100}, {110} and {111} if fixed are displayed as black solid, dotted and dash-dotted lines respectively. In (b) the markers of the final orientations obtained in (a) are represented too.

Acknowledgments

We thank Simon Cox for stimulating discussions and Ken Brakke for developing and maintaining the Surface Evolver. TD acknowledges Supercomputing Wales facilities. CR acknowledges support by the French government, through the National Research Agency (ANR-20-CE30-0019) and through the UCA JEDI Investments in the Future project of the National Research Agency (ANR-15-IDEX-01).

Conflicts of interest

There are no conflicts to declare.